PROGRAM linear1
!
! Author:  Br. David Carlson
!
! Date:  December 31, 1999
!
! Last revised:  May 2, 2004
!
! This program uses the sgefs library routine to solve a system of N linear
! equations in N variables.

IMPLICIT NONE

INTEGER, PARAMETER::N = 3
INTEGER, PARAMETER::LDA = 3
INTEGER k
REAL, DIMENSION(LDA, N)::A, C
REAL, DIMENSION(N)::V, WORK
REAL::RCOND

INTEGER::IND
INTEGER, DIMENSION(N)::IWORK

! Could also initialize using reshape.
! Indices are always row, column.
A(1, 1) = 10.0
A(2, 1) = -3.0
A(3, 1) = 5.0
A(1, 2) = -7.0
A(2, 2) = 2.0
A(3, 2) = -1.0
A(1, 3) = 0.0
A(2, 3) = 6.0
A(3, 3) = 5.0

! Save a copy of matrix A:
C = A

V(1) = 7.0
V(2) = 4.0
V(3) = 6.0

WRITE (*, 20)
WRITE (*, 30)

DO k = 1, n
   WRITE (*, 40) k, V(k)
END DO

20 FORMAT (1X, /, 'Program to solve a system of linear equations.')
30 FORMAT (1X, 'The product of the matrix A and the solution vector should be:')
40 FORMAT (1X, 'b', I1, ' = ', E16.8)

CALL SGEFS(A, LDA, N, V, 1, IND, WORK, IWORK, RCOND) 
!  Error Messages Printed:
!    IND=-1  fatal    N is greater than LDA. 
!    IND=-2  fatal    N is less than 1.  
!    IND=-3  fatal    ITASK is less than 1 or greater than 4. 
!    IND=-4  fatal    The matrix A is computationally singular. 
!                     A solution has not been computed. 
!    IND=-10 warning  The solution has no apparent significance.
!                     The solution may be inaccurate or the matrix 
!                     A may be poorly scaled. 
 
WRITE (*, 50)

DO k = 1, n
   WRITE (*, 60) k, V(k)
END DO

WRITE (*, 70) IND
WRITE (*, 80) 1.0 / RCOND

50 FORMAT (1X, /, 'Approximate solution:')
60 FORMAT (1X, 'x', I1, ' = ', E16.8)
70 FORMAT (1X, 'Rough number of digits of accuracy: ', I8)
80 FORMAT (1X, 'Estimated condition number of matrix A: ', E16.5)

CALL check(n, C, V)

END PROGRAM


! Given:  n    Number of rows (and columns) in matrix a.
!         c    n by n matrix with real coefficients
!         b    Column vector with n real coefficients (approximate solution)
! Task:   To find and print the product of c by b.
! Return: Nothing.
SUBROUTINE check(n, c, b)

IMPLICIT NONE

INTEGER, INTENT(IN)::n 
REAL, DIMENSION(n, n), INTENT(IN)::c
REAL, DIMENSION(n), INTENT(IN)::b

INTEGER::row, col
REAL::temp

WRITE (*, 90)

DO row = 1, n
   temp = 0.0
   DO col = 1, n
      temp = temp + c(row, col) * b(col)
   END DO
   WRITE (*, 100) row, temp
END DO

90 FORMAT (1X, /, 'Product of matrix and approx. solution gives:')
100 FORMAT (1X, 'b', I1, ' = ', E16.8)

END SUBROUTINE