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