C     SOLUTION OF THE MATRIX EQUATION AX + XB = C.
C
C     THIS IS A SAMPLE MAIN PROGRAM WHICH SOLVES RANDOMLY GENERATED
C     PROBLEMS.
C
C     A IS REDUCED TO SCHUR FORM AND B TO UPPER-HESSENBERG FORM.
C     THE RESULTING BLOCK-SCHUR SYSTEM OF EQUATIONS IS
C     SOLVED AND THE SOLUTION TRANSFORMED BACK TO THAT OF THE
C     ORIGINAL EQUATION.
C
C     APRIL 9, 1978.  WRITTEN BY STEPHEN NASH
C
C
      REAL XX(10),YY(10),ZZ(10)
      REAL*8 EPS,ORT(10),A(10,10),B(10,10),C(10,10),V(10,10)
     *      ,AA(10,10),BB(10,10),TT(10,10),DABS,T(10,10),DMAX,WORK(5310)
c      DATA EPS/Z3410000000000000/
C
C     WITHIN THE MAIN PROGRAM, A,B AND C ARE STORED NORMALLY.
C     AT THE END OF THE COMPUTATION C HOLDS THE SOLUTION X. THE ARRAY
C     WORK IS USED FOR TEMPORARY STORAGE WITHIN THE SUBROUTINES.
C     V IS USED TO STORE THE TRANSFORMATION MATRIX FOR A.
C     THE TRANSFORMATION MATRIX FOR B IS STORED IN ORT AND B.
C     C AND X ARE M * N MATRICES. A IS M * M, B IS N * N.
C
      INTEGER IPR(200)
      eps = d1mach(4)
      NDIM = 10
      MDIM = 10
      MODE = 1
      ISEED = 14365725
      MAXIIJ = 2
      DO 90 IIJJ = 1,MAXIIJ
C
C     N AND M ARE RANDOMLY GENERATED IN THE RANGE U10,50e.
C     AFTERWARDS, THE MATRICES A,B, AND X (STORED IN T) ARE
C     ALSO GENERATED RANDOMLY WITH ENTRIES IN THE INTERVAL U0,1).
C     ISEED IS A SEED FOR THE RANDOM NUMBER GENERATOR.
C     IF IIJJ=1 WE ARE GENERATING A NEW PROBLEM.  WHEN
C     IIJJ=2, THE RIGHT HAND SIDE IS CHANGED, BUT THE
C     COEFFICIENT MATRICES ARE LEFT UNCHANGED.
C
      CALL RANDK(ISEED,XX(1),0)
      CALL RANDK(ISEED,XX(2),0)
      M = 10
      N = 5
      IF (IIJJ .EQ. 2) MODE = 2
      NM = MAX0(N,M)
      WRITE(6,10) M,N
10    FORMAT('-SOLUTION OF',I3,' BY',I3,' MATRIX EQUATION')
      DO 20 I = 1,NM
         DO 15 K = 1,NM
            CALL RANDK(ISEED,XX(K),0)
            CALL RANDK(ISEED,YY(K),0)
            CALL RANDK(ISEED,ZZ(K),0)
15       CONTINUE
         DO 20 J = 1,NM
            IF (IIJJ .EQ. 2) GO TO 17
            A(I,J) = XX(J)
            B(I,J) = YY(J)
            AA(I,J) = A(I,J)
            BB(I,J) = B(I,J)
17          T(I,J) = ZZ(J)
20    CONTINUE
      DO 30 I = 1,NM
      DO 30 J = 1,NM
         TT(I,J) = T(I,J)
30    CONTINUE
C
C     GENERATE C FROM A, B, AND X
C
      CALL SLVCHK(AA,BB,T,C,NDIM,N,MDIM,M)
C
C     NOW SOLVE THE MATRIX EQUATION
C
      CALL AXXBC(MODE,NDIM,N,MDIM,M,A,V,B,C,ORT,WORK,IPR,EPS,IERR)
      IF (IERR.NE.0) GO TO 60
C
C     CHECK THE SIZE OF THE ERROR.
C
      DMAX = 0.D0
      DO 40 I = 1,M
      DO 40 J = 1,N
         IF (DMAX.LT.DABS(C(I,J)-TT(I,J))) DMAX = DABS(C(I,J)-TT(I,J))
40    CONTINUE
      WRITE(6,50) DMAX
50    FORMAT(' MAXIMAL ERROR = ', D14.6)
      GO TO 90
60    WRITE(6,70)
70    FORMAT('0SINGULAR SYSTEM OF EQUATIONS')
90    CONTINUE
      STOP
      END
C
C
C
C
C
      SUBROUTINE RANDK (IY, YFL, INDEX)
C
C     THIS IS A UNIFORM RANDOM NUMBER GENERATOR WRITTEN BY G. E.
C  FORSYTHE IN SPRING 1969, FOLLOWING D. KNUTH, THE ART OF COMPUTER
C  PROGRAMMING, VOL. 2, PP. 155-156.  IT IS MUCH SUPERIOR TO RANDU,
C  THE RANDOM NUMBER GENERATOR FOUND IN IBM'S SCIENTIFIC SUBROUTINE
C  PACKAGE.
C
C     BEFORE THE FIRST CALL OF RANDK, IY SHOULD BE SET OUTSIDE RANDK
C  TO AN ARBITRARY INTEGER VALUE.  (IN WATFOR THIS IS ESSENTIAL.)
C  FOR PROGRAM CHECKOUT, THE INITIAL VALUE OF IY SHOULD BE A FIXED
C  INTEGER.  FOR RANDOM NUMBERS DIFFERENT ON EVERY RUN (AND HENCE
C  NOT REPRODUCIBLE), DECLARE INTEGER CLOCK1 AND THEN INITIALIZE
C  IY TO CLOCK1(4).
C
C     IF RANDK IS CALLED WITH AN INTEGER INDEX = 1, THEN THE OUTPUT
C  VALUE OF IY IS A PSEUDORANDOM INTEGER UNIFORMLY DISTRIBUTED IN THE
C  RANGE 0 <= IY < 2**31.
C
C     IF RANDK IS CALLED WITH INDEX = 0, THEN NOT ONLY IS IY PRODUCED,
C  BUT ALSO (AT SOME EXTRA COST IN TIME) A FLOATING NUMBER YFL, UNI-
C  FORMLY DISTRIBUTED IN THE INTERVAL 0.0 <= YFL < 1.0.
C
      IY = IY*314159269 + 453806245
4     IF (IY .GE. 0) GO TO 6
C
C  CAUTION:  THE STATEMENT LABEL 4 IS ESSENTIAL IN ORDER TO PREVENT
C  CERTAIN COMPILERS (E.G., FORTRAN H WITH OPT 0) FROM PERFORMING
C  UNWANTED "OPTIMIZATIONS."  IT SHOULD NOT BE REMOVED.
C
5        IY = IY + 2147483647 + 1
C        STATEMENT 5 ADDS 2**31 TO NEGATIVE VALUES OF IY
C
6     IF (INDEX) 7, 7, 8
C
7        YFL = IY
         YFL = YFL*.4656613E-9
C
8     RETURN
      END
C
C
C
C
C
      SUBROUTINE SLVCHK(A,B,X,C,NDIM,N,MDIM,M)
C
      INTEGER I,J,K,M,MDIM,N,NDIM
      REAL*8 A(MDIM,M),B(NDIM,N),X(MDIM,N),C(MDIM,N)
C
C     THIS SUBROUTINE FORMS AX + XB IN ORDER TO CHECK THE SOLUTION
C     TO THE PROBLEM.  THE RESULT IS STORED IN THE MATRIX C.
C     PARAMETERS:
C     A,B,X - INPUT MATRICES
C     C     - OUTPUT MATRIX
C     NDIM  - DECLARED ROW DIMENSION OF THE MATRIX B
C     N     - ACTUAL ROW DIMENSION OF B AND COLUMN DIMENSION
C             OF B, C, AND X
C     MDIM  - DECLARED ROW DIMENSION OF THE MATRICES A, X, AND C
C     M     - ACTUAL COLUMN DIMESION OF A AND ROW DIMENSION OF
C             A, X, AND C.
C
      DO 10 I = 1,M
         DO 10 J = 1,N
            C(I,J) = 0.D0
            DO 10 K = 1,M
               C(I,J) = C(I,J) + A(I,K) * X(K,J)
10    CONTINUE
      DO 20 I = 1,M
         DO 20 J = 1,N
            DO 20 K = 1,N
               C(I,J) = C(I,J) + X(I,K) * B(K,J)
20    CONTINUE
      RETURN
      END
C
C------------------------------------------------------------------------
C
      SUBROUTINE AXXBC(MODE,NDIM,N,MDIM,M,A,V,B,C,ORT,D,IPR,EPS,IERR)
C
      INTEGER I,IERR,IND,IPR(1),J,M,MDIM,MODE,N,NDIM
      REAL*8 A(MDIM,M),B(NDIM,N),V(NDIM,N),
     *     C(MDIM,N),ORT(1),EPS,REPS,D(1)
C
C     THIS ROUTINE SOLVES THE MATRIX EQUATION AX + XB = C WHERE
C
C        A - M * M MATRIX
C        X - M * N MATRIX
C        B - N * N MATRIX
C        C - M * N MATRIX.
C
C     THE PARAMETERS FOR THIS SUBROUTINE ARE:
C
C     MODE - INDICATES WHETHER A NEW MATRIX EQUATION IS BEING SOLVED,
C            OR THAT JUST THE RIGHT HAND SIDE HAS BEEN CHANGED.
C            MODE=1 IMPLIES THAT THE EQUATION IS NEW
C            MODE=2 IMPLIES THAT ONLY THE RIGHT HAND SIDE IS NEW
C     NDIM - DECLARED ROW DIMENSION OF THE MATRIX B
C     N    - ROW DIMENSION OF B; COLUMN DIMENSION OF B AND C
C     MDIM - DECLARED ROW DIMENSION OF A, C, AND V
C     M    - ROW DIMENSION OF A, C, AND V; COLUMN DIMENSION OF A AND V
C     A    - M * M MATRIX (DOUBLE PRECISION)
C     B    - N * N MATRIX (DOUBLE PRECISION)
C     V    - N * N MATRIX (DOUBLE PRECISION) USED TO STORE THE MATRIX
C            WHICH TRANSFORMS B(T) TO UPPER SCHUR FORM
C     C    - M * N MATRIX WHICH HOLDS C ON INPUT AND X ON OUTPUT
C     ORT  - DOUBLE PRECISION VECTOR USED BY SUBROUTINE RG
C            OF LENGTH AT LEAST MAX(N,M)
C     D    - DOUBLE PRECISION VECTOR FOR INTERNAL USE OF LENGTH
C            AT LEAST 2M**2 + 7M
C     IPR  - INTEGER VECTOR FOR INTERNAL USE OF LENGTH AT LEAST 4N
C     EPS  - ERROR TOLERANCE. EPS SHOULD EQUAL THE SMALLEST NUMBER
C            WHICH SATISFIES FL(1 + EPS) > 1.
C     IERR - INTEGER VARIABLE USED TO SET ERROR CODES.
C            IERR=0 => PROPER RETURN.
C            IERR=J => J-TH EIGENVALUE OF B HAS NOT BEEN DETERMINED
C                      AFTER 30 QR ITERATIONS.
C            IERR=-J => A SINGULAR MATRIX WAS ENCOUNTERED WHEN
C                       SOLVING FOR THE J-TH COLUMN OF X.
C
C     METHOD:
C
C     FIRST OF ALL, B(TRANSPOSE) IS FORMED AND IS USED INSTEAD OF
C     B THROUGHOUT THE PROGRAM. (IN THE FOLLOWING DISCUSSION, B(T) IS
C     USED TO REPRESENT B(TRANSPOSE)).
C     B(T) IS TRANSFORMED TO REAL UPPER SCHUR FORM AND A TO UPPER
C     HESSENBERG FORM USING ORTHOGONAL TRANSFORMATIONS.
C     THE RIGHT HAND SIDE C IS MULTIPLIED BY THESE TRANSFORMATION
C     MATRICES AND THE SOLUTION OF THIS TRANSFORMED SYSTEM IS
C     COMPUTED.  THIS SOLUTION IS THEN MULTIPLIED BY THE TRANS-
C     FORMATION MATRICES IN ORDER TO OBTAIN THE SOLUTION TO THE
C     ORIGINAL PROBLEM.
C
C     NOTE:
C
C     IN ORDER TO GET OPTIMAL EFFICIENCY FROM THIS ROUTINE, N SHOULD
C     BE <= M.  IF THIS IS NOT TRUE FOR YOUR PROBLEM, TRANSPOSE THE
C     EQUATION AND SOLVE THE RESULTING PROBLEM.
C
C     -------------------------------------------------------------
C
C     FORM B(TRANSPOSE) FOR INTERNAL USE
C
      IF (MODE .EQ. 2) GO TO 7
      DO 5 I = 1,N
         DO 5 J = I,N
            D(1) = B(I,J)
            B(I,J) = B(J,I)
            B(J,I) = D(1)
5     CONTINUE
C
C     A AND B(T) WILL BE TRANSFORMED INTO THE APPROPRIATE FORMS USING
C     MODIFICATIONS OF EISPACK ROUTINES (MODIFIED BY S. NASH).
C
C     WARNING:  DO NOT CHANGE THE ORDER OF THE NEXT TWO STATEMENTS
C
      CALL RG(NDIM,N,B,1,V,ORT,EPS,IERR)
      CALL RG(MDIM,M,A,0,A,ORT,EPS,IERR)
C
C     NOW, TRANSFORM THE RIGHT HAND SIDE.
C
7     CALL TRANSF(A,ORT,1,C,V,0,M,N,MDIM,NDIM,D)
C
C     NOW SOLVE THE SYSTEM OF EQUATIONS (BLOCK BACK SUBSTITUTION)
C     THE SYSTEM IS BLOCK UPPER HESSENBERG WITH UPPER HESSENBERG
C     MATRICES ON THE DIAGONAL (A + B(I,I) * I) AND MULTIPLES OF
C     THE IDENTITY OFF THE DIAGONAL.
C     TESTS ARE MADE TO SEE WHETHER A SINGLE OR A DOUBLE BLOCK
C     MUST BE HANDLED AT THE CURRENT STAGE AND APPROPRIATE
C     SUBROUTINES ARE CALLED.
C
C     SET RELATIVE ERROR TOLERANCE
C
      REPS = EPS*N*N*M*M
      IND = N - 1
      IF (IND.EQ.0) GO TO 40
10    IF (DABS(B(IND + 1,IND)).LE.REPS)  GO TO 20
      CALL N2SOLV(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR)
      IF (IERR.NE.0) RETURN
      GO TO 30
20    CALL NSOLVE(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR)
      IF (IERR.NE.0) RETURN
30    IF (IND.GT.0) GO TO 10
40    IF (IND.EQ.0) CALL NSOLVE(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR)
C
C     THE SOLUTION OF THE TRANSFORMED SYSTEM OF EQUATIONS
C     HAS NOW BEEN OBTAINED.  IT IS TRANSFORMED BACK INTO
C     THE SOLUTION OF THE ORIGINAL SYSTEM OF EQUATIONS.
C
      CALL TRANSF(A,ORT,0,C,V,1,M,N,MDIM,NDIM,D)
      RETURN
      END
C
C
C
C
C
      SUBROUTINE RG(NDIM,N,A,MATZ,Z,ORT,EPS,IERR)
C
      INTEGER N,NDIM,IERR,MATZ
      REAL*8 A(NDIM,N),Z(NDIM,N),ORT(N),EPS
C
C     THIS SUBROUTINE REDUCES A TO EITHER UPPER-HESSENBERG
C     OR SCHUR FORM DEPENDING ON WHETHER MATZ IS ZERO OR NOT.
C     THE ROUTINES USED ARE FROM EISPACK EXCEPT THAT HQR2 HAS
C     BEEN MODIFIED SO AS NOT TO COMPUTE THE EIGENVECTORS OF
C     A AND ONLY REDUCES A TO SCHUR FORM. ORTHOGONAL TRANSFORMATIONS
C     ARE USED THROUGHOUT.
C     PARAMETERS:
C     NDIM - DECLARED ROW DIMENSION OF THE MATRIX A
C     N    - ACTUAL DIMENSION OF THE MATRIX A
C     A    - MATRIX TO BE TRANSFORMED (ON INPUT) AND TRANSFORMED
C            MATRIX (ON OUTPUT)
C     MATZ - INTEGER VARIABLE; MATZ = 0 => TRANSFORM A TO UPPER-
C            HESSENBERG FORM; OTHERWISE COMPUTE THE SCHUR FORM OF A
C     Z    - STORES THE TRANSFORMATION MATRIX IF A IS TRANSFORMED TO
C          - SCHUR FORM (IF MATZ=0, Z CAN BE THE SAME AS A).
C     ORT  - VECTOR USED FOR INTERNAL STORAGE; IF MATZ=0, ORT HOLDS
C            INFORMATION ABOUT THE TRANSFORMATION MATRIX
C     EPS  - ERROR TOLERANCE, AS IN SUBROUTINE AXXBC.
C     IERR - ERROR CODE; SEE ROUTINE AXXBC FOR DETAILS
C
      IERR=0
      CALL ORTHES(NDIM,N,A,ORT)
      IF (MATZ.EQ.0) RETURN
      CALL ORTRAN(NDIM,N,A,ORT,Z)
      CALL HQR2(NDIM,N,A,Z,EPS,IERR)
10    RETURN
      END
C
C
C
C
C
      SUBROUTINE ORTHES(NDIM,N,A,ORT)
C
      INTEGER I,J,M,N,II,JJ,LA,MP,NDIM,KP1
      REAL*8 A(NDIM,N),ORT(N)
      REAL*8 F,G,H,SCALE
      REAL*8 DSQRT,DABS,DSIGN
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES,
C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
C
C     GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE
C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
C     1 THROUGH N TO UPPER HESSENBERG FORM BY
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT:
C
C        NDIM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT;
C
C        N IS THE ORDER OF THE MATRIX;
C
C        A CONTAINS THE INPUT MATRIX.
C
C     ON OUTPUT:
C
C        A CONTAINS THE HESSENBERG MATRIX.  INFORMATION ABOUT
C          THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION
C          IS STORED IN THE REMAINING TRIANGLE UNDER THE
C          HESSENBERG MATRIX;
C
C        ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
C          ONLY ELEMENTS 1 THROUGH N ARE USED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     -----------------------------------------------------------------
C
      LA = N - 1
      KP1 = 2
      IF (LA .LT. KP1) GO TO 200
C
      DO 180 M = KP1, LA
         H = 0.0D0
         ORT(M) = 0.0D0
         SCALE = 0.0D0
C
C     SCALE COLUMN (ALGOL TOL THEN NOT NEEDED)
C
         DO 90 I = M, N
90       SCALE = SCALE + DABS(A(I,M-1))
C
         IF (SCALE .EQ. 0.0D0) GO TO 180
         MP = M + N
C
C        FOR I=N STEP -1 UNTIL M DO --
C
         DO 100 II = M, N
            I = MP - II
            ORT(I) = A(I,M-1) / SCALE
            H = H + ORT(I) * ORT(I)
100      CONTINUE
C
         G = -DSIGN(DSQRT(H),ORT(M))
         H = H - ORT(M) * G
         ORT(M) = ORT(M) - G
C
C        FORM (I-(U*UT)/H) * A
C
         DO 130 J = M, N
            F = 0.0D0
C
C           FOR I=N STEP -1 UNTIL M DO --
C
            DO 110 II = M, N
               I = MP - II
               F = F + ORT(I) * A(I,J)
110         CONTINUE
C
            F = F / H
C
            DO 120 I = M, N
120         A(I,J) = A(I,J) - F * ORT(I)
C
130      CONTINUE
C
C        FORM (I-(U*UT)/H)*A*(I-(U*UT)/H)
C
         DO 160 I = 1, N
            F = 0.0D0
C
C           FOR J=N STEP -1 UNTIL M DO --
C
            DO 140 JJ = M, N
               J = MP - JJ
               F = F + ORT(J) * A(I,J)
140         CONTINUE
C
            F = F / H
C
            DO 150 J = M, N
150         A(I,J) = A(I,J) - F * ORT(J)
C
160      CONTINUE
C
         ORT(M) = SCALE * ORT(M)
         A(M,M-1) = SCALE * G
180   CONTINUE
200   RETURN
      END
C
C
C
C
C
      SUBROUTINE ORTRAN(NDIM,N,A,ORT,Z)
C
      INTEGER I,J,N,KL,MM,MP,NDIM,MP1
      REAL*8 A(NDIM,N),ORT(N),Z(NDIM,N)
      REAL*8 G
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTRANS,
C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C
C     THIS SUBROUTINE ACCUMULATES THE ORTHOGONAL SIMILARITY
C     TRANSFORMATIONS USED IN THE REDUCTION OF A REAL GENERAL
C     MATRIX TO UPPER HESSENBERG FORM BY  ORTHES.
C
C     ON INPUT:
C
C        NDIM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT;
C
C        N IS THE ORDER OF THE MATRIX;
C
C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
C          FORMATIONS USED IN THE REDUCTION BY  ORTHES
C          IN ITS STRICT LOWER TRIANGLE;
C
C        ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS-
C          FORMATIONS USED IN THE REDUCTION BY  ORTHES.
C          ONLY ELEMENTS 1 THROUGH N ARE USED.
C
C     ON OUTPUT:
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C          REDUCTION BY  ORTHES;
C
C        ORT HAS BEEN ALTERED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     -----------------------------------------------------------
C
C     INITIALIZE Z TO IDENTITY MATRIX
C
      DO 80 I = 1, N
C
         DO 60 J = 1, N
60       Z(I,J) = 0.0D0
C
         Z(I,I) = 1.0D0
80    CONTINUE
C
      KL = N - 2
      IF (KL .LT. 1) GO TO 200
C
C     FOR MP=N-1 STEP -1 UNTIL 2 DO --
C
      DO 140 MM = 1, KL
         MP = N - MM
         IF (A(MP,MP-1) .EQ. 0.0D0) GO TO 140
         MP1 = MP + 1
C
         DO 100 I = MP1, N
100      ORT(I) = A(I,MP-1)
C
         DO 130 J = MP, N
            G = 0.0D0
C
            DO 110 I = MP, N
110         G = G + ORT(I) * Z(I,J)
C
C           DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES.
C           DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW
C
            G = (G / ORT(MP)) / A(MP,MP-1)
C
            DO 120 I = MP, N
120         Z(I,J) = Z(I,J) + G * ORT(I)
C
130      CONTINUE
C
140   CONTINUE
C
200   RETURN
      END
C
C
C
C
C
      SUBROUTINE HQR2(NDIM,N,H,Z,EPS,IERR)
C
      INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NDIM,NN,
     X        ITS,MP2,ENM2,IERR
      REAL*8 H(NDIM,N),Z(NDIM,N)
      REAL*8 P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,EPS
      REAL*8 DSQRT,DABS,DSIGN
      INTEGER MIN0
      LOGICAL NOTLAS
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2,
C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C     IT HAS BEEN FURTHER MODIFIED BY STEPHEN NASH SO THAT IT ONLY
C     COMPUTES THE SCHUR FORM OF THE UPPER-HESSENBERG MATRIX H.
C
C     ON INPUT:
C
C        NDIM MUST BE SET TO THE ROW DIMENSION OF TWO - DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT;
C
C        N IS THE ORDER OF THE MATRIX;
C
C        H CONTAINS THE UPPER HESSENBERG MATRIX;
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ORTRAN
C          AFTER THE REDUCTION BY ORTHES.
C
C     ON OUTPUT:
C
C        H HAS BEEN DESTROYED;
C
C        Z CONTAINS THE TRANSFORMATION MATRIX.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     -----------------------------------------------------------------
C
C     EPS IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C     THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
      IERR = 0
      NORM = 0.0D0
      K = 1
C
C     COMPUTE MATRIX NORM
C
      DO 50 I = 1, N
C
         DO 40 J = K, N
40       NORM = NORM + DABS(H(I,J))
C
         K = I
50    CONTINUE
C
      EN = N
      T = 0.0D0
C
C     SEARCH FOR NEXT EIGENVALUES
C
60    IF (EN .LT. 1) GO TO 340
      ITS = 0
      NA = EN - 1
      ENM2 = NA - 1
C
C     LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C     FOR L=EN STEP -1 UNTIL 1 DO --
C
70    DO 80 LL = 1, EN
         L = EN + 1 - LL
         IF (L .EQ. 1) GO TO 100
         S = DABS(H(L-1,L-1)) + DABS(H(L,L))
         IF (S .EQ. 0.0D0) S = NORM
         IF (DABS(H(L,L-1)) .LE. EPS * S) GO TO 100
80    CONTINUE
C
C     FORM SHIFT
C
100   X = H(EN,EN)
      IF (L .EQ. EN) GO TO 270
      Y = H(NA,NA)
      W = H(EN,NA) * H(NA,EN)
      IF (L .EQ. NA) GO TO 280
      IF (ITS .EQ. 30) GO TO 1000
      IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
C
C     FORM EXCEPTIONAL SHIFT
C
      T = T + X
C
      DO 120 I = 1, EN
120   H(I,I) = H(I,I) - X
C
      S = DABS(H(EN,NA)) + DABS(H(NA,ENM2))
      X = 0.75D0 * S
      Y = X
      W = -0.4375D0 * S * S
130   ITS = ITS + 1
C
C     LOOK FOR TWO CONSECUTIVE SMALL
C     SUB-DIAGONAL ELEMENTS.
C     FOR M=EN-2 STEP -1 UNTIL L DO --
C
      DO 140 MM = L, ENM2
         M = ENM2 + L - MM
         ZZ = H(M,M)
         R = X - ZZ
         S = Y - ZZ
         P = (R * S - W) / H(M+1,M) + H(M,M+1)
         Q = H(M+1,M+1) - ZZ - R - S
         R = H(M+2,M+1)
         S = DABS(P) + DABS(Q) + DABS(R)
         P = P / S
         Q = Q / S
         R = R / S
         IF (M .EQ. L) GO TO 150
         IF (DABS(H(M,M-1)) * (DABS(Q) + DABS(R)) .LE. EPS * DABS(P)
     X    * (DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))) GO TO 150
140   CONTINUE
C
150   MP2 = M + 2
C
      DO 160 I = MP2, EN
         H(I,I-2) = 0.0D0
         IF (I .EQ. MP2) GO TO 160
         H(I,I-3) = 0.0D0
160   CONTINUE
C
C     DOUBLE QR STEP INVOLVING ROWS L TO EN AND
C     COLUMNS M TO EN
C
      DO 260 K = M, NA
         NOTLAS = K .NE. NA
         IF (K .EQ. M) GO TO 170
         P = H(K,K-1)
         Q = H(K+1,K-1)
         R = 0.0D0
         IF (NOTLAS) R = H(K+2,K-1)
         X = DABS(P) + DABS(Q) + DABS(R)
         IF (X .EQ. 0.0D0) GO TO 260
         P = P / X
         Q = Q / X
         R = R / X
170      S = DSIGN(DSQRT(P*P+Q*Q+R*R),P)
         IF (K .EQ. M) GO TO 180
         H(K,K-1) = -S * X
         GO TO 190
180      IF (L .NE. M) H(K,K-1) = -H(K,K-1)
190      P = P + S
         X = P / S
         Y = Q / S
         ZZ = R / S
         Q = Q / P
         R = R / P
C
C        ROW MODIFICATION
C
         DO 210 J = K, N
            P = H(K,J) + Q * H(K+1,J)
            IF (.NOT. NOTLAS) GO TO 200
            P = P + R * H(K+2,J)
            H(K+2,J) = H(K+2,J) - P * ZZ
200         H(K+1,J) = H(K+1,J) - P * Y
            H(K,J) = H(K,J) - P * X
210      CONTINUE
C
         J = MIN0(EN,K+3)
C
C        COLUMN MODIFICATION
C
         DO 230 I = 1, J
            P = X * H(I,K) + Y * H(I,K+1)
            IF (.NOT. NOTLAS) GO TO 220
            P = P + ZZ * H(I,K+2)
            H(I,K+2) = H(I,K+2) - P * R
220         H(I,K+1) = H(I,K+1) - P * Q
            H(I,K) = H(I,K) - P
230      CONTINUE
C
C        ACCUMULATE TRANSFORMATIONS
C
         DO 250 I = 1, N
            P = X * Z(I,K) + Y * Z(I,K+1)
            IF (.NOT. NOTLAS) GO TO 240
            P = P + ZZ * Z(I,K+2)
            Z(I,K+2) = Z(I,K+2) - P * R
240         Z(I,K+1) = Z(I,K+1) - P * Q
            Z(I,K) = Z(I,K) - P
250      CONTINUE
C
260   CONTINUE
C
      GO TO 70
C
C     ONE ROOT FOUND
C
270   H(EN,EN) = X + T
      EN = NA
      GO TO 60
C
C     TWO ROOTS FOUND
C
280   P = (Y - X) / 2.0D0
      Q = P * P + W
      ZZ = DSQRT(DABS(Q))
      H(EN,EN) = X + T
      X = H(EN,EN)
      H(NA,NA) = Y + T
      IF (Q .LT. 0.0D0) GO TO 320
C
C     REAL PAIR
C
      ZZ = P + DSIGN(ZZ,P)
      X = H(EN,NA)
      S = DABS(X) + DABS(ZZ)
      P = X / S
      Q = ZZ / S
      R = DSQRT(P*P+Q*Q)
      P = P / R
      Q = Q / R
C
C     ROW MODIFICATION
C
      DO 290 J = NA, N
         ZZ = H(NA,J)
         H(NA,J) = Q * ZZ + P * H(EN,J)
         H(EN,J) = Q * H(EN,J) - P * ZZ
290   CONTINUE
C
C     COLUMN MODIFICATION
C
      DO 300 I = 1, EN
         ZZ = H(I,NA)
         H(I,NA) = Q * ZZ + P * H(I,EN)
         H(I,EN) = Q * H(I,EN) - P * ZZ
300   CONTINUE
C
C     ACCUMULATE TRANSFORMATIONS
C
      DO 310 I = 1, N
         ZZ = Z(I,NA)
         Z(I,NA) = Q * ZZ + P * Z(I,EN)
         Z(I,EN) = Q * Z(I,EN) - P * ZZ
310   CONTINUE
C
      GO TO 330
C
C     COMPLEX PAIR
C
320   CONTINUE
330   EN = ENM2
      GO TO 60
C
C     ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
C     VECTORS OF UPPER TRIANGULAR FORM
C
1000  IERR = EN
340   RETURN
      END
C
C
C
C
C
      SUBROUTINE TRANSF(A,ORT,IT1,C,V,IT2,M,N,MDIM,NDIM,D)
C
      INTEGER I,IT1,IT2,J,K,K1,K2,KK,M,MDIM,N,N2
      REAL*8 V(NDIM,N),C(MDIM,N),A(MDIM,M),ORT(M),D(1),G
C
C     SUBROUTINE TO FORM C(NEW) = U * C * V.  THE INTEGER VARIABLES
C     IT1 AND IT2 INDICATE WHETHER U(TRANSPOSE) OR V(TRANSPOSE)
C     ARE REQUIRED INSTEAD OF U OR V.
C     IF ITX = 0 THEN THE ACTUAL MATRIX IS USED, IF ITX = 1, THEN
C     ITS TRANSPOSE IS USED.
C     IN ORDER TO SAVE STORAGE, THE MATRIX U (THE TRANSFORMATION
C     MATRIX FOR THE MATRIX A) IS STORED IN THE LOWER PART OF A
C     AND IN THE VECTOR ORT IN CODED FORM. (SEE SUBROUTINE ORTHES)
C     THE PARAMETERS FOR THIS SUBROUTINE ARE AS FOR THE SUBROUTINE
C     AXXBC EXCEPT THAT HERE ORT TAKES THE PLACE OF ORT.
C
C     --------------------------------------------------------------
C
C     FORM U * C AND STORE IN C (AS IN SUBROUTINE ORTRAN)
C
      M2 = M - 2
C
C     CHECK FOR SMALL MATRIX
C
      IF (M2.LE.0) GO TO 45
      DO 40 KK = 1,M2
         K = M2 - KK + 1
         IF (IT1.EQ.1) K = KK
         K1 = K + 1
         IF (A(K1,K).EQ.0.D0) GO TO 40
         D(K1) = ORT(K1)
         K2 = K + 2
         DO 10 I = K2,M
            D(I) = A(I,K)
10       CONTINUE
         DO 30 J = 1,N
            G = 0.D0
            DO 20 I = K1,M
               G = G + D(I) * C(I,J)
20          CONTINUE
            G = G / ORT(K1) / A(K1,K)
            DO 30 I = K1,M
            C(I,J) = C(I,J) + G * D(I)
30       CONTINUE
40    CONTINUE
C
C     FORM C(NEW) = (U * C) * V
C
45    DO 60 I = 1,M
         DO 50 J = 1,N
            D(J) = 0.D0
            DO 50 K = 1,N
               IF (IT2.EQ.0) D(J) = D(J) + C(I,K) * V(K,J)
               IF (IT2.EQ.1) D(J) = D(J) + C(I,K) * V(J,K)
50          CONTINUE
            DO 60 J = 1,N
               C(I,J) = D(J)
60    CONTINUE
      RETURN
      END
C
C
C
C
C
      SUBROUTINE NSOLVE(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR)
C
      INTEGER I,I1,IERR,IND,IPR(1),IROW1,J,M,M1,MDIM,N,NDIM,MFIN
      REAL*8 A(MDIM,M),B(NDIM,N),C(MDIM,N),D(1),REPS
C
C     THIS SUBROUTINE SETS UP AND SOLVES (WITH THE HELP OF
C     THE SUBROUTINE HESOLV) A SINGLE SYSTEM OF EQUATIONS
C
C      ( A + B(IND+1,IND+1) * I )(X(IND+1)) = (C(IND+1))
C
C     (X(I),C(I) - COLUMN I OF THE MATRICES X AND C)
C     AND STORES THE RESULT IN THE APPROPRIATE ROW OF THE
C     MATRIX C.  THE UPPER-HESSENBERG MATRIX IS STORED IN THE
C     VECTOR D IN ORDER TO SAVE STORAGE (ENTRIES BELOW THE SUB-
C     DIAGONAL ARE IGNORED).
C     THE RIGHT HAND SIDE IS ALSO STORED IN THE VECTOR D IN ORDER
C     TO REDUCE THE NUMBER OF VECTORS USED IN THE SUBROUTINES.
C     THE PARAMETERS ARE AS IN THE SUBROUTINE AXXBC.
C
C     ----------------------------------------------------------
C
C     PERFORM BLOCK BACK-SUBSTITUTION IF NECESSARY
C
      IF (IND.LT.N - 1) CALL BACKSB(C,B,IND,N,M,MDIM,NDIM)
C
C     SET UP THE SYSTEM OF EQUATIONS
C
      MFIN = (M * (M + 1)) / 2 + M
      DO 20 I = 1,M
C
C        FIND INDEX OF BEGINNING OF ROW I WITHIN THE VECTOR D
C
         M1 = IROW1(I,M)
         I1 = I - 1
         IF (I.EQ.1) I1 = 1
         DO 10 J = I1,M
            D(M1+J-I1+1) = A(I,J)
10       CONTINUE
         J = 2
         IF (I.EQ.1) J = 1
         D(M1+J) = D(M1+J) + B(IND+1,IND+1)
C
C        STORE THE RIGHT HAND SIDE
C
         D(MFIN+I) = C(I,IND+1)
20    CONTINUE
C
C     SOLVE THE SYSTEM OF EQUATIONS AND STORE THE RESULT IN C.
C
      CALL HESOLV(D,IPR,M,REPS,IERR)
      IF (IERR.NE.0) GO TO 40
      DO 30 I = 1,M
         C(I,IND+1) = D(IPR(I))
30    CONTINUE
      IND = IND - 1
      RETURN
40    IERR = -IND - 1
      RETURN
      END
C
C
C
C
C
      SUBROUTINE HESOLV(D,IPR,M,REPS,IERR)
C
      INTEGER I,I1,I2,IERR,IPR(1),IROW1,J,J1,K,N,N1,MFIN
      REAL*8 D(1),DABS,MULT,REPS
C
C     THIS SUBROUTINE SOLVES AN UPPER-HESSENBERG SYSTEM OF EQUATIONS
C     WHICH IS STORED IN PACKED FORM IN THE VECTOR D.  THE METHOD USED
C     IS GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING.  THE PARAMETERS
C     IN THE SUBROUTINE ARE:
C     D      -  MATRIX OF COEFFICIENTS, RIGHT HAND SIDE (ON INPUT),
C               AND SOLUTION (ON OUTPUT)
C     IPR    -  INTEGER VECTOR WHICH RECORDS INTERCHANGES
C     M      -  SIZE OF SYSTEM
C     REPS   -  ERROR TOLERANCE
C     IERR   -  ERROR CODE (IERR=0 => NORMAL RETURN)
C                          (IERR=-1 => SINGULAR MATRIX)
C
C     ----------------------------------------------------------------
C
C     INITIALIZE INTERCHANGE VECTORS AND PARAMETERS.
C
      IERR = 0
      MFIN = (M * (M + 1)) / 2 + M
      DO 10 I = 1,M
         IPR(M+I) = IROW1(I,M)
         IPR(I) = I + MFIN
10    CONTINUE
      M1 = M - 1
C
C     REDUCE THE MATRIX TO UPPER TRIANGULAR FORM
C
      IF (M.EQ.1) GO TO 35
      DO 30 I = 1,M1
         IF (DABS(D(IPR(M+I)+1)).GT.DABS(D(IPR(M+I+1)+1))) GO TO 20
C
C        INTERCHANGE ROWS IF NECESSARY
C
         K = IPR(M+I)
         IPR(M+I) = IPR(M+I+1)
         IPR(M+I+1) = K
         K = IPR(I)
         IPR(I) = IPR(I+1)
         IPR(I+1) = K
C
C        CHECK FOR COMPUTATIONALLY SINGULAR MATRIX
C
20       IF (DABS(D(IPR(M+I)+1)).LT.REPS) GO TO 60
         IPR(M+I+1) = IPR(M+I+1) + 1
C
C        ELIMINATE SUBDIAGONAL ELEMENT OF A
C
         MULT = D(IPR(M+I+1)) / D(IPR(M+I)+1)
         D(IPR(I+1)) = D(IPR(I+1)) - MULT * D(IPR(I))
         I1 = I + 1
         DO 30 J = I1,M
            D(IPR(M+I+1)+J-I) = D(IPR(M+I+1)+J-I) -
     *                            MULT*D(IPR(M+I)+J+1-I)
30    CONTINUE
35    IF (DABS(D(IPR(M+M)+1)).LT.REPS) GO TO 60
C
C     PERFORM BACK - SUBSTITUTION
C
      D(IPR(M)) = D(IPR(M)) / D(IPR(M+M)+1)
      IF (M1.EQ.0) RETURN
      DO 50 I1 = 1,M1
         I = M - I1
         I2 = I + 1
         MULT = 0.D0
         DO 40 J1 = I2,M
            J = J1 - I2 + 2
            MULT = MULT + D(IPR(J1)) * D(IPR(M+I)+J)
40       CONTINUE
         D(IPR(I)) = (D(IPR(I)) - MULT) / D(IPR(M+I)+1)
50    CONTINUE
      RETURN
60    IERR = -1
      RETURN
      END
C
C
C
C
      SUBROUTINE BACKSB(C,B,IND,N,M,MDIM,NDIM)
C
      INTEGER I,IND,IND1,IND2,J,M,MDIM,N,NDIM
      REAL*8 B(NDIM,N),C(MDIM,N)
C
C     BLOCK BACK - SUBSTITUTION FOR ONE 'ROW'.
C     PARAMETERS ARE AS IN SUBROUTINE AXXBC.
C
      IND1 = IND + 1
      IND2 = IND + 2
      DO 10 I = IND2,N
         DO 10 J = 1,M
            C(J,IND1) = C(J,IND1) - B(IND1,I) * C(J,I)
10    CONTINUE
      RETURN
      END
C
C
C
C
C
      INTEGER FUNCTION IROW1(I,M)
C
      INTEGER I,M
C
C     THIS FUNCTION COMPUTES THE INDEX (FOR THE VECTOR D) ONE BEFORE
C     THE BEGINNING OF ROW I IN THE M*M MATRIX.
C
      IROW1 = (I-1) * M - (I-2) * (I-3) / 2
      IF (I.EQ.1) IROW1 = 0
      RETURN
      END
C
C
C
C
C
      SUBROUTINE N2SOLV(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR)
C
      INTEGER I,I1,IERR,IND,IPR(1),IROW2,J,J1,J2,K,LROW2,M,M1,
     *   MDIM,N,NDIM,MFIN
      REAL*8 A(MDIM,M),B(NDIM,N),C(MDIM,N),D(1),REPS
C
C     THIS SUBROUTINE SETS UP A DOUBLE SYSTEM OF EQUATIONS.
C     IT SOLVES (USING SUBROUTINE H2SOLV) THE SET OF EQUATIONS
C
C     ( A+B(IND,IND)*I        B(IND,IND+1) )( X(IND) )   ( C(IND) )
C     (                                    )(        ) = (        )
C     ( B(IND+1,IND)*I  A+B(IND+1,IND+1)*I )(X(IND+1))   (C(IND+1))
C
C     BY INTERCHANGING ROWS AND COLUMNS (TO INSURE THAT THERE ARE
C     ZEROES BELOW THE SUB-SUB-DIAGONAL) AND THEN USING GAUSSIAN
C     ELIMINATION.  IN ORDER TO TAKE ADVANTAGE OF THE NUMBER OF
C     ZEROES IN THIS MATRIX, THE MATRIX IS STORED IN THE VECTOR
C     D IN PACKED FORM. (ENTRIES BELOW THE SUB-SUB-DIAGONAL ARE
C     IGNORED) THE RIGHT HAND SIDE AND SOLUTION ARE ALSO STORED IN D.
C     THE PARAMETERS ARE AS IN THE SUBROUTINE AXXBC.
C
C     ------------------------------------------------------------
C
C     PERFORM BLOCK BACK-SUBSTITUTION IF NECESSARY
C
      IF (IND.LT.N-1) CALL BACKS2(C,B,IND,N,M,MDIM,NDIM)
C
C     SET UP THE SYSTEM OF EQUATIONS FOR H2SOLV
C
      M2 = 2 * M
      MFIN = (M2 * (M2 + 1)) / 2 + 4 * M
      DO 20 I = 1,M
C
C        FIND BEGINNING AND LENGTH OF ROW 2I-1 IN THE VECTOR D.
C
         M1 = IROW2(2*I-1,M)
         K = LROW2(2*I-1,M)
         I1 = I - 1
         IF (I.EQ.1) I1 = 1
         DO 10 J = I1,M
            J1 = 2 * (J - I1 + 1)
            J2 = 1
            IF (M1.EQ.0) J2 = 0
            D(M1+J1-1) = A(I,J)
            D(M1+J1) = 0.D0
            D(M1+K+J1-J2) = A(I,J)
            D(M1+K+J1-1-J2) = 0.D0
10       CONTINUE
         J1 = 3
         IF (I.EQ.1) J1 = 1
         D(J1+M1) = D(J1+M1) + B(IND,IND)
         D(J1+M1+1) = D(J1+M1+1) + B(IND,IND+1)
         IF (I.NE.1) J1 = 2
         D(J1+M1+K) = D(J1+M1+K) + B(IND+1,IND)
         D(J1+M1+K+1) = D(J1+M1+K+1) + B(IND+1,IND+1)
C
C        STORE RIGHT HAND SIDE
C
         D(2*I+MFIN) = C(I,IND+1)
         D(2*I-1+MFIN) = C(I,IND)
20    CONTINUE
C
C     SOLVE THE SYSTEM OF EQUATIONS AND STORE THE RESULT BACK IN C
C
      CALL H2SOLV(D,IPR,M,REPS,IERR)
      IF (IERR.NE.0) GO TO 40
      DO 30 I = 1,M
         C(I,IND) = D(IPR(2*I-1))
         C(I,IND+1) = D(IPR(2*I))
30    CONTINUE
      IND = IND - 2
      RETURN
40    IERR = -IND - 1
      RETURN
      END
C
C
C
C
C
      SUBROUTINE H2SOLV(D,IPR,M,REPS,IERR)
C
      INTEGER I,I1,I2,IERR,IPR(1),IROW2,J,J1,K,K1,L,N,N2,N21,MFIN
      REAL*8 REPS,D(1),MAX,DABS
C
C     THIS SUBROUTINE SOLVES A SYSTEM OF 2M * 2M EQUATIONS
C     WHERE THE MATRIX HAS ZEROES BELOW THE SUB - SUB - DIAGONAL
C     AND WHERE THE MATRIX IS STORED AS A VECTOR IN PACKED FORM.
C     THE METHOD USED IS GAUSSIAN ELIMINATION.
C     PARAMETERS:
C     D      -  MATRIX OF COEFFICIENTS, RIGHT HAND SIDE (ON INPUT),
C               AND SOLUTION (ON OUTPUT).
C     IPR    -  VECTOR USED TO RECORD ROW PERMUTATIONS
C     M      -  SIZE OF SYSTEM
C     REPS   -  ERROR TOLERANCE
C     IERR   -  ERROR CODE (IERR=0 => NORMAL RETURN)
C                          (IERR=-1 => SINGULAR MATRIX)
C
C     --------------------------------------------------------------
C
C     INITIALIZE PIVOT VECTORS AND SIZE PARAMETERS.
C
      IERR = 0
      M2 = 2 * M
      MFIN = (M2 * (M2 + 1)) / 2 + 4 * M
      DO 10 I = 1,M2
         IPR(M2+I) = IROW2(I,M)
         IPR(I) = I + MFIN
10    CONTINUE
      M21 = M2 - 1
C
C     TRANSFORM MATRIX TO TRIANGULAR FORM
C
      DO 40 I = 1,M21
         I1 = 2
         IF (I.EQ.M21) I1 = 1
         L = 0
C
C        INTERCHANGE ROWS IF NECESSARY
C
         MAX = DABS(D(IPR(M2+I) + 1))
         DO 20 J = 1,I1
            IF (DABS(D(IPR(M2+J+I) + 1)).LE.MAX) GO TO 20
            MAX = DABS(D(IPR(M2+J+I) + 1))
            L = J
20       CONTINUE
C
C        CHECK FOR COMPUTATIONALLY SINGULAR MATRIX
C
         IF (MAX.LE.REPS) GO TO 80
         IF (L.EQ.0) GO TO 30
         K = IPR(M2+I)
         IPR(M2+I) = IPR(M2+L+I)
         IPR(M2+L+I) = K
         K = IPR(I)
         IPR(I) = IPR(I+L)
         IPR(I+L) = K
30       CONTINUE
C
C        ADJUST POINTERS TO BEGINNINGS OF ROWS
C
         IPR(M2+I+1) = IPR(M2+I+1) + 1
         IF (I.NE.M21) IPR(M2+I+2) = IPR(M2+I+2) + 1
         IP1 = I + 1
C
C        ELIMINATE SUBDIAGONAL ELEMENTS IN THE MATRIX
C
         DO 40 J = 1,I1
            MAX = D(IPR(M2+I+J)) / D(IPR(M2+I)+1)
            D(IPR(I+J)) = D(IPR(I+J)) - MAX * D(IPR(I))
            DO 40 K1 = IP1,M2
               K = K1 - I
               D(IPR(M2+I+J)+K) = D(IPR(M2+I+J)+K) -
     *                                    MAX * D(IPR(M2+I)+1+K)
40    CONTINUE
      IF (DABS(D(IPR(M2+M2)+1)).LE.REPS) GO TO 80
C
C     PERFORM BACK SUBSTITUTION
C
      D(IPR(M2)) = D(IPR(M2)) / D(IPR(M2+M2)+1)
      DO 60 I1 = 1,M21
         I = M2 - I1
         I2 = I + 1
         MAX = 0.D0
         DO 50 J1 = I2,M2
            J = J1 - I2 + 2
            MAX = MAX + D(IPR(J1)) * D(IPR(M2+I)+J)
50       CONTINUE
         D(IPR(I)) = (D(IPR(I)) - MAX) / D(IPR(M2+I)+1)
60    CONTINUE
70    RETURN
80    IERR = -1
      GO TO 70
      END
C
C
C
C
C
      SUBROUTINE BACKS2(C,B,IND,N,M,MDIM,NDIM)
C
      INTEGER I,IND,IND1,IND2,J,M,MDIM,N,NDIM
      REAL*8 B(NDIM,N),C(MDIM,N)
C
C     BLOCK BACK - SUBSTITUTION FOR TWO 'ROWS'.
C     PARAMETERS ARE AS IN SUBROUTINE AXXBC.
C
      IND1 = IND + 1
      IND2 = IND + 2
      DO 10 I = IND2,N
         DO 10 J = 1,M
            C(J,IND1) = C(J,IND1) - B(IND1,I) * C(J,I)
            C(J,IND) = C(J,IND) - B(IND,I) * C(J,I)
10    CONTINUE
      RETURN
      END
C
C
C
C
C
      INTEGER FUNCTION IROW2(I,M)
C
      INTEGER I,K,M
C
C     THIS FUNCTION FINDS THE INDEX (FOR THE VECTOR D) ONE
C     BEFORE THE BEGINNING OF ROW I IN THE 2N*2N MATRIX.
C
      IROW2 = 2 * (I - 1) * M
      K = (I - 4) * (I - 3) / 2
      IF (I.LE.2) K = 0
      IROW2 = IROW2 - K
      RETURN
      END
C
C
C
C
C
      INTEGER FUNCTION LROW2(I,M)
C
      INTEGER I,K,M
C
C     THIS FUNCTION FINDS THE LENGTH OF ROW I IN THE
C     DOUBLE MATRIX (WHEN IT IS STORED IN PACKED FORM IN THE
C     VECTOR D)
C
      K = I - 3
      IF (I.LE.2) K = 0
      LROW2 = 2 * M - K
      IF (I.LE.2) LROW2 = 2 * M
      RETURN
      END