C IMSL ROUTINE NAME - LINV2F C C----------------------------------------------------------------------- C C COMPUTER - IBM77/DOUBLE C C LATEST REVISION - JUNE 1, 1982 C C PURPOSE - INVERSION OF A MATRIX - FULL STORAGE MODE - C HIGH ACCURACY SOLUTION C C USAGE - CALL LINV2F (A,N,IA,AINV,IDGT,WKAREA,IER) C C ARGUMENTS A - INPUT MATRIX OF DIMENSION N BY N CONTAINING C THE MATRIX TO BE INVERTED. C N - ORDER OF A. (INPUT) C IA - ROW DIMENSION OF MATRICES A AND AINV EXACTLY C AS SPECIFIED IN THE DIMENSION STATEMENT IN C THE CALLING PROGRAM. (INPUT) C AINV - OUTPUT MATRIX OF DIMENSION N BY N CONTAINING C THE INVERSE OF A. A AND AINV MUST OCCUPY C SEPARATE CORE LOCATIONS. C IDGT - INPUT OPTION. C IF IDGT IS GREATER THAN 0, THE ELEMENTS OF A C ARE ASSUMED TO BE CORRECT TO IDGT DECIMAL C DIGITS AND THE ROUTINE PERFORMS AN ACCURACY C TEST. C IF IDGT EQUALS 0, THE ACCURACY TEST IS C BYPASSED. C ON OUTPUT, IDGT CONTAINS THE APPROXIMATE C NUMBER OF DIGITS IN THE ANSWER WHICH C WERE UNCHANGED AFTER IMPROVEMENT. C WKAREA - WORK AREA OF DIMENSION GREATER THAN OR EQUAL C TO N**2+3N. C IER - ERROR PARAMETER. (OUTPUT) C TERMINAL ERROR C IER=129 INDICATES THAT THE MATRIX IS C ALGORITHMICALLY SINGULAR. (SEE THE C CHAPTER L PRELUDE). C IER=131 INDICATES THAT THE MATRIX IS TOO C ILL-CONDITIONED FOR ITERATIVE IMPROVEMENT C TO BE EFFECTIVE. C WARNING ERROR C IER=34 INDICATES THAT THE ACCURACY TEST C FAILED. THE COMPUTED SOLUTION MAY BE IN C ERROR BY MORE THAN CAN BE ACCOUNTED FOR C BY THE UNCERTAINTY OF THE DATA. THIS C WARNING CAN BE PRODUCED ONLY IF IDGT IS C GREATER THAN 0 ON INPUT. SEE CHAPTER L C PRELUDE FOR FURTHER DISCUSSION. C C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 C - SINGLE/H36,H48,H60 C C REQD. IMSL ROUTINES - SINGLE/LEQT2F,LUDATN,LUELMN,LUREFN,UERTST, C UGETIO C - DOUBLE/LEQT2F,LUDATN,LUELMN,LUREFN,UERTST, C UGETIO,VXADD,VXMUL,VXSTO 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 LINV2F (A,N,IA,AINV,IDGT,WKAREA,IER) C DOUBLE PRECISION A(IA,N),AINV(IA,N),WKAREA(1),ZERO,ONE DATA ONE/1.0D0/,ZERO/0.0D0/ C FIRST EXECUTABLE STATEMENT C INITIALIZE IER IER=0 C SET AINV TO THE N X N C IDENTITY MATRIX DO 10 I = 1,N DO 5 J = 1,N AINV(I,J) = ZERO 5 CONTINUE AINV(I,I) = ONE 10 CONTINUE C COMPUTE THE INVERSE OF A CALL LEQT2F (A,N,N,IA,AINV,IDGT,WKAREA,IER) IF (IER.EQ.0) GO TO 9005 9000 CONTINUE CALL UERTST (IER,'LINV2F') 9005 RETURN END C IMSL ROUTINE NAME - LEQT2F C C----------------------------------------------------------------------- C C COMPUTER - IBM77/DOUBLE C C LATEST REVISION - JUNE 1, 1982 C C PURPOSE - LINEAR EQUATION SOLUTION - FULL STORAGE C MODE - HIGH ACCURACY SOLUTION C C USAGE - CALL LEQT2F (A,M,N,IA,B,IDGT,WKAREA,IER) C C ARGUMENTS A - INPUT MATRIX OF DIMENSION N BY N CONTAINING C THE COEFFICIENT MATRIX OF THE EQUATION C AX = B. C M - NUMBER OF RIGHT-HAND SIDES. (INPUT) C N - ORDER OF A AND NUMBER OF ROWS IN B. (INPUT) C IA - ROW DIMENSION OF A AND B EXACTLY AS SPECIFIED C IN THE DIMENSION STATEMENT IN THE CALLING C PROGRAM. (INPUT) C B - INPUT MATRIX OF DIMENSION N BY M CONTAINING C THE RIGHT-HAND SIDES OF THE EQUATION AX = B. C ON OUTPUT, THE N BY M MATRIX OF SOLUTIONS C REPLACES B. C IDGT - INPUT OPTION. C IF IDGT IS GREATER THAN 0, THE ELEMENTS OF C A AND B ARE ASSUMED TO BE CORRECT TO IDGT C DECIMAL DIGITS AND THE ROUTINE PERFORMS C AN ACCURACY TEST. C IF IDGT EQUALS 0, THE ACCURACY TEST IS C BYPASSED. C ON OUTPUT, IDGT CONTAINS THE APPROXIMATE C NUMBER OF DIGITS IN THE ANSWER WHICH C WERE UNCHANGED AFTER IMPROVEMENT. C WKAREA - WORK AREA OF DIMENSION GREATER THAN OR EQUAL C TO N**2+3N. C IER - ERROR PARAMETER. (OUTPUT) C WARNING ERROR C IER = 34 INDICATES THAT THE ACCURACY TEST C FAILED. THE COMPUTED SOLUTION MAY BE IN C ERROR BY MORE THAN CAN BE ACCOUNTED FOR C BY THE UNCERTAINTY OF THE DATA. THIS C WARNING CAN BE PRODUCED ONLY IF IDGT IS C GREATER THAN 0 ON INPUT. (SEE THE C CHAPTER L PRELUDE FOR FURTHER DISCUSSION.) C TERMINAL ERROR C IER = 129 INDICATES THAT THE MATRIX IS C ALGORITHMICALLY SINGULAR. (SEE THE C CHAPTER L PRELUDE). C IER = 131 INDICATES THAT THE MATRIX IS TOO C ILL-CONDITIONED FOR ITERATIVE IMPROVEMENT C TO BE EFFECTIVE. C C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 C - SINGLE/H36,H48,H60 C C REQD. IMSL ROUTINES - SINGLE/LUDATN,LUELMN,LUREFN,UERTST,UGETIO C - DOUBLE/LUDATN,LUELMN,LUREFN,UERTST,UGETIO, C VXADD,VXMUL,VXSTO 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 LEQT2F (A,M,N,IA,B,IDGT,WKAREA,IER) C DIMENSION A(IA,1),B(IA,1),WKAREA(1) DOUBLE PRECISION A,B,WKAREA,D1,D2,WA C FIRST EXECUTABLE STATEMENT C INITIALIZE IER IER=0 JER=0 J = N*N+1 K = J+N MM = K+N KK = 0 MM1 = MM-1 JJ=1 DO 5 L=1,N DO 5 I=1,N WKAREA(JJ)=A(I,L) JJ=JJ+1 5 CONTINUE C DECOMPOSE A CALL LUDATN (WKAREA,N,N,A,IA,IDGT,D1,D2,WKAREA(J),WKAREA(K), * WA,IER) IF (IER.GT.128) GO TO 25 IF (IDGT .EQ. 0 .OR. IER .NE. 0) KK = 1 DO 15 I = 1,M C PERFORMS THE ELIMINATION PART OF C AX = B CALL LUELMN (A,IA,N,B(1,I),WKAREA(J),WKAREA(MM)) C REFINEMENT OF SOLUTION TO AX = B IF (KK .NE. 0) * CALL LUREFN (WKAREA,N,N,A,IA,B(1,I),IDGT,WKAREA(J),WKAREA(MM), * WKAREA(K),WKAREA(K),JER) DO 10 II=1,N B(II,I) = WKAREA(MM1+II) 10 CONTINUE IF (JER.NE.0) GO TO 20 15 CONTINUE GO TO 25 20 IER = 131 25 JJ=1 DO 30 J = 1,N DO 30 I = 1,N A(I,J)=WKAREA(JJ) JJ=JJ+1 30 CONTINUE IF (IER .EQ. 0) GO TO 9005 9000 CONTINUE CALL UERTST (IER,'LEQT2F') 9005 RETURN END C IMSL ROUTINE NAME - LUDATN C C----------------------------------------------------------------------- C C COMPUTER - IBM77/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,1),LU(ILU,1),APVT(1),EQUIL(1) 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 Q = RN+P 71 IF (Q .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 Q = WA+10.D0**(-IDGT) 106 IF (Q .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,'LUDATN') 9005 RETURN END C IMSL ROUTINE NAME - LUELMN C C----------------------------------------------------------------------- C C COMPUTER - IBM77/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,1),B(1),APVT(1),X(1) 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 IMSL ROUTINE NAME - LUREFN C C----------------------------------------------------------------------- C C COMPUTER - IBM77/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 - SINGLE/LUELMN,UERTST,UGETIO C - DOUBLE/LUELMN,UERTST,UGETIO,VXADD,VXMUL, C VXSTO 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 LUREFN (A,IA,N,UL,IUL,B,IDGT,APVT,X,RES,DX,IER) C DIMENSION A(IA,1),UL(IUL,1),B(1),X(1),RES(1),DX(1) DIMENSION APVT(1) DOUBLE PRECISION DNORM DIMENSION ACCXT(2) DOUBLE PRECISION A,ACCXT,B,UL,X,RES,DX,ZERO,XNORM,DXNORM,APVT DATA ITMAX/75/,ZERO/0.D0/ C FIRST EXECUTABLE STATEMENT IER=0 XNORM = ZERO DO 10 I=1,N XNORM = DMAX1(XNORM,DABS(X(I))) 10 CONTINUE IF (XNORM .NE. ZERO) GO TO 20 IDGT = 50 GO TO 9005 20 DO 45 ITER=1,ITMAX DO 30 I=1,N ACCXT(1) = 0.0D0 ACCXT(2) = 0.0D0 CALL VXADD(B(I),ACCXT) DO 25 J=1,N CALL VXMUL(-A(I,J),X(J),ACCXT) 25 CONTINUE CALL VXSTO(ACCXT,RES(I)) 30 CONTINUE CALL LUELMN (UL,IUL,N,RES,APVT,DX) DXNORM = ZERO XNORM = ZERO DO 35 I=1,N X(I) = X(I) + DX(I) DXNORM = DMAX1(DXNORM,DABS(DX(I))) XNORM = DMAX1(XNORM,DABS(X(I))) 35 CONTINUE IF (ITER .NE. 1) GO TO 40 IDGT = 50 IF (DXNORM .NE. ZERO) IDGT = -DLOG10(DXNORM/XNORM) 40 DNORM = XNORM+DXNORM 41 IF (DNORM .EQ. XNORM) GO TO 9005 45 CONTINUE C ITERATION DID NOT CONVERGE IER = 129 9000 CONTINUE CALL UERTST(IER,'LUREFN') 9005 RETURN END C IMSL ROUTINE NAME - VXADD C C----------------------------------------------------------------------- C C COMPUTER - IBM77/DOUBLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - EXTENDED PRECISION ADD C C USAGE - CALL VXADD (A,ACC) C C ARGUMENTS A - DOUBLE PRECISION NUMBER TO BE ADDED TO THE C ACCUMULATOR. (INPUT) C ACC - ACCUMULATOR. (INPUT AND OUTPUT) C ACC IS A DOUBLE PRECISION VECTOR OF LENGTH C 2. ON OUTPUT, ACC CONTAINS THE SUM OF C INPUT ACC AND A. C C PRECISION/HARDWARE - DOUBLE/H32 C - NOT AVAILABLE/H36,H48,H60 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 REMARKS VXADD ADDS THE DOUBLE PRECISION NUMBER A TO THE C EXTENDED PRECISION ACCUMULATOR, ACC. THE SUBROUTINE C ASSUMES THAT AN EXTENDED PRECISION NUMBER IS ALREADY IN C THE ACCUMULATOR. THEREFORE, BEFORE THE FIRST CALL TO C VXADD, ACC(1) AND ACC(2) MUST BE SET TO ZERO. 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 VXADD(A,ACC) C C SPECIFICATIONS FOR ARGUMENTS DOUBLE PRECISION A,ACC(2) C SPECIFICATIONS FOR LOCAL VARIABLES DOUBLE PRECISION X,Y,Z,ZZ C FIRST EXECUTABLE STATEMENT X = ACC(1) Y = A IF (DABS(ACC(1)).GE.DABS(A)) GO TO 1 X = A Y = ACC(1) C COMPUTE Z+ZZ = ACC(1)+A EXACTLY 1 Z = X+Y ZZ = (X-Z)+Y C COMPUTE ZZ+ACC(2) USING DOUBLE C PRECISION ARITHMETIC ZZ = ZZ+ACC(2) C COMPUTE ACC(1)+ACC(2) = Z+ZZ EXACTLY ACC(1) = Z+ZZ ACC(2) = (Z-ACC(1))+ZZ RETURN END C IMSL ROUTINE NAME - VXMUL C C----------------------------------------------------------------------- C C COMPUTER - IBM77/DOUBLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - EXTENDED PRECISION MULTIPLY C C USAGE - CALL VXMUL (A,B,ACC) C C ARGUMENTS A - INPUT DOUBLE PRECISION NUMBER C B - INPUT DOUBLE PRECISION NUMBER C ACC - ACCUMULATOR. (INPUT AND OUTPUT) C ACC IS A DOUBLE PRECISION VECTOR OF LENGTH C 2. ON OUTPUT, ACC CONTAINS THE SUM OF C INPUT ACC AND A*B. C C PRECISION/HARDWARE - DOUBLE/H32 C - NOT AVAILABLE/H36,H48,H60 C C REQD. IMSL ROUTINES - VXADD C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C REMARKS VXMUL ADDS THE PRODUCT A*B TO THE EXTENDED PRECISION C ACCUMULATOR, ACC. THE SUBROUTINE ASSUMES THAT AN C EXTENDED PRECISION NUMBER IS ALREADY IN THE C ACCUMULATOR. THEREFORE, BEFORE THE FIRST CALL TO C VXMUL, ACC(1) AND ACC(2) MUST BE SET TO ZERO. 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 VXMUL (A,B,ACC) C SPECIFICATIONS FOR ARGUMENTS DOUBLE PRECISION A,B,ACC(2) C SPECIFICATIONS FOR LOCAL VARIABLES DOUBLE PRECISION X,HA,TA,HB,TB INTEGER IX(2),I LOGICAL*1 LX(8),LI(4) EQUIVALENCE (X,LX(1),IX(1)),(I,LI(1)) DATA I/0/ C SPLIT A = HA+TA C B = HB+TB C FIRST EXECUTABLE STATEMENT X = A LI(4) = LX(5) IX(2) = 0 I = (I/16)*16 LX(5) = LI(4) HA=X TA=A-HA X = B LI(4) = LX(5) IX(2) = 0 I = (I/16)*16 LX(5) = LI(4) HB = X TB = B-HB C COMPUTE HA*HB,HA*TB,TA*HB, AND TA*TB C AND CALL VXADD TO ACCUMULATE THE C SUM X = TA*TB CALL VXADD(X,ACC) X = HA*TB CALL VXADD(X,ACC) X = TA*HB CALL VXADD(X,ACC) X = HA*HB CALL VXADD(X,ACC) RETURN END C IMSL ROUTINE NAME - VXSTO C C----------------------------------------------------------------------- C C COMPUTER - IBM77/DOUBLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - DOUBLE PRECISION STORE. C C USAGE - CALL VXSTO(ACC,D) C C ARGUMENTS ACC - ACCUMULATOR. (INPUT) C ACC IS A DOUBLE PRECISION VECTOR OF LENGTH C 2. ACC IS ASSUMED TO BE THE RESULT OF C CALLING VXADD OR VXMUL TO PERFORM EXTENDED C PRECISION OPERATIONS. C D - DOUBLE PRECISION SCALAR. (OUTPUT) C ON OUTPUT, D CONTAINS A DOUBLE PRECISION C APPROXIMATION TO THE VALUE OF THE EXTENDED C PRECISION ACCUMULATOR. C C PRECISION/HARDWARE - DOUBLE/H32 C - NOT AVAILABLE/H36,H48,H60 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 - 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 VXSTO (ACC,D) C SPECIFICATIONS FOR ARGUMENTS DOUBLE PRECISION ACC(2),D C FIRST EXECUTABLE STATEMENT D = ACC(1)+ACC(2) RETURN END *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CINV(A,B,C,D,N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NP=50) DOUBLE PRECISION A(NP,NP), B(NP,NP), C(NP,NP), D(NP,NP), $ CJ(NP), DJ(NP) DO 10 I=1,N DO 15 J=1,N C(I,J)=0.0D0 D(I,J)=0.0D0 15 CONTINUE C(I,I)=1.0D0 10 CONTINUE CALL CLUDCMP(A,B,N,NP,INDX,DD) C DO 11 I=1,N C DO 11 J=1,N C WRITE(9,*) A(I,J) C 11 CONTINUE C DO 12 I=1,N C DO 12 J=1,N C WRITE(9,*) B(I,J) C 12 CONTINUE DO 20 J=1,N DO 25 I=1,N CJ(I)=C(I,J) DJ(I)=D(I,J) 25 CONTINUE CALL CLUBKSB(A,B,N,NP,INDX,CJ,DJ) DO 30 I=1,N C(I,J)=CJ(I) D(I,J)=DJ(I) 30 CONTINUE 20 CONTINUE RETURN END *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CLUDCMP(A,B,N,NP,INDX,DD) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NMAX=100,TINY=1.0E-10) DOUBLE PRECISION A(NP,NP), B(NP,NP), INDX(N), VV(NMAX) DD=1. DO 10 I=1,N AAMAX=0. DO 20 J=1,N ABSA=SQRT(A(I,J)**2+B(I,J)**2) IF (ABSA.GT.AAMAX) AAMAX=ABSA 20 CONTINUE IF (AAMAX.EQ.0.) WRITE(9,*) 'SINGULAR MATRIX' VV(I)=1./AAMAX 10 CONTINUE DO 30 J=1,N DO 40 I=1,J-1 SUM1=A(I,J) SUM2=B(I,J) DO 50 K=1,I-1 SUM1=SUM1-(A(I,K)*A(K,J)-B(I,K)*B(K,J)) SUM2=SUM2-(A(I,K)*B(K,J)+B(I,K)*A(K,J)) 50 CONTINUE A(I,J)=SUM1 B(I,J)=SUM2 40 CONTINUE AAMAX=0. DO 60 I=J,N SUM1=A(I,J) SUM2=B(I,J) DO 70 K=1,J-1 SUM1=SUM1-(A(I,K)*A(K,J)-B(I,K)*B(K,J)) SUM2=SUM2-(A(I,K)*B(K,J)+B(I,K)*A(K,J)) 70 CONTINUE A(I,J)=SUM1 B(I,J)=SUM2 DUM=VV(I)*SQRT(SUM1**2+SUM2**2) IF (DUM.GE.AAMAX) THEN IMAX=I AAMAX=DUM ENDIF 60 CONTINUE IF (J.NE.IMAX) THEN DO 80 K=1,N DUM1=A(IMAX,K) DUM2=B(IMAX,K) A(IMAX,K)=A(J,K) B(IMAX,K)=B(J,K) A(J,K)=DUM1 B(J,K)=DUM2 80 CONTINUE DD=-DD VV(IMAX)=VV(J) ENDIF INDX(J)=IMAX IF ((A(J,J).EQ.0.).AND.(B(J,J).EQ.0.)) A(J,J)=TINY WRITE(*,*) 'A(J,J)', A(J,J) IF (J.NE.N) THEN ABSA2=A(J,J)**2+B(J,J)**2 WRITE(*,*) 'ABSA2', ABSA2 DUM1=A(J,J)/ABSA2 DUM2=-B(J,J)/ABSA2 DO 90 I=J+1,N TEM=A(I,J) A(I,J)=A(I,J)*DUM1-B(I,J)*DUM2 B(I,J)=B(I,J)*DUM1+TEM*DUM2 90 CONTINUE ENDIF 30 CONTINUE RETURN END *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CLUBKSB(A,B,N,NP,INDX,F,G) IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION A(NP,NP), B(NP,NP), INDX(N), F(NP), G(NP) II=0 DO 10 I=1,N LL=INDX(I) SUM1=F(LL) SUM2=G(LL) F(LL)=F(I) G(LL)=G(I) IF (II.NE.0) THEN DO 20 J=II,I-1 SUM1=SUM1-(A(I,J)*F(J)-B(I,J)*G(J)) SUM2=SUM2-(B(I,J)*F(J)+A(I,J)*G(J)) 20 CONTINUE ELSEIF ((SUM1.NE.0.).OR.(SUM2.NE.0.)) THEN II=I ENDIF F(I)=SUM1 G(I)=SUM2 10 CONTINUE DO 30 I=N,1,-1 SUM1=F(I) SUM2=G(I) IF (I.LT.N) THEN DO 40 J=I+1,N SUM1=SUM1-(A(I,J)*F(J)-B(I,J)*G(J)) SUM2=SUM2-(B(I,J)*F(J)+A(I,J)*G(J)) 40 CONTINUE ENDIF DEN=A(I,I)**2+B(I,I)**2 F(I)=(SUM1*A(I,I)+SUM2*B(I,I))/DEN G(I)=(SUM2*A(I,I)-SUM1*B(I,I))/DEN 30 CONTINUE RETURN END