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