SUBROUTINE DGEFS(A,LDA,N,V,ITASK,IND,WORK,IWORK,RCOND) C***BEGIN PROLOGUE DGEFS C***DATE WRITTEN 800326 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D2A1 C***AUTHOR VOORHEES, E., (LANL) C***PURPOSE DGEFS solves a GENERAL double precision C NXN system of linear equations. C***DESCRIPTION C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Subroutine DGEFS solves a general NxN system of double C precision linear equations using LINPACK subroutines DGECO C and DGESL. That is, if A is an NxN double precision matrix C and if X and B are double precision N-vectors, then DGEFS C 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.GT.1) 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 DGEFS C in this case. C C Argument Description *** C C A DOUBLE PRECISION(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 DOUBLE PRECISION(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 .GT. 1, the equation is solved using the existing C factored matrix A and IWORK. C If ITASK .LT. 1, then terminal error message IND=-3 is C 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 DOUBLE PRECISION(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 C Error Messages Printed *** C C IND=-1 terminal N is greater than LDA. C IND=-2 terminal N is less than 1. C IND=-3 terminal ITASK is less than 1. C IND=-4 terminal 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 Note- The above terminal(*fatal*) error messages are C designed to be handled by XERRWV in which C LEVEL=1 (recoverable) and IFLAG=2 . LEVEL=0 C for warning error messages from XERROR. Unless C the user provides otherwise, an error message C will be printed followed by an abort. C***REFERENCES SUBROUTINE DGEFS WAS DEVELOPED BY GROUP C-3, LOS ALAMOS C SCIENTIFIC LABORATORY, LOS ALAMOS, NM 87545. C THE LINPACK SUBROUTINES USED BY DGEFS 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 D1MACH,DGECO,DGESL,XERROR,XERRWV C***END PROLOGUE DGEFS C INTEGER LDA,N,ITASK,IND,IWORK(N) DOUBLE PRECISION A(LDA,N),V(N),WORK(N),D1MACH DOUBLE PRECISION RCOND C***FIRST EXECUTABLE STATEMENT DGEFS IF (LDA.LT.N) GO TO 101 IF (N.LE.0) GO TO 102 IF (ITASK.LT.1) GO TO 103 IF (ITASK.GT.1) GO TO 20 C C FACTOR MATRIX A INTO LU CALL DGECO(A,LDA,N,IWORK,RCOND,WORK) C C CHECK FOR COMPUTATIONALLY SINGULAR MATRIX IF (RCOND.EQ.0.0D0) GO TO 104 C C COMPUTE IND (ESTIMATE OF NO. OF SIGNIFICANT DIGITS) IND=-IDINT(DLOG10(D1MACH(4)/RCOND)) C C CHECK FOR IND GREATER THAN ZERO IF (IND.GT.0) GO TO 20 IND=-10 CALL XERROR( 'DGEFS ERROR (IND=-10) -- SOLUTION MAY HAVE NO SIGNIF 1ICANCE',58,-10,0) C C SOLVE AFTER FACTORING 20 CALL DGESL(A,LDA,N,IWORK,V,0) RETURN C C IF LDA.LT.N, IND=-1, TERMINAL XERRWV MESSAGE 101 IND=-1 CALL XERRWV( 'DGEFS ERROR (IND=-1) -- LDA=I1 IS LESS THAN N=I2', 148,-1,1,2,LDA,N,0,0,0) RETURN C C IF N.LT.1, IND=-2, TERMINAL XERRWV MESSAGE 102 IND=-2 CALL XERRWV( 'DGEFS ERROR (IND=-2) -- N=I1 IS LESS THAN 1', 143,-2,1,1,N,0,0,0,0) RETURN C C IF ITASK.LT.1, IND=-3, TERMINAL XERRWV MESSAGE 103 IND=-3 CALL XERRWV( 'DGEFS ERROR (IND=-3) -- ITASK=I1 IS LESS THAN 1', 147,-3,1,1,ITASK,0,0,0,0) RETURN C C IF SINGULAR MATRIX, IND=-4, TERMINAL XERRWV MESSAGE 104 IND=-4 CALL XERRWV( 'DGEFS ERROR (IND=-4) -- SINGULAR MATRIX A - NO SOLUT 1ION',55,-4,1,0,0,0,0,0,0) RETURN C END SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) C***BEGIN PROLOGUE DGEFA C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SGEFA-S DGEFA-D CGEFA-C), C GENERAL MATRIX,LINEAR ALGEBRA,MATRIX,MATRIX FACTORIZATION C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Factors a double precision matrix by Gaussian elimination. C***DESCRIPTION C C DGEFA factors a double precision matrix by Gaussian elimination. C C DGEFA is usually called by DGECO, but it can be called C directly with a saving in time if RCOND is not needed. C (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) . C C On Entry C C A DOUBLE PRECISION(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 DGESL or DGEDI will divide by zero C if called. Use RCOND in DGECO 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 DAXPY,DSCAL,IDAMAX C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DSCAL,IDAMAX C***END PROLOGUE DGEFA INTEGER LDA,N,IPVT(1),INFO DOUBLE PRECISION A(LDA,1) C DOUBLE PRECISION T INTEGER IDAMAX,J,K,KP1,L,NM1 C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C C***FIRST EXECUTABLE STATEMENT DGEFA 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 = IDAMAX(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.0D0) 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.0D0/A(K,K) CALL DSCAL(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 DAXPY(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.0D0) INFO = N RETURN END SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) C***BEGIN PROLOGUE DGESL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SGESL-S DGESL-D CGESL-C), C LINEAR ALGEBRA,MATRIX,SOLVE C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Solves the double precision system A*X=B or TRANS(A)*X=B C using the factors computed by DGECO or DGEFA. C***DESCRIPTION C C DGESL solves the double precision system C A * X = B or TRANS(A) * X = B C using the factors computed by DGECO or DGEFA. C C On Entry C C A DOUBLE PRECISION(LDA, N) C the output from DGECO or DGEFA. 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 DGECO or DGEFA. C C B DOUBLE PRECISION(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 DGECO has set RCOND .GT. 0.0 C or DGEFA has set INFO .EQ. 0 . C C To compute INVERSE(A) * C where C is a matrix C with P columns C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND is too small) GO TO ... C DO 10 J = 1, P C CALL DGESL(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 DAXPY,DDOT C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DDOT C***END PROLOGUE DGESL INTEGER LDA,N,IPVT(1),JOB DOUBLE PRECISION A(LDA,1),B(1) C DOUBLE PRECISION DDOT,T INTEGER K,KB,L,NM1 C***FIRST EXECUTABLE STATEMENT DGESL 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 DAXPY(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 DAXPY(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 = DDOT(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) + DDOT(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 SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z) C***BEGIN PROLOGUE DGECO C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 861211 (YYMMDD) C***CATEGORY NO. D2A1 C***KEYWORDS LIBRARY=SLATEC(LINPACK), C TYPE=DOUBLE PRECISION(SGECO-S DGECO-D CGECO-C), C CONDITION NUMBER,GENERAL MATRIX,LINEAR ALGEBRA,MATRIX, C MATRIX FACTORIZATION C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) C***PURPOSE Factors a double precision matrix by Gaussian elimination C and estimates the condition of the matrix. C***DESCRIPTION C C DGECO factors a double precision matrix by Gaussian elimination C and estimates the condition of the matrix. C C If RCOND is not needed, DGEFA is slightly faster. C To solve A*X = B , follow DGECO by DGESL. C To compute INVERSE(A)*C , follow DGECO by DGESL. C To compute DETERMINANT(A) , follow DGECO by DGEDI. C To compute INVERSE(A) , follow DGECO by DGEDI. C C On Entry C C A DOUBLE PRECISION(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 DOUBLE PRECISION 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 DOUBLE PRECISION(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 DGEFA C BLAS DAXPY,DDOT,DSCAL,DASUM C Fortran DABS,DMAX1,DSIGN C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DASUM,DAXPY,DDOT,DGEFA,DSCAL C***END PROLOGUE DGECO INTEGER LDA,N,IPVT(1) DOUBLE PRECISION A(LDA,1),Z(1) DOUBLE PRECISION RCOND C DOUBLE PRECISION DDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L C C COMPUTE 1-NORM OF A C C***FIRST EXECUTABLE STATEMENT DGECO ANORM = 0.0D0 DO 10 J = 1, N ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1)) 10 CONTINUE C C FACTOR C CALL DGEFA(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.0D0 DO 20 J = 1, N Z(J) = 0.0D0 20 CONTINUE DO 100 K = 1, N IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30 S = DABS(A(K,K))/DABS(EK-Z(K)) CALL DSCAL(N,S,Z,1) EK = S*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = DABS(WK) SM = DABS(WKM) IF (A(K,K) .EQ. 0.0D0) GO TO 40 WK = WK/A(K,K) WKM = WKM/A(K,K) GO TO 50 40 CONTINUE WK = 1.0D0 WKM = 1.0D0 50 CONTINUE KP1 = K + 1 IF (KP1 .GT. N) GO TO 90 DO 60 J = KP1, N SM = SM + DABS(Z(J)+WKM*A(K,J)) Z(J) = Z(J) + WK*A(K,J) S = S + DABS(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.0D0/DASUM(N,Z,1) CALL DSCAL(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) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 S = 1.0D0/DABS(Z(K)) CALL DSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) C YNORM = 1.0D0 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 DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 S = 1.0D0/DABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE U*Z = V C DO 160 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150 S = DABS(A(K,K))/DABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 T = -Z(K) CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END