C
C     LINALG
C
C     **** VERSION OF AUGUST 15, 1983
C
C     LINALG IS A SET OF FOUR SUBROUTINES AND THREE FUNCTIONS
C     THAT PERFORM LINEAR ALGEBRA TASKS.
C
C     THE SUBROUTINES ARE QRGEN, QRHESS, SOLVE, AND QRUPDAT.
C         THEY RESPECTIVELY OBTAIN THE QR DECOMPOSITION OF A
C         GENERAL MATRIX, OF AN UPPER HESSENBERG MATRIX, SOLVE
C         A LINEAR SYSTEM WHEN SUCH A FACTORIZATION OF THE
C         COEFFICIENT MATRIX IS KNOWN, AND UPDATE AN UPPER
C         TRIANGULAR MATRIX BY ROTATIONS WHEN A RANK ONE MATRIX
C         IS ADDED TO IT.
C
C     THE FUNCTIONS ARE DOTPRD, V2NORM AND VSNORM, WHICH RESPECTIVELY
C         COMPUTE THE INNER PRODUCT OF TWO VECTORS, THE 2-NORM OF
C         A SINGLE VECTOR AND THE SUP-NORM OF A SINGLE VECTOR.
C
      SUBROUTINE QRGEN(M,N,A,LDA,PIVOT,IPVT,SIGMA,WA)
      INTEGER M,N,LDA
      LOGICAL PIVOT
      INTEGER IPVT(N)
      DOUBLE PRECISION A(LDA,N),SIGMA(N),WA(N)
C     **********
C
C     SUBROUTINE QRGEN (MODIFIED FROM MORE'S PROGRAM QR)
C
C     THIS SUBROUTINE USES HOUSEHOLDER TRANSFORMATIONS AND (OPTIONAL)
C     COLUMN PIVOTING TO FIND A QR DECOMPOSITION OF THE M BY N MATRIX
C     A. THUS QRGEN DETERMINES AN ORTHOGONAL MATRIX Q, A PERMUTATION P,
C     AND AN UPPER TRAPEZOIDAL MATRIX R SUCH THAT A*P = Q*R.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE QRGEN(M,N,A,LDA,PIVOT,IPVT,SIGMA,WA)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER NOT GREATER THAN LDA.
C
C       N IS A POSITIVE INTEGER NOT GREATER THAN M.
C
C       A IS AN M BY N ARRAY. THIS SUBROUTINE WILL COMPUTE THE
C         QR DECOMPOSITION OF THE INPUT A. ON OUTPUT THE UPPER
C         TRAPEZOIDAL PART OF A CONTAINS R.
C
C       LDA IS THE DECLARED ROW DIMENSION OF A
C
C       PIVOT IS A LOGICAL INPUT VARIABLE. IF PIVOT IS SET .TRUE.,
C         THEN COLUMN PIVOTING IS ENFORCED AMONG THE FIRST LDA - 1
C         COLUMNS. IF PIVOT IS SET .FALSE.,
C         THEN THERE IS NO PIVOTING.
C
C       IPVT IS AN INTEGER ARRAY OF LENGTH N WHICH DEFINES THE PERMUTATI
C         P SUCH THAT A*P = Q*R. THE J-TH COLUMN OF P IS THE IVPT(J)-TH
C         COLUMN OF THE IDENTITY. ALTERNATIVELY, A(.,IPVT(.) = Q*R.
C
C       SIGMA AND WA ARE SCRATCH ARRAYS OF LENGTH N.
C
C     FUNCTIONS USED -- V2NORM, DMAX1, DSQRT
C
C     ********
      INTEGER I,J,JMAX,K,KP1,NN
      DOUBLE PRECISION AKNORM,ONE,P05,SUM,TEMP,ZERO
      DOUBLE PRECISION DMAX1,DSQRT,V2NORM
      DATA ZERO,ONE,P05 /0.0D0,1.0D0,0.05D0/
C
C     BEGIN PROGRAM
C
C     INITIALIZE THE PERMUTATION ARRAYS AND THE ARRAYS SIGMA AND WA
C
      IF (.NOT. PIVOT) GO TO 15
      DO 10 J=1,N
         SIGMA(J)=V2NORM(M,A(1,J))
         WA(J)=SIGMA(J)
         IPVT(J)=J
   10    CONTINUE
      NN=LDA-1
      IF (NN.GT.N) NN=N
C
C     REDUCTION OF A TO UPPER TRAPEZOIDAL FORM
C
   15 DO 100 K=1,N
         IF (.NOT. PIVOT) GO TO 40
C
C        BRING COLUMN OF LARGEST NORM INTO PIVOT POSITION
C
         AKNORM=ZERO
         JMAX=K
         DO 20 J=K,NN
            IF (SIGMA(J).LE.AKNORM) GO TO 20
            AKNORM=SIGMA(J)
            JMAX=J
   20       CONTINUE
         IF (JMAX.EQ.K) GO TO 40
         DO 30 I=1,M
            TEMP=A(I,K)
            A(I,K)=A(I,JMAX)
            A(I,JMAX)=TEMP
   30       CONTINUE
         SIGMA(JMAX)=SIGMA(K)
         WA(JMAX)=WA(K)
         J=IPVT(JMAX)
         IPVT(JMAX)=IPVT(K)
         IPVT(K)=J
   40    CONTINUE
C
C        COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN K
C
         AKNORM=V2NORM(M-K+1,A(K,K))
         IF (AKNORM.EQ.ZERO) GO TO 90
         IF (A(K,K).LT.ZERO) AKNORM=-AKNORM
         DO 50 I=K,M
            A(I,K)=A(I,K)/AKNORM
   50       CONTINUE
         A(K,K)=ONE+A(K,K)
C
C        APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS
C        AND UPDATE THE NORMS.
C
         IF (K.EQ.N) GO TO 90
         KP1=K+1
         DO 80 J=KP1,N
            SUM=ZERO
            DO 60 I=K,M
               SUM=SUM+A(I,K)*A(I,J)
   60          CONTINUE
            TEMP=SUM/A(K,K)
            DO 70 I=K,M
               A(I,J)=A(I,J)-TEMP*A(I,K)
   70          CONTINUE
            IF (.NOT. PIVOT) GO TO 80
            IF (SIGMA(J).EQ.ZERO) GO TO 80
            TEMP=A(K,J)/SIGMA(J)
            SIGMA(J)=SIGMA(J)*DSQRT(DMAX1(ZERO,ONE-TEMP*TEMP))
            TEMP=ONE+P05*(SIGMA(J)/WA(J))**2
            IF (TEMP.NE.ONE) GO TO 80
            SIGMA(J)=V2NORM(M-K,A(KP1,J))
           WA(J)=SIGMA(J)
   80       CONTINUE
   90    CONTINUE
         A(K,K)=-AKNORM
         IF (K.EQ.M) GO TO 100
         DO 95 I=KP1,M
            A(I,K)=ZERO
   95       CONTINUE
  100    CONTINUE
      RETURN
C
C     LAST CARD OF QRGEN FOLLOWS
C
      END
      SUBROUTINE QRHESS(R,LDR,M,P)
      INTEGER LDR,M,P
      DOUBLE PRECISION R(LDR,M)
C
C     QRHESS USES JACOBI ROTATIONS TO ZERO OUT THE SUBDIAGONAL
C     OF THE M BY M UPPER HESSENBERG MATRIX R.
C
C     LDR IS THE DECLARED ROW DIMENSION OF R
C
C     P IS THE COLUMN INDEX OF THE FIRST NONZERO SUBDIAGONAL ENTRY
C
C     FUNCTIONS USED -- DOTPRD, V2NORM
C
      INTEGER I,IP1,J,MM1
      DOUBLE PRECISION X(2),Y(2),TEMP,DOTPRD,V2NORM,ZERO
      DATA ZERO /0.0D0/
C
C     BEGIN PROGRAM
C
      IF (P.GE.M) RETURN
      MM1=M-1
      DO 100 I=P,MM1
C
C        COMPUTE THE ROTATION FOR COLUMN I
C
         TEMP=V2NORM(2,R(I,I))
         IF (TEMP.EQ.ZERO) GO TO 100
         IP1=I+1
         X(1)=R(I,I)/TEMP
         X(2)=R(IP1,I)/TEMP
         R(I,I)=TEMP
         R(IP1,I)=ZERO
         DO 50 J=IP1,M
C
C           APPLY ROTATION TO REMAINING COLUMNS
C
            Y(2)=-R(I,J)
            R(I,J)=DOTPRD(2,X,R(I,J))
            Y(1)=R(IP1,J)
            R(IP1,J)=DOTPRD(2,X,Y)
   50       CONTINUE
  100    CONTINUE
      RETURN
C
C     LAST CARD OF QRHESS FOLLOWS
C
      END
      SUBROUTINE SOLVE(A,LDA,M,LPERM,R,LDR,B,X,QTB,RES,IFREF)
      INTEGER LDA,LDR,M,LPERM(M)
      DOUBLE PRECISION A(LDA,M),R(LDR,M),B(M),X(M),QTB(M),RES(M)
      LOGICAL IFREF
C
C     SOLVE OBTAINS A SOLUTION TO A * P * X = B FOR A CERTAIN
C     PERMUTATION MATRIX P.
C
C     A IS AN M BY M MATRIX WITH DECLARED ROW DIMENSION LDA.
C     LPERM IS AN M-VECTOR GIVING A PERMUTATION OF (1,2,..,M).
C     LET P BE THE PERMUTATION MATRIX WITH P(L,J)=1 IFF LPERM(J)=L.
C     THEN A * P IS SOME ORTHOGONAL MATRIX Q TIMES THE UPPER
C     TRIANGULAR MATRIX R. ON RETURN, X SOLVES THE SYSTEM ABOVE
C     AND QTB IS Q(TRANSPOSE) * B = R * X. IF QTB IS NOT REQUIRED,
C     IT CAN SHARE STORAGE WITH X. IF IFREF = .TRUE., ONE STEP OF
C     ITERATIVE REFINEMENT IS PERFORMED USING THE RESIDUAL RES.
C     IF IFREF = .FALSE. OR IF B CAN BE ALTERED THEN RES CAN
C     SHARE STORAGE WITH B.
C
C     FUNCTIONS USED -- DOTPRD, VSNORM
C
      INTEGER I,J,JJ,JM1,L
      DOUBLE PRECISION DOTPRD,VSNORM,TEMP,TEMPX,TQTB1,TQTB2,ZERO
      DATA ZERO /0.0D0/
C
C     BEGIN PROGRAM
C
C     COMPUTE P(TRANSPOSE) * A(TRANSPOSE) * B  AND WRITE IN X
C
      DO 10 J=1,M
         L=LPERM(J)
         X(J)=DOTPRD(M,A(1,L),B)
   10    CONTINUE
C
C     SOLVE R(TRANSPOSE) * Y = X  AND WRITE THE SOLUTION IN X
C
      X(1)=X(1)/R(1,1)
      IF (M.EQ.1) GO TO 30
      DO 20 J=2,M
         X(J)=(X(J)-DOTPRD(J-1,R(1,J),X(1)))/R(J,J)
   20    CONTINUE
C
C     COPY X INTO QTB
C
   30 DO 40 J=1,M
         QTB(J)=X(J)
   40    CONTINUE
C
C     SOLVE R * Y = X  AND WRITE THE SOLUTION IN X
C
      DO 70 JJ=1,M
         J=M+1-JJ
         TEMP=X(J)/R(J,J)
         X(J)=TEMP
         IF (JJ.EQ.M) GO TO 60
         JM1=J-1
         DO 50 I=1,JM1
            X(I)=X(I)-TEMP*R(I,J)
   50       CONTINUE
   60    CONTINUE
   70    CONTINUE
      IF (.NOT. IFREF) GO TO 300
C
C     COMPUTE RESIDUAL RES.
C
      TEMPX=VSNORM(M,B)
      DO 80 I=1,M
         RES(I)=B(I)
   80    CONTINUE
      DO 100 J=1,M
         L=LPERM(J)
         TEMP=X(J)
         DO 90 I=1,M
            RES(I)=RES(I)-A(I,L)*TEMP
   90       CONTINUE
  100    CONTINUE
C
C     NO REFINEMENT IF RESIDUAL SMALL
C
      TEMP=VSNORM(M,RES)
      IF ((TEMPX+TEMP).EQ.TEMPX) GO TO 300
C
C     REFINEMENT. FIRST STORE X, QTB IN UNUSED PARTS
C     OF FIRST, SECOND COLUMNS OF R.
C
      TEMPX=X(1)
      TQTB1=QTB(1)
      IF (M.EQ.1) GO TO 120
      R(2,1)=X(2)
      TQTB2=QTB(2)
      IF (M.EQ.2) GO TO 120
      DO 110 I=3,M
         R(I,1)=X(I)
         R(I,2)=QTB(I)
  110    CONTINUE
  120 CONTINUE
C
C     COMPUTE P(TRANSPOSE) * A(TRANSPOSE) * RES  AND WRITE IN X
C
      DO 130 J=1,M
         L=LPERM(J)
         X(J)=DOTPRD(M,A(1,L),RES)
  130    CONTINUE
C
C     SOLVE R(TRANSPOSE) * Y = X  AND WRITE THE SOLUTION IN X
C
      X(1)=X(1)/R(1,1)
      IF (M.EQ.1) GO TO 150
      DO 140 J=2,M
         X(J)=(X(J)-DOTPRD(J-1,R(1,J),X(1)))/R(J,J)
  140    CONTINUE
C
C     COPY X INTO QTB
C
  150 CONTINUE
      DO 160 J=1,M
         QTB(J)=X(J)
  160    CONTINUE
C
C     SOLVE R * Y = X  AND WRITE THE SOLUTION IN X
C
      DO 190 JJ=1,M
         J=M+1-JJ
         TEMP=X(J)/R(J,J)
         X(J)=TEMP
         IF (JJ.EQ.M) GO TO 180
         JM1=J-1
         DO 170 I=1,JM1
            X(I)=X(I)-TEMP*R(I,J)
  170       CONTINUE
  180    CONTINUE
  190    CONTINUE
C
C     REFINE X, QTB
C
      TEMP=QTB(1)
      X(1)=X(1)+TEMPX
      QTB(1)=TEMP+TQTB1
      IF (M.EQ.1) GO TO 210
      TEMP=QTB(2)
      X(2)=X(2)+R(2,1)
      R(2,1)=ZERO
      QTB(2)=TEMP+TQTB2
      IF (M.EQ.2) GO TO 210
      DO 200 I=3,M
         TEMP=QTB(I)
         X(I)=X(I)+R(I,1)
         QTB(I)=TEMP+R(I,2)
         R(I,1)=ZERO
         R(I,2)=ZERO
  200    CONTINUE
  210 CONTINUE
  300 CONTINUE
      RETURN
C
C     LAST CARD OF SOLVE FOLLOWS
C
      END
      SUBROUTINE QRUPDT (LDR, N, R, U, V)
      INTEGER N, LDR
      DOUBLE PRECISION R(LDR, N), U(N), V(N)
C
C     QRUPDT UPDATES R OF QR FACTORIZATION AFTER A RANK ONE CHANGE
C
C     LDR       DECLARED ROW SIZE OF  R
C     N         NO. COLUMNS OF R, SIZE OF U, V
C     R         ON INPUT CONTAINS UPPER TRIANGULAR MATRIX TO BE UPDATED
C               AND UPON OUTPUT CONTAINS UPDATED MATRIX
C     U,V       VECTORS WHOSE OUTER PRODUCT IS TO BE ADDED TO R
C
C
C     UPDATES  R + U * V (TRNS)  BY APPLYING JACOBI ROTATIONS TO ZERO
C        THE LOWER PORTION OF  U * V (TRNS).  THE PRODUCT  U * V (TRNS)
C        DOES NOT HAVE TO BE EXPLICITLY FORMED.  THE JACOBI ROTATIONS
C        ARE, HOWEVER, EXPLICITLY APPLIED TO R.  AFTER THIS PORTION
C        OF THE ALGORITHM IS COMPLETED,  THE  L2-NORM (U) * V (TRNS) IS
C        ADDED TO R.  ANOTHER SEQUENCE OF JACOBI ROTATIONS IS THEN
C        APPLIED TO ZERO THE SUBDIAGONAL OF R, WHICH WAS CREATED DURING
C        THE APPLICATION OF THE FIRST SEQUENCE OF JACOBI ROTATIONS.
C
C     FUNCTIONS USED -- DOTPRD,V2NORM,DSQRT
C
      INTEGER I,IP1,J,K,NM1
      DOUBLE PRECISION X(2),Y(2),TEMP,TEMP2,DSQRT,DOTPRD,V2NORM,ZERO,UL
      DATA ZERO /0.0D0/
C
C     BEGIN PROGRAM
C
      TEMP=U(1)
      IF (N.EQ.1) GO TO 105
C
C     REDUCE U TO L2-NORM(U) * E1
C
      NM1=N-1
      I=N
      UL=U(N)
      TEMP2=UL*UL
      DO 100 K=1,NM1
         IP1=I
         I=I-1
         TEMP2=TEMP2+U(I)*U(I)
         IF (TEMP2.EQ.ZERO) GO TO 90
         TEMP=DSQRT(TEMP2)
         X(1)=U(I)/TEMP
         X(2)=UL/TEMP
         UL=TEMP
         DO 50 J=I,N
            Y(1)=R(IP1,J)
            Y(2)=-R(I,J)
            R(I,J)=DOTPRD(2,X,R(I,J))
            R(IP1,J)=DOTPRD(2,X,Y)
   50       CONTINUE
   90    CONTINUE
  100    CONTINUE
C
C     NOW ADD IN RANK ONE CHANGE TO FIRST ROW OF R
C
      TEMP = V2NORM(N, U)
  105 DO 110 J=1,N
         R(1,J) = R(1,J) + TEMP * V(J)
  110    CONTINUE
      IF (N.EQ.1) GO TO 300
C
C     APPLY JACOBI ROTATIONS TO ZERO OUT SUBDIAGONAL
C
      DO 200 I=1,NM1
         IP1=I+1
         TEMP=V2NORM(2,R(I,I))
         IF (TEMP .EQ. ZERO) GO TO 190
         X(1)=R(I,I)/TEMP
         X(2)=R(IP1,I)/TEMP
         R(I,I)=TEMP
         R(IP1,I)=ZERO
         DO 180 J=IP1,N
            Y(1)=R(IP1,J)
            Y(2)=-R(I,J)
            R(I,J)=DOTPRD(2,X,R(I,J))
            R(IP1,J)=DOTPRD(2,X,Y)
  180       CONTINUE
  190    CONTINUE
  200    CONTINUE
  300 RETURN
C
C     LAST CARD OF QRUPDT FOLLOWS
C
      END
      DOUBLE PRECISION FUNCTION DOTPRD(N,X,Y)
C
C     COMPUTES THE INNER PRODUCT OF THE N-VECTORS X AND Y
C
      INTEGER N
      DOUBLE PRECISION X(N),Y(N)
      INTEGER I
      DOTPRD=0.D0
      DO 10 I=1,N
         DOTPRD=DOTPRD+X(I)*Y(I)
   10    CONTINUE
      RETURN
C
C     LAST CARD OF FUNCTION DOTPRD FOLLOWS
C
      END
      DOUBLE PRECISION FUNCTION V2NORM(N,X)
C
C     COMPUTES THE 2-NORM OF THE N-VECTOR X
C     AVOIDS OVERFLOW AND UNDERFLOW ON MACHINES WITH
C     UNDERFLOW LIMIT .LT. MACHEPS**2.
C
      INTEGER N
      DOUBLE PRECISION X(N)
      INTEGER I
      DOUBLE PRECISION XMAX,SUM,TEMP,ZERO
      DOUBLE PRECISION DABS,DMAX1,DSQRT,DD,EE
      DATA ZERO /0.0D0/
      XMAX=ZERO
      V2NORM=ZERO
      DO 10 I=1,N
         DD=DABS(X(I))
         EE=XMAX
         XMAX=DMAX1(EE,DD)
   10    CONTINUE
      IF (XMAX.LE.ZERO) GO TO 40
      SUM=ZERO
      DO 30 I=1,N
         IF ((XMAX+DABS(X(I))).EQ.XMAX) GO TO 20
            TEMP=X(I)/XMAX
            SUM=SUM+TEMP*TEMP
   20       CONTINUE
   30    CONTINUE
      V2NORM=XMAX*DSQRT(SUM)
   40 CONTINUE
      RETURN
C
C     LAST CARD OF FUNCTION V2NORM FOLLOWS
C
      END
      DOUBLE PRECISION FUNCTION VSNORM(N,X)
      INTEGER N
      DOUBLE PRECISION X(N)
C
C     RETURNS THE SUP-NORM OF THE N-VECTOR X
C
      INTEGER J
      DOUBLE PRECISION DABS,DMAX1,DD,EE
      VSNORM=0.0D0
      DO 10 J=1,N
         DD=DABS(X(J))
         EE=VSNORM
         VSNORM=DMAX1(EE,DD)
   10    CONTINUE
      RETURN
      END