C ALGORITHM 646 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, C SEPT., 1986, P. 278. C TEST DRIVER FOR PDFIND USING FULL AND BAND SYMMETRIC MATRICES. C C AUGUST 1985 C C.R. CRAWFORD C 39 MACPHERSON AVENUE C TORONTO, ONTARIO C CANADA M5R 1W7 C C PARAMETER (MAXROW=50,MAXCOL=50,PRUNIT=6) PARAMETER (FULL=1,BAND=2) PARAMETER (DEF=0,NDEF=-1,UABORT=1,FAIL=2) C********************************************************************** C C COMMON STORAGE FOR MATRIX DATA. C C A,B & C ARE ARRAYS FOR THE STORAGE OF THREE REAL SYMMETRIC C MATRICES. C C THE VECTORS X AND Z MUST BE AS LONG AS THE ORDER OF THE C MATRICES. C C IPVT IS AN ARRAY REQUIRED BY SOME LINPAC ROUTINES. C C TYPE = 1 (FULL) OR 2 (BAND) TO INDICATE MATRIX STORAGE. C INTEGER LDA,N,M,IPVT(MAXCOL),TYPE DOUBLE PRECISION X(MAXCOL),Z(MAXCOL),A(MAXROW,MAXCOL), . B(MAXROW,MAXCOL),C(MAXROW,MAXCOL) COMMON /MATRIX/A,B,C,X,Z,LDA,N,M,IPVT,TYPE C********************************************************************** C C PARAMETERS IN COMMON FOR TESTING AND CONTROL OF CHLSKY. C C LUP IS LOGICAL UNIT USED FOR PRINTING. C LIMIT = MAXIMUM NUMBER OF CALLS TO CHLSKY ALLOWED. C NITER = CURRENT NUMBER OF CALLS TO CHLSKY. INTEGER LUP,LIMIT,NITER COMMON /TEST/LUP,LIMIT,NITER C C********************************************************************** DOUBLE PRECISION S(2) EXTERNAL CHLSKY * LDA = MAXROW LUP = PRUNIT WRITE (LUP,9001) C C TESTS FOR THE BAND-SYMMETRIC CASE. C CYCLE THROUGH 7 EXAMPLES. C TYPE = BAND DO 10 JKL = 1,7 WRITE (LUP,9081) JKL CALL GETMAT(A,B,LDA,N,M,LUP,'BAND',JKL) CALL PRTMAT(A,B,LDA,N,M,LUP,'BAND') WRITE (LUP,9011) NITER = 0 IF (JKL.EQ.3) THEN LIMIT = 5 * ELSE IF (JKL.EQ.4) THEN LIMIT = 30 * ELSE LIMIT = 10 END IF C THE ITERATION LIMIT WILL BE EXCEEDED IN THE 3RD EXAMPLE. C CONVERGENCE FAILS FOR THE LAST EXAMPLE ON A VAX/11/780 C UNDER VMS. C C********************************************************************** C C FIND THE LINEAR COMBINATION. C CALL PDFIND(CHLSKY,S,INFO) C C********************************************************************** C C DISPLAY THE RESULTS. C IF (INFO.EQ.NDEF) THEN WRITE (LUP,9031) * ELSE IF (INFO.EQ.DEF) THEN WRITE (LUP,9041) WRITE (LUP,9021) S(1),S(2) * ELSE IF (INFO.EQ.UABORT) THEN WRITE (LUP,9051) NITER WRITE (LUP,9021) S(1),S(2) * ELSE IF (INFO.EQ.FAIL) THEN WRITE (LUP,9061) S(1) END IF * 10 CONTINUE C********************************************************************** C TESTS FOR SYMMETRIC MATRICES C TYPE = FULL C C CYCLE THROUGH 7 EXAMPLES. C DO 20 JKL = 1,7 WRITE (LUP,9091) JKL CALL GETMAT(A,B,LDA,N,M,LUP,'FULL',JKL) CALL PRTMAT(A,B,LDA,N,M,LUP,'FULL') WRITE (LUP,9011) NITER = 0 IF (JKL.EQ.3) THEN LIMIT = 5 * ELSE IF (JKL.EQ.4) THEN LIMIT = 30 * ELSE LIMIT = 10 END IF C THE ITERATION LIMIT WILL BE EXCEEDED IN THE 3RD EXAMPLE. C CONVERGENCE DOES NOT FAIL FOR THE 4TH EXAMPLE ON A VAX/11/780 C UNDER VMS THOUGH THE SAME PAIR FAILS AS A BAND-SYMMETRIC PAIR. C********************************************************************** C C FIND THE LINEAR COMBINATION. C CALL PDFIND(CHLSKY,S,INFO) C C********************************************************************** C C DISPLAY THE RESULTS. C IF (INFO.EQ.NDEF) THEN WRITE (LUP,9031) * ELSE IF (INFO.EQ.DEF) THEN WRITE (LUP,9041) WRITE (LUP,9021) S(1),S(2) * ELSE IF (INFO.EQ.UABORT) THEN WRITE (LUP,9051) NITER WRITE (LUP,9021) S(1),S(2) * ELSE IF (INFO.EQ.FAIL) THEN WRITE (LUP,9061) S(1) END IF C********************************************************************** 20 CONTINUE STOP * 9001 FORMAT (/' INTERMEDIATE RESULTS FROM CHLSKY.'/ . ' (PX,PY) ARE THE NEW VALUES FOUND BY CHLSKY.'/ . ' (SX,SY) ARE THE COEFFICIENTS OF THE LINEAR COMBINATION.'/ . ' AP = PX*SX + PY*SY?5 EXCEPT FOR ROUNDING ERROR, AP IS THE '/ . ' NEGATIVE DIAGONAL WHICH CAUSED THE FAILIURE OF CHOLESKY.'/ . ' KPD IS THE ORDER OF THE POSITIVE DEFINITE SUB-MATRIX.'/, . ' THE SAME MATRICES ARE USED IN THE BAND-SYMMETRIC AND FULL-'/, . ' SYMMETRIC EXAMPLES SO THE RESULTS SHOULD BE NEARLY IDENTICAL.' . /' FLOATING-POINT UNDERFLOWS MAY OCCUR IN DDOT?5 IN ALL CASES'/ . ' THE RESULT SHOULD BE SET TO ZERO.'//) 9011 FORMAT (/11X,'PX',19X,'PY',19X,'SX',19X,'SY',19X,'AP',11X, . 'KPD TYPE') 9021 FORMAT (' ALPHA = ',1PD21.12,'. BETA = ',1PD21.12) 9031 FORMAT (/' THE PAIR IS NOT DEFINITE.') 9041 FORMAT (/' THE PAIR IS DEFINITE.') 9051 FORMAT (/' TERMINATED BY CHLSKY AFTER ',I3,' ITERATIONS.') 9061 FORMAT (/' THE PAIR IS NOT CLEARLY DEFINITE. 1+COS(QR) = ',E9.3) 9071 FORMAT (///) 9081 FORMAT (/50X,' BAND-SYMMETRIC EXAMPLE ',I2) 9091 FORMAT (/50X,' SYMMETRIC EXAMPLE ',I2) END SUBROUTINE PDFIND(CHLSKY,S,INFO) EXTERNAL CHLSKY * DOUBLE PRECISION S(2) INTEGER INFO C C GIVEN REAL SYMMETRIC MATRICES A AND B FIND REAL NUMBERS S(1) AND C S(2) SUCH THAT C C C = S(1)*A + S(2)*B C C IS A POSITIVE DEFINITE MATRIX WHERE S(1)**2 + S(2)**2 = 1. C C INPUT: CHLSKY(K,U,V) - C C A SUBROUTINE WITH THREE PARAMETERS: AN INTEGER K AND C TWO DOUBLE PRECISION PARAMETERS U AND V. CHLSKY C TESTS IF THE MATRIX C C = U*A + V*B C IS POSITIVE DEFINITE AND REPORTS THE RESULTS THROUGH C THE PARAMETERS U, V AND K. C 1) IF THE C IS POSITIVE DEFINITE, K = 0 AND C U AND V ARE UNCHANGED. C 2) IF C IS NOT POSITIVE DEFINITE, K = -1 AND C U = X*A*X, C V = X*B*X, C WHERE X IS A VECTOR SUCH THAT C X*C*X .LE. 0. C NOTE - THE ABOVE INEQUALITY IMPLIES THAT IN THIS C CASE THE RETURNED VALUES OF U AND V MUST SATISFY: C U0*U + V0*V .LE. 0 TO WORKING PRECISION, C WHERE U0 AND V0 ARE THE VALUES OF U AND V ON ENTRY C TO CHLSKY. C 3) IF THE USER WISHES TO TERMINATE THE ITERATION, C K = 1 AND U AND V ARE UNCHANGED. C PARAMETER (PD=0,NPD=-1,UABORT=1) C C OUTPUT OF PDFIND: C INFO - FLAG INDICATING THE RESULTS; C SEE PARAMETER DEFINITIONS BELOW: C = 0; (A,B) IS A DEFINITE PAIR. C C = -1; (A,B) IS NOT A DEFINITE PAIR, I.E. C THERE EXISTS A VECTOR X SUCH THAT C X*A*X = X*B*X = 0. C C = 1; THE USER TERMINATED THE ITERATION C THROUGH THE SUBROUTINE CHLSKY. C C = 2; ITERATION FAILED TO CONVERGE. S(1) C CONTAINS 1+COSINE OF THE ANGLE C BETWEEN Q AND R. C C S(1),S(2) = THE SOLUTION IN CASE INFO = 0. C = THE LATEST ESTIMATE IN CASE INFO = 1. C S(1) = 1+COSINE OF THE ANGLE BETWEEN Q AND R C IN CASE INFO = 2. C PARAMETER (DEF=0,NDEF=-1,FAIL=2) C C AUGUST 1985 C C.R. CRAWFORD C 39 MACPHERSON AVENUE C TORONTO, ONTARIO C CANADA M5R 1W7 C C********************************************************************** DOUBLE PRECISION Q(2),R(2),NORM,SR,SQ,QR,CR,CQ,PREVQR INTEGER FLAG LOGICAL QUIT,UPDATE C********************************************************************** NFAC = 0 QR = 1.0D0 S(1) = 0.0D0 S(2) = 1.0D0 QUIT = .FALSE. UPDATE = .FALSE. C C LOOP UNTIL QUIT. C 10 CONTINUE PREVQR = QR C C IS THE COMBINATION POS. DEF. ? C CALL CHLSKY(FLAG,S(1),S(2)) NFAC = NFAC + 1 NORM = DMAX1(DABS(S(1)),DABS(S(2))) IF (NORM.GT.0.0D0) THEN NORM = NORM*DSQRT((S(1)/NORM)**2+ (S(2)/NORM)**2) S(1) = S(1)/NORM S(2) = S(2)/NORM END IF * IF (FLAG.EQ.PD) THEN INFO = DEF QUIT = .TRUE. * ELSE IF (FLAG.EQ.UABORT) THEN INFO = UABORT QUIT = .TRUE. * ELSE IF (NORM.EQ.0.0D0) THEN INFO = NDEF QUIT = .TRUE. * ELSE IF (NFAC.EQ.1) THEN C C AFTER THE FIRST TRY, SET Q AND R. C Q(1) = S(1) Q(2) = S(2) S(1) = Q(1) S(2) = Q(2) * ELSE IF (NFAC.EQ.2) THEN C C AFTER SECOND TRY KEEP Q AND SET R = S. C R(1) = S(1) R(2) = S(2) * ELSE C C AFTER THIRD OR MORE, COMPARE Q AND R TO S. C EXPRESS S AS A LINEAR COMB. OF Q AND R; C S = POSITIVE*(CQ*Q + CR*R). C SR = S(1)*R(1) + S(2)*R(2) SQ = S(1)*Q(1) + S(2)*Q(2) QR = Q(1)*R(1) + Q(2)*R(2) CQ = SQ - SR*QR CR = SR - SQ*QR C C OF THE 'QUADRANTS' IN THE PLANE FORMED BY Q, -Q, C R AND -R, WHERE IS S? C IF ((CQ.LE.0.D0) .AND. (CR.LE.0.D0)) THEN C C IF -Q & -R, THE PAIR IS INDEFINITE. C INFO = NDEF QUIT = .TRUE. * ELSE IF (CR.LE.0.D0) THEN C C IF +Q & -R, S REPLACES Q. C Q(1) = S(1) Q(2) = S(2) * ELSE IF (CQ.LE.0.D0) THEN C C IF -Q & +R, S REPLACES R. C R(1) = S(1) R(2) = S(2) END IF * END IF * IF ( .NOT. QUIT .AND. (NFAC.GT.1)) THEN C C SET S = MIDPOINT, GET THE NORM, C CHECK FOR ZERO NORM, AND NORMALIZE. C S(1) = Q(1) + R(1) S(2) = Q(2) + R(2) NORM = DMAX1(DABS(S(1)),DABS(S(2))) IF (NORM.GT.0.0D0) THEN NORM = NORM*DSQRT((S(1)/NORM)**2+ (S(2)/NORM)**2) S(1) = S(1)/NORM S(2) = S(2)/NORM * ELSE INFO = NDEF QUIT = .TRUE. END IF * END IF C C CHECK FOR NO-CONVERGENCE. C IF ((.NOT.QUIT) .AND. (NFAC.GT.2) .AND. (PREVQR.LE.QR)) THEN C C QR IS NO LONGER DECREASING TO -1. C INFO = FAIL S(1) = 1.D0 + QR QUIT = .TRUE. END IF * IF ( .NOT. QUIT) GO TO 10 C C END OF LOOP C RETURN END SUBROUTINE CHLSKY(FLAG,U,V) DOUBLE PRECISION U,V INTEGER FLAG C*********************************************************************** C AUXILIARY ROUTINE USED WITH PDFIND. C A SUBROUTINE WITH THREE PARAMETERS: AN INTEGER FLAG AND C TWO DOUBLE PRECISION PARAMETERS U AND V. CHLSKY C TESTS IF THE MATRIX C C = U*A + V*B C IS POSITIVE DEFINITE AND REPORTS THE RESULTS THROUGH C THE PARAMETERS U, V AND FLAG. SEE PARAMETER DEFINITIONS C BELOW. C 1) IF THE C IS POSITIVE DEFINITE, FLAG = 0 AND C U AND V ARE UNCHANGED. C 2) IF C IS NOT POSITIVE DEFINITE, FLAG = -1 AND C U = X*A*X, C V = X*B*X, C WHERE X IS A VECTOR SUCH THAT C X*C*X .LE. 0. C NOTE - THE ABOVE INEQUALITY IMPLIES THAT IN THIS C CASE THE RETURNED VALUES OF U AND V MUST SATISFY: C U0*U + V0*V .LE. 0 C TO WORKING PRECISION, WHERE U0 AND V0 ARE THE C VALUES OF U AND V ON ENTRY TO CHLSKY. C 3) IF THE USER WISHES TO TERMINATE THE ITERATION, C FLAG = 1 AND U AND V ARE UNCHANGED. C PARAMETER (PD=0,NPD=-1,UABORT=1) C C THIS IMPLEMENTATION IS FOR REAL SYMMETRIC MATRICES AND BAND- C SYMMETRIC MATRICES. IT CALLS ROUTINES FROM LINPACK(*) TO C ATTEMPT A CHOLESKY DECOMPOSITION AND SOLVE A SYSTEM OF EQUATIONS C TO COMPUTE THE RESULTS REQUIRED BY PDFIND IN CASE THE C FACTORIZATION FAILS. C (*) DONGARRA,BUNCH,MOLER AND STEWART, LINPACK USERS' GUIDE, C SIAM PUBLICATIONS, PHILADELPHIA, 1979. C C ROUTINES CALLED: C C DDOT - A LINPACK ROUTINE WHICH COMPUTES THE DOT PRODUCT C OF TWO VECTORS. C C DPBFA - A LINPACK ROUTINE WHICH ATTEMPTS A CHOLESKY C FACTORIZATION OF A REAL BANDED SYMMETRIC MATRIX. C THE LINPACK DOCUMENTATION DOES NOT SPECIFY THE C CONTENTS OF THE ARRAY IN CASE THE MATRIX IS C INDEFINITE. THUS THIS ROUTINE IS DEPENDENT ON A C PARTICULAR IMPLEMENTATION OF LINPACK. C C DGBSL - A LINPACK ROUTINE WHICH SOLVES A BANDED SYSTEM OF C EQUATIONS. FOR THIS APPLICATION WE HAVE A BANDED C UPPER-TRIANGULAR SYSTEM WHICH IS PRESENTED TO DGBSL C HAVING ZERO LOWER BANDWIDTH. C C DSBIP - A ROUTINE WHICH TAKES THE INNER PRODUCT OF A VECTOR C AND ITSELF USING A REAL SYMMETRIC BANDED MATRIX C STORED AS IN LINPACK. C C DPOFA - A LINPACK ROUTINE WHICH ATTEMPTS A CHOLESKY C FACTORIZATION OF A REAL SYMMETRIC MATRIX. C THE REMARKS ABOUT DPBFA ALSO APPLY TO THIS ROUTINE. C C DTRSL - A LINPACK ROUTINE WHICH SOLVES A TRIANGULAR SYSTEM C OF EQUATIONS. C C DGEIP - A ROUTINE WHICH TAKES THE INNER PRODUCT OF A VECTOR C WITH ITSELF USING A REAL MATRIX. C C C. R. CRAWFORD, TORONTO, ONTARIO. C AUGUST 1985 C*********************************************************************** PARAMETER (MAXROW=50,MAXCOL=50,PRUNIT=6) PARAMETER (FULL=1,BAND=2) C*********************************************************************** C C COMMON STORAGE FOR MATRIX DATA. C C A,B & C ARE ARRAYS FOR THE STORAGE OF THREE REAL SYMMETRIC C MATRICES. C C THE VECTORS X AND Z MUST BE AS LONG AS THE ORDER OF THE C MATRICES. C C IPVT IS AN ARRAY REQUIRED BY SOME LINPAC ROUTINES. C C TYPE = 1 (FULL) OR 2 (BAND) TO INDICATE MATRIX STORAGE. C INTEGER LDA,N,M,IPVT(MAXCOL),TYPE DOUBLE PRECISION X(MAXCOL),Z(MAXCOL),A(MAXROW,MAXCOL), . B(MAXROW,MAXCOL),C(MAXROW,MAXCOL) COMMON /MATRIX/A,B,C,X,Z,LDA,N,M,IPVT,TYPE C********************************************************************* C C PARAMETERS IN COMMON FOR TESTING AND CONTROL OF CHLSKY. C THESE QUANTITIES ARE SET IN THE DRIVER PROGRAM. C LUP IS LOGICAL UNIT USED FOR PRINTING. C LIMIT = MAXIMUM NUMBER OF CALLS TO CHLSKY ALLOWED. C NITER = CURRENT NUMBER OF CALLS TO CHLSKY. INTEGER LUP,LIMIT,NITER COMMON /TEST/LUP,LIMIT,NITER C C*********************************************************************** C C FUNCTIONS AND LOCAL VARIABLES. C DOUBLE PRECISION DSBIP,DGEIP,DDOT DOUBLE PRECISION U1,V1,NORM INTEGER INFO C*********************************************************************** IF (NITER.GT.LIMIT) THEN FLAG = UABORT RETURN * END IF * NITER = NITER + 1 C C SAVE U AND V FOR INTERMEDIATE PRINT-OUT. C U1 = U V1 = V C C SET UP THE LINEAR COMBINATION MATRIX U*A + V*B C IF (TYPE.EQ.BAND) THEN M1 = M + 1 DO 20 I = 1,N DO 10 J = 1,M1 IF (I.GT.M1-J) C(J,I) = U*A(J,I) + V*B(J,I) 10 CONTINUE 20 CONTINUE C ELSE IF (TYPE.EQ.FULL) THEN N1 = N - 1 DO 40 I = 1,N1 C(I,I) = U*A(I,I) + V*B(I,I) I1 = I + 1 DO 30 J = I1,N C(I,J) = U*A(I,J) + V*B(I,J) C(J,I) = C(I,J) 30 CONTINUE 40 CONTINUE C(N,N) = U*A(N,N) + V*B(N,N) END IF C********************************************************************* C C ATTEMPT A CHOLESKY FACTORIZATION. C IF (TYPE.EQ.BAND) THEN CALL DPBFA(C,LDA,N,M,INFO) * ELSE IF (TYPE.EQ.FULL) THEN CALL DPOFA(C,LDA,N,INFO) END IF * IF (INFO.NE.0) THEN KPD = INFO - 1 IF (KPD.GT.0) THEN C C THE KPD-BY-KPD PRINCIPAL MINOR IS POSITIVE DEFINITE, C AND ITS CHOLESKY FACTOR HAS BEEN COMPUTED. C MOVE C(1:KPD,KPD+1) (FULL STORAGE) INTO X(1:KPD). C DO 50 I = 1,KPD IF (TYPE.EQ.BAND) THEN KK = M + I - KPD IF (KK.GE.1) THEN X(I) = C(KK,KPD+1) * ELSE X(I) = 0.0D0 END IF * ELSE IF (TYPE.EQ.FULL) THEN X(I) = C(I,KPD+1) END IF * 50 CONTINUE C C********************************************************************* C C SOLVE THE KPD-BY-KPD TRIANGULAR SYSTEM. THE SOLUTION CANNOT C FAIL USING THE SUBMATRIX LEFT BY DP?FA. C IF (TYPE.EQ.BAND) THEN DO 60 I = 1,KPD IPVT(I) = I 60 CONTINUE CALL DGBSL(C,LDA,KPD,0,M,IPVT,X,0) * ELSE IF (TYPE.EQ.FULL) THEN KK = 0 CALL DTRSL(C,LDA,KPD,X,1,KK) END IF C C********************************************************************* END IF C C FILL IN THE REST OF X, NORMALIZE IT, C AND COMPUTE THE INNER PRODUCTS WITH A AND B. C IF KPD = 0, THEN X(1) = -1 C X(KPD+1) = -1.0D0 NORM = DSQRT(DDOT(KPD+1,X,1,X,1)) DO 70 I = 1,KPD + 1 X(I) = X(I)/NORM 70 CONTINUE IF (TYPE.EQ.BAND) THEN U = DSBIP(A,LDA,KPD+1,M,X,Z) V = DSBIP(B,LDA,KPD+1,M,X,Z) * ELSE IF (TYPE.EQ.FULL) THEN U = DGEIP(A,LDA,KPD+1,X,Z) V = DGEIP(B,LDA,KPD+1,X,Z) END IF C********************************************************************* C C DISPLAY INTERMEDIATE RESULTS C AP = U*U1 + V*V1 WRITE (LUP,9001) U,V,U1,V1,AP,KPD,TYPE C********************************************************************* END IF C IF (INFO.EQ.0) THEN C C IS POSITIVE DEFINITE. FLAG = PD * ELSE C CHOLESKY FAILED; U AND V HAVE BEEN RECOMPUTED. FLAG = NPD END IF * RETURN * 9001 FORMAT (1X,5 (1PD21.12),1X,I3,2X,I2) END SUBROUTINE GETMAT(A,B,LDA,N,M,LUP,TYPE,KPAR) INTEGER LDA,N,M,LUP,KPAR CHARACTER *4 TYPE DOUBLE PRECISION A(LDA,1),B(LDA,1) C********************************************************************** C C ROUTINE TO LOAD TEST EXAMPLES. C C SET A = TRANS(X)*D1*X C AND B = TRANS(X)*D2*X C C WHERE X IS AN UPPER-TRIANGULAR MATRIX WITH (HALF)BAND-WIDTH C M AND ENTRIES ALL ONES EXCEPT THE LOWER-RIGHT 5X5 WHICH ARE ALL C EPS, AND C C D1 = DIAG(SIN(V1), SIN(V2), ..., SIN(VN)) C D2 = DIAG(COS(V1), COS(V2), ..., COS(VN)) C VI = (T(I) + R)*PIBY2. C C THE MATRIX STORAGE FORMAT IS DETERMINED BY TYPE. C TYPE = FULL - N0RMAL FULL MATRIX STORAGE. C BAND - BAND STORAGE AS PER LINPAC, TYPE 'SB' C C C.R. CRAWFORD, TORONTO, ONTARIO, AUG. 1985. C******************************************************************** DOUBLE PRECISION T(20),R,PIBY2,SA,SB,FAC,EPS PIBY2 = DATAN2(1.0D0,0.0D0) C SET INITIAL T(I) = 2**(-I+2) SO THAT THE V(I) NEARLY C FILL A SEMI-CIRCLE. T(1) = 2.0D0 DO 10 I = 2,20 T(I) = T(I-1)/2.0D0 10 CONTINUE EPS = 1.0D0 R = 0.0D0 IF (KPAR.EQ.1) THEN C C THIS EXAMPLE IS A DEFINITE PAIR.(+-) C M = 3 N = 5 R = 1.0D0 * ELSE IF (KPAR.EQ.2) THEN C C THIS EXAMPLE IS A DEFINITE PAIR.(-+) C M = 3 N = 5 T(1) = T(5) T(1) = 2.0D0 * ELSE IF (KPAR.EQ.3) THEN C C THIS PAIR IS DEFINITE BUT WILL FAIL IF C THE ITERATION LIMIT SET BY THE DRIVER IS TOO SMALL. C M = 4 N = 10 R = 1.0D0 * ELSE IF (KPAR.EQ.4) THEN C C THIS PAIR IS SUCH THAT ON VAX/8600 UNDER VMS C PDFIND DOES NOT CONVERGE IN BAND FORM BUT DOES CONVERGE C IN FULL SYMMETRIC FORM. C M = 5 N = 20 EPS = 0.5D-6 R = 1.1 * ELSE IF (KPAR.EQ.5) THEN C C FIRST DIAGONALS WILL BE SET TO ZERO. C SO THAT CHLSKY WILL RETURN (0,0). C M = 3 N = 9 * ELSE IF (KPAR.EQ.6) THEN C C THIS EXAMPLE IS A NON-DEFINITE PAIR. C M = 4 N = 10 R = 1.0D0 T(10) = -1.0D0 * ELSE IF (KPAR.EQ.7) THEN C C THIS EXAMPLE IS NON-DEFINITE WITH Q = -R. C M = 0 N = 10 T(2) = 0.0 END IF C C SET UP THE MATRICES BASED ON T(*) AND R. C IF (TYPE.EQ.'BAND') THEN C C USE BAND-MATRIX STORAGE A' LA LINPAC. C SEE PAGE 4.3 OF LINPACK USER'S GUIDE. C DO 40 J = 1,N I1 = MAX0(1,J-M) DO 30 I = I1,J SA = 0.0D0 SB = 0.0D0 DO 20 II = I1,I FAC = 1.0D0 IF (II.GT.N-5) THEN IF (I.GT.N-5) FAC = EPS IF (J.GT.N-5) FAC = FAC*EPS END IF * SA = SA + FAC*DSIN((T(II)+R)*PIBY2) SB = SB + FAC*DCOS((T(II)+R)*PIBY2) 20 CONTINUE K = I - J + M + 1 A(K,J) = SA B(K,J) = SB 30 CONTINUE 40 CONTINUE IF (KPAR.EQ.5) THEN C C CAUSE CHLSKY TO RETURN U = 0; V = 0. C A(M+1,1) = 0.0D0 B(M+1,1) = 0.0D0 * ELSE IF (KPAR.EQ.7) THEN C C FORCE Q = -R C A(M+1,1) = 0.0D0 A(M+1,2) = 0.0D0 END IF * ELSE IF (TYPE.EQ.'FULL') THEN C C USE FULL MATRIX STORAGE. C DO 70 J = 1,N I1 = MAX0(1,J-M) DO 60 I = 1,J SA = 0.0D0 SB = 0.0D0 IF (I.GE.I1) THEN DO 50 II = I1,I SA = SA + DSIN((T(II)+R)*PIBY2) SB = SB + DCOS((T(II)+R)*PIBY2) 50 CONTINUE END IF * A(I,J) = SA B(I,J) = SB A(J,I) = SA B(J,I) = SB 60 CONTINUE 70 CONTINUE IF (KPAR.EQ.5) THEN C C CAUSE CHLSKY TO RETURN U = 0; V = 0. C A(1,1) = 0.0D0 B(1,1) = 0.0D0 END IF * END IF END SUBROUTINE PRTMAT(A,B,LDA,N,M,LUP,TYPE) INTEGER LDA,N,M,LUP CHARACTER *4 TYPE DOUBLE PRECISION A(LDA,1),B(LDA,1) C********************************************************************** C C ROUTINE TO PRINT PAIRS OF TEST MATRICES. C C C.R. CRAWFORD, TORONTO, ONTARIO, AUG. 1985. C********************************************************************** WRITE (LUP,9011) N,M IF (TYPE.EQ.'BAND') THEN C C USE BAND-MATRIX STORAGE A' LA LINPAC. C PRINT THE LOWER TRIANGLE BY ROWS STOPPING AT THE DIAGONAL. C SEE PAGE 4.3 OF LINPACK USER'S GUIDE. C WRITE (LUP,*) ' MATRIX A ' DO 10 J = 1,N I1 = MAX0(1,J-M) WRITE (LUP,9001) (A(I-J+M+1,J),I=I1,J) 10 CONTINUE WRITE (LUP,9021) WRITE (LUP,*) ' MATRIX B ' DO 20 J = 1,N I1 = MAX0(1,J-M) WRITE (LUP,9001) (B(I-J+M+1,J),I=I1,J) 20 CONTINUE WRITE (LUP,9021) * ELSE C C USE FULL MATRIX STORAGE. C WRITE (LUP,*) ' MATRIX A ' DO 30 I = 1,N WRITE (LUP,9001) (A(I,J),J=1,I) 30 CONTINUE WRITE (LUP,9021) WRITE (LUP,*) ' MATRIX B ' DO 40 I = 1,N WRITE (LUP,9001) (B(I,J),J=1,I) 40 CONTINUE WRITE (LUP,9021) END IF * RETURN * 9001 FORMAT (10E13.5) 9011 FORMAT (' N = ',I5,'?5 M = ',I5) 9021 FORMAT (/) END DOUBLE PRECISION FUNCTION DGEIP(A,LDA,N,X,Z) C C COMPUTE THE INNER-PRODUCT TRANS(X)*A*X WHERE A IS A C REAL MATRIX. C C INPUT: C A THE MATRIX STORED AS LINPAC TYPE "GE". C LDA FIRST DIMENSION OF THE ARRAY A. C N ORDER OF THE MATRIX AND VECTOR. C X VECTOR IN THE INNER PRODUCT. C Z SCRATCH VECTOR. C C ROUTINES CALLED: DDOT C C APRIL 1983, C.R. CRAWFORD. C********************************************************************** DOUBLE PRECISION A(LDA,N),X(N),Z(N) DOUBLE PRECISION DDOT DO 10 I = 1,N Z(I) = DDOT(N,X,1,A(I,1),LDA) 10 CONTINUE DGEIP = DDOT(N,X,1,Z,1) RETURN END DOUBLE PRECISION FUNCTION DSBIP(ABD,LDA,N,M,X,Z) C C COMPUTE THE INNER-PRODUCT TRANS(X)*A*X WHERE A IS A C REAL SYMMETRIC BAND MATRIX. C C INPUT: C ABD THE MATRIX STORED AS LINPAC TYPE "SB". C LDA FIRST DIMENSION OF THE ARRAY A. C N ORDER OF THE MATRIX AND VECTOR. C M HALF BANDWIDTH OF THE MATRIX A, THAT IS C A(I,J) = 0 IF ABS(I-J) > M. C X VECTOR IN THE INNER PRODUCT. C Z SCRATCH VECTOR. C ROUTINES CALLED: DDOT C C APRIL 1983 - C.R. CRAWFORD C********************************************************************** DOUBLE PRECISION ABD(LDA,N),X(N),Z(N) DOUBLE PRECISION DDOT DO 10 I = 1,N L1 = MIN0(M+1,I) L2 = MIN0(M,N-I) K1 = I - L1 + 1 KK = M + 2 - L1 Z(I) = DDOT(L1,X(K1),1,ABD(KK,I),1) IF (I.NE.N) Z(I) = Z(I) + DDOT(L2,X(I+1),1,ABD(M,I+1),LDA-1) 10 CONTINUE DSBIP = DDOT(N,X,1,Z,1) RETURN END DOUBLE PRECISION FUNCTION DNRM2(N,DX,INCX) INTEGER NEXT DOUBLE PRECISION DX(1),CUTLO,CUTHI,HITEST,SUM,XMAX,ZERO,ONE DATA ZERO,ONE/0.0D0,1.0D0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO.(UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO,CUTHI/8.232D-11,1.304D19/ C IF (N.GT.0) GO TO 10 DNRM2 = ZERO GO TO 140 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N*INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT, (30,40,70,80) * 30 IF (DABS(DX(I)).GT.CUTLO) GO TO 110 ASSIGN 40 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 40 IF (DX(I).EQ.ZERO) GO TO 130 IF (DABS(DX(I)).GT.CUTLO) GO TO 110 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 60 C C PREPARE FOR PHASE 4. C 50 I = J ASSIGN 80 TO NEXT SUM = (SUM/DX(I))/DX(I) 60 XMAX = DABS(DX(I)) GO TO 90 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF (DABS(DX(I)).GT.CUTLO) GO TO 100 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 80 IF (DABS(DX(I)).LE.XMAX) GO TO 90 SUM = ONE + SUM* (XMAX/DX(I))**2 XMAX = DABS(DX(I)) GO TO 130 C 90 SUM = SUM + (DX(I)/XMAX)**2 GO TO 130 C C C PREPARE FOR PHASE 3. C 100 SUM = (SUM*XMAX)*XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 110 HITEST = CUTHI/FLOAT(N) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 120 J = I,NN,INCX IF (DABS(DX(J)).GE.HITEST) GO TO 50 SUM = SUM + DX(J)**2 120 CONTINUE DNRM2 = DSQRT(SUM) GO TO 140 C 130 CONTINUE I = I + INCX IF (I.LE.NN) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX*DSQRT(SUM) 140 CONTINUE RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,IXIY,M,MP1,N C IF (N.LE.0) RETURN IF (DA.EQ.0.0D0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF (M.EQ.0) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF (N.LT.4) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I+1) = DY(I+1) + DA*DX(I+1) DY(I+2) = DY(I+2) + DA*DX(I+2) DY(I+3) = DY(I+3) + DA*DX(I+3) 50 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF (M.EQ.0) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF (N.LT.5) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) + . DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) 50 CONTINUE 60 DDOT = DTEMP RETURN END SUBROUTINE DPBFA(ABD,LDA,N,M,INFO) INTEGER LDA,N,M,INFO DOUBLE PRECISION ABD(LDA,1) C C DPBFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX STORED IN BAND FORM. C C DPBFA IS USUALLY CALLED BY DPBCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE MATRIX TO BE FACTORED. THE COLUMNS OF THE UPPER C TRIANGLE ARE STORED IN THE COLUMNS OF ABD AND THE C DIAGONALS OF THE UPPER TRIANGLE ARE STORED IN THE C ROWS OF ABD . SEE THE COMMENTS BELOW FOR DETAILS. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C LDA MUST BE .GE. M + 1 . C C N INTEGER C THE ORDER OF THE MATRIX A . C C M INTEGER C THE NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C 0 .LE. M .LT. N . C C ON RETURN C C ABD AN UPPER TRIANGULAR MATRIX R , STORED IN BAND C FORM, SO THAT A = TRANS(R)*R . C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K IF THE LEADING MINOR OF ORDER K IS NOT C POSITIVE DEFINITE. C C BAND STORAGE C C IF A IS A SYMMETRIC POSITIVE DEFINITE BAND MATRIX, C THE FOLLOWING PROGRAM SEGMENT WILL SET UP THE INPUT. C C M = (BAND WIDTH ABOVE DIAGONAL) C DO 20 J = 1, N C I1 = MAX0(1, J-M) C DO 10 I = I1, J C K = I-J+M+1 C ABD(K,J) = A(I,J) C 10 CONTINUE C 20 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 DDOT C FORTRAN MAX0,DSQRT C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,T DOUBLE PRECISION S INTEGER IK,J,JK,K,MU C BEGIN BLOCK WITH ...EXITS TO 40 C C DO 30 J = 1,N INFO = J S = 0.0D0 IK = M + 1 JK = MAX0(J-M,1) MU = MAX0(M+2-J,1) IF (M.LT.MU) GO TO 20 DO 10 K = MU,M T = ABD(K,J) - DDOT(K-MU,ABD(IK,JK),1,ABD(MU,J),1) T = T/ABD(M+1,JK) ABD(K,J) = T S = S + T*T IK = IK - 1 JK = JK + 1 10 CONTINUE 20 CONTINUE S = ABD(M+1,J) - S C ......EXIT IF (S.LE.0.0D0) GO TO 40 ABD(M+1,J) = DSQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END SUBROUTINE DGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB) INTEGER LDA,N,ML,MU,IPVT(1),JOB DOUBLE PRECISION ABD(LDA,1),B(1) C C DGBSL SOLVES THE DOUBLE PRECISION BAND SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY DGBCO OR DGBFA. C C ON ENTRY C C ABD DOUBLE PRECISION(LDA, N) C THE OUTPUT FROM DGBCO OR DGBFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM DGBCO OR DGBFA. 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 DGBCO HAS SET RCOND .GT. 0.0 C OR DGBFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL DGBSL(ABD,LDA,N,ML,MU,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 FORTRAN MIN0 C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,T INTEGER K,KB,L,LA,LB,LM,M,NM1 C M = MU + ML + 1 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 (ML.EQ.0) GO TO 30 IF (NM1.LT.1) GO TO 30 DO 20 K = 1,NM1 LM = MIN0(ML,N-K) L = IPVT(K) T = B(L) IF (L.EQ.K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL DAXPY(LM,T,ABD(M+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)/ABD(M,K) LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = -B(K) CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),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 LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = DDOT(LM,ABD(LA,K),1,B(LB),1) B(K) = (B(K)-T)/ABD(M,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (ML.EQ.0) GO TO 90 IF (NM1.LT.1) GO TO 90 DO 80 KB = 1,NM1 K = N - KB LM = MIN0(ML,N-K) B(K) = B(K) + DDOT(LM,ABD(M+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 DTRSL(T,LDT,N,B,JOB,INFO) INTEGER LDT,N,JOB,INFO DOUBLE PRECISION T(LDT,1),B(1) C C C DTRSL SOLVES SYSTEMS OF THE FORM C C T * X = B C OR C TRANS(T) * X = B C C WHERE T IS A TRIANGULAR MATRIX OF ORDER N. HERE TRANS(T) C DENOTES THE TRANSPOSE OF THE MATRIX T. C C ON ENTRY C C T DOUBLE PRECISION(LDT,N) C T CONTAINS THE MATRIX OF THE SYSTEM. THE ZERO C ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND C THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE C USED TO STORE OTHER INFORMATION. C C LDT INTEGER C LDT IS THE LEADING DIMENSION OF THE ARRAY T. C C N INTEGER C N IS THE ORDER OF THE SYSTEM. C C B DOUBLE PRECISION(N). C B CONTAINS THE RIGHT HAND SIDE OF THE SYSTEM. C C JOB INTEGER C JOB SPECIFIES WHAT KIND OF SYSTEM IS TO BE SOLVED. C IF JOB IS C C 00 SOLVE T*X=B, T LOWER TRIANGULAR, C 01 SOLVE T*X=B, T UPPER TRIANGULAR, C 10 SOLVE TRANS(T)*X=B, T LOWER TRIANGULAR, C 11 SOLVE TRANS(T)*X=B, T UPPER TRIANGULAR. C C ON RETURN C C B B CONTAINS THE SOLUTION, IF INFO .EQ. 0. C OTHERWISE B IS UNALTERED. C C INFO INTEGER C INFO CONTAINS ZERO IF THE SYSTEM IS NONSINGULAR. C OTHERWISE INFO CONTAINS THE INDEX OF C THE FIRST ZERO DIAGONAL ELEMENT OF T. C C LINPACK. THIS VERSION DATED 08/14/78 . C G. W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS DAXPY,DDOT C FORTRAN MOD C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,TEMP INTEGER CASE,J,JJ C C BEGIN BLOCK PERMITTING ...EXITS TO 150 C C CHECK FOR ZERO DIAGONAL ELEMENTS. C DO 10 INFO = 1,N C ......EXIT IF (T(INFO,INFO).EQ.0.0D0) GO TO 150 10 CONTINUE INFO = 0 C C DETERMINE THE TASK AND GO TO IT. C CASE = 1 IF (MOD(JOB,10).NE.0) CASE = 2 IF (MOD(JOB,100)/10.NE.0) CASE = CASE + 2 GO TO (20,50,80,110),CASE C C SOLVE T*X=B FOR T LOWER TRIANGULAR C 20 CONTINUE B(1) = B(1)/T(1,1) IF (N.LT.2) GO TO 40 DO 30 J = 2,N TEMP = -B(J-1) CALL DAXPY(N-J+1,TEMP,T(J,J-1),1,B(J),1) B(J) = B(J)/T(J,J) 30 CONTINUE 40 CONTINUE GO TO 140 C C SOLVE T*X=B FOR T UPPER TRIANGULAR. C 50 CONTINUE B(N) = B(N)/T(N,N) IF (N.LT.2) GO TO 70 DO 60 JJ = 2,N J = N - JJ + 1 TEMP = -B(J+1) CALL DAXPY(J,TEMP,T(1,J+1),1,B(1),1) B(J) = B(J)/T(J,J) 60 CONTINUE 70 CONTINUE GO TO 140 C C SOLVE TRANS(T)*X=B FOR T LOWER TRIANGULAR. C 80 CONTINUE B(N) = B(N)/T(N,N) IF (N.LT.2) GO TO 100 DO 90 JJ = 2,N J = N - JJ + 1 B(J) = B(J) - DDOT(JJ-1,T(J+1,J),1,B(J+1),1) B(J) = B(J)/T(J,J) 90 CONTINUE 100 CONTINUE GO TO 140 C C SOLVE TRANS(T)*X=B FOR T UPPER TRIANGULAR. C 110 CONTINUE B(1) = B(1)/T(1,1) IF (N.LT.2) GO TO 130 DO 120 J = 2,N B(J) = B(J) - DDOT(J-1,T(1,J),1,B(1),1) B(J) = B(J)/T(J,J) 120 CONTINUE 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END SUBROUTINE DPOFA(A,LDA,N,INFO) INTEGER LDA,N,INFO DOUBLE PRECISION A(LDA,1) C C DPOFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE C MATRIX. C C DPOFA IS USUALLY CALLED BY DPOCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR DPOCO) = (1 + 18/N)*(TIME FOR DPOFA) . C C ON ENTRY C C A DOUBLE PRECISION(LDA, N) C THE SYMMETRIC MATRIX TO BE FACTORED. ONLY THE C DIAGONAL AND UPPER TRIANGLE ARE USED. 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 R SO THAT A = TRANS(R)*R C WHERE TRANS(R) IS THE TRANSPOSE. C THE STRICT LOWER TRIANGLE IS UNALTERED. C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR C OF ORDER K IS NOT POSITIVE DEFINITE. 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 DDOT C FORTRAN DSQRT C C INTERNAL VARIABLES C DOUBLE PRECISION DDOT,T DOUBLE PRECISION S INTEGER J,JM1,K C BEGIN BLOCK WITH ...EXITS TO 40 C C DO 30 J = 1,N INFO = J S = 0.0D0 JM1 = J - 1 IF (JM1.LT.1) GO TO 20 DO 10 K = 1,JM1 T = A(K,J) - DDOT(K-1,A(1,K),1,A(1,J),1) T = T/A(K,K) A(K,J) = T S = S + T*T 10 CONTINUE 20 CONTINUE S = A(J,J) - S C ......EXIT IF (S.LE.0.0D0) GO TO 40 A(J,J) = DSQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END