PROGRAM example_3_1
  USE precisions; USE lu_dcmps, ONLY: lu_gauss, lu_solve
  IMPLICIT NONE

  INTEGER :: i, j, vp(3)
  REAL(KIND=wp), DIMENSION(3,3) :: A = RESHAPE(&
    (/ 33.0, -24.0, -8.0, 16.0, -10.0, -4.0, 72.0, -57.0, -17.0 /),(/3,3/)),&
                                   INVA6 = RESHAPE(&
    (/ -58.0, 48.0, 16.0, -16.0, 15.0, 4.0, -192.0, 153.0, 54.0 /),(/3,3/))
  REAL(KIND=wp), DIMENSION(3) :: b = (/ 359.0, -281.0, -85.0 /),&
                                 xbar = (/ -1.0, 2.0, 5.0 /)
  REAL(KIND=wp) :: AA(3,3), X(3,3), r(3), xtilde(3)

  PRINT*
  PRINT*, 'A e 6*INVA'
  PRINT*
  DO i =1, 3
    PRINT('(6f10.1)'), (A(i,j),j=1,3), (INVA6(i,j),j=1,3)
  END DO

  AA = A; X = 0.0_wp;
  DO i = 1, 3
    X(i,i) = 1.0_wp
  END DO

  CALL lu_gauss(AA,vp)
  DO i=1,3
    CALL lu_solve(AA,X(:,i),vp)
  END DO
  X = 6.0_wp*X
  PRINT*,
  PRINT*, '6*INVA calculada por lu_gauss'
  PRINT*
  DO i =1, 3
    PRINT('(3f10.1)'), (X(i,j), j=1,3)
  END DO
  PRINT*; PRINT*, 'lado direito b'; PRINT*; PRINT('(3f10.1)'), b
  xtilde = b; CALL lu_solve(AA,xtilde,vp)
  PRINT*; PRINT*, 'solucao exata'; PRINT*; PRINT('(3f10.1)'), xbar
  PRINT*; PRINT*, 'solucao calculada'; PRINT*
  DO i = 1, 3
    PRINT('(ES25.14)'), xtilde(i)
  END DO
  PRINT*; PRINT*, 'permutacoes'; PRINT*
  PRINT('(3i3)'), vp
!  PRINT*; PRINT*, 'residuo'; PRINT*
!  r = b - MATMUL(A+0.0_hp,X(:,1)+0.0_hp); PRINT('(3e10.2)'), r  
END PROGRAM example_3_1