C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/exminv.for C C====Short test driver for IMSL Matrix inverse routine DLINRG. C To run this program issue: WATFOR77 EXMINV FOR (NOWARN NOEXT C Resultant matrix will display on the screen. C INTEGER I, J, LDA, LDAINV, N, IER, IJOB DOUBLE PRECISION A(500,500), AORI(500,500) DOUBLE PRECISION MIN, MAX, AVE, SUM, D1, D2 DOUBLE PRECISION B(500), WKAREA(1000) DATA LDA/500/, LDAINV/500/, N/3/ C---Sample Matrix whose output should be: CA-Inverse= (IBM 370 mainframe): C 12.448985208229890 -4.714287271305429 -6.734697936924458 C -1.428571817826357 0.000000000000000 1.428571817826357 C -2.040817438687704 1.428571817826357 0.612245620861346 CMIN,AVE,MAX= 0.000000000000000 4.625929269271485D-18 4.163336342344337 D-17 C CA-Inverse= (PC Salford Fortran): C 12.4489781877 -4.71428532503 -6.73469286271 C -1.42857133126 0.000000000000E+0000 1.42857133126 C -2.04081604849 1.42857133126 0.612244717234 CMIN,AVE,MAX= 0.000000000000E+0000 1.056253849817E-16 8.881784197001E-16 C DO I=1,N DO J=1,N A(I,J)=1.0 END DO END DO A(1,2)=3.3 A(2,2)=3.0 A(3,2)=4.0 A(1,3)=3.3 A(2,3)=4.0 A(3,3)=3.3 IJOB=1 D1=1 C---Make a copy of A DO I=1,N DO J=1,N AORI(I,J)=A(I,J) END DO END DO C---Invert A. CALL LINV3F (A,B,IJOB,N,LDA,D1,D2,WKAREA,IER) C---Display the results. IF(IER.NE.0) THEN WRITE(6,*) 'Conditioning Problem; IER = ',IER ENDIF WRITE(6,*) 'A-Inverse=' DO I=1,N WRITE(6,*) (A(I,J),J=1,N) END DO C---Check results by computing the MIN, MAX, and AVERAGE Norms C of the residual matrix: AINV*A-I MIN=1.D75 MAX=-1.D-75 AVE=0. DO I=1,N DO J=1,N SUM=0. DO K=1,N SUM=SUM+A(I,K)*AORI(K,J) END DO IF (I.EQ.J) SUM=SUM-1.D0 SUM=ABS(SUM) IF(MIN.GT.SUM) MIN=SUM IF(MAX.LT.SUM) MAX=SUM AVE=AVE+SUM END DO END DO AVE=AVE/N**2 WRITE(*,*) 'MIN,AVE,MAX=',MIN,AVE,MAX C PAUSE STOP END C IMSL ROUTINE NAME - LINV3F C C----------------------------------------------------------------------- C C COMPUTER - IBM/DOUBLE C C LATEST REVISION - JUNE 1, 1982 C C PURPOSE - IN PLACE INVERSE, EQUATION SOLUTION, AND/OR C DETERMINANT EVALUATION - FULL STORAGE MODE C C USAGE - CALL LINV3F (A,B,IJOB,N,IA,D1,D2,WKAREA,IER) C C ARGUMENTS A - INPUT/OUTPUT MATRIX OF DIMENSION N BY N. SEE C PARAMETER IJOB. C B - INPUT/OUTPUT VECTOR OF LENGTH N WHEN IJOB = C 2 OR 3. OTHERWISE, B IS NOT USED. C ON INPUT, B CONTAINS THE RIGHT HAND SIDE OF C OF THE EQUATION AX = B. C ON OUTPUT, THE SOLUTION X REPLACES B. C IJOB - INPUT OPTION PARAMETER. IJOB = I IMPLIES: C I = 1, INVERT MATRIX A. A IS REPLACED BY C ITS INVERSE. C I = 2, SOLVE THE EQUATION AX = B. A IS C REPLACED BY THE LU DECOMPOSITION OF A C ROWWISE PERMUTATION OF A, WHERE U IS C UPPER TRIANGULAR AND L IS LOWER C TRIANGULAR WITH UNIT DIAGONAL. C THE UNIT DIAGONAL OF L IS NOT STORED. C I = 3, SOLVE AX = B AND INVERT MATRIX A. C A IS REPLACED BY ITS INVERSE. C I = 4, COMPUTE THE DETERMINANT OF A. C A IS REPLACED BY THE LU DECOMPOSITION C OF A ROWWISE PERMUTATION OF A. C N - ORDER OF A. (INPUT) C IA - ROW DIMENSION OF MATRIX A EXACTLY AS C SPECIFIED IN THE DIMENSION STATEMENT IN THE C CALLING PROGRAM. (INPUT) C D1 - INPUT/OUTPUT. IF THE D1 AND D2 COMPONENTS C D2 OF DETERMINANT(A) = D1*2**D2 ARE DESIRED, C INPUT D1.GE.0. OTHERWISE, INPUT D1.LT.0. C D2 IS NEVER INPUT. C WKAREA - WORK AREA OF LENGTH AT LEAST 2*N FOR IJOB=1 C OR IJOB=3. C WORK AREA OF LENGTH AT LEAST N FOR IJOB=2 C OR IJOB=4. C IER - ERROR PARAMETER. (OUTPUT) C WARNING WITH FIX C IER = 65 INDICATES THAT IJOB WAS LESS THAN C 1 OR GREATER THAN 4. IJOB IS ASSUMED TO C BE 4. C TERMINAL ERROR C IER = 130 INDICATES THAT MATRIX A IS C ALGORITHMICALLY SINGULAR. (SEE THE C CHAPTER L PRELUDE). C C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 C - SINGLE/H36,H48,H60 C C REQD. IMSL ROUTINES - LUDATN,LUELMN,UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE LINV3F (A,B,IJOB,N,IA,D1,D2,WKAREA,IER) C DOUBLE PRECISION A(IA,N),B(N),WKAREA(2*N),C1,C2,D1,D2,WA,ZERO, * ONE,SUM,C DATA ZERO/0.0D0/,ONE/1.0D0/ C FIRST EXECUTABLE STATEMENT C LU DECOMPOSITION OF A CALL LUDATN (A,IA,N,A,IA,0,C1,C2,WKAREA,WKAREA,WA,IER) IF (D1 .LT. ZERO .AND. IJOB .GE. 1 .AND. IJOB .LT. 4) GO TO 5 D1 = C1 D2 = C2 5 IF (IER .GE. 128) GO TO 60 IF (IJOB .LE. 0 .OR. IJOB .GT. 4) GO TO 55 C SOLVE AX = B IF (IJOB .EQ. 2 .OR . IJOB .EQ. 3) CALL LUELMN (A,IA,N,B,WKAREA,B) IF (IJOB .NE. 1 .AND. IJOB .NE. 3) GO TO 9005 C MATRIX INVERSION A(N,N) = ONE/A(N,N) NM1 = N-1 IF (NM1 .LT. 1) GO TO 9005 DO 40 II=1,NM1 L = N-II M = L+1 DO 15 I=M,N SUM = ZERO DO 10 K=M,N SUM = SUM-A(I,K)*A(K,L) 10 CONTINUE WKAREA(N+I) = SUM 15 CONTINUE DO 20 I=M,N A(I,L) = WKAREA(N+I) 20 CONTINUE DO 30 J=L,N SUM = ZERO IF (J .EQ. L) SUM = ONE DO 25 K=M,N SUM = SUM-A(L,K)*A(K,J) 25 CONTINUE WKAREA(N+J) = SUM/A(L,L) 30 CONTINUE DO 35 J=L,N A(L,J) = WKAREA(N+J) 35 CONTINUE 40 CONTINUE C PERMUTE COLUMNS OF A INVERSE DO 50 I=1,N J = N-I+1 JP = WKAREA(J) IF (J .EQ. JP) GO TO 50 DO 45 K=1,N C = A(K,JP) A(K,JP) = A(K,J) A(K,J) = C 45 CONTINUE 50 CONTINUE GO TO 9005 55 CONTINUE C WARNING WITH FIX - IJOB WAS SET C INCORRECTLY IER = 65 GO TO 9000 C TERMINAL ERROR - MATRIX A IS C ALGORITHMICALLY SINGULAR 60 IER = 130 9000 CONTINUE CALL UERTST(IER,6HLINV3F) 9005 RETURN END C IMSL ROUTINE NAME - LUDATN C C----------------------------------------------------------------------- C C COMPUTER - IBM/DOUBLE C C LATEST REVISION - JUNE 1, 1982 C C PURPOSE - NUCLEUS CALLED ONLY BY IMSL SUBROUTINE LEQT2F C C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 C - SINGLE/H36,H48,H60 C C REQD. IMSL ROUTINES - UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE LUDATN (A,IA,N,LU,ILU,IDGT,D1,D2,APVT,EQUIL,WA,IER) C DIMENSION A(IA,N),LU(ILU,N),APVT(N),EQUIL(N) DOUBLE PRECISION A,LU,D1,D2,EQUIL,WA,ZERO,ONE,FOUR,SIXTN,SIXTH, * RN,WREL,BIGA,BIG,P,SUM,AI,WI,T,TEST,Q,APVT DATA ZERO,ONE,FOUR,SIXTN,SIXTH/0.D0,1.D0,4.D0, * 16.D0,.0625D0/ C FIRST EXECUTABLE STATEMENT C INITIALIZATION IER = 0 RN = N WREL = ZERO D1 = ONE D2 = ZERO BIGA = ZERO DO 10 I=1,N BIG = ZERO DO 5 J=1,N P = A(I,J) LU(I,J) = P P = DABS(P) IF (P .GT. BIG) BIG = P 5 CONTINUE IF (BIG .GT. BIGA) BIGA = BIG IF (BIG .EQ. ZERO) GO TO 110 EQUIL(I) = ONE/BIG 10 CONTINUE DO 105 J=1,N JM1 = J-1 IF (JM1 .LT. 1) GO TO 40 C COMPUTE U(I,J), I=1,...,J-1 DO 35 I=1,JM1 SUM = LU(I,J) IM1 = I-1 IF (IDGT .EQ. 0) GO TO 25 C WITH ACCURACY TEST AI = DABS(SUM) WI = ZERO IF (IM1 .LT. 1) GO TO 20 DO 15 K=1,IM1 T = LU(I,K)*LU(K,J) SUM = SUM-T WI = WI+DABS(T) 15 CONTINUE LU(I,J) = SUM 20 WI = WI+DABS(SUM) IF (AI .EQ. ZERO) AI = BIGA TEST = WI/AI IF (TEST .GT. WREL) WREL = TEST GO TO 35 C WITHOUT ACCURACY 25 IF (IM1 .LT. 1) GO TO 35 DO 30 K=1,IM1 SUM = SUM-LU(I,K)*LU(K,J) 30 CONTINUE LU(I,J) = SUM 35 CONTINUE 40 P = ZERO C COMPUTE U(J,J) AND L(I,J), I=J+1,..., DO 70 I=J,N SUM = LU(I,J) IF (IDGT .EQ. 0) GO TO 55 C WITH ACCURACY TEST AI = DABS(SUM) WI = ZERO IF (JM1 .LT. 1) GO TO 50 DO 45 K=1,JM1 T = LU(I,K)*LU(K,J) SUM = SUM-T WI = WI+DABS(T) 45 CONTINUE LU(I,J) = SUM 50 WI = WI+DABS(SUM) IF (AI .EQ. ZERO) AI = BIGA TEST = WI/AI IF (TEST .GT. WREL) WREL = TEST GO TO 65 C WITHOUT ACCURACY TEST 55 IF (JM1 .LT. 1) GO TO 65 DO 60 K=1,JM1 SUM = SUM-LU(I,K)*LU(K,J) 60 CONTINUE LU(I,J) = SUM 65 Q = EQUIL(I)*DABS(SUM) IF (P .GE. Q) GO TO 70 P = Q IMAX = I 70 CONTINUE C TEST FOR ALGORITHMIC SINGULARITY IF (RN+P .EQ. RN) GO TO 110 IF (J .EQ. IMAX) GO TO 80 C INTERCHANGE ROWS J AND IMAX D1 = -D1 DO 75 K=1,N P = LU(IMAX,K) LU(IMAX,K) = LU(J,K) LU(J,K) = P 75 CONTINUE EQUIL(IMAX) = EQUIL(J) 80 APVT(J) = IMAX D1 = D1*LU(J,J) 85 IF (DABS(D1) .LE. ONE) GO TO 90 D1 = D1*SIXTH D2 = D2+FOUR GO TO 85 90 IF (DABS(D1) .GE. SIXTH) GO TO 95 D1 = D1*SIXTN D2 = D2-FOUR GO TO 90 95 CONTINUE JP1 = J+1 IF (JP1 .GT. N) GO TO 105 C DIVIDE BY PIVOT ELEMENT U(J,J) P = LU(J,J) DO 100 I=JP1,N LU(I,J) = LU(I,J)/P 100 CONTINUE 105 CONTINUE C PERFORM ACCURACY TEST IF (IDGT .EQ. 0) GO TO 9005 P = 3*N+3 WA = P*WREL IF (WA+10.D0**(-IDGT) .NE. WA) GO TO 9005 IER = 34 GO TO 9000 C ALGORITHMIC SINGULARITY 110 IER = 129 D1 = ZERO D2 = ZERO 9000 CONTINUE C PRINT ERROR CALL UERTST(IER,6HLUDATN) 9005 RETURN END C IMSL ROUTINE NAME - LUELMN C C----------------------------------------------------------------------- C C COMPUTER - IBM/DOUBLE C C LATEST REVISION - JUNE 1, 1982 C C PURPOSE - NUCLEUS CALLED ONLY BY IMSL SUBROUTINE LEQT2F C C REQD. IMSL ROUTINES - NONE REQUIRED C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE LUELMN (A,IA,N,B,APVT,X) C DIMENSION A(IA,N),B(N),APVT(N),X(N) DOUBLE PRECISION A,B,X,SUM,APVT C FIRST EXECUTABLE STATEMENT C SOLVE LY = B FOR Y DO 5 I=1,N 5 X(I) = B(I) IW = 0 DO 20 I=1,N IP = APVT(I) SUM = X(IP) X(IP) = X(I) IF (IW .EQ. 0) GO TO 15 IM1 = I-1 DO 10 J=IW,IM1 SUM = SUM-A(I,J)*X(J) 10 CONTINUE GO TO 20 15 IF (SUM .NE. 0.D0) IW = I 20 X(I) = SUM C SOLVE UX = Y FOR X DO 30 IB=1,N I = N+1-IB IP1 = I+1 SUM = X(I) IF (IP1 .GT. N) GO TO 30 DO 25 J=IP1,N SUM = SUM-A(I,J)*X(J) 25 CONTINUE 30 X(I) = SUM/A(I,I) RETURN END C to enable Watfor77 to run callers of these. C********************************************************************** C ROUTINE: UERTST FORTRAN C PURPOSE: Replacement for IMSL version 9.2 UERTST and UGETIO C to enable Watfor77 to run callers of these. C PACKAGE: IMSL 9.2 C CALLS: none C CALLED BY: IMSL Version 2.9 C AUTHOR: H. D. Knoble hdkLESS at SPAM psu dot edu 10/22/91 C REVISIONS: C********************************************************************** SUBROUTINE UERSET (LEVEL,LEVOLD) C To use this abbreviated version of UERTST in Watfor77, after your C last Fortran END statement add the Watfor77 compiler directive C (beginning in column 1): C$INCLUDE UERTST FORTRAN * C C To use this abbreviated version of UERTST in VS Fortran, after your C last Fortran END statement add the VS Fortran compiler directive C (beginning in column 7): INCLUDE 'UERTST FORTRAN *' C SPECIFICATIONS FOR ARGUMENTS INTEGER LEVEL,LEVOLD C FIRST EXECUTABLE STATEMENT LEVOLD = LEVEL CALL UERTST (LEVOLD,6HUERSET) RETURN END SUBROUTINE UERTST (IER,NAME) C==== Dummy version to debug in WATFOR77 with. INTEGER IER CHARACTER*6 NAME WRITE(6,*) '$$UERTST NAME=',NAME,', IER=',IER RETURN END