!    USE MKL95_LAPACK
!    USE mkl95_blas
    IMPLICIT   NONE
    
    INTEGER*4, parameter :: N     = 16      ! size of problem
    integer*4, parameter :: LIMIT = 11      ! number of itterations
!
    INTEGER    I, J, k, ERROR_CODE

    real*8 A(N,N), B(N,1),  AA(N,N), X(N,1) , XX(N,LIMIT), E(N,1)
    real*8 ErrorMax, s, determinant, v
!
    real*16 acum    
!
    INTERFACE
     SUBROUTINE MATRIX_SOLVE (MATRIX_IN, RHS, determinant, ERROR_CODE)
!
      real*8,   dimension(:,:), intent (in)    :: matrix_in
      real*8,   dimension(:,:), intent (inout) :: rhs
      real*8,                   intent (out)   :: determinant
      integer*4,                intent (out)   :: error_code
!
     END SUBROUTINE MATRIX_SOLVE
    END INTERFACE
!
!   BEGIN Read Input File 
    OPEN(111,FILE='TDReview_FEM-K_REDUCE.csv',STATUS='UNKNOWN', ACTION='READ')
    REWIND(111)
    READ(111,*)	((A(i,j), j=1,N), i=1,N)      
    CLOSE(111)
    
    OPEN(112,FILE='TDReview_FEM-LOAD_REDUCE.csv',STATUS='UNKNOWN', ACTION='READ')
    REWIND(112)
    READ(112,*)	((B(i,j), j=1,1), i=1,N)      
    CLOSE(112)
!   END Read Input File 

! check A for symmetry
    s = 0
    v = 0
    do i = 1,n
      do j = i+1,n
        s = max (s, abs(a(i,j)-a(j,i)) )
        v = max (v, abs(a(i,j)) )
      end do
    end do
    write (*,*) 'Max symmetry error =',s
    write (*,*) 'Max coefficient    =',v
!    
    XX=0D0
    DO i=1,LIMIT
!
        AA     = A                   ! Reasign AA which is modified when calling GESV 
!
!        X(:,1) = SUM(XX,DIM=2)
        do j = 1,n
          acum = 0
          do k = 1,i
            acum = acum + xx(j,k)
          end do
          X(j,1) = acum
        end do
!
!        E = B - MATMUL(A,X)
        v = 0
        do j = 1,n
          acum = B(j,1)
          do k = 1,n
            acum = acum - A(j,k)*X(k,1)
            v = max ( v, abs(A(j,k)*X(k,1)) )
          end do
          E(j,1) = acum
        end do

        ErrorMax = MAXVAL(ABS(E))
        s = 0
        do k = 1,N
          s = s + ABS (E(k,1))
        end do
        WRITE (*,*) 'LOOP',I
        write (*,*) ' Max error term  :', ErrorMax
        write (*,*) ' Avg error term  :', s/n
        write (*,*) ' Max coefficient :', v
        IF (ErrorMax<1.D-10) EXIT
!
!        CALL GESV (AA, E)
        call MATRIX_SOLVE (AA, E, determinant, ERROR_CODE)
        write (*,*) ' Solve error =',error_code
        write (*,*) ' determinant =',determinant
        XX(:,i) = E(:,1)
    END DO

!    OPEN(111,FILE='TDReview_FEM-DISP_VEC.csv',STATUS='UNKNOWN', ACTION='WRITE')

    DO i=1,N
        WRITE (*,"(1000(1x,F55.19,','))") X(i,:)
    END DO
    

    END

     SUBROUTINE MATRIX_SOLVE (MATRIX_IN, RHS, determinant, ERROR_CODE)
!
!  solves for N simultaneous equations for M right hand sides
!  uses partial pivoting to select the next row to reduce for the
!  next active column
!
!  Inputs
!   MATRIX_IN   equations are stored in rows
!   RHS         the right hand sides
!
!  Results
!   RHS         the solution vectors
!   determinant the determinant of MATRIX_IN
!   ERROR_CODE  error flag for solution
!               0   for solution
!               1-n for zero row in matrix
!               999 for zero diagonal during solution
!
!   rank        the number of simultaneous equations provided
!   cases       the number of right hand sides provided
!
      real*8,   dimension(:,:), intent (in)    :: matrix_in
      real*8,   dimension(:,:), intent (inout) :: rhs
      real*8,                   intent (out)   :: determinant
      integer*4,                intent (out)   :: error_code
!
      real*16,  parameter                   :: accuracy = 1.0D-20
      real*16,  dimension(:,:), allocatable :: matrix
      real*16                               :: C
      INTEGER*4 rank, cases, I, J, K, L
      INTRINSIC ABS, SIZE
!
      rank  = size (matrix_in,1)
      cases = size (rhs,2)
!
      ERROR_CODE   = -1
      determinant  = 1.
      if (rank < 1) return
!
!    NORMALISE THE ROWS WITH MAXIMUM ROW VALUE
!
      allocate (matrix(rank,rank))
      DO I = 1,rank
         C = 0.
         DO J = 1,rank
            IF (ABS(MATRIX_IN(i,j)) > ABS(C)) C = MATRIX_IN(i,j)
         END DO
         determinant = determinant * C
         IF (ABS(C) < accuracy) then
            WRITE(98,903) I
  903       FORMAT(' **** ZERO ROW ',I0,' IN EQUATION SOLUTION ****')
            ERROR_CODE  = I
            determinant = 0.
            RETURN
         end if
         DO J = 1,rank
            MATRIX(i,j) = MATRIX_IN(i,j) / C
         END DO
         DO L = 1,cases
            RHS(I,L) = RHS(I,L) / C
         END DO
      END DO
!
!    NOW DO FOWARD REDUCTION
!
      DO I = 1,rank
!
!       SELECT FOWARD PIVOT FROM COLUMN "I"
!
         C = 0.
         DO J = I,rank
            IF (ABS(MATRIX(j,i)) <= C) CYCLE
            K = J
            C = ABS(MATRIX(j,i))
         END DO
!
         IF (C < accuracy) then
            WRITE (98,901)
  901       FORMAT(' **** ZERO DIVISION IN EQUATION SOLUTION ****')
            ERROR_CODE  = 999
            determinant = 0.
            RETURN
         end if
         IF (K .NE. I) THEN
!
!          SWAP THE ROWS
!
            DO J = 1,rank
               C = MATRIX(i,j)
               MATRIX(i,j) = MATRIX(k,j)
               MATRIX(k,j) = C
            END DO
            DO L = 1,cases
               C = RHS(I,L)
               RHS(I,L) = RHS(k,L)
               RHS(k,L) = C
            END DO
            determinant = -determinant          ! CHANGE SIGN OF DETERMINANT
         END IF
!
!       NOW REDUCE THIS ROW
!
         determinant = determinant * MATRIX(I,I)
!
         DO J = I+1,rank
            C = MATRIX(j,i)/MATRIX(I,I)
            DO K = I+1,rank
               MATRIX(j,k) = MATRIX(j,k) - MATRIX(i,k) * C
            END DO
            DO L = 1,cases
               RHS(J,L) = RHS(J,L) - RHS(I,L) * C
            END DO
         END DO
!
      END DO
!
!    NOW BACKWARD SUBSTITUTION
!
      DO L = 1,cases
         DO I = rank,1,-1
            C = 0.
            DO J = I+1,rank
               C = C + MATRIX(i,j) * RHS(J,L)
            END DO
            RHS(I,L) = (RHS(I,L)-C) / MATRIX(I,I)
         END DO
      END DO
!
      ERROR_CODE = 0
      RETURN
!
      END SUBROUTINE MATRIX_SOLVE
  