C
C     SLE.F   VER 1.0  NOVEMBER 24, 2011
C
C     BY TOM IRVINE
C
C
C     SOLVE A SYSTEM OF LINEAR EQUATIONS
C     OF THE FORM:
C
C         A*X=B
C
C23456789
      EXTERNAL DGESV
      REAL :: START, FINISH
      DOUBLE PRECISION, ALLOCATABLE :: A(:,:),B(:,:),X(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: AOLD(:,:),BOLD(:,:)
      DOUBLE PRECISION :: ALPHA,BETA
      INTEGER DeAllocateStatus
      INTEGER :: I,J,N,INFO,NRHS
      INTEGER :: LDA,LDB
      INTEGER, ALLOCATABLE :: IPIV(:)
      CHARACTER :: A_filename*20,B_filename*20

C
      PRINT '("")'
      PRINT '( "Enter the filename name of matrix A")'
      READ(*,*)A_filename
      PRINT '( "Enter the filename of matrix B")'
      READ(*,*)B_filename
C
      OPEN(UNIT=10,FILE=A_filename,STATUS='OLD',ACTION='READ')
      OPEN(UNIT=11,FILE=B_filename,STATUS='OLD',ACTION='READ')
C
      PRINT '("")'
      PRINT '( "Enter the number of rows for matrix A")'
      READ(*,*)N

      PRINT '( "Enter the number of columns for matrix B")'
      READ(*,*)NRHS
C
      LDA=N
      LDB=N
      ALLOCATE(A(LDA,N))
      ALLOCATE(AOLD(LDA,N))
      ALLOCATE(B(LDB,NRHS))
      ALLOCATE(BOLD(LDB,NRHS))
      ALLOCATE(X(LDB,NRHS))
      ALLOCATE(IPIV(N))
C
      DO I=1,N
            READ(10,*) (A(I,J),J=1,N)
      ENDDO
      DO I=1,N
            READ(11,*) (B(I,J),J=1,NRHS)
      ENDDO
C
      AOLD=A
      BOLD=B
C
      CALL CPU_TIME(START)
C
      CALL DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
C
C     Check for the exact singularity.
C
      CALL CPU_TIME(FINISH)
C
      IF( INFO.GT.0 ) THEN
        WRITE(*,*)'The diagonal element of the triangular factor of A,'
        WRITE(*,*)'U(',INFO,',',INFO,') is zero, so that'
        WRITE(*,*)'A is singular; the solution could not be computed.'
        STOP
      END IF
C
      X=B
C
      PRINT '("")'
      PRINT '("Time = ",F6.3," seconds")',FINISH-START
C
      PRINT '("")'
      PRINT '( "Display Results?  1=yes 2=no ")'
      READ(*,*)IDR
C
      IF(IDR.EQ.1)THEN
         PRINT '("")'
         PRINT '("A")'
         DO 100 I=1,N
            WRITE(*,33) (AOLD(I,J),J=1,N)
100      CONTINUE
C
         PRINT '("")'
         PRINT '("B")'
         DO 300 I=1,N
            WRITE(*,33)(BOLD(I,J),J=1,NRHS)
300      CONTINUE
C
         PRINT '("")'
         PRINT '("X")'
         DO 200 I=1,N
            WRITE(*,33) (X(I,J),J=1,NRHS)
200      CONTINUE
      ENDIF
33    FORMAT(1000(E10.4,1X))
C
      DEALLOCATE (A, STAT = DeAllocateStatus)
      DEALLOCATE (B, STAT = DeAllocateStatus)
      DEALLOCATE (X, STAT = DeAllocateStatus)
      DEALLOCATE (AOLD, STAT = DeAllocateStatus)
      DEALLOCATE (BOLD, STAT = DeAllocateStatus)
C
      STOP
      END