SUBROUTINE LINEQ (NA,N,A,B,COND,IPVT,WORK) C C *****PARAMETERS: INTEGER NA,N,IPVT(N) DOUBLE PRECISION A(NA,N),B(N),COND,WORK(N) C C *****SUBROUTINES CALLED: C DGECOM,DGESLM C C ------------------------------------------------------------------ C C *****PURPOSE: C THIS SUBROUTINE SOLVES THE SYSTEM OF LINEAR EQUATIONS C A*X = B C WHERE B IS AN N-VECTOR. SUBROUTINES DGECOM (LU-FACTORIZATION) C AND DGESLM (FORWARD ELIMINATION AND BACK SUBSTITUTION) ARE C EMPLOYED. AN ESTIMATE OF THE CONDITION OF A IS RETURNED. C SHOULD A BE SINGULAR TO WORKING ACCURACY, COND IS SET TO 1.0D+20. C C *****PARAMETER DESCRIPTION: C C ON INPUT: C C NA ROW DIMENSION OF THE ARRAY CONTAINING A AS C DECLARED IN THE CALLING PROGRAM DIMENSION C STATEMENT; C C N ORDER OF THE MATRIX A; C C A N X N COEFFICIENT MATRIX; C C B RIGHT HAND SIDE VECTOR OF LENGTH N; C C ON OUTPUT: C C B SOLUTION VECTOR X = (A-INVERSE)*B; C C COND AN ESTIMATE OF THE CONDITION OF A. C =1/RCOND WHERE RCOND IS THE INVERSE C OF THE CONDITION ESTIMATE (SEE THE LINPACK C USER'S GUIDE FOR DETAILS); C C IPVT PIVOT VECTOR OF LENGTH N (SEE DGECOM C DOCUMENTATION); C C WORK A REAL SCRATCH VECTOR OF LENGTH N. C C *****APPLICATIONS AND USAGE RESTRICTIONS: C THE VALUE OF COND SHOULD ALWAYS BE CHECKED BY THE CALLING C PROGRAM. SHOULD A BE NEAR-SINGULAR (OR SINGULAR TO WORKING C ACCURACY) THE DATA SHOULD BE INVESTIGATED FOR POSSIBLE C ERRORS. IF THERE ARE NONE AND THE PROBLEM IS APPARENTLY C WELL-POSED AND/OR MEANINGFUL, SINGULAR VALUE ANALYSIS MAY C THEN BE A MORE RELIABLE SOLUTION TECHNIQUE (E.G., EISPACK C SUBROUTINES SVD OR MINFIT). C C *****ALGORITHM NOTES: C THE CONTENTS OF A ARE MODIFIED BY THIS SUBROUTINE. SHOULD THE C ORIGINAL COEFFICIENTS OF A BE NEEDED SUBSEQUENTLY, THE C CONTENTS OF A SHOULD BE SAVED PRIOR TO THE CALL TO LINEQ. C C *****HISTORY: C WRITTEN BY ALAN J. LAUB (DEP'T. OF ELEC. ENGRG. - SYSTEMS, C UNIVERSITY OF SOUTHERN CALIFORNIA, LOS ANGELES, CA 90007, C PH.: (213)-743-5535), MAY 1980. C MOST RECENT VERSION: MAY 6, 1980. C C ------------------------------------------------------------------ C CALL DGECOM (A,NA,N,IPVT,COND,WORK) IF ((1.0D0 + COND) .GT. 1.0D0) GO TO 20 COND = 1.0D+20 RETURN 20 CONTINUE COND = 1.0D0/COND CALL DGESLM (A,NA,N,IPVT,B) RETURN END