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