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

  INTEGER :: i, j, vp(4)
  REAL(KIND=wp), DIMENSION(4,4) :: A = RESHAPE(&
    (/ 1.0, -2.0, 3.0, 1.0, -2.0, 1.0, -2.0, -1.0,&
     & 3.0, -2.0, 1.0, 5.0, 1.0, -1.0, 5.0, 3.0 /),(/4,4/)),&
                                   INVA52 = RESHAPE(&
    (/ -15.0, -38.0, -1.0, -6.0, -38.0, -20.0, -6.0, 16.0,&
     & -1.0, -6.0, -7.0, 10.0, -6.0, 16.0, 10.0, 8.0 /),(/4,4/))
  REAL(KIND=wp), DIMENSION(4) :: b = (/ 3.0, -4.0, 7.0, 8.0 /),&
                              xbar = (/ 1.0, 1.0, 1.0, 1.0 /)
  REAL(KIND=wp) :: AA(4,4), X(4,4), r(4), xtilde(4)

  PRINT*; PRINT*, 'A e 52*INVA'; PRINT*
  DO i =1, 4
    PRINT('(8f10.1)'), (A(i,j),j=1,4), (INVA52(i,j),j=1,4)
  END DO

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

  CALL lu_gauss(AA,vp)
  DO i=1,4
    CALL lu_solve(AA,X(:,i),vp)
  END DO
  X = 5.2e1_wp*X
  PRINT*; PRINT*, '52*INVA calculada'; PRINT*
  DO i =1, 4
    PRINT('(4f10.1)'), (X(i,j), j=1,4)
  END DO
  PRINT*; PRINT*, 'lado direito b'; PRINT*; PRINT('(4f10.1)'), b
  xtilde = b; CALL lu_solve(AA,xtilde,vp)
  PRINT*; PRINT*, 'solucao exata'; PRINT*; PRINT('(4f10.1)'), xbar
  PRINT*; PRINT*, 'solucao calculada'; PRINT*
  DO i = 1, 4
    PRINT('(ES30.14)'), xtilde(i)
  END DO
  PRINT*; PRINT*, 'permutacoes'; PRINT*
  PRINT('(4i3)'), vp
!  PRINT*; PRINT*, 'residuo'; PRINT*
!  r = b - MATMUL(A+0.0_hp,X(:,1)+0.0_hp); PRINT('(4e10.2)'), r  
END PROGRAM example_3_2