SUBROUTINE SGECO(A,LDA,N,IPVT,RCOND,Z) C***BEGIN PROLOGUE SGECO C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS CONDITION,FACTOR,LINEAR ALGEBRA,LINPACK,MATRIX C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Factors a real matrix by Gaussian elimination and estimates C the condition number of the matrix. C***DESCRIPTION C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C SGECO factors a real matrix by Gaussian elimination C and estimates the condition of the matrix. C C If RCOND is not needed, SGEFA is slightly faster. C To solve A*X = B , follow SGECO by SGESL. C To compute INVERSE(A)*C , follow SGECO by SGESL. C To compute DETERMINANT(A) , follow SGECO by SGEDI. C To compute INVERSE(A) , follow SGECO by SGEDI. C C On Entry C C A REAL(LDA, N) C the matrix to be factored. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C On Return C C A an upper triangular matrix and the multipliers C which were used to obtain it. C The factorization can be written A = L*U , where C L is a product of permutation and unit lower C triangular matrices and U is upper triangular. C C IPVT INTEGER(N) C an integer vector of pivot indices. C C RCOND REAL C an estimate of the reciprocal condition of A . C For the system A*X = B , relative perturbations C in A and B of size EPSILON may cause C relative perturbations in X of size EPSILON/RCOND . C If RCOND is so small that the logical expression C 1.0 + RCOND .EQ. 1.0 C is true, then A may be singular to working C precision. In particular, RCOND is zero if C exact singularity is detected or the estimate C underflows. C C Z REAL(N) C a work vector whose contents are usually unimportant. C If A is close to a singular matrix, then Z is C an approximate null vector in the sense that C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C LINPACK SGEFA C BLAS SAXPY,SDOT,SSCAL,SASUM C Fortran ABS,AMAX1,SIGN C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED SASUM,SAXPY,SDOT,SGEFA,SSCAL C***END PROLOGUE SGECO INTEGER LDA,N,IPVT(*) REAL A(LDA,*),Z(*) REAL RCOND C REAL SDOT,EK,T,WK,WKM REAL ANORM,S,SASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L C C COMPUTE 1-NORM OF A C C***FIRST EXECUTABLE STATEMENT SGECO ANORM = 0.0E0 DO 10 J = 1, N ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1)) 10 CONTINUE C C FACTOR C CALL SGEFA(A,LDA,N,IPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID C OVERFLOW. C C SOLVE TRANS(U)*W = E C EK = 1.0E0 DO 20 J = 1, N Z(J) = 0.0E0 20 CONTINUE DO 100 K = 1, N IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K)) IF (ABS(EK-Z(K)) .LE. ABS(A(K,K))) GO TO 30 S = ABS(A(K,K))/ABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = ABS(WK) SM = ABS(WKM) IF (A(K,K) .EQ. 0.0E0) GO TO 40 WK = WK/A(K,K) WKM = WKM/A(K,K) GO TO 50 40 CONTINUE WK = 1.0E0 WKM = 1.0E0 50 CONTINUE KP1 = K + 1 IF (KP1 .GT. N) GO TO 90 DO 60 J = KP1, N SM = SM + ABS(Z(J)+WKM*A(K,J)) Z(J) = Z(J) + WK*A(K,J) S = S + ABS(Z(J)) 60 CONTINUE IF (S .GE. SM) GO TO 80 T = WKM - WK WK = WKM DO 70 J = KP1, N Z(J) = Z(J) + T*A(K,J) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE TRANS(L)*Y = W C DO 120 KB = 1, N K = N + 1 - KB IF (K .LT. N) Z(K) = Z(K) + SDOT(N-K,A(K+1,K),1,Z(K+1),1) IF (ABS(Z(K)) .LE. 1.0E0) GO TO 110 S = 1.0E0/ABS(Z(K)) CALL SSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0E0 C C SOLVE L*V = Y C DO 140 K = 1, N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T IF (K .LT. N) CALL SAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (ABS(Z(K)) .LE. 1.0E0) GO TO 130 S = 1.0E0/ABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE U*Z = V C DO 160 KB = 1, N K = N + 1 - KB IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 150 S = ABS(A(K,K))/ABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (A(K,K) .NE. 0.0E0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0E0) Z(K) = 1.0E0 T = -Z(K) CALL SAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0 RETURN END SUBROUTINE SGEFA(A,LDA,N,IPVT,INFO) C***BEGIN PROLOGUE SGEFA C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS FACTOR,LINEAR ALGEBRA,LINPACK,MATRIX C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Factors a real matrix by Gaussian elimination. C***DESCRIPTION C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C SGEFA factors a real matrix by Gaussian elimination. C C SGEFA is usually called by SGECO, but it can be called C directly with a saving in time if RCOND is not needed. C (Time for SGECO) = (1 + 9/N)*(Time for SGEFA) . C C On Entry C C A REAL(LDA, N) C the matrix to be factored. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C On Return C C A an upper triangular matrix and the multipliers C which were used to obtain it. C The factorization can be written A = L*U , where C L is a product of permutation and unit lower C triangular matrices and U is upper triangular. C C IPVT INTEGER(N) C an integer vector of pivot indices. C C INFO INTEGER C = 0 normal value. C = K if U(K,K) .EQ. 0.0 . This is not an error C condition for this subroutine, but it does C indicate that SGESL or SGEDI will divide by zero C if called. Use RCOND in SGECO for a reliable C indication of singularity. C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C BLAS SAXPY,SSCAL,ISAMAX C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED ISAMAX,SAXPY,SSCAL C***END PROLOGUE SGEFA INTEGER LDA,N,IPVT(*),INFO REAL A(LDA,*) C REAL T INTEGER ISAMAX,J,K,KP1,L,NM1 C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C C***FIRST EXECUTABLE STATEMENT SGEFA INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = ISAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (A(L,K) .EQ. 0.0E0) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0E0/A(K,K) CALL SSCAL(N-K,T,A(K+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL SAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. 0.0E0) INFO = N RETURN END SUBROUTINE SGEFS(A,LDA,N,V,ITASK,IND,WORK,IWORK,RCOND) C***BEGIN PROLOGUE SGEFS C***DATE WRITTEN 800317 (YYMMDD) C***REVISION DATE 870916 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS GENERAL SYSTEM OF LINEAR EQUATIONS,LINEAR EQUATIONS C***AUTHOR VOORHEES, E., (LOS ALAMOS NATIONAL LABORATORY) C***PURPOSE SGEFS solves a GENERAL single precision real C NXN system of linear equations. C***DESCRIPTION C C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Subroutine SGEFS solves a general NxN system of single C precision linear equations using LINPACK subroutines SGECO C and SGESL. That is, if A is an NxN real matrix and if X C and B are real N-vectors, then SGEFS solves the equation C C A*X=B. C C The matrix A is first factored into upper and lower tri- C angular matrices U and L using partial pivoting. These C factors and the pivoting information are used to find the C solution vector X. An approximate condition number is C calculated to provide a rough estimate of the number of C digits of accuracy in the computed solution. C C If the equation A*X=B is to be solved for more than one vector C B, the factoring of A does not need to be performed again and C the option to only solve (ITASK .EQ. 2) will be faster for C the succeeding solutions. In this case, the contents of A, C LDA, N and IWORK must not have been altered by the user follow- C ing factorization (ITASK=1). IND will not be changed by SGEFS C in this case. Other settings of ITASK are used to solve linear C systems involving the transpose of A. C C Argument Description *** C C A REAL(LDA,N) C on entry, the doubly subscripted array with dimension C (LDA,N) which contains the coefficient matrix. C on return, an upper triangular matrix U and the C multipliers necessary to construct a matrix L C so that A=L*U. C LDA INTEGER C the leading dimension of the array A. LDA must be great- C er than or equal to N. (terminal error message IND=-1) C N INTEGER C the order of the matrix A. The first N elements of C the array A are the elements of the first column of C the matrix A. N must be greater than or equal to 1. C (terminal error message IND=-2) C V REAL(N) C on entry, the singly subscripted array(vector) of di- C mension N which contains the right hand side B of a C system of simultaneous linear equations A*X=B. C on return, V contains the solution vector, X . C ITASK INTEGER C If ITASK=1, the matrix A is factored and then the C linear equation is solved. C If ITASK=2, the equation is solved using the existing C factored matrix A and IWORK. C If ITASK=3, the matrix is factored and A'x=b is solved C If ITASK=4, the transposed equation is solved using the C existing factored matrix A and IWORK. C If ITASK .LT. 1 or ITASK .GT. 4, then the terminal error C message IND=-3 is printed. C IND INTEGER C GT. 0 IND is a rough estimate of the number of digits C of accuracy in the solution, X. C LT. 0 see error message corresponding to IND below. C WORK REAL(N) C a singly subscripted array of dimension at least N. C IWORK INTEGER(N) C a singly subscripted array of dimension at least N. C RCOND REAL C estimate of 1.0/cond(A) C C Error Messages Printed *** C C IND=-1 fatal N is greater than LDA. C IND=-2 fatal N is less than 1. C IND=-3 fatal ITASK is less than 1 or greater than 4. C IND=-4 fatal The matrix A is computationally singular. C A solution has not been computed. C IND=-10 warning The solution has no apparent significance. C The solution may be inaccurate or the matrix C A may be poorly scaled. C C***REFERENCES SUBROUTINE SGEFS WAS DEVELOPED BY GROUP C-3, LOS ALAMOS C SCIENTIFIC LABORATORY, LOS ALAMOS, NM 87545. C THE LINPACK SUBROUTINES USED BY SGEFS ARE DESCRIBED IN C DETAIL IN THE *LINPACK USERS GUIDE* PUBLISHED BY C THE SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS C (SIAM) DATED 1979. C***ROUTINES CALLED R1MACH,SGECO,SGESL,XERROR C***END PROLOGUE SGEFS C INTEGER LDA,N,ITASK,IND,IWORK(*) REAL A(LDA,*),V(*),WORK(*),R1MACH REAL RCOND CHARACTER MSG*54 C***FIRST EXECUTABLE STATEMENT SGEFS IF (LDA.LT.N) GO TO 101 IF (N.LE.0) GO TO 102 IF (ITASK.LT.1) GO TO 103 IF (ITASK.GT.4) GO TO 103 IF (ITASK.EQ.2 .OR. ITASK.GT.3) GO TO 20 C C FACTOR MATRIX A INTO LU CALL SGECO(A,LDA,N,IWORK,RCOND,WORK) C C CHECK FOR COMPUTATIONALLY SINGULAR MATRIX IF (RCOND.EQ.0.0) GO TO 104 C C COMPUTE IND (ESTIMATE OF NO. OF SIGNIFICANT DIGITS) IND=-INT(ALOG10(R1MACH(4)/RCOND)) C C CHECK FOR IND GREATER THAN ZERO IF (IND.GT.0) GO TO 20 IND=-10 CALL XERROR( 'SGEFS ERROR (IND=-10) -- SOLUTION MAY HAVE NO SIGNIF 1ICANCE',58,-10,0) C C SOLVE AFTER FACTORING 20 JOB=0 IF (ITASK.GT.2) JOB=1 CALL SGESL(A,LDA,N,IWORK,V,JOB) RETURN C C IF LDA.LT.N, IND=-1, FATAL XERROR MESSAGE 101 IND=-1 WRITE(MSG, '( * ''SGEFS ERROR (IND=-1) -- LDA='', I5, '' IS LESS THAN N='', * I5 )' ) LDA, N CALL XERROR(MSG(1:54), 54, -1, 0) RETURN C C IF N.LT.1, IND=-2, FATAL XERROR MESSAGE 102 IND=-2 WRITE(MSG, '( * ''SGEFS ERROR (IND=-2) -- N='', I5, '' IS LESS THAN 1.'') ')N CALL XERROR(MSG(1:47), 47, -2, 0) RETURN C C IF ITASK.LT.1, IND=-3, FATAL XERROR MESSAGE 103 IND=-3 WRITE(MSG, '( * ''SGEFS ERROR (IND=-3) -- ITASK='', I5, '' IS LT 1 OR GT 4.'') * ') ITASK CALL XERROR(MSG(1:52), 52, -3, 0) RETURN C C IF SINGULAR MATRIX, IND=-4, FATAL XERROR MESSAGE 104 IND=-4 CALL XERROR( 'SGEFS ERROR (IND=-4) -- SINGULAR MATRIX A - NO SOLUT 1ION',55,-4,0) RETURN C END SUBROUTINE SGESL(A,LDA,N,IPVT,B,JOB) C***BEGIN PROLOGUE SGESL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS LINEAR ALGEBRA,LINPACK,MATRIX,SOLVE C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Solves the real system A*X=B or TRANS(A)*X=B C using the factors of SGECO or SGEFA C***DESCRIPTION C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C SGESL solves the real system C A * X = B or TRANS(A) * X = B C using the factors computed by SGECO or SGEFA. C C On Entry C C A REAL(LDA, N) C the output from SGECO or SGEFA. C C LDA INTEGER C the leading dimension of the array A . C C N INTEGER C the order of the matrix A . C C IPVT INTEGER(N) C the pivot vector from SGECO or SGEFA. C C B REAL(N) C the right hand side vector. C C JOB INTEGER C = 0 to solve A*X = B , C = nonzero to solve TRANS(A)*X = B where C TRANS(A) is the transpose. C C On Return C C B the solution vector X . C C Error Condition C C A division by zero will occur if the input factor contains a C zero on the diagonal. Technically, this indicates singularity, C but it is often caused by improper arguments or improper C setting of LDA . It will not occur if the subroutines are C called correctly and if SGECO has set RCOND .GT. 0.0 C or SGEFA has set INFO .EQ. 0 . C C To compute INVERSE(A) * C where C is a matrix C with P columns C CALL SGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND is too small) GO TO ... C DO 10 J = 1, P C CALL SGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. This version dated 08/14/78 . C Cleve Moler, University of New Mexico, Argonne National Lab. C C Subroutines and Functions C C BLAS SAXPY,SDOT C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED SAXPY,SDOT C***END PROLOGUE SGESL INTEGER LDA,N,IPVT(*),JOB REAL A(LDA,*),B(*) C REAL SDOT,T INTEGER K,KB,L,NM1 C***FIRST EXECUTABLE STATEMENT SGESL NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL SAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL SAXPY(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N T = SDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + SDOT(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END