C     ALGORITHM 603, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 3,
C     SEP., 1983, P. 376-380.
C***********************************************************************MAN   10
C                                                                       MAN   20
C   THIS DRIVER RUNS A TEST PROBLEM USING C O L R O W.                  MAN   30
C   THE ARRAYS TOP, AR, BOT, AND B ARE SET UP IN DATA                   MAN   40
C   STATEMENTS, AS ARE THE PARAMETERS NRWTOP, NRWBLK, NCLBLK,           MAN   50
C   NBLOKS, NRWBOT, WHICH DEFINE THE STRUCTURE OF THE COEFF-            MAN   60
C   ICIENT MATRIX.                                                      MAN   70
C                                                                       MAN   80
C***********************************************************************MAN   90
C                                                                       MAN  100
      DOUBLE PRECISION TOP, AR, BOT, B, X                               MAN  110
      DOUBLE PRECISION ERROR, ERR                                       MAN  120
      DIMENSION TOP(2,4), AR(4,8,2), BOT(2,4), B(12), X(12)             MAN  130
      INTEGER PIVOT(12)                                                 MAN  140
      DATA N, NRWTOP, NOVRLP, NRWBLK, NCLBLK, NBLOKS, NRWBOT            MAN  150
     * /12,2,4,4,8,2,2/                                                 MAN  160
      DATA TOP(1,1), TOP(1,2), TOP(1,3), TOP(1,4), TOP(2,1), TOP(2,2),  MAN  170
     * TOP(2,3), TOP(2,4) /0.0D0,-0.98D0,-0.79D0,-0.15D0,-1.00D0,0.25D0,MAN  180
     * -0.87D0,0.35D0/                                                  MAN  190
      DATA AR(1,1,1), AR(1,2,1), AR(1,3,1), AR(1,4,1), AR(1,5,1),       MAN  200
     * AR(1,6,1), AR(1,7,1), AR(1,8,1) /0.78D0,0.31D0,-0.85D0,0.89D0,   MAN  210
     * -0.69D0,-0.98D0,-0.76D0,-0.82D0/                                 MAN  220
      DATA AR(2,1,1), AR(2,2,1), AR(2,3,1), AR(2,4,1), AR(2,5,1),       MAN  230
     * AR(2,6,1), AR(2,7,1), AR(2,8,1) /0.12D0,-0.01D0,0.75D0,0.32D0,   MAN  240
     * -1.00D0,-0.53D0,-0.83D0,-0.98D0/                                 MAN  250
      DATA AR(3,1,1), AR(3,2,1), AR(3,3,1), AR(3,4,1), AR(3,5,1),       MAN  260
     * AR(3,6,1), AR(3,7,1), AR(3,8,1) /-0.58D0,0.04D0,0.87D0,0.38D0,   MAN  270
     * -1.00D0,-0.21D0,-0.93D0,-0.84D0/                                 MAN  280
      DATA AR(4,1,1), AR(4,2,1), AR(4,3,1), AR(4,4,1), AR(4,5,1),       MAN  290
     * AR(4,6,1), AR(4,7,1), AR(4,8,1) /-0.21D0,-0.91D0,-0.09D0,-0.62D0,MAN  300
     * -1.99D0,-1.12D0,-1.21D0,0.07D0/                                  MAN  310
      DATA AR(1,1,2), AR(1,2,2), AR(1,3,2), AR(1,4,2), AR(1,5,2),       MAN  320
     * AR(1,6,2), AR(1,7,2), AR(1,8,2) /0.78D0,-0.93D0,-0.76D0,0.48D0,  MAN  330
     * -0.87D0,-0.14D0,-1.00D0,-0.59D0/                                 MAN  340
      DATA AR(2,1,2), AR(2,2,2), AR(2,3,2), AR(2,4,2), AR(2,5,2),       MAN  350
     * AR(2,6,2), AR(2,7,2), AR(2,8,2) /-0.99D0,0.21D0,-0.73D0,-0.48D0, MAN  360
     * -0.93D0,-0.91D0,0.10D0,-0.89D0/                                  MAN  370
      DATA AR(3,1,2), AR(3,2,2), AR(3,3,2), AR(3,4,2), AR(3,5,2),       MAN  380
     * AR(3,6,2), AR(3,7,2), AR(3,8,2) /-0.68D0,-0.09D0,-0.58D0,-0.21D0,MAN  390
     * 0.85D0,-0.39D0,0.79D0,-0.71D0/                                   MAN  400
      DATA AR(4,1,2), AR(4,2,2), AR(4,3,2), AR(4,4,2), AR(4,5,2),       MAN  410
     * AR(4,6,2), AR(4,7,2), AR(4,8,2) /0.39D0,-0.99D0,-0.12D0,-0.75D0, MAN  420
     * -0.68D0,-0.99D0,0.50D0,-0.88D0/                                  MAN  430
      DATA BOT(1,1), BOT(1,2), BOT(1,3), BOT(1,4), BOT(2,1), BOT(2,2),  MAN  440
     * BOT(2,3), BOT(2,4) /0.71D0,-0.64D0,0.0D0,0.48D0,0.08D0,100.0D0,  MAN  450
     * 50.00D0,15.00D0/                                                 MAN  460
      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), MAN  470
     * B(11), B(12) /-1.92D0,-1.27D0,-2.12D0,-2.16D0,-2.27D0,-6.08D0,   MAN  480
     * -3.03D0,-4.62D0,-1.02D0,-3.52D0,.55D0,165.08D0/                  MAN  490
C                                                                       MAN  500
C***********************************************************************MAN  510
C                                                                       MAN  520
C   THE INPUT MATRIX IS GIVEN BY:                                       MAN  530
C                                                                       MAN  540
C  0.0  -0.98 -0.79 -0.15                                               MAN  550
C -1.00  0.25 -0.87  0.35                                               MAN  560
C  0.78  0.31 -0.85  0.89 -0.69 -0.98 -0.76 -0.82                       MAN  570
C  0.12 -0.01  0.75  0.32 -1.00 -0.53 -0.83 -0.98                       MAN  580
C -0.58  0.04  0.87  0.38 -1.00 -0.21 -0.93 -0.84                       MAN  590
C -0.21 -0.91 -0.09 -0.62 -1.99 -1.12 -1.21  0.07                       MAN  600
C                          0.78 -0.93 -0.76  0.48 -0.87 -0.14 -1.00 -0.5MAN  610
C                         -0.99  0.21 -0.73 -0.48 -0.93 -0.91  0.10 -0.8MAN  620
C                         -0.68 -0.09 -0.58 -0.21  0.85 -0.39  0.79 -0.7MAN  630
C                          0.39 -0.99 -0.12 -0.75 -0.68 -0.99  0.50 -0.8MAN  640
C                                                  0.71 -0.64  0.0   0.4MAN  650
C                                                  0.08 100.0 50.00 15.0MAN  660
C                                                                       MAN  670
C       THE RIGHT HAND SIDE IS GIVEN BY:                                MAN  680
C                                                                       MAN  690
C         B = (-1.92,-1.27,-2.12,-2.16,-2.27,-6.08,-3.03,-4.62,         MAN  700
C              -1.02,-3.52,0.55,165.08)                                 MAN  710
C                                                                       MAN  720
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY;                        MAN  730
C                                                                       MAN  740
C          X = (1,1,1,1,1,1,1,1,1,1,1,1)                                MAN  750
C                                                                       MAN  760
C***********************************************************************MAN  770
C                                                                       MAN  780
      CALL COLROW(N, TOP, NRWTOP, NOVRLP, AR, NRWBLK, NCLBLK, NBLOKS,   MAN  790
     * BOT, NRWBOT, PIVOT, B, X, IFLAG)                                 MAN  800
      IF (IFLAG.NE.0) GO TO 20                                          MAN  810
      ERROR = 0.D0                                                      MAN  820
      DO 10 I=1,N                                                       MAN  830
        ERR = 1.D0 - X(I)                                               MAN  840
        ERROR = DMAX1(ERROR,DABS(ERR))                                  MAN  850
        WRITE (6,99998) X(I), ERR                                       MAN  860
   10 CONTINUE                                                          MAN  870
      WRITE (6,99999) ERROR                                             MAN  880
99999 FORMAT (12H MAX ERROR =, D15.7)                                   MAN  890
99998 FORMAT (1H , F15.7, D15.7)                                        MAN  900
      RETURN                                                            MAN  910
   20 CONTINUE                                                          MAN  920
      WRITE (6,99997) IFLAG                                             MAN  930
99997 FORMAT (9H IFLAG = , I3)                                          MAN  940
      RETURN                                                            MAN  950
      END                                                               MAN  960
      SUBROUTINE COLROW(N, TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK,       COL   10
     * NCLBLK, NBLOKS, BOTBLK, NRWBOT, PIVOT, B, X, IFLAG)
C
C***************************************************************
C
C  THIS PROGRAM SOLVES THE LINEAR SYSTEM  A*X = B  WHERE  A IS
C  AN ALMOST BLOCK DIAGONAL MATRIX OF THE FORM
C
C               TOPBLK
C               ARRAY(1)
C                     ARRAY(2)
C                          .
C                             .
C                                .
C                                   .
C                                    ARRAY(NBLOKS)
C                                           BOTBLK
C
C  WHERE
C           TOPBLK IS  NRWTOP  BY NOVRLP
C           ARRAY(K), K=1,NBLOKS, ARE NRWBLK BY NRWBLK+NOVRLP
C           BOTBLK IS NRWBOT BY NOVRLP,
C  AND
C           NOVRLP = NRWTOP + NRWBOT
C  WITH
C           NOVRLP.LE.NRWBLK .
C
C  THE LINEAR SYSTEM IS OF ORDER  N = NBLOKS*NRWBLK + NOVRLP.
C
C  THE METHOD IMPLEMENTED IS BASED ON GAUSS ELIMINATION WITH
C  ALTERNATE ROW AND COLUMN ELIMINATION WITH PARTIAL PIVOTING,
C  WHICH PRODUCES A STABLE DECOMPOSITION OF THE MATRIX  A
C  WITHOUT INTRODUCING FILL-IN.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE, REMOVE
C  ALL DOUBLE PRECISION STATEMENTS.  THERE IS ONE SUCH STATEMENT
C  IN C O L R O W, THREE IN C R D C M P, AND TWO IN C R S O L V.
C  IN ADDITION, REFERENCES TO BUILT-IN FUNCTIONS DABS AND DMAX1
C  MUST BE REPLACED BY ABS AND AMAX1, RESPECTIVELY.  DABS OCCURS
C  NINE TIMES, IN C R D C M P.  DMAX1 OCCURS FOUR TIMES, IN
C  C R D C M P.  FINALLY, ZERO IS INITIALISED TO 0.D0 IN A
C  DATA STATEMENT IN C R D C M P.  THIS MUST BE REPLACED BY:
C               DATA ZERO/0.0/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR (IF IFLAG = 0)
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  AUXILIARY PROGRAMS  *****
C
C       CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,IFLAG)
C            - DECOMPOSES THE MATRIX  A  USING MODIFIED
C              ALTERNATE ROW AND COLUMN ELIMINATON WITH
C              PARTIAL PIVOTING, AND IS USED FOR THIS
C              PURPOSE IN C O L R O W.
C              THE ARGUMENTS ARE AS IN C O L R O W.
C
C       CRSLVE(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,B,X)
C            - SOLVES THE SYSTEM A*X = B ONCE A IS DECOMPOSED.
C              THE ARGUMENTS ARE ALLAS IN C O L R O W.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       THE SUBROUTINE  C O L R O W  AUTOMATICALLY SOLVES THE
C  INPUT SYSTEM WHEN IFLAG=0.  C O L R O W  IS CALLED ONLY ONCE
C  FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT
C  HAND SIDES CAN BE OBTAINED BY ONE CALL TO  C O L R O W  AND
C  P-1 CALLS TO CRSLVE ONLY. SINCE THE ARRAYS TOPBLK,ARRAY,
C  BOTBLK AND PIVOT CONTAIN THE DECOMPOSITION OF THE GIVEN
C  COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM
C  C O L R O W , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE
C  CALLS TO CRSLVE WITH THE SAME LEFT HAND SIDES. FOR THE
C  SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT
C  MATRIX, THE ARRAYS TOPBLK,ARRAY,BOTBLK MUST BE COPIED
C  BEFORE A CALL TO  C O L R O W .
C
C***********************************************************************
C
C             *****  SAMPLE CALLING PROGRAM  *****
C
C      THE FOLLOWING PROGRAM WILL EXERCISE COLROW, IN THE
C      CASE WHEN THE COEFFICIENT MATRIX IS NON-SINGULAR.
C
C       DOUBLE PRECISION TOP,AR,BOT,B,X
C       DOUBLE PRECISION ERROR,ERR
C       DIMENSION TOP(2,4),AR(4,8,2),BOT(2,4),B(12),X(12)
C       INTEGER PIVOT(12)
C       DATA N,NRWTOP,NOVRLP,NRWBLK,NCLBLK,NBLOKS,NRWBOT/12,2,4,4,8,2,2/
C       DATA TOP(1,1),TOP(1,2),TOP(1,3),TOP(1,4),
C    *       TOP(2,1),TOP(2,2),TOP(2,3),TOP(2,4)/
C    *0.0 D0,-0.98D0,-0.79D0,-0.15D0,
C    *-1.00D0, 0.25D0,-0.87D0, 0.35D0/
C       DATA AR(1,1,1),AR(1,2,1),AR(1,3,1),AR(1,4,1),
C    *       AR(1,5,1),AR(1,6,1),AR(1,7,1),AR(1,8,1)/
C    *0.78D0, 0.31D0,-0.85D0, 0.89D0,-0.69D0,-0.98D0,-0.76D0,-0.82D0/
C       DATA AR(2,1,1),AR(2,2,1),AR(2,3,1),AR(2,4,1),
C    *       AR(2,5,1),AR(2,6,1),AR(2,7,1),AR(2,8,1)/
C    *0.12D0,-0.01D0, 0.75D0, 0.32D0,-1.00D0,-0.53D0,-0.83D0,-0.98D0/
C       DATA AR(3,1,1),AR(3,2,1),AR(3,3,1),AR(3,4,1),
C    *       AR(3,5,1),AR(3,6,1),AR(3,7,1),AR(3,8,1)/
C    *-0.58D0, 0.04D0, 0.87D0, 0.38D0,-1.00D0,-0.21D0,-0.93D0,-0.84D0/
C       DATA AR(4,1,1),AR(4,2,1),AR(4,3,1),AR(4,4,1),
C    *       AR(4,5,1),AR(4,6,1),AR(4,7,1),AR(4,8,1)/
C    *-0.21D0,-0.91D0,-0.09D0,-0.62D0,-1.99D0,-1.12D0,-1.21D0, 0.07D0/
C       DATA AR(1,1,2),AR(1,2,2),AR(1,3,2),AR(1,4,2),
C    *       AR(1,5,2),AR(1,6,2),AR(1,7,2),AR(1,8,2)/
C    *0.78D0,-0.93D0,-0.76D0, 0.48D0,-0.87D0,-0.14D0,-1.00D0,-0.59D0/
C       DATA AR(2,1,2),AR(2,2,2),AR(2,3,2),AR(2,4,2),
C    *       AR(2,5,2),AR(2,6,2),AR(2,7,2),AR(2,8,2)/
C    *-0.99D0, 0.21D0,-0.73D0,-0.48D0,-0.93D0,-0.91D0, 0.10D0,-0.89D0/
C       DATA AR(3,1,2),AR(3,2,2),AR(3,3,2),AR(3,4,2),
C    *       AR(3,5,2),AR(3,6,2),AR(3,7,2),AR(3,8,2)/
C    *-0.68D0,-0.09D0,-0.58D0,-0.21D0, 0.85D0,-0.39D0, 0.79D0,-0.71D0/
C       DATA AR(4,1,2),AR(4,2,2),AR(4,3,2),AR(4,4,2),
C    *       AR(4,5,2),AR(4,6,2),AR(4,7,2),AR(4,8,2)/
C    *0.39D0,-0.99D0,-0.12D0,-0.75D0,-0.68D0,-0.99D0, 0.50D0,-0.88D0/
C       DATA BOT(1,1),BOT(1,2),BOT(1,3),BOT(1,4),
C    *       BOT(2,1),BOT(2,2),BOT(2,3),BOT(2,4)/
C    *0.71D0,-0.64D0, 0.0 D0, 0.48D0,
C    *0.08D0,100.0D0,50.00D0,15.00D0/
C       DATA B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9),B(10),B(11),
C    *       B(12)/
C    *-1.92D0,-1.27D0,-2.12D0,-2.16D0,-2.27D0,-6.08D0,-3.03D0,
C    *-4.62D0,-1.02D0,-3.52D0,.55D0,165.08D0/
C
C***********************************************************************
C
C   THE INPUT MATRIX IS GIVEN BY:
C
C  0.0  -0.98 -0.79 -0.15
C -1.00  0.25 -0.87  0.35
C  0.78  0.31 -0.85  0.89 -0.69 -0.98 -0.76 -0.82
C  0.12 -0.01  0.75  0.32 -1.00 -0.53 -0.83 -0.98
C -0.58  0.04  0.87  0.38 -1.00 -0.21 -0.93 -0.84
C -0.21 -0.91 -0.09 -0.62 -1.99 -1.12 -1.21  0.07
C                          0.78 -0.93 -0.76  0.48 -0.87 -0.14 -1.00 -0.5
C                         -0.99  0.21 -0.73 -0.48 -0.93 -0.91  0.10 -0.8
C                         -0.68 -0.09 -0.58 -0.21  0.85 -0.39  0.79 -0.7
C                          0.39 -0.99 -0.12 -0.75 -0.68 -0.99  0.50 -0.8
C                                                  0.71 -0.64  0.0   0.4
C                                                  0.08 100.0 50.00 15.0
C
C       THE RIGHT HAND SIDE IS GIVEN BY:
C
C         B = (-1.92,-1.27,-2.12,-2.16,-2.27,-6.08,-3.03,-4.62,
C              -1.02,-3.52,0.55,165.08)
C
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY;
C
C          X = (1,1,1,1,1,1,1,1,1,1,1,1)
C
C***********************************************************************
C
C       CALL COLROW(N,TOP,NRWTOP,NOVRLP,AR,NRWBLK,NCLBLK,NBLOKS,
C    *              BOT,NRWBOT,PIVOT,B,X,IFLAG)
C       IF(IFLAG.NE.0)GO TO 1000
C       ERROR = 0.D0
C       DO 10 I=1,N
C          ERR = 1.D0 - X(I)
C          ERROR = DMAX1(ERROR,DABS(ERR))
C          WRITE(6,100)X(I),ERR
C  10   CONTINUE
C       WRITE(6,200)ERROR
C 200   FORMAT(12H MAX ERROR = ,D15.7)
C 100   FORMAT(1H ,F15.7,D15.7)
C       RETURN
C1000   CONTINUE
C       WRITE(6,300)IFLAG
C 300   FORMAT(9H IFLAG =  ,I3)
C       RETURN
C       END
C
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, ARRAY, BOTBLK, B, X
      INTEGER PIVOT(1)
      DIMENSION TOPBLK(NRWTOP,1), ARRAY(NRWBLK,NCLBLK,1),
     * BOTBLK(NRWBOT,1), B(1), X(1)
      CALL CRDCMP(N, TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK, NCLBLK,
     * NBLOKS, BOTBLK, NRWBOT, PIVOT, IFLAG)
      IF (IFLAG.NE.0) RETURN
      CALL CRSLVE(TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK, NCLBLK,
     * NBLOKS, BOTBLK, NRWBOT, PIVOT, B, X)
      RETURN
      END
      SUBROUTINE CRDCMP(N, TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK,       CRD   10
     * NCLBLK, NBLOKS, BOTBLK, NRWBOT, PIVOT, IFLAG)
C
C***************************************************************
C
C  C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A
C  USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH
C  PARTIAL PIVOTING.  THE MATRIX  A  IS STORED IN THE ARRAYS
C  TOPBLK, ARRAY, AND BOTBLK.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A TO BE DECOMPOSED
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, ARRAY, BOTBLK
      DOUBLE PRECISION ROWMAX, ROWPIV, ROWMLT, COLMAX, COLPIV
      DOUBLE PRECISION SWAP, COLMLT, PIVMAX, ZERO, TEMPIV
      INTEGER PIVOT(1)
      DIMENSION TOPBLK(NRWTOP,1), ARRAY(NRWBLK,NCLBLK,1),
     * BOTBLK(NRWBOT,1)
      DATA ZERO /0.0D0/
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
      IFLAG = 0
      PIVMAX = ZERO
      NRWTP1 = NRWTOP + 1
      NROWEL = NRWBLK - NRWTOP
      NRWEL1 = NROWEL + 1
      NVRLP0 = NOVRLP - 1
C
C***************************************************************
C
C          ****  CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C               IF PARAMETERS ARE INVALID THEN TERMINATE AT 10;
C                                         ELSE CONTINUE AT 100.
C
C***************************************************************
C
      IF (N.NE.NBLOKS*NRWBLK+NOVRLP) GO TO 10
      IF (NOVRLP.NE.NRWTOP+NRWBOT) GO TO 10
      IF (NCLBLK.NE.NOVRLP+NRWBLK) GO TO 10
      IF (NOVRLP.GT.NRWBLK) GO TO 10
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      GO TO 20
   10 CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE INVALID.  SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = 1
      RETURN
   20 CONTINUE
C
C***************************************************************
C
C               ****  FIRST, IN TOPBLK....
C
C***************************************************************
C
C          ***  APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN
C                 PIVOTING ....
C
C***************************************************************
C
      DO 110 I=1,NRWTOP
        IPLUS1 = I + 1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IPVT = I
        COLMAX = DABS(TOPBLK(I,I))
        DO 30 J=IPLUS1,NOVRLP
          TEMPIV = DABS(TOPBLK(I,J))
          IF (TEMPIV.LE.COLMAX) GO TO 30
          IPVT = J
          COLMAX = TEMPIV
   30   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410
        PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        PIVOT(I) = IPVT
        IF (IPVT.EQ.I) GO TO 60
        DO 40 L=I,NRWTOP
          SWAP = TOPBLK(L,IPVT)
          TOPBLK(L,IPVT) = TOPBLK(L,I)
          TOPBLK(L,I) = SWAP
   40   CONTINUE
        DO 50 L=1,NRWBLK
          SWAP = ARRAY(L,IPVT,1)
          ARRAY(L,IPVT,1) = ARRAY(L,I,1)
          ARRAY(L,I,1) = SWAP
   50   CONTINUE
   60   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        COLPIV = TOPBLK(I,I)
        DO 100 J=IPLUS1,NOVRLP
          COLMLT = TOPBLK(I,J)/COLPIV
          TOPBLK(I,J) = COLMLT
          IF (IPLUS1.GT.NRWTOP) GO TO 80
          DO 70 L=IPLUS1,NRWTOP
            TOPBLK(L,J) = TOPBLK(L,J) - COLMLT*TOPBLK(L,I)
   70     CONTINUE
   80     CONTINUE
          DO 90 L=1,NRWBLK
            ARRAY(L,J,1) = ARRAY(L,J,1) - COLMLT*ARRAY(L,I,1)
   90     CONTINUE
  100   CONTINUE
  110 CONTINUE
C
C***************************************************************
C
C          ****  IN EACH BLOCK ARRAY(,,K)....
C
C***************************************************************
C
      INCR = 0
      DO 320 K=1,NBLOKS
        KPLUS1 = K + 1
C
C          *****************************************************
C
C          ***  FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH
C                       ROW PIVOTING....
C
C          *****************************************************
C
        DO 180 J=NRWTP1,NRWBLK
          JPLUS1 = J + 1
          JMINN = J - NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IPVT = JMINN
          ROWMAX = DABS(ARRAY(JMINN,J,K))
          LOOP = JMINN + 1
          DO 120 I=LOOP,NRWBLK
            TEMPIV = DABS(ARRAY(I,J,K))
            IF (TEMPIV.LE.ROWMAX) GO TO 120
            IPVT = I
            ROWMAX = TEMPIV
  120     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410
          PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          INCRJ = INCR + J
          PIVOT(INCRJ) = INCR + IPVT + NRWTOP
          IF (IPVT.EQ.JMINN) GO TO 140
          DO 130 L=J,NCLBLK
            SWAP = ARRAY(IPVT,L,K)
            ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K)
            ARRAY(JMINN,L,K) = SWAP
  130     CONTINUE
  140     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          ROWPIV = ARRAY(JMINN,J,K)
          DO 150 I=LOOP,NRWBLK
            ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV
  150     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          DO 170 L=JPLUS1,NCLBLK
            ROWMLT = ARRAY(JMINN,L,K)
            DO 160 I=LOOP,NRWBLK
              ARRAY(I,L,K) = ARRAY(I,L,K) - ROWMLT*ARRAY(I,J,K)
  160       CONTINUE
  170     CONTINUE
  180   CONTINUE
C
C          *****************************************************
C
C          ***  NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH
C                      COLUMN PIVOTING....
C
C          *****************************************************
C
        DO 310 I=NRWEL1,NRWBLK
          IPLUSN = I + NRWTOP
          IPLUS1 = I + 1
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IPVT = IPLUSN
          COLMAX = DABS(ARRAY(I,IPVT,K))
          LOOP = IPLUSN + 1
          DO 190 J=LOOP,NCLBLK
            TEMPIV = DABS(ARRAY(I,J,K))
            IF (TEMPIV.LE.COLMAX) GO TO 190
            IPVT = J
            COLMAX = TEMPIV
  190     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410
          PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          INCRN = INCR + IPLUSN
          PIVOT(INCRN) = INCR + IPVT
          IRWBLK = IPLUSN - NRWBLK
          IF (IPVT.EQ.IPLUSN) GO TO 240
          DO 200 L=I,NRWBLK
            SWAP = ARRAY(L,IPVT,K)
            ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K)
            ARRAY(L,IPLUSN,K) = SWAP
  200     CONTINUE
          IPVBLK = IPVT - NRWBLK
          IF (K.EQ.NBLOKS) GO TO 220
          DO 210 L=1,NRWBLK
            SWAP = ARRAY(L,IPVBLK,KPLUS1)
            ARRAY(L,IPVBLK,KPLUS1) = ARRAY(L,IRWBLK,KPLUS1)
            ARRAY(L,IRWBLK,KPLUS1) = SWAP
  210     CONTINUE
          GO TO 240
  220     CONTINUE
          DO 230 L=1,NRWBOT
            SWAP = BOTBLK(L,IPVBLK)
            BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK)
            BOTBLK(L,IRWBLK) = SWAP
  230     CONTINUE
  240     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          COLPIV = ARRAY(I,IPLUSN,K)
          DO 300 J=LOOP,NCLBLK
            COLMLT = ARRAY(I,J,K)/COLPIV
            ARRAY(I,J,K) = COLMLT
            IF (I.EQ.NRWBLK) GO TO 260
            DO 250 L=IPLUS1,NRWBLK
              ARRAY(L,J,K) = ARRAY(L,J,K) - COLMLT*ARRAY(L,IPLUSN,K)
  250       CONTINUE
  260       CONTINUE
            JRWBLK = J - NRWBLK
            IF (K.EQ.NBLOKS) GO TO 280
            DO 270 L=1,NRWBLK
              ARRAY(L,JRWBLK,KPLUS1) = ARRAY(L,JRWBLK,KPLUS1) -
     *         COLMLT*ARRAY(L,IRWBLK,KPLUS1)
  270       CONTINUE
            GO TO 300
  280       CONTINUE
            DO 290 L=1,NRWBOT
              BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK) -
     *         COLMLT*BOTBLK(L,IRWBLK)
  290       CONTINUE
  300     CONTINUE
  310   CONTINUE
        INCR = INCR + NRWBLK
  320 CONTINUE
C
C***************************************************************
C
C          ****  FINALLY, IN BOTBLK....
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          ***  APPLY NRWBOT ROW ELIMINATIONS WITH ROW
C                  PIVOTING....
C
C               IF BOT HAS JUST ONE ROW GO TO 500
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWBOT.EQ.1) GO TO 400
      DO 390 J=NRWTP1,NVRLP0
        JPLUS1 = J + 1
        JMINN = J - NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IPVT = JMINN
        ROWMAX = DABS(BOTBLK(JMINN,J))
        LOOP = JMINN + 1
        DO 330 I=LOOP,NRWBOT
          TEMPIV = DABS(BOTBLK(I,J))
          IF (TEMPIV.LE.ROWMAX) GO TO 330
          IPVT = I
          ROWMAX = TEMPIV
  330   CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410
        PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCRJ = INCR + J
        PIVOT(INCRJ) = INCR + IPVT + NRWTOP
        IF (IPVT.EQ.JMINN) GO TO 350
        DO 340 L=J,NOVRLP
          SWAP = BOTBLK(IPVT,L)
          BOTBLK(IPVT,L) = BOTBLK(JMINN,L)
          BOTBLK(JMINN,L) = SWAP
  340   CONTINUE
  350   CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        ROWPIV = BOTBLK(JMINN,J)
        DO 360 I=LOOP,NRWBOT
          BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV
  360   CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 380 L=JPLUS1,NOVRLP
          ROWMLT = BOTBLK(JMINN,L)
          DO 370 I=LOOP,NRWBOT
            BOTBLK(I,L) = BOTBLK(I,L) - ROWMLT*BOTBLK(I,J)
  370     CONTINUE
  380   CONTINUE
  390 CONTINUE
  400 CONTINUE
C
C***************************************************************
C
C          DONE PROVIDED THE LAST ELEMENT IS NOT ZERO
C
C***************************************************************
C
      IF (PIVMAX+DABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN
C
C***************************************************************
C
C       ****  MATRIX IS SINGULAR - SET IFLAG = - 1.
C                                  TERMINATE AT 1000.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
  410 CONTINUE
      IFLAG = -1
      RETURN
      END
      SUBROUTINE CRSLVE(TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK, NCLBLK,  CRS   10
     * NBLOKS, BOTBLK, NRWBOT, PIVOT, B, X)
C
C***************************************************************
C
C  C R S L V E  SOLVES THE LINEAR SYSTEM
C                       A*X = B
C  USING THE DECOMPOSITION ALREADY GENERATED IN  C R D C M P.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY  ...
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         THE PIVOT VECTOR FROM  C R D C M P
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR
C
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, ARRAY, BOTBLK, X, B
      DOUBLE PRECISION DOTPRD, XJ, XINCRJ, BINCRJ, SWAP
      INTEGER PIVOT(1)
      DIMENSION TOPBLK(NRWTOP,1), ARRAY(NRWBLK,NCLBLK,1),
     * BOTBLK(NRWBOT,1), B(1), X(1)
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
      NRWTP1 = NRWTOP + 1
      NRWBK1 = NRWBLK + 1
      NVRLP1 = NOVRLP + 1
      NRWTP0 = NRWTOP - 1
      NRWBT1 = NRWBOT + 1
      NROWEL = NRWBLK - NRWTOP
      NRWEL1 = NROWEL + 1
      NVRLP0 = NOVRLP - 1
      NBLKS1 = NBLOKS + 1
      NBKTOP = NRWBLK + NRWTOP
C
C***************************************************************
C
C               ****  FORWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST, IN TOPBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 30 J=1,NRWTOP
        X(J) = B(J)/TOPBLK(J,J)
        IF (J.EQ.NRWTOP) GO TO 20
        XJ = -X(J)
        LOOP = J + 1
        DO 10 I=LOOP,NRWTOP
          B(I) = B(I) + TOPBLK(I,J)*XJ
   10   CONTINUE
   20   CONTINUE
   30 CONTINUE
C
C       ********************************************************
C
C          ***  IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INCR = 0
      DO 120 K=1,NBLOKS
        INCRTP = INCR + NRWTOP
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 50 J=1,NRWTOP
          INCRJ = INCR + J
          XINCRJ = -X(INCRJ)
          DO 40 I=1,NRWBLK
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
   40     CONTINUE
   50   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 80 J=NRWTP1,NRWBLK
          INCRJ = INCR + J
          JPIVOT = PIVOT(INCRJ)
          IF (JPIVOT.EQ.INCRJ) GO TO 60
          SWAP = B(INCRJ)
          B(INCRJ) = B(JPIVOT)
          B(JPIVOT) = SWAP
   60     CONTINUE
          BINCRJ = -B(INCRJ)
          LOOP = J - NRWTP0
          DO 70 I=LOOP,NRWBLK
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*BINCRJ
   70     CONTINUE
   80   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 110 J=NRWBK1,NBKTOP
          INCRJ = INCR + J
          JRWTOP = J - NRWTOP
          X(INCRJ) = B(INCRJ)/ARRAY(JRWTOP,J,K)
          IF (J.EQ.NBKTOP) GO TO 100
          XINCRJ = -X(INCRJ)
          LOOP = J - NRWTP0
          DO 90 I=LOOP,NRWBLK
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
   90     CONTINUE
  100     CONTINUE
  110   CONTINUE
        INCR = INCR + NRWBLK
  120 CONTINUE
C
C       ********************************************************
C
C          ***  FINALLY, IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INCRTP = INCR + NRWTOP
      DO 140 J=1,NRWTOP
        INCRJ = INCR + J
        XINCRJ = -X(INCRJ)
        DO 130 I=1,NRWBOT
          INCRI = INCRTP + I
          B(INCRI) = B(INCRI) + BOTBLK(I,J)*XINCRJ
  130   CONTINUE
  140 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWBOT.EQ.1) GO TO 180
      DO 170 J=NRWTP1,NVRLP0
        INCRJ = INCR + J
        JPIVOT = PIVOT(INCRJ)
        IF (JPIVOT.EQ.INCRJ) GO TO 150
        SWAP = B(INCRJ)
        B(INCRJ) = B(JPIVOT)
        B(JPIVOT) = SWAP
  150   CONTINUE
        BINCRJ = -B(INCRJ)
        LOOP = J - NRWTP0
        DO 160 I=LOOP,NRWBOT
          INCRI = INCRTP + I
          B(INCRI) = B(INCRI) + BOTBLK(I,J)*BINCRJ
  160   CONTINUE
  170 CONTINUE
  180 CONTINUE
C
C***************************************************************
C
C               ****  BACKWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 210 LL=1,NRWBOT
        J = NVRLP1 - LL
        INCRJ = INCR + J
        NRWBTL = NRWBT1 - LL
        X(INCRJ) = B(INCRJ)/BOTBLK(NRWBTL,J)
        IF (LL.EQ.NRWBOT) GO TO 200
        XINCRJ = -X(INCRJ)
        LOOP = NRWBOT - LL
        DO 190 I=1,LOOP
          INCRI = INCRTP + I
          B(INCRI) = B(INCRI) + BOTBLK(I,J)*XINCRJ
  190   CONTINUE
  200   CONTINUE
  210 CONTINUE
C
C       ********************************************************
C
C          ***  THEN IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 300 L=1,NBLOKS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        K = NBLKS1 - L
        INCR = INCR - NRWBLK
        DO 240 L1=NRWEL1,NRWBLK
          I = NRWBLK + NRWEL1 - L1
          IPLUSN = I + NRWTOP
          LOOP = IPLUSN + 1
          INCRN = INCR + IPLUSN
          DOTPRD = X(INCRN)
          DO 220 J=LOOP,NCLBLK
            INCRJ = INCR + J
            DOTPRD = DOTPRD - ARRAY(I,J,K)*X(INCRJ)
  220     CONTINUE
          X(INCRN) = DOTPRD
          IPVTN = PIVOT(INCRN)
          IF (INCRN.EQ.IPVTN) GO TO 230
          SWAP = X(INCRN)
          X(INCRN) = X(IPVTN)
          X(IPVTN) = SWAP
  230     CONTINUE
  240   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCRTP = INCR + NRWTOP
        DO 260 J=NRWBK1,NCLBLK
          INCRJ = INCR + J
          XINCRJ = -X(INCRJ)
          DO 250 I=1,NROWEL
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
  250     CONTINUE
  260   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 290 LL=1,NROWEL
          J = NRWBK1 - LL
          INCRJ = INCR + J
          NRWELL = NRWEL1 - LL
          X(INCRJ) = B(INCRJ)/ARRAY(NRWELL,J,K)
          IF (LL.EQ.NROWEL) GO TO 280
          XINCRJ = -X(INCRJ)
          LOOP = NROWEL - LL
          DO 270 I=1,LOOP
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
  270     CONTINUE
  280     CONTINUE
  290   CONTINUE
  300 CONTINUE
C
C       ********************************************************
C
C          ***  IN TOPBLK FINISH WITH....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 330 L=1,NRWTOP
        I = NRWTP1 - L
        LOOP = I + 1
        DOTPRD = X(I)
        DO 310 J=LOOP,NOVRLP
          DOTPRD = DOTPRD - TOPBLK(I,J)*X(J)
  310   CONTINUE
        X(I) = DOTPRD
        IPVTI = PIVOT(I)
        IF (I.EQ.IPVTI) GO TO 320
        SWAP = X(I)
        X(I) = X(IPVTI)
        X(IPVTI) = SWAP
  320   CONTINUE
  330 CONTINUE
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN   10
C                                                                       MAN   20
C       THIS DRIVER RUNS A TEST PROBLEM USING THE ROUTINE               MAN   30
C       A R C E C O. THE ARRAYS MTR (FOR THE STRUCTURE OF THE           MAN   40
C       LINEAR SYSTEM), AR (THE COEFFICIENT MATRIX), AND                MAN   50
C       B (THE RIGHT HAND SIDE) ARE SET UP IN DATA STATE-               MAN   60
C       MENTS.                                                          MAN   70
C                                                                       MAN   80
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN   90
      DOUBLE PRECISION AR, X, B, ERROR, ERR                             MAN  100
      DIMENSION AR(114), MTR(3,5), B(18), X(18)                         MAN  110
      INTEGER PIVOT(18)                                                 MAN  120
      DATA N /18/, NMBLKS /5/                                           MAN  130
      DATA MTR(1,1), MTR(2,1), MTR(3,1), MTR(1,2), MTR(2,2), MTR(3,2),  MAN  140
     * MTR(1,3), MTR(2,3), MTR(3,3), MTR(1,4), MTR(2,4), MTR(3,4),      MAN  150
     * MTR(1,5), MTR(2,5), MTR(3,5) /2,4,3,4,7,4,5,8,2,3,6,3,4,5,0/     MAN  160
      DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7), AR(8),      MAN  170
     * AR(9), AR(10), AR(11), AR(12), AR(13), AR(14), AR(15), AR(16),   MAN  180
     * AR(17), AR(18), AR(19), AR(20), AR(21), AR(22), AR(23), AR(24)   MAN  190
     * /-1.00D0,-1.00D0,-0.98D0,0.25D0,-0.79D0,-0.87D0,-0.15D0,0.35D0,  MAN  200
     * 0.78D0,-0.82D0,-0.83D0,-0.21D0,0.31D0,0.12D0,-0.98D0,-0.93D0,    MAN  210
     * -0.85D0,-0.01D0,-0.58D0,-0.84D0,0.89D0,0.75D0,0.04D0,0.37D0/     MAN  220
      DATA AR(25), AR(26), AR(27), AR(28), AR(29), AR(30), AR(31),      MAN  230
     * AR(32), AR(33), AR(34), AR(35), AR(36), AR(37), AR(38), AR(39),  MAN  240
     * AR(40), AR(41), AR(42), AR(43), AR(44), AR(45), AR(46), AR(47),  MAN  250
     * AR(48) /-0.69D0,0.32D0,0.87D0,-0.94D0,-0.98D0,-1.00D0,0.38D0,    MAN  260
     * -0.96D0,-0.76D0,-0.53D0,-1.00D0,-1.00D0,-0.99D0,-0.87D0,-0.93D0, MAN  270
     * 0.85D0,0.17D0,-0.91D0,-0.14D0,-0.91D0,-0.39D0,-1.37D0,-0.28D0,   MAN  280
     * -1.00D0/                                                         MAN  290
      DATA AR(49), AR(50), AR(51), AR(52), AR(53), AR(54), AR(55),      MAN  300
     * AR(56), AR(57), AR(58), AR(59), AR(60), AR(61), AR(62), AR(63),  MAN  310
     * AR(64), AR(65), AR(66), AR(67), AR(68), AR(69), AR(70), AR(71),  MAN  320
     * AR(72) /0.10D0,0.79D0,1.29D0,0.90D0,-0.59D0,-0.89D0,-0.71D0,     MAN  330
     * -1.59D0,0.78D0,-0.99D0,-0.68D0,0.39D0,1.10D0,-0.93D0,0.21D0,     MAN  340
     * -0.09D0,-0.99D0,-1.63D0,-0.76D0,-0.73D0,-0.58D0,-0.12D0,-1.01D0, MAN  350
     * 0.48D0/                                                          MAN  360
      DATA AR(73), AR(74), AR(75), AR(76), AR(77), AR(78), AR(79),      MAN  370
     * AR(80), AR(81), AR(82), AR(83), AR(84), AR(85), AR(86), AR(87),  MAN  380
     * AR(88), AR(89), AR(90), AR(91), AR(92), AR(93), AR(94), AR(95),  MAN  390
     * AR(96) /-0.48D0,-0.21D0,-0.75D0,-0.27D0,0.08D0,-0.67D0,-0.24D0,  MAN  400
     * 0.61D0,0.56D0,-0.41D0,0.54D0,-0.99D0,0.40D0,-0.41D0,0.16D0,      MAN  410
     * -0.93D0,0.16D0,-0.16D0,0.70D0,-0.46D0,0.98D0,0.43D0,0.71D0,      MAN  420
     * -0.47D0/                                                         MAN  430
      DATA AR(97), AR(98), AR(99), AR(100), AR(101), AR(102), AR(103),  MAN  440
     * AR(104), AR(105), AR(106), AR(107), AR(108), AR(109), AR(110),   MAN  450
     * AR(111), AR(112), AR(113), AR(114) /-0.25D0,0.89D0,-0.97D0,      MAN  460
     * -0.98D0,-0.92D0,-0.94D0,-0.60D0,-0.73D0,-0.52D0,-0.54D0,-0.30D0, MAN  470
     * 0.07D0,-0.46D0,-1.00D0,0.18D0,0.04D0,-0.58D0,-0.36D0/            MAN  480
      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), MAN  490
     * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18)           MAN  500
     * /-2.92D0,-1.27D0,-1.30D0,-1.17D0,-2.10D0,-4.51D0,-1.71D0,-4.59D0,MAN  510
     * -4.19D0,-0.93D0,-3.31D0,0.52D0,-0.12D0,-0.05D0,-0.98D0,-2.07D0,  MAN  520
     * -2.73D0,-1.95D0/                                                 MAN  530
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN  540
C                                                                       MAN  550
C       THE INPUT MATRIX HAS 5 BLOCKS, AND HAS THE FOLLOWING            MAN  560
C       STRUCTURE, WHERE THE INTEGER ENTRIES ARE THE INDICES            MAN  570
C       OF THE ELEMENTS IN THE ARRAY AR, GIVEN MODULO 100.              MAN  580
C                                                                       MAN  590
C                                                                       MAN  600
C                                                                       MAN  610
C       1  3  5  7                                                      MAN  620
C       2  4  6  8                                                      MAN  630
C          9 13 17 21 25 29 33                                          MAN  640
C         10 14 18 22 26 30 34                                          MAN  650
C         11 15 19 23 27 31 35                                          MAN  660
C         12 16 20 24 28 32 36                                          MAN  670
C                  37 42 47 52 57 62 67 72                              MAN  680
C                  38 43 48 53 58 63 68 73                              MAN  690
C                  39 44 49 54 59 64 69 74                              MAN  700
C                  40 45 50 55 60 65 70 75                              MAN  710
C                  41 46 51 56 61 66 71 76                              MAN  720
C                                    77 80 83 86 89 92                  MAN  730
C                                    78 81 84 87 90 93                  MAN  740
C                                    79 82 85 88 91 94                  MAN  750
C                                             95 99  3  7 11            MAN  760
C                                             96  0  4  8 12            MAN  770
C                                             97  1  5  9 13            MAN  780
C                                             98  2  6 10 14            MAN  790
C                                                                       MAN  800
C                                                                       MAN  810
C       THE RIGHT HAND SIDE IS GIVEN BY :                               MAN  820
C                                                                       MAN  830
C       B =                                                             MAN  840
C                                                                       MAN  850
C       ( -2.92,-1.27,-1.30,-1.17,-2.10,-4.51,-1.71,-4.59,-4.19,        MAN  860
C       -0.93,-3.31, 0.52,-0.12,-0.05,-0.98,-2.07,-2.73,-1.95)          MAN  870
C                                                                       MAN  880
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY :                       MAN  890
C                                                                       MAN  900
C       X = (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)                       MAN  910
C                                                                       MAN  920
C                                                                       MAN  930
C       THE STRUCTURE MATRIX, MTR(1:3,1:5), IS GIVEN BY :               MAN  940
C                                                                       MAN  950
C       MTR =                                                           MAN  960
C                                                                       MAN  970
C                 2   4   5   3   4                                     MAN  980
C                 4   7   8   6   5                                     MAN  990
C                 3   4   2   3   0                                     MAN 1000
C                                                                       MAN 1010
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN 1020
      CALL ARCECO(N, AR, MTR, NMBLKS, PIVOT, B, X, IFLAG)               MAN 1030
      IF (IFLAG.NE.0) GO TO 20                                          MAN 1040
      ERROR = 0.D0                                                      MAN 1050
      DO 10 I=1,18                                                      MAN 1060
        ERR = 1.D0 - X(I)                                               MAN 1070
        ERROR = DMAX1(ERROR,DABS(ERR))                                  MAN 1080
        WRITE (6,99998) X(I), ERR                                       MAN 1090
   10 CONTINUE                                                          MAN 1100
      WRITE (6,99999) ERROR                                             MAN 1110
99999 FORMAT (12H MAX ERROR =, D15.7)                                   MAN 1120
99998 FORMAT (1H , F15.7, D15.7)                                        MAN 1130
      RETURN                                                            MAN 1140
   20 CONTINUE                                                          MAN 1150
      WRITE (6,99997) IFLAG                                             MAN 1160
99997 FORMAT (9H IFLAG = , I3)                                          MAN 1170
      RETURN                                                            MAN 1180
      END                                                               MAN 1190
      SUBROUTINE ARCECO(N, ARRAY, MTRSTR, NMBLKS, PIVOT, B, X, IFLAG)   ARC   10
C
C***************************************************************
C
C  THIS PROGRAM SOLVES THE LINEAR SYSTEM  A*X = B  WHERE  A IS
C  AN ALMOST BLOCK DIAGONAL MATRIX.  THE METHOD IMPLEMENTED IS
C  BASED ON GAUSS ELIMINATION WITH ALTERNATE ROW AND COLUMN
C  ELIMINATION WITH PARTIAL PIVOTING, WHICH PRODUCES A STABLE
C  DECOMPOSITION OF THE MATRIX  A  WITHOUT INTRODUCING FILL-IN.
C
C***************************************************************
C
C    TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE,
C    REMOVE ALL DOUBLE PRECISION STATEMENTS. THERE IS ONE
C    SUCH STATEMENT IN EACH OF :-  ARCECO,ARCEDC,ARCEPR,
C    ARCESL,ARCEFE,ARCEFS,ARCEFM,ARCEBM,ARCEBS,ARCEBE; AND
C    TWO IN ARCEPC.
C    IN ADDITION, REFERENCES TO THE BUILT-IN FUNCTIONS DABS
C    AND DMAX1 MUST BE REPLACED BY ABS, AND AMAX1, RESPECTIVELY,
C        DABS OCCURS TWICE IN ARCEPR (JUST BEFORE THE LOOP
C    DO 20, AND INSIDE THAT LOOP), AND TWICE IN ARCEPC (AGAIN
C    JUST BEFORE THE LOOP DO 20, AND INSIDE THAT LOOP).
C        DMAX1 OCCURS TWICE IN THE PACKAGE, ONCE IN ARCEPR AFTER
C    THE TEST FOR SINGULARTY) AND ONCE IN ARCEPC, ALSO AFTER
C    THE TEST FOR SINGULARITY.
C    FINALLY, ZERO IS INITIALISED TO 0.D0 IN A DATA STATEMENT
C    IN ARCEDC. THIS MUST BE REPLACED BY :
C        DATA ZERO/0.0/
C
C***************************************************************
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C                N     - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM, WHERE
C                           N = SUM(MTRSTR(1,K),K=1,NMBLKS)
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         WHERE
C                          NUMELS = SUM(MTRSTR(1,K)*MTRSTR(2,K);
C                               K=1,NMBLKS).
C                         CONTAINS THE ENTRIES OF THE ALMOST
C                         BLOCK DIAGONAL MATRIX  A  WHOSE BLOCK
C                         STRUCTURE IS GIVEN BY THE INTEGER ARRAY
C                         MTRSTR. THE ELEMENTS OF A ARE STORED BY
C                         COLUMNS, IN BLOCKS CORRESPONDING TO THE
C                         GIVEN STRUCTURE.
C
C               MTRSTR - INTEGER(3,NMBLKS)
C                         DESCRIBES THE BLOCK STRUCTURE OF  A:
C                           MTRSTR(1,K)  = NUMBER OF ROWS IN
C                                            BLOCK K.
C                           MTRSTR(2,K)  = NUMBER OF COLUMNS IN
C                                            BLOCK K.
C                           MTRSTR(3,K)  = NUMBER OF COLUMNS
C                                          OVERLAPPED BY BLOCK K
C                                          AND BLOCK (K+1).
C                        MTRSTR MUST SATISFY SOME RESTRICTIONS.
C                         IN ORDER THAT A BE SQUARE, WE NEED
C                         SUM(MTRSTR(1,K(,K=1,NMBLKS) = N =
C                         SUM((MTRSTR(2,K)-MTRSTR(3,K)),K=1,NMBLKS).
C                        IN ADDITION, TO ENSURE THAT THREE SUCCESS-
C                         IVE BLOCKS DO NOT HAVE COLUMNS IN COMMON,
C                         MTRSTR MUST SATISFY
C                         MTRSTR(3,K-1)+MTRSTR(3,K) <= MTRSTR(2,K),
C                         FOR K = 2,NMBLKS.
C                        FINALLY,  A R C E C O,  SETS
C                         MTRSTR(3,NMBLKS) = 0, IN ARCECD.
C
C               NMBLKS - INTEGER
C                         TOTAL NUMBER OF BLOCKS IN  A
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         CONTAINS THE MODIFIED ALTERNATE ROW
C                         AND COLUMN DECOMPOSITION OF  A (IF
C                         IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR (IF IFLAG = 0)
C
C                IFLAG - INTEGER
C                         =  1,IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
C               *****  AUXILIARY PROGRAMS  *****
C
C            ARCEDC(ARRAY,MTRSTR,NMBLKS,PIVOT,IFLAG)
C               - DECOMPOSES THE MATRIX  A  USING MODIFIED
C                 ALTERNATE ROW AND COLUMN ELIMINATION
C                 WITH PARTIAL PIVOTING, AND IS USED FOR
C                 THIS PURPOSE IN  A R C E C O.
C                 THE ARGUMENTS ARE ALL AS IN A R C E C O.
C
C            ARCESL(ARRAY,MTRSTR,NMBLKS,PIVOT,B,X)
C               - SOLVES THE SYSTEM  A*X = B  ONCE  A  IS
C                 DECOMPOSED.
C                 THE ARGUMENTS ARE ALL AS IN A R C E C O .
C
C***************************************************************
C
C               *****  BLOCK STRUCTURE OF  A  *****
C
C  THE NMBLKS BLOCKS OF A ARE STORED CONSECUTIVELY IN THE ONE
C  DIMENSIONAL MATRIX  ARRAY, THE ENTRIES OF  A  BEING STORED
C  AS FOLLOWS:
C
C       IN ARRAY(1)     THE (1,1) ENTRY OF THE TOP BLOCK,
C
C       IN ARRAY(INDEX) THE (1,1) ENTRY OF THE ITH BLOCK WHERE
C                       INDEX = 1 + SUM(MTRSTR(1,J)*MTRSTR(2,J);
C                               J=1,I-1), I=2,NMBLKS.
C
C***************************************************************
C
C       THE SUBROUTINE  A R C E C O  AUTOMATICALLY SOLVES THE
C  INPUT SYSTEM WHEN IFLAG=0.  A R C E C O  IS CALLED ONLY ONCE
C  FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT
C  HAND SIDES CAN BE OBTAINED BY ONE CALL TO  A R C E C O  AND
C  P-1 CALLS TO ARCESL ONLY. SINCE THE ARRAYS ARRAY AND
C  PIVOT CONTAIN, RESPECTIVELY, THE DECOMPOSITION OF THE GIVEN
C  COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM
C  A R C E C O , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE
C  CALLS TO ARCESL WITH THE SAME RIGHT HAND SIDES. FOR THE
C  SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT
C  MATRIX, THE ARRAY ARRAY MUST BE COPIED BEFORE A CALL
C  TO  A R C E C O .
C
C**********************************************************************
C
C               *****  SAMPLE CALLING PROGRAM  *****
C
C       THE FOLLOWING PROGRAM WILL EXERCISE  A R C E C O , IN
C       THE CASE WHEN THE COEFFICIENT MATRIX IS NONSINGULAR.
C
C       DOUBLE PRECISION AR,X,B,ERROR,ERR
C       DIMENSION AR(114),MTR(3,5),B(18),X(18)
C       INTEGER PIVOT(18)
C       DATA MTR(1,1),MTR(2,1),MTR(3,1),MTR(1,2),MTR(2,2),MTR(3,2),
C    *       MTR(1,3),MTR(2,3),MTR(3,3),MTR(1,4),MTR(2,4),MTR(3,4),
C    *       MTR(1,5),MTR(2,5),MTR(3,5)/
C    *  2,4,3,4,7,4,5,8,2,3,6,3,4,5,0/
C       DATA AR(1),AR(2),AR(3),AR(4),AR(5),AR(6),AR(7),AR(8),AR(9),
C    *       AR(10),AR(11),AR(12),AR(13),AR(14),AR(15),AR(16),AR(17),
C    *       AR(18),AR(19),AR(20),AR(21),AR(22),AR(23),AR(24)/
C    *-1.00D0,-1.00D0,-0.98D0, 0.25D0,-0.79D0,-0.87D0,-0.15D0, 0.35D0,
C    * 0.78D0,-0.82D0,-0.83D0,-0.21D0, 0.31D0, 0.12D0,-0.98D0,-0.93D0,
C    *-0.85D0,-0.01D0,-0.58D0,-0.84D0, 0.89D0, 0.75D0, 0.04D0, 0.37D0/
C       DATA AR(25),AR(26),AR(27),AR(28),AR(29),AR(30),AR(31),AR(32),
C    *       AR(33),AR(34),AR(35),AR(36),AR(37),AR(38),AR(39),AR(40),
C    *       AR(41),AR(42),AR(43),AR(44),AR(45),AR(46),AR(47),AR(48)/
C    *-0.69D0, 0.32D0, 0.87D0,-0.94D0,-0.98D0,-1.00D0, 0.38D0,-0.96D0,
C    *-0.76D0,-0.53D0,-1.00D0,-1.00D0,-0.99D0,-0.87D0,-0.93D0, 0.85D0,
C    * 0.17D0,-0.91D0,-0.14D0,-0.91D0,-0.39D0,-1.37D0,-0.28D0,-1.00D0/
C       DATA AR(49),AR(50),AR(51),AR(52),AR(53),AR(54),AR(55),AR(56),
C    *       AR(57),AR(58),AR(59),AR(60),AR(61),AR(62),AR(63),AR(64),
C    *       AR(65),AR(66),AR(67),AR(68),AR(69),AR(70),AR(71),AR(72)/
C    * 0.10D0, 0.79D0, 1.29D0, 0.90D0,-0.59D0,-0.89D0,-0.71D0,-1.59D0,
C    * 0.78D0,-0.99D0,-0.68D0, 0.39D0, 1.10D0,-0.93D0, 0.21D0,-0.09D0,
C    *-0.99D0,-1.63D0,-0.76D0,-0.73D0,-0.58D0,-0.12D0,-1.01D0, 0.48D0/
C       DATA AR(73),AR(74),AR(75),AR(76),AR(77),AR(78),AR(79),AR(80),
C    *       AR(81),AR(82),AR(83),AR(84),AR(85),AR(86),AR(87),AR(88),
C    *       AR(89),AR(90),AR(91),AR(92),AR(93),AR(94),AR(95),AR(96)/
C    *-0.48D0,-0.21D0,-0.75D0,-0.27D0, 0.08D0,-0.67D0,-0.24D0, 0.61D0,
C    * 0.56D0,-0.41D0, 0.54D0,-0.99D0, 0.40D0,-0.41D0, 0.16D0,-0.93D0,
C    * 0.16D0,-0.16D0, 0.70D0,-0.46D0, 0.98D0, 0.43D0, 0.71D0,-0.47D0/
C       DATA AR(97),AR(98),AR(99),AR(100),AR(101),AR(102),AR(103),
C    *       AR(104),AR(105),AR(106),AR(107),AR(108),AR(109),AR(110),
C    *       AR(111),AR(112),AR(113),AR(114)/
C    *-0.25D0, 0.89D0,-0.97D0,-0.98D0,-0.92D0,-0.94D0,-0.60D0,-0.73D0,
C    *-0.52D0,-0.54D0,-0.30D0, 0.07D0,-0.46D0,-1.00D0, 0.18D0, 0.04D0,
C    *-0.58D0,-0.36D0/
C       DATA B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9),B(10),B(11),
C    *       B(12),B(13),B(14),B(15),B(16),B(17),B(18)/
C    *-2.92D0,-1.27D0,-1.30D0,-1.17D0,-2.10D0,-4.51D0,-1.71D0,-4.59D0,
C    *-4.19D0,-0.93D0,-3.31D0, 0.52D0,-0.12D0,-0.05D0,-0.98D0,-2.07D0,
C    *-2.73D0,-1.95D0/
C       DATA N/18/,NMBLKS/5/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       THE MATRIX AR HAS 5 BLOCKS, WITH THE FOLLOWING STRUCTURE, THE
C       ENTRIES BEING THE INDICES OF THE ELEMENTS OF AR, MODULO 100.
C
C
C       1  3  5  7
C       2  4  6  8
C          9 13 17 21 25 29 33
C         10 14 18 22 26 30 34
C         11 15 19 23 27 31 35
C         12 16 20 24 28 32 36
C                  37 42 47 52 57 62 67 72
C                  38 43 48 53 58 63 68 73
C                  39 44 49 54 59 64 69 74
C                  40 45 50 55 60 65 70 75
C                  41 46 51 56 61 66 71 76
C                                    77 80 83 86 89 92
C                                    78 81 84 87 90 93
C                                    79 82 85 88 91 94
C                                             95 99  3  7 11
C                                             96  0  4  8 12
C                                             97  1  5  9 13
C                                             98  2  6 10 14
C
C
C       THE RIGHT HAND SIDE IS GIVEN BY :
C
C       B = ( -2.92,-1.27,-1.30,-1.17,-2.10,-4.51,-1.71,-4.59,-4.19,
C       -0.93,-3.31, 0.52,-0.12,-0.05,-0.98,-2.07,-2.73,-1.95)
C
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY :
C
C       X = (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
C
C       THE STRUCTURE MATRIX, MTR(1:3,1:5), IS GIVEN BY :
C
C                 2   4   5   3   4
C       MTR =     4   7   8   6   5
C                 3   4   2   3   0
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       CALL ARCECO(N,AR,MTR,NMBLKS,PIVOT,B,X,IFLAG)
C       IF(IFLAG.NE.0)GO TO 1000
C       ERROR = 0.D0
C       DO 10 I=1,18
C          ERR = 1.D0 - X(I)
C          ERROR = DMAX1(ERROR,DABS(ERR))
C          WRITE(5,100)X(I),ERR
C  10   CONTINUE
C       WRITE(5,200)ERROR
C 200   FORMAT(12H MAX ERROR = ,D15.7)
C 100   FORMAT(1H ,F15.7,D15.7)
C       RETURN
C1000   CONTINUE
C       WRITE(5,300)IFLAG
C 300   FORMAT(9H IFLAG =  ,I3)
C       RETURN
C       END
      DOUBLE PRECISION ARRAY, B, X
      INTEGER MTRSTR(3,1), PIVOT(1)
      DIMENSION ARRAY(1), B(1), X(1)
      CALL ARCEDC(N, ARRAY, MTRSTR, NMBLKS, PIVOT, IFLAG)
      IF (IFLAG.NE.0) RETURN
      CALL ARCESL(ARRAY, MTRSTR, NMBLKS, PIVOT, B, X)
      RETURN
      END
      SUBROUTINE ARCEDC(N, ARRAY, MTRSTR, NMBLKS, PIVOT, IFLAG)         ARC   10
C
C***************************************************************
C
C  A R C E D C SUPERVISES THE MODIFIED ALTERNATE ROW AND COLUMN
C  DECOMPOSITION WITH PARTIAL PIVOTING OF THE ALMOST BLOCK
C  DIAGONAL MATRIX  A   STORED IN THE ARRAYS A R R A Y  AND
C  M T R S T R .
C
C***************************************************************
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C                N     - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM, WHERE
C                           N = SUM(MTRSTR(1,K),K=1,NMBLKS)
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         WHERE
C                          NUMELS = SUM(MTRSTR(1,K)*MTRSTR(2,K);
C                               K=1,NMBLKS).
C                         CONTAINS THE ENTRIES OF THE ALMOST
C                         BLOCK DIAGONAL MATRIX  A  WHOSE BLOCK
C                         STRUCTURE IS GIVEN BY THE INTEGER ARRAY
C                         MTRSTR. THE ELEMENTS OF A ARE STORED BY
C                         COLUMNS, IN BLOCKS CORRESPONDING TO THE
C                         GIVEN STRUCTURE.
C               MTRSTR - INTEGER(3,NMBLKS)
C                         DESCRIBES THE BLOCK STRUCTURE OF  A:
C                           MTRSTR(1,K)  = NUMBER OF ROWS IN
C                                            BLOCK K.
C                           MTRSTR(2,K)  = NUMBER OF COLUMNS IN
C                                            BLOCK K.
C                           MTRSTR(3,K)  = NUMBER OF COLUMNS
C                                          OVERLAPPED BY BLOCK K
C                                          AND BLOCK (K+1).
C                        MTRSTR MUST SATISFY SOME RESTRICTIONS.
C                         IN ORDER THAT A BE SQUARE, WE NEED
C                         SUM(MTRSTR(1,K(,K=1,NMBLKS) = N =
C                         SUM((MTRSTR(2,K)-MTRSTR(3,K)),K=1,NMBLKS).
C                        IN ADDITION, TO ENSURE THAT THREE SUCCESS-
C                         IVE BLOCKS DO NOT HAVE COLUMNS IN COMMON,
C                         MTRSTR MUST SATISFY
C                         MTRSTR(3,K-1)+MTRSTR(3,K) <= MTRSTR(2,K),
C                         FOR K = 2,NMBLKS.
C                        FINALLY,  A R C E C O,  SETS
C                         MTRSTR(3,NMBLKS) = 0, IN ARCECD.
C
C               NMBLKS - INTEGER
C                         TOTAL NUMBER OF BLOCKS
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         CONTAINS THE MODIFIED ALTERNATE ROW
C                         AND COLUMN DECOMPOSITION OF  A (IF
C                         IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                IFLAG - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
C               *****  AUXILIARY PROGRAMS  *****
C
C       ARCEPR(BLOCK,NRWBLK,NCLBLK,NRWPIV,PIVOT,PIVMAX,IFLAG)
C        CARRIES OUT THE ROW ELIMINATIONS
C
C       ARCEPC(TOPBLK,NRWTOP,NOVRLP,BOTBLK,NRWBOT,NCLPIV,
C            PIVOT,PIVMAX,IFLAG)
C        CARRIES OUT THE COLUMN ELIMINATIONS
C
C***************************************************************
C
      DOUBLE PRECISION ARRAY, PIVMAX, ZERO
      INTEGER PIVOT(1)
      DIMENSION ARRAY(1), MTRSTR(3,1)
      DATA ZERO /0.D0/
C
C***************************************************************
C
C       ****  CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C            IF PARAMETERS ARE INVALID THEN TERMINATE AT 7;
C                                      ELSE CONTINUE AT 8.
C
C***************************************************************
C
C
      MTRSTR(3,NMBLKS) = 0
      DO 10 K=2,NMBLKS
        IF (MTRSTR(3,K-1)+MTRSTR(3,K).GT.MTRSTR(2,K)) GO TO 30
   10 CONTINUE
      ISUM1 = 0
      ISUM2 = 0
      DO 20 K=1,NMBLKS
        ISUM1 = ISUM1 + MTRSTR(1,K)
        ISUM2 = ISUM2 + MTRSTR(2,K) - MTRSTR(3,K)
   20 CONTINUE
      IF (ISUM1.NE.ISUM2) GO TO 30
      IF (ISUM1.NE.N) GO TO 30
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE ACCEPTABLE - CONTINUE AT 8
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      GO TO 40
   30 CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C        PARAMETERS ARE INVALID.  SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = 1
      RETURN
   40 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C  INTERNAL PARAMETERS:
C       C
C       C       INDEX1: POINTER TO THE ELEMENT IN THE COLUMN
C       C               WHERE ROW PIVOTING STARTS.
C       C
C       C       INDEX2: POINTER TO THE ELEMENT IN THE COLUMN
C       C               WHERE COLUMN PIVOTING STARTS.
C       C
C       C       INDEX3: POINTER TO 1ST ELEMENT IN 1ST COLUMN
C       C               OF NEXT BLOCK.
C       C
C       C       INDPIV: POINTER TO 1ST ELEMENT OF BLOCK OF PIVOT
C       C
C       C       NRWBLK: NUMBER OF ROWS IN BLOCK.
C       C
C       C       NRWBK2: NUMBER OF ROWS IN NEXT BLOCK.
C       C
C       C       NRWPIV: NUMBER OF ROW ELIMINATIONS.
C       C
C       C       NCLBLK: NUMBER OF COLUMNS IN BLOCK TO BE
C       C               ROW PIVOTED.
C       C
C       C       NCLPIV: NUMBER OF COLUMN ELIMINATIONS.
C       C
C       C       NOVRLP: NUMBER OF COLUMNS OVERLAPPED BY THE
C       C               CURRENT BLOCK AND THE NEXT BLOCK.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      PIVMAX = ZERO
      IFLAG = 0
      INDEX1 = 1
      INDPIV = 1
      NRWBLK = MTRSTR(1,1)
      NCLBLK = MTRSTR(2,1)
      NOVRLP = MTRSTR(3,1)
      NRWPIV = NCLBLK - NOVRLP
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       CALL ARCEPR TO PERFORM NRWPIV ROW ELIMINATIONS
C       C               ON TOP BLOCK.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWPIV.GT.0) CALL ARCEPR(ARRAY(INDEX1), NRWBLK, NCLBLK,
     * NRWPIV, PIVOT(INDPIV), PIVMAX, IFLAG)
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       IF MATRIX SINGULAR RETURN.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (IFLAG.LT.0) RETURN
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       NOW DO DECOMPOSITION PROCEEDING ONE BLOCK AT A
C       C               TIME.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 70 K=2,NMBLKS
        INDPIV = INDPIV + NRWPIV
        INDEX2 = INDEX1 + NRWBLK*NRWPIV
        INDEX3 = INDEX2 + NRWBLK*NOVRLP
        NCLPIV = NRWBLK - NRWPIV
        NRWBK2 = MTRSTR(1,K)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C  CALL ARCEPC TO PERFORM NCLPIV COLUMN ELIMINATIONS.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.EQ.0) GO TO 50
        CALL ARCEPC(ARRAY(INDEX2), NRWBLK, NOVRLP, ARRAY(INDEX3),
     *   NRWBK2, NCLPIV, PIVOT(INDPIV), PIVMAX, IFLAG)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C             C
C             C    IF MATRIX IS SINGULAR RETURN.
C             C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (IFLAG.LT.0) RETURN
   50   CONTINUE
        NRWBLK = NRWBK2
        INDEX1 = INDEX3 + NRWBLK*NCLPIV
        NCLBLK = MTRSTR(2,K) - NCLPIV
        NOVRLP = MTRSTR(3,K)
        NRWPIV = NCLBLK - NOVRLP
        INDPIV = INDPIV + NCLPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEPR TO PERFORM NRWPIV ROW ELIMINATIONS.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NRWPIV.EQ.0) GO TO 60
        CALL ARCEPR(ARRAY(INDEX1), NRWBLK, NCLBLK, NRWPIV,
     *   PIVOT(INDPIV), PIVMAX, IFLAG)
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C             C
C             C    IF MATRIX IS SINGULAR RETURN.
C             C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (IFLAG.LT.0) RETURN
   60   CONTINUE
   70 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEPR(BLOCK, NRWBLK, NCLBLK, NRWPIV, PIVOT, PIVMAX,   ARC   10
     * IFLAG)
C
C***************************************************************
C
C  A R C E P R  PERFORMS  NRWPIV  ROW ELIMINATIONS ON THE MATRIX
C               BLOCK
C
C***************************************************************
C
      INTEGER PIVOT(NRWBLK)
      DOUBLE PRECISION BLOCK, ROWMAX, PIVMAX, TEMPIV, ROWPIV, SWAP
      DIMENSION BLOCK(NRWBLK,NCLBLK)
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C        PERFORM NRWPIV ROW ELIMINATIONS...
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 90 J=1,NRWPIV
        JPLUS1 = J + 1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    DETERMINE ROW PIVOT AND PIVOT INDEX
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        MAX = J
        ROWMAX = DABS(BLOCK(J,J))
        IF (J.EQ.NRWBLK) GO TO 30
        DO 20 I1=JPLUS1,NRWBLK
          TEMPIV = DABS(BLOCK(I1,J))
          IF (TEMPIV.LE.ROWMAX) GO TO 10
          ROWMAX = TEMPIV
          MAX = I1
   10     CONTINUE
   20   CONTINUE
   30   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    TEST FOR SINGULARITY:
C          C
C          C             IF SINGULAR THEN TERMINATE AT 90;
C          C                         ELSE CONTINUE.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 100
        PIVMAX = DMAX1(PIVMAX,ROWMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    IF NECESSARY INTERCHANGE ROWS
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        PIVOT(J) = MAX
        IF (J.EQ.MAX) GO TO 50
        DO 40 J1=J,NCLBLK
          SWAP = BLOCK(MAX,J1)
          BLOCK(MAX,J1) = BLOCK(J,J1)
          BLOCK(J,J1) = SWAP
   40   CONTINUE
   50   CONTINUE
        IF (J.EQ.NRWBLK) RETURN
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    COMPUTE THE MULTIPLIERS
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        ROWPIV = BLOCK(J,J)
        DO 60 I1=JPLUS1,NRWBLK
          BLOCK(I1,J) = BLOCK(I1,J)/ROWPIV
   60   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    PERFORM ROW ELIMINATIONS WITH COLUMN INDEXING
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 80 J1=JPLUS1,NCLBLK
          DO 70 L1=JPLUS1,NRWBLK
            BLOCK(L1,J1) = BLOCK(L1,J1) - BLOCK(L1,J)*BLOCK(J,J1)
   70     CONTINUE
   80   CONTINUE
   90 CONTINUE
      RETURN
  100 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       MATRIX IS SINGULAR - SET IFLAG = -1.
C       C                            TERMINATE AT 90.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = -1
      RETURN
      END
      SUBROUTINE ARCEPC(TOPBLK, NRWTOP, NOVRLP, BOTBLK, NRWBOT, NCLPIV, ARC   10
     * PIVOT, PIVMAX, IFLAG)
C
C***************************************************************
C
C     A R C E P C PERFORMS NCLPIV COLUMN ELIMINATIONS ON THE
C                 MATRICES TOPBLK AND BOTBLK
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, BOTBLK, COLMAX, PIVMAX, COLMLT
      DOUBLE PRECISION TEMPIV, SWAP
      INTEGER PIVOT(NRWTOP)
      DIMENSION TOPBLK(NRWTOP,NOVRLP), BOTBLK(NRWBOT,NOVRLP)
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       PERFORM THE COLUMN ELIMINATIONS ON A LOOP.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 110 J=1,NCLPIV
        I = NRWTOP - NCLPIV + J
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    DETERMINE COLUMN PIVOT AND PIVOT INDEX
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        MAX = J
        COLMAX = DABS(TOPBLK(I,J))
        IF (J.EQ.NOVRLP) GO TO 30
        JPLUS1 = J + 1
        DO 20 J1=JPLUS1,NOVRLP
          TEMPIV = DABS(TOPBLK(I,J1))
          IF (TEMPIV.LE.COLMAX) GO TO 10
          COLMAX = TEMPIV
          MAX = J1
   10     CONTINUE
   20   CONTINUE
   30   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    TEST FOR SINGULARITY:
C          C
C          C             IF SINGULAR THEN TERMINATE AT 110;
C          C                         ELSE CONTINUE.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 120
        PIVMAX = DMAX1(PIVMAX,COLMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    IF NECESSARY INTERCHANGE COLUMNS
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        PIVOT(J) = MAX
        IF (J.EQ.MAX) GO TO 60
        DO 40 I1=I,NRWTOP
          SWAP = TOPBLK(I1,J)
          TOPBLK(I1,J) = TOPBLK(I1,MAX)
          TOPBLK(I1,MAX) = SWAP
   40   CONTINUE
        DO 50 I2=1,NRWBOT
          SWAP = BOTBLK(I2,J)
          BOTBLK(I2,J) = BOTBLK(I2,MAX)
          BOTBLK(I2,MAX) = SWAP
   50   CONTINUE
   60   CONTINUE
        IF (J.EQ.NOVRLP) RETURN
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    COMPUTE MULTIPLIERS AND PERFORM COLUMN
C          C            ELIMINATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 100 J1=JPLUS1,NOVRLP
          COLMLT = TOPBLK(I,J1)/TOPBLK(I,J)
          TOPBLK(I,J1) = COLMLT
          IF (I.EQ.NRWTOP) GO TO 80
          IPLUS1 = I + 1
          DO 70 L1=IPLUS1,NRWTOP
            TOPBLK(L1,J1) = TOPBLK(L1,J1) - COLMLT*TOPBLK(L1,J)
   70     CONTINUE
   80     CONTINUE
          DO 90 L1=1,NRWBOT
            BOTBLK(L1,J1) = BOTBLK(L1,J1) - COLMLT*BOTBLK(L1,J)
   90     CONTINUE
  100   CONTINUE
  110 CONTINUE
      RETURN
  120 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       MATRIX IS SINGULAR - SET IFLAG = -1.
C       C                            TERMINATE AT 110.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = -1
      RETURN
      END
      SUBROUTINE ARCESL(ARRAY, MTRSTR, NMBLKS, PIVOT, B, X)             ARC   10
C
C***************************************************************
C
C  A R C E S L  SUPERVISES THE SOLUTION OF THE LINEAR SYSTEM
C                          A*X = B
C  USING THE DECOMPOSITION OF THE MATRIX  A  ALREADY GENERATED
C  IN  A R C E D C.  IT INVOLVES TWO LOOPS, THE FORWARD LOOP,
C  CONSISTING OF FORWARD SOLUTION, FORWARD MODIFICATION, AND
C  FORWARD ELIMINATION, AND THE BACKWARD LOOP, CONSISTING OF
C  BACKWARD SOLUTION, BACKWARD MODIFICATION, AND BACKWARD
C  ELIMINATION.
C
C***************************************************************
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         WHERE
C                          NUMELS = SUM(MTRSTR(1,K)*MTRSTR(2,K);
C                               K=1,NMBLKS).
C                         OUTPUT FROM A R C E D C
C
C               MTRSTR - INTEGER(3,NMBLKS)
C                         DESCRIBES THE BLOCK STRUCTURE OF  A:
C                           MTRSTR(1,K)  = NUMBER OF ROWS IN
C                                            BLOCK K.
C                           MTRSTR(2,K)  = NUMBER OF COLUMNS IN
C                                            BLOCK K.
C                           MTRSTR(3,K)  = NUMBER OF COLUMNS
C                                          OVERLAPPED BY BLOCK K
C                                          AND BLOCK (K+1).
C
C                        THE LINEAR SYSTEM IS OF ORDER:
C                           N  = SUM(MTRSTR(1,K),K=1,NMBLKS)
C
C               NMBLKS - INTEGER
C                         TOTAL NUMBER OF BLOCKS IN  A
C
C                PIVOT - INTEGER(N)
C                         OUTPUT FROM A R C E D C
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR
C
C***************************************************************
C
C               *****  AUXILIARY PROGRAMS  *****
C
C
C       ARCEFS - PERFORMS FORWARD SOLUTION STEP
C
C       ARCEFM - PERFORMS FORWARD MODIFICATION STEP
C
C       ARCEFE - PERFORMS FORWARD ELIMINATION STEP
C
C       ARCEBS - PERFORMS BACKWARD SOLUTION STEP
C
C       ARCEBM - PERFORMS BACKWARD MODIFICATION STEP
C
C       ARCEBE - PERFORMS BACKWARD ELIMINATION STEP
C
C***************************************************************
C
      DOUBLE PRECISION ARRAY, B, X
      INTEGER PIVOT(1)
      DIMENSION ARRAY(1), MTRSTR(3,1), B(1), X(1)
      INDPIV = 1
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C  INTERNAL PARAMETERS:
C       C
C       C       INDEXA: POINTER TO 1ST ELEMENT OF BLOCK OF  A  .
C       C
C       C       INDEXB: POINTER TO 1ST ELEMENT OF BLOCK OF B.
C       C
C       C       INDPIV,NRWBLK,NRWPIV,NCLBLK,NCLPIV,NOVRLP
C       C       ARE AS IN ARCEDC.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INDEXA = 1
      NRWBLK = MTRSTR(1,1)
      NCLBLK = MTRSTR(2,1)
      NOVRLP = MTRSTR(3,1)
      NRWPIV = NCLBLK - NOVRLP
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       CALL ARCEFE TO PERFORM FORWARD ELIMINATION.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWPIV.GT.0) CALL ARCEFE(ARRAY(INDEXA), NRWBLK, NRWPIV,
     * PIVOT(INDPIV), B(INDPIV))
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       FORWARD LOOP
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 10 K=2,NMBLKS
        INDEXA = INDEXA + NRWBLK*NRWPIV
        NCLPIV = NRWBLK - NRWPIV
        INDPIV = INDPIV + NRWPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEFS TO PERFORM FORWARD SOLUTION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.GT.0) CALL ARCEFS(ARRAY(INDEXA), NRWBLK, NCLPIV,
     *   NOVRLP, B(INDPIV), X(INDPIV))
        INDEXA = INDEXA + NOVRLP*NRWBLK
        NRWBLK = MTRSTR(1,K)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEFM TO PERFORM FORWARD MODIFICATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.GT.0) CALL ARCEFM(ARRAY(INDEXA), NRWBLK, NCLPIV,
     *   B(INDPIV), X(INDPIV))
        INDEXA = INDEXA + NRWBLK*NCLPIV
        NCLBLK = MTRSTR(2,K) - NCLPIV
        NOVRLP = MTRSTR(3,K)
        NRWPIV = NCLBLK - NOVRLP
        INDPIV = INDPIV + NCLPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEFE TO PERFORM FORWARD ELIMINATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NRWPIV.GT.0) CALL ARCEFE(ARRAY(INDEXA), NRWBLK, NRWPIV,
     *   PIVOT(INDPIV), B(INDPIV))
   10 CONTINUE
      INDEXB = INDPIV + NRWPIV - 1
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       BACKWARD LOOP
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 30 LL=2,NMBLKS
        K = NMBLKS - LL + 1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEBM TO PERFORM BACKWARD MODIFICATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NRWPIV.EQ.0) GO TO 20
        IF (NRWPIV.NE.NCLBLK) CALL ARCEBM(ARRAY(INDEXA), NRWBLK,
     *   NCLBLK, NRWPIV, B(INDPIV), X(INDPIV))
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEBS TO PERFORM BACKWARD SOLUTION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        CALL ARCEBS(ARRAY(INDEXA), NRWBLK, NCLBLK, NRWPIV, B(INDPIV),
     *   X(INDPIV))
   20   CONTINUE
        INDEXA = INDEXA - NRWBLK*NCLPIV
        NRWBLK = MTRSTR(1,K)
        NOVRLP = MTRSTR(3,K)
        INDEXA = INDEXA - NRWBLK*NOVRLP
        INDPIV = INDPIV - NCLPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEBE TO PERFORM BACKWARD ELIMINATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.GT.0) CALL ARCEBE(ARRAY(INDEXA), NRWBLK, NCLPIV,
     *   NOVRLP, PIVOT(INDPIV), X(INDPIV))
        NRWPIV = NRWBLK - NCLPIV
        NCLBLK = NOVRLP + NRWPIV
        INDEXA = INDEXA - NRWBLK*NRWPIV
        INDPIV = INDPIV - NRWPIV
        NCLPIV = MTRSTR(2,K) - NCLBLK
   30 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       IF ROW ELIMINATIONS WERE DONE IN TOPBLOCK ,CALL
C       C       ARCEBS TO PERFORM BACKWARD SOLUTION
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWPIV.EQ.0) RETURN
      IF (NRWPIV.NE.NCLBLK) CALL ARCEBM(ARRAY(INDEXA), NRWBLK, NCLBLK,
     * NRWPIV, B(INDPIV), X(INDPIV))
      CALL ARCEBS(ARRAY(INDEXA), NRWBLK, NCLBLK, NRWPIV, B(INDPIV),
     * X(INDPIV))
      RETURN
      END
      SUBROUTINE ARCEFE(BLOCK, NRWBLK, NRWPIV, PIVOT, B)                ARC   10
C
C***************************************************************
C
C  A R C E F E  PERFORMS THE FORWARD ELIMINATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, BI, SWAP
      INTEGER PIVOT(NRWPIV), PIVOTI
      DIMENSION BLOCK(NRWBLK,NRWPIV), B(1)
      DO 30 I=1,NRWPIV
        PIVOTI = PIVOT(I)
        IF (PIVOTI.EQ.I) GO TO 10
        SWAP = B(I)
        B(I) = B(PIVOTI)
        B(PIVOTI) = SWAP
   10   CONTINUE
        IF (I.EQ.NRWBLK) RETURN
        BI = B(I)
        IPLUS1 = I + 1
        DO 20 L=IPLUS1,NRWBLK
          B(L) = B(L) - BLOCK(L,I)*BI
   20   CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEFS(BLOCK, NRWBLK, NCLPIV, NOVRLP, B, X)            ARC   10
C
C***************************************************************
C
C  A R C E F S PERFORMS THE FORWARD SOLUTION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NOVRLP), B(1), X(1)
      DO 20 J=1,NCLPIV
        I = NRWBLK - NCLPIV + J
        X(J) = B(J)/BLOCK(I,J)
        IF (I.EQ.NRWBLK) RETURN
        IPLUS1 = I + 1
        LONG = NRWBLK - I
        JPLUS1 = J + 1
        XJ = X(J)
        DO 10 L=1,LONG
          IPLUSL = I + L
          JPLUSL = J + L
          B(JPLUSL) = B(JPLUSL) - BLOCK(IPLUSL,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEFM(BLOCK, NRWBLK, NCLPIV, B, X)                    ARC   10
C
C***************************************************************
C
C  A R C E F M  PERFORMS THE FORWARD MODIFICATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NCLPIV), B(1), X(1)
      NCLPV1 = NCLPIV + 1
      DO 20 J=1,NCLPIV
        XJ = X(J)
        DO 10 L=1,NRWBLK
          NCLPVL = NCLPIV + L
          B(NCLPVL) = B(NCLPVL) - BLOCK(L,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEBM(BLOCK, NRWBLK, NCLBLK, NRWPIV, B, X)            ARC   10
C
C***************************************************************
C
C  A R C E B M  PERFORMS THE BACKWARD MODIFICATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NCLBLK), B(1), X(1)
      NRWPV1 = NRWPIV + 1
      DO 20 J=NRWPV1,NCLBLK
        XJ = X(J)
        DO 10 L=1,NRWPIV
          B(L) = B(L) - BLOCK(L,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEBS(BLOCK, NRWBLK, NCLBLK, NRWPIV, B, X)            ARC   10
C
C***************************************************************
C
C  A R C E B S PERFORMS THE BACKWARD SOLUTION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NCLBLK), B(1), X(1)
      DO 20 NJ=1,NRWPIV
        J = NRWPIV - NJ + 1
        X(J) = B(J)/BLOCK(J,J)
        IF (J.EQ.1) RETURN
        JMIN1 = J - 1
        XJ = X(J)
        DO 10 L=1,JMIN1
          B(L) = B(L) - BLOCK(L,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEBE(BLOCK, NRWBLK, NCLPIV, NOVRLP, PIVOT, X)        ARC   10
C
C***************************************************************
C
C  A R C E B E PERFORMS THE BACKWARD ELIMINATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, X, DOTPRD, SWAP
      INTEGER PIVOT(NRWBLK), PIVOTJ
      DIMENSION BLOCK(NRWBLK,NOVRLP), X(1)
      DO 40 NJ=1,NCLPIV
        J = NCLPIV + 1 - NJ
        I = NRWBLK + 1 - NJ
        DOTPRD = X(J)
        IF (J.EQ.NOVRLP) GO TO 20
        JPLUS1 = J + 1
        DO 10 J1=JPLUS1,NOVRLP
          DOTPRD = DOTPRD - X(J1)*BLOCK(I,J1)
   10   CONTINUE
   20   CONTINUE
        X(J) = DOTPRD
        PIVOTJ = PIVOT(J)
        IF (PIVOTJ.EQ.J) GO TO 30
        SWAP = X(PIVOTJ)
        X(PIVOTJ) = X(J)
        X(J) = SWAP
   30   CONTINUE
   40 CONTINUE
      RETURN
      END
C***********************************************************************MAN   10
C                                                                       MAN   20
C   THIS DRIVER RUNS A TEST PROBLEM USING C O L R O W.                  MAN   30
C   THE ARRAYS TOP, AR, BOT, AND B ARE SET UP IN DATA                   MAN   40
C   STATEMENTS, AS ARE THE PARAMETERS NRWTOP, NRWBLK, NCLBLK,           MAN   50
C   NBLOKS, NRWBOT, WHICH DEFINE THE STRUCTURE OF THE COEFF-            MAN   60
C   ICIENT MATRIX.                                                      MAN   70
C                                                                       MAN   80
C***********************************************************************MAN   90
C                                                                       MAN  100
      DOUBLE PRECISION TOP, AR, BOT, B, X                               MAN  110
      DOUBLE PRECISION ERROR, ERR                                       MAN  120
      DIMENSION TOP(2,4), AR(4,8,2), BOT(2,4), B(12), X(12)             MAN  130
      INTEGER PIVOT(12)                                                 MAN  140
      DATA N, NRWTOP, NOVRLP, NRWBLK, NCLBLK, NBLOKS, NRWBOT            MAN  150
     * /12,2,4,4,8,2,2/                                                 MAN  160
      DATA TOP(1,1), TOP(1,2), TOP(1,3), TOP(1,4), TOP(2,1), TOP(2,2),  MAN  170
     * TOP(2,3), TOP(2,4) /0.0D0,-0.98D0,-0.79D0,-0.15D0,-1.00D0,0.25D0,MAN  180
     * -0.87D0,0.35D0/                                                  MAN  190
      DATA AR(1,1,1), AR(1,2,1), AR(1,3,1), AR(1,4,1), AR(1,5,1),       MAN  200
     * AR(1,6,1), AR(1,7,1), AR(1,8,1) /0.78D0,0.31D0,-0.85D0,0.89D0,   MAN  210
     * -0.69D0,-0.98D0,-0.76D0,-0.82D0/                                 MAN  220
      DATA AR(2,1,1), AR(2,2,1), AR(2,3,1), AR(2,4,1), AR(2,5,1),       MAN  230
     * AR(2,6,1), AR(2,7,1), AR(2,8,1) /0.12D0,-0.01D0,0.75D0,0.32D0,   MAN  240
     * -1.00D0,-0.53D0,-0.83D0,-0.98D0/                                 MAN  250
      DATA AR(3,1,1), AR(3,2,1), AR(3,3,1), AR(3,4,1), AR(3,5,1),       MAN  260
     * AR(3,6,1), AR(3,7,1), AR(3,8,1) /-0.58D0,0.04D0,0.87D0,0.38D0,   MAN  270
     * -1.00D0,-0.21D0,-0.93D0,-0.84D0/                                 MAN  280
      DATA AR(4,1,1), AR(4,2,1), AR(4,3,1), AR(4,4,1), AR(4,5,1),       MAN  290
     * AR(4,6,1), AR(4,7,1), AR(4,8,1) /-0.21D0,-0.91D0,-0.09D0,-0.62D0,MAN  300
     * -1.99D0,-1.12D0,-1.21D0,0.07D0/                                  MAN  310
      DATA AR(1,1,2), AR(1,2,2), AR(1,3,2), AR(1,4,2), AR(1,5,2),       MAN  320
     * AR(1,6,2), AR(1,7,2), AR(1,8,2) /0.78D0,-0.93D0,-0.76D0,0.48D0,  MAN  330
     * -0.87D0,-0.14D0,-1.00D0,-0.59D0/                                 MAN  340
      DATA AR(2,1,2), AR(2,2,2), AR(2,3,2), AR(2,4,2), AR(2,5,2),       MAN  350
     * AR(2,6,2), AR(2,7,2), AR(2,8,2) /-0.99D0,0.21D0,-0.73D0,-0.48D0, MAN  360
     * -0.93D0,-0.91D0,0.10D0,-0.89D0/                                  MAN  370
      DATA AR(3,1,2), AR(3,2,2), AR(3,3,2), AR(3,4,2), AR(3,5,2),       MAN  380
     * AR(3,6,2), AR(3,7,2), AR(3,8,2) /-0.68D0,-0.09D0,-0.58D0,-0.21D0,MAN  390
     * 0.85D0,-0.39D0,0.79D0,-0.71D0/                                   MAN  400
      DATA AR(4,1,2), AR(4,2,2), AR(4,3,2), AR(4,4,2), AR(4,5,2),       MAN  410
     * AR(4,6,2), AR(4,7,2), AR(4,8,2) /0.39D0,-0.99D0,-0.12D0,-0.75D0, MAN  420
     * -0.68D0,-0.99D0,0.50D0,-0.88D0/                                  MAN  430
      DATA BOT(1,1), BOT(1,2), BOT(1,3), BOT(1,4), BOT(2,1), BOT(2,2),  MAN  440
     * BOT(2,3), BOT(2,4) /0.71D0,-0.64D0,0.0D0,0.48D0,0.08D0,100.0D0,  MAN  450
     * 50.00D0,15.00D0/                                                 MAN  460
      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), MAN  470
     * B(11), B(12) /-1.92D0,-1.27D0,-2.12D0,-2.16D0,-2.27D0,-6.08D0,   MAN  480
     * -3.03D0,-4.62D0,-1.02D0,-3.52D0,.55D0,165.08D0/                  MAN  490
C                                                                       MAN  500
C***********************************************************************MAN  510
C                                                                       MAN  520
C   THE INPUT MATRIX IS GIVEN BY:                                       MAN  530
C                                                                       MAN  540
C  0.0  -0.98 -0.79 -0.15                                               MAN  550
C -1.00  0.25 -0.87  0.35                                               MAN  560
C  0.78  0.31 -0.85  0.89 -0.69 -0.98 -0.76 -0.82                       MAN  570
C  0.12 -0.01  0.75  0.32 -1.00 -0.53 -0.83 -0.98                       MAN  580
C -0.58  0.04  0.87  0.38 -1.00 -0.21 -0.93 -0.84                       MAN  590
C -0.21 -0.91 -0.09 -0.62 -1.99 -1.12 -1.21  0.07                       MAN  600
C                          0.78 -0.93 -0.76  0.48 -0.87 -0.14 -1.00 -0.5MAN  610
C                         -0.99  0.21 -0.73 -0.48 -0.93 -0.91  0.10 -0.8MAN  620
C                         -0.68 -0.09 -0.58 -0.21  0.85 -0.39  0.79 -0.7MAN  630
C                          0.39 -0.99 -0.12 -0.75 -0.68 -0.99  0.50 -0.8MAN  640
C                                                  0.71 -0.64  0.0   0.4MAN  650
C                                                  0.08 100.0 50.00 15.0MAN  660
C                                                                       MAN  670
C       THE RIGHT HAND SIDE IS GIVEN BY:                                MAN  680
C                                                                       MAN  690
C         B = (-1.92,-1.27,-2.12,-2.16,-2.27,-6.08,-3.03,-4.62,         MAN  700
C              -1.02,-3.52,0.55,165.08)                                 MAN  710
C                                                                       MAN  720
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY;                        MAN  730
C                                                                       MAN  740
C          X = (1,1,1,1,1,1,1,1,1,1,1,1)                                MAN  750
C                                                                       MAN  760
C***********************************************************************MAN  770
C                                                                       MAN  780
      CALL COLROW(N, TOP, NRWTOP, NOVRLP, AR, NRWBLK, NCLBLK, NBLOKS,   MAN  790
     * BOT, NRWBOT, PIVOT, B, X, IFLAG)                                 MAN  800
      IF (IFLAG.NE.0) GO TO 20                                          MAN  810
      ERROR = 0.D0                                                      MAN  820
      DO 10 I=1,N                                                       MAN  830
        ERR = 1.D0 - X(I)                                               MAN  840
        ERROR = DMAX1(ERROR,DABS(ERR))                                  MAN  850
        WRITE (6,99998) X(I), ERR                                       MAN  860
   10 CONTINUE                                                          MAN  870
      WRITE (6,99999) ERROR                                             MAN  880
99999 FORMAT (12H MAX ERROR =, D15.7)                                   MAN  890
99998 FORMAT (1H , F15.7, D15.7)                                        MAN  900
      RETURN                                                            MAN  910
   20 CONTINUE                                                          MAN  920
      WRITE (6,99997) IFLAG                                             MAN  930
99997 FORMAT (9H IFLAG = , I3)                                          MAN  940
      RETURN                                                            MAN  950
      END                                                               MAN  960
      SUBROUTINE COLROW(N, TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK,       COL   10
     * NCLBLK, NBLOKS, BOTBLK, NRWBOT, PIVOT, B, X, IFLAG)
C
C***************************************************************
C
C  THIS PROGRAM SOLVES THE LINEAR SYSTEM  A*X = B  WHERE  A IS
C  AN ALMOST BLOCK DIAGONAL MATRIX OF THE FORM
C
C               TOPBLK
C               ARRAY(1)
C                     ARRAY(2)
C                          .
C                             .
C                                .
C                                   .
C                                    ARRAY(NBLOKS)
C                                           BOTBLK
C
C  WHERE
C           TOPBLK IS  NRWTOP  BY NOVRLP
C           ARRAY(K), K=1,NBLOKS, ARE NRWBLK BY NRWBLK+NOVRLP
C           BOTBLK IS NRWBOT BY NOVRLP,
C  AND
C           NOVRLP = NRWTOP + NRWBOT
C  WITH
C           NOVRLP.LE.NRWBLK .
C
C  THE LINEAR SYSTEM IS OF ORDER  N = NBLOKS*NRWBLK + NOVRLP.
C
C  THE METHOD IMPLEMENTED IS BASED ON GAUSS ELIMINATION WITH
C  ALTERNATE ROW AND COLUMN ELIMINATION WITH PARTIAL PIVOTING,
C  WHICH PRODUCES A STABLE DECOMPOSITION OF THE MATRIX  A
C  WITHOUT INTRODUCING FILL-IN.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE, REMOVE
C  ALL DOUBLE PRECISION STATEMENTS.  THERE IS ONE SUCH STATEMENT
C  IN C O L R O W, THREE IN C R D C M P, AND TWO IN C R S O L V.
C  IN ADDITION, REFERENCES TO BUILT-IN FUNCTIONS DABS AND DMAX1
C  MUST BE REPLACED BY ABS AND AMAX1, RESPECTIVELY.  DABS OCCURS
C  NINE TIMES, IN C R D C M P.  DMAX1 OCCURS FOUR TIMES, IN
C  C R D C M P.  FINALLY, ZERO IS INITIALISED TO 0.D0 IN A
C  DATA STATEMENT IN C R D C M P.  THIS MUST BE REPLACED BY:
C               DATA ZERO/0.0/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR (IF IFLAG = 0)
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  AUXILIARY PROGRAMS  *****
C
C       CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,IFLAG)
C            - DECOMPOSES THE MATRIX  A  USING MODIFIED
C              ALTERNATE ROW AND COLUMN ELIMINATON WITH
C              PARTIAL PIVOTING, AND IS USED FOR THIS
C              PURPOSE IN C O L R O W.
C              THE ARGUMENTS ARE AS IN C O L R O W.
C
C       CRSLVE(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,B,X)
C            - SOLVES THE SYSTEM A*X = B ONCE A IS DECOMPOSED.
C              THE ARGUMENTS ARE ALLAS IN C O L R O W.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       THE SUBROUTINE  C O L R O W  AUTOMATICALLY SOLVES THE
C  INPUT SYSTEM WHEN IFLAG=0.  C O L R O W  IS CALLED ONLY ONCE
C  FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT
C  HAND SIDES CAN BE OBTAINED BY ONE CALL TO  C O L R O W  AND
C  P-1 CALLS TO CRSLVE ONLY. SINCE THE ARRAYS TOPBLK,ARRAY,
C  BOTBLK AND PIVOT CONTAIN THE DECOMPOSITION OF THE GIVEN
C  COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM
C  C O L R O W , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE
C  CALLS TO CRSLVE WITH THE SAME LEFT HAND SIDES. FOR THE
C  SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT
C  MATRIX, THE ARRAYS TOPBLK,ARRAY,BOTBLK MUST BE COPIED
C  BEFORE A CALL TO  C O L R O W .
C
C***********************************************************************
C
C             *****  SAMPLE CALLING PROGRAM  *****
C
C      THE FOLLOWING PROGRAM WILL EXERCISE COLROW, IN THE
C      CASE WHEN THE COEFFICIENT MATRIX IS NON-SINGULAR.
C
C       DOUBLE PRECISION TOP,AR,BOT,B,X
C       DOUBLE PRECISION ERROR,ERR
C       DIMENSION TOP(2,4),AR(4,8,2),BOT(2,4),B(12),X(12)
C       INTEGER PIVOT(12)
C       DATA N,NRWTOP,NOVRLP,NRWBLK,NCLBLK,NBLOKS,NRWBOT/12,2,4,4,8,2,2/
C       DATA TOP(1,1),TOP(1,2),TOP(1,3),TOP(1,4),
C    *       TOP(2,1),TOP(2,2),TOP(2,3),TOP(2,4)/
C    *0.0 D0,-0.98D0,-0.79D0,-0.15D0,
C    *-1.00D0, 0.25D0,-0.87D0, 0.35D0/
C       DATA AR(1,1,1),AR(1,2,1),AR(1,3,1),AR(1,4,1),
C    *       AR(1,5,1),AR(1,6,1),AR(1,7,1),AR(1,8,1)/
C    *0.78D0, 0.31D0,-0.85D0, 0.89D0,-0.69D0,-0.98D0,-0.76D0,-0.82D0/
C       DATA AR(2,1,1),AR(2,2,1),AR(2,3,1),AR(2,4,1),
C    *       AR(2,5,1),AR(2,6,1),AR(2,7,1),AR(2,8,1)/
C    *0.12D0,-0.01D0, 0.75D0, 0.32D0,-1.00D0,-0.53D0,-0.83D0,-0.98D0/
C       DATA AR(3,1,1),AR(3,2,1),AR(3,3,1),AR(3,4,1),
C    *       AR(3,5,1),AR(3,6,1),AR(3,7,1),AR(3,8,1)/
C    *-0.58D0, 0.04D0, 0.87D0, 0.38D0,-1.00D0,-0.21D0,-0.93D0,-0.84D0/
C       DATA AR(4,1,1),AR(4,2,1),AR(4,3,1),AR(4,4,1),
C    *       AR(4,5,1),AR(4,6,1),AR(4,7,1),AR(4,8,1)/
C    *-0.21D0,-0.91D0,-0.09D0,-0.62D0,-1.99D0,-1.12D0,-1.21D0, 0.07D0/
C       DATA AR(1,1,2),AR(1,2,2),AR(1,3,2),AR(1,4,2),
C    *       AR(1,5,2),AR(1,6,2),AR(1,7,2),AR(1,8,2)/
C    *0.78D0,-0.93D0,-0.76D0, 0.48D0,-0.87D0,-0.14D0,-1.00D0,-0.59D0/
C       DATA AR(2,1,2),AR(2,2,2),AR(2,3,2),AR(2,4,2),
C    *       AR(2,5,2),AR(2,6,2),AR(2,7,2),AR(2,8,2)/
C    *-0.99D0, 0.21D0,-0.73D0,-0.48D0,-0.93D0,-0.91D0, 0.10D0,-0.89D0/
C       DATA AR(3,1,2),AR(3,2,2),AR(3,3,2),AR(3,4,2),
C    *       AR(3,5,2),AR(3,6,2),AR(3,7,2),AR(3,8,2)/
C    *-0.68D0,-0.09D0,-0.58D0,-0.21D0, 0.85D0,-0.39D0, 0.79D0,-0.71D0/
C       DATA AR(4,1,2),AR(4,2,2),AR(4,3,2),AR(4,4,2),
C    *       AR(4,5,2),AR(4,6,2),AR(4,7,2),AR(4,8,2)/
C    *0.39D0,-0.99D0,-0.12D0,-0.75D0,-0.68D0,-0.99D0, 0.50D0,-0.88D0/
C       DATA BOT(1,1),BOT(1,2),BOT(1,3),BOT(1,4),
C    *       BOT(2,1),BOT(2,2),BOT(2,3),BOT(2,4)/
C    *0.71D0,-0.64D0, 0.0 D0, 0.48D0,
C    *0.08D0,100.0D0,50.00D0,15.00D0/
C       DATA B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9),B(10),B(11),
C    *       B(12)/
C    *-1.92D0,-1.27D0,-2.12D0,-2.16D0,-2.27D0,-6.08D0,-3.03D0,
C    *-4.62D0,-1.02D0,-3.52D0,.55D0,165.08D0/
C
C***********************************************************************
C
C   THE INPUT MATRIX IS GIVEN BY:
C
C  0.0  -0.98 -0.79 -0.15
C -1.00  0.25 -0.87  0.35
C  0.78  0.31 -0.85  0.89 -0.69 -0.98 -0.76 -0.82
C  0.12 -0.01  0.75  0.32 -1.00 -0.53 -0.83 -0.98
C -0.58  0.04  0.87  0.38 -1.00 -0.21 -0.93 -0.84
C -0.21 -0.91 -0.09 -0.62 -1.99 -1.12 -1.21  0.07
C                          0.78 -0.93 -0.76  0.48 -0.87 -0.14 -1.00 -0.5
C                         -0.99  0.21 -0.73 -0.48 -0.93 -0.91  0.10 -0.8
C                         -0.68 -0.09 -0.58 -0.21  0.85 -0.39  0.79 -0.7
C                          0.39 -0.99 -0.12 -0.75 -0.68 -0.99  0.50 -0.8
C                                                  0.71 -0.64  0.0   0.4
C                                                  0.08 100.0 50.00 15.0
C
C       THE RIGHT HAND SIDE IS GIVEN BY:
C
C         B = (-1.92,-1.27,-2.12,-2.16,-2.27,-6.08,-3.03,-4.62,
C              -1.02,-3.52,0.55,165.08)
C
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY;
C
C          X = (1,1,1,1,1,1,1,1,1,1,1,1)
C
C***********************************************************************
C
C       CALL COLROW(N,TOP,NRWTOP,NOVRLP,AR,NRWBLK,NCLBLK,NBLOKS,
C    *              BOT,NRWBOT,PIVOT,B,X,IFLAG)
C       IF(IFLAG.NE.0)GO TO 1000
C       ERROR = 0.D0
C       DO 10 I=1,N
C          ERR = 1.D0 - X(I)
C          ERROR = DMAX1(ERROR,DABS(ERR))
C          WRITE(6,100)X(I),ERR
C  10   CONTINUE
C       WRITE(6,200)ERROR
C 200   FORMAT(12H MAX ERROR = ,D15.7)
C 100   FORMAT(1H ,F15.7,D15.7)
C       RETURN
C1000   CONTINUE
C       WRITE(6,300)IFLAG
C 300   FORMAT(9H IFLAG =  ,I3)
C       RETURN
C       END
C
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, ARRAY, BOTBLK, B, X
      INTEGER PIVOT(1)
      DIMENSION TOPBLK(NRWTOP,1), ARRAY(NRWBLK,NCLBLK,1),
     * BOTBLK(NRWBOT,1), B(1), X(1)
      CALL CRDCMP(N, TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK, NCLBLK,
     * NBLOKS, BOTBLK, NRWBOT, PIVOT, IFLAG)
      IF (IFLAG.NE.0) RETURN
      CALL CRSLVE(TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK, NCLBLK,
     * NBLOKS, BOTBLK, NRWBOT, PIVOT, B, X)
      RETURN
      END
      SUBROUTINE CRDCMP(N, TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK,       CRD   10
     * NCLBLK, NBLOKS, BOTBLK, NRWBOT, PIVOT, IFLAG)
C
C***************************************************************
C
C  C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A
C  USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH
C  PARTIAL PIVOTING.  THE MATRIX  A  IS STORED IN THE ARRAYS
C  TOPBLK, ARRAY, AND BOTBLK.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A TO BE DECOMPOSED
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, ARRAY, BOTBLK
      DOUBLE PRECISION ROWMAX, ROWPIV, ROWMLT, COLMAX, COLPIV
      DOUBLE PRECISION SWAP, COLMLT, PIVMAX, ZERO, TEMPIV
      INTEGER PIVOT(1)
      DIMENSION TOPBLK(NRWTOP,1), ARRAY(NRWBLK,NCLBLK,1),
     * BOTBLK(NRWBOT,1)
      DATA ZERO /0.0D0/
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
      IFLAG = 0
      PIVMAX = ZERO
      NRWTP1 = NRWTOP + 1
      NROWEL = NRWBLK - NRWTOP
      NRWEL1 = NROWEL + 1
      NVRLP0 = NOVRLP - 1
C
C***************************************************************
C
C          ****  CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C               IF PARAMETERS ARE INVALID THEN TERMINATE AT 10;
C                                         ELSE CONTINUE AT 100.
C
C***************************************************************
C
      IF (N.NE.NBLOKS*NRWBLK+NOVRLP) GO TO 10
      IF (NOVRLP.NE.NRWTOP+NRWBOT) GO TO 10
      IF (NCLBLK.NE.NOVRLP+NRWBLK) GO TO 10
      IF (NOVRLP.GT.NRWBLK) GO TO 10
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      GO TO 20
   10 CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE INVALID.  SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = 1
      RETURN
   20 CONTINUE
C
C***************************************************************
C
C               ****  FIRST, IN TOPBLK....
C
C***************************************************************
C
C          ***  APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN
C                 PIVOTING ....
C
C***************************************************************
C
      DO 110 I=1,NRWTOP
        IPLUS1 = I + 1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IPVT = I
        COLMAX = DABS(TOPBLK(I,I))
        DO 30 J=IPLUS1,NOVRLP
          TEMPIV = DABS(TOPBLK(I,J))
          IF (TEMPIV.LE.COLMAX) GO TO 30
          IPVT = J
          COLMAX = TEMPIV
   30   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410
        PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        PIVOT(I) = IPVT
        IF (IPVT.EQ.I) GO TO 60
        DO 40 L=I,NRWTOP
          SWAP = TOPBLK(L,IPVT)
          TOPBLK(L,IPVT) = TOPBLK(L,I)
          TOPBLK(L,I) = SWAP
   40   CONTINUE
        DO 50 L=1,NRWBLK
          SWAP = ARRAY(L,IPVT,1)
          ARRAY(L,IPVT,1) = ARRAY(L,I,1)
          ARRAY(L,I,1) = SWAP
   50   CONTINUE
   60   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        COLPIV = TOPBLK(I,I)
        DO 100 J=IPLUS1,NOVRLP
          COLMLT = TOPBLK(I,J)/COLPIV
          TOPBLK(I,J) = COLMLT
          IF (IPLUS1.GT.NRWTOP) GO TO 80
          DO 70 L=IPLUS1,NRWTOP
            TOPBLK(L,J) = TOPBLK(L,J) - COLMLT*TOPBLK(L,I)
   70     CONTINUE
   80     CONTINUE
          DO 90 L=1,NRWBLK
            ARRAY(L,J,1) = ARRAY(L,J,1) - COLMLT*ARRAY(L,I,1)
   90     CONTINUE
  100   CONTINUE
  110 CONTINUE
C
C***************************************************************
C
C          ****  IN EACH BLOCK ARRAY(,,K)....
C
C***************************************************************
C
      INCR = 0
      DO 320 K=1,NBLOKS
        KPLUS1 = K + 1
C
C          *****************************************************
C
C          ***  FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH
C                       ROW PIVOTING....
C
C          *****************************************************
C
        DO 180 J=NRWTP1,NRWBLK
          JPLUS1 = J + 1
          JMINN = J - NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IPVT = JMINN
          ROWMAX = DABS(ARRAY(JMINN,J,K))
          LOOP = JMINN + 1
          DO 120 I=LOOP,NRWBLK
            TEMPIV = DABS(ARRAY(I,J,K))
            IF (TEMPIV.LE.ROWMAX) GO TO 120
            IPVT = I
            ROWMAX = TEMPIV
  120     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410
          PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          INCRJ = INCR + J
          PIVOT(INCRJ) = INCR + IPVT + NRWTOP
          IF (IPVT.EQ.JMINN) GO TO 140
          DO 130 L=J,NCLBLK
            SWAP = ARRAY(IPVT,L,K)
            ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K)
            ARRAY(JMINN,L,K) = SWAP
  130     CONTINUE
  140     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          ROWPIV = ARRAY(JMINN,J,K)
          DO 150 I=LOOP,NRWBLK
            ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV
  150     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          DO 170 L=JPLUS1,NCLBLK
            ROWMLT = ARRAY(JMINN,L,K)
            DO 160 I=LOOP,NRWBLK
              ARRAY(I,L,K) = ARRAY(I,L,K) - ROWMLT*ARRAY(I,J,K)
  160       CONTINUE
  170     CONTINUE
  180   CONTINUE
C
C          *****************************************************
C
C          ***  NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH
C                      COLUMN PIVOTING....
C
C          *****************************************************
C
        DO 310 I=NRWEL1,NRWBLK
          IPLUSN = I + NRWTOP
          IPLUS1 = I + 1
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IPVT = IPLUSN
          COLMAX = DABS(ARRAY(I,IPVT,K))
          LOOP = IPLUSN + 1
          DO 190 J=LOOP,NCLBLK
            TEMPIV = DABS(ARRAY(I,J,K))
            IF (TEMPIV.LE.COLMAX) GO TO 190
            IPVT = J
            COLMAX = TEMPIV
  190     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410
          PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          INCRN = INCR + IPLUSN
          PIVOT(INCRN) = INCR + IPVT
          IRWBLK = IPLUSN - NRWBLK
          IF (IPVT.EQ.IPLUSN) GO TO 240
          DO 200 L=I,NRWBLK
            SWAP = ARRAY(L,IPVT,K)
            ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K)
            ARRAY(L,IPLUSN,K) = SWAP
  200     CONTINUE
          IPVBLK = IPVT - NRWBLK
          IF (K.EQ.NBLOKS) GO TO 220
          DO 210 L=1,NRWBLK
            SWAP = ARRAY(L,IPVBLK,KPLUS1)
            ARRAY(L,IPVBLK,KPLUS1) = ARRAY(L,IRWBLK,KPLUS1)
            ARRAY(L,IRWBLK,KPLUS1) = SWAP
  210     CONTINUE
          GO TO 240
  220     CONTINUE
          DO 230 L=1,NRWBOT
            SWAP = BOTBLK(L,IPVBLK)
            BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK)
            BOTBLK(L,IRWBLK) = SWAP
  230     CONTINUE
  240     CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
          COLPIV = ARRAY(I,IPLUSN,K)
          DO 300 J=LOOP,NCLBLK
            COLMLT = ARRAY(I,J,K)/COLPIV
            ARRAY(I,J,K) = COLMLT
            IF (I.EQ.NRWBLK) GO TO 260
            DO 250 L=IPLUS1,NRWBLK
              ARRAY(L,J,K) = ARRAY(L,J,K) - COLMLT*ARRAY(L,IPLUSN,K)
  250       CONTINUE
  260       CONTINUE
            JRWBLK = J - NRWBLK
            IF (K.EQ.NBLOKS) GO TO 280
            DO 270 L=1,NRWBLK
              ARRAY(L,JRWBLK,KPLUS1) = ARRAY(L,JRWBLK,KPLUS1) -
     *         COLMLT*ARRAY(L,IRWBLK,KPLUS1)
  270       CONTINUE
            GO TO 300
  280       CONTINUE
            DO 290 L=1,NRWBOT
              BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK) -
     *         COLMLT*BOTBLK(L,IRWBLK)
  290       CONTINUE
  300     CONTINUE
  310   CONTINUE
        INCR = INCR + NRWBLK
  320 CONTINUE
C
C***************************************************************
C
C          ****  FINALLY, IN BOTBLK....
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          ***  APPLY NRWBOT ROW ELIMINATIONS WITH ROW
C                  PIVOTING....
C
C               IF BOT HAS JUST ONE ROW GO TO 500
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWBOT.EQ.1) GO TO 400
      DO 390 J=NRWTP1,NVRLP0
        JPLUS1 = J + 1
        JMINN = J - NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IPVT = JMINN
        ROWMAX = DABS(BOTBLK(JMINN,J))
        LOOP = JMINN + 1
        DO 330 I=LOOP,NRWBOT
          TEMPIV = DABS(BOTBLK(I,J))
          IF (TEMPIV.LE.ROWMAX) GO TO 330
          IPVT = I
          ROWMAX = TEMPIV
  330   CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410
        PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCRJ = INCR + J
        PIVOT(INCRJ) = INCR + IPVT + NRWTOP
        IF (IPVT.EQ.JMINN) GO TO 350
        DO 340 L=J,NOVRLP
          SWAP = BOTBLK(IPVT,L)
          BOTBLK(IPVT,L) = BOTBLK(JMINN,L)
          BOTBLK(JMINN,L) = SWAP
  340   CONTINUE
  350   CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        ROWPIV = BOTBLK(JMINN,J)
        DO 360 I=LOOP,NRWBOT
          BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV
  360   CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 380 L=JPLUS1,NOVRLP
          ROWMLT = BOTBLK(JMINN,L)
          DO 370 I=LOOP,NRWBOT
            BOTBLK(I,L) = BOTBLK(I,L) - ROWMLT*BOTBLK(I,J)
  370     CONTINUE
  380   CONTINUE
  390 CONTINUE
  400 CONTINUE
C
C***************************************************************
C
C          DONE PROVIDED THE LAST ELEMENT IS NOT ZERO
C
C***************************************************************
C
      IF (PIVMAX+DABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN
C
C***************************************************************
C
C       ****  MATRIX IS SINGULAR - SET IFLAG = - 1.
C                                  TERMINATE AT 1000.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
  410 CONTINUE
      IFLAG = -1
      RETURN
      END
      SUBROUTINE CRSLVE(TOPBLK, NRWTOP, NOVRLP, ARRAY, NRWBLK, NCLBLK,  CRS   10
     * NBLOKS, BOTBLK, NRWBOT, PIVOT, B, X)
C
C***************************************************************
C
C  C R S L V E  SOLVES THE LINEAR SYSTEM
C                       A*X = B
C  USING THE DECOMPOSITION ALREADY GENERATED IN  C R D C M P.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY  ...
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         THE PIVOT VECTOR FROM  C R D C M P
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR
C
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, ARRAY, BOTBLK, X, B
      DOUBLE PRECISION DOTPRD, XJ, XINCRJ, BINCRJ, SWAP
      INTEGER PIVOT(1)
      DIMENSION TOPBLK(NRWTOP,1), ARRAY(NRWBLK,NCLBLK,1),
     * BOTBLK(NRWBOT,1), B(1), X(1)
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
      NRWTP1 = NRWTOP + 1
      NRWBK1 = NRWBLK + 1
      NVRLP1 = NOVRLP + 1
      NRWTP0 = NRWTOP - 1
      NRWBT1 = NRWBOT + 1
      NROWEL = NRWBLK - NRWTOP
      NRWEL1 = NROWEL + 1
      NVRLP0 = NOVRLP - 1
      NBLKS1 = NBLOKS + 1
      NBKTOP = NRWBLK + NRWTOP
C
C***************************************************************
C
C               ****  FORWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST, IN TOPBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 30 J=1,NRWTOP
        X(J) = B(J)/TOPBLK(J,J)
        IF (J.EQ.NRWTOP) GO TO 20
        XJ = -X(J)
        LOOP = J + 1
        DO 10 I=LOOP,NRWTOP
          B(I) = B(I) + TOPBLK(I,J)*XJ
   10   CONTINUE
   20   CONTINUE
   30 CONTINUE
C
C       ********************************************************
C
C          ***  IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INCR = 0
      DO 120 K=1,NBLOKS
        INCRTP = INCR + NRWTOP
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 50 J=1,NRWTOP
          INCRJ = INCR + J
          XINCRJ = -X(INCRJ)
          DO 40 I=1,NRWBLK
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
   40     CONTINUE
   50   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 80 J=NRWTP1,NRWBLK
          INCRJ = INCR + J
          JPIVOT = PIVOT(INCRJ)
          IF (JPIVOT.EQ.INCRJ) GO TO 60
          SWAP = B(INCRJ)
          B(INCRJ) = B(JPIVOT)
          B(JPIVOT) = SWAP
   60     CONTINUE
          BINCRJ = -B(INCRJ)
          LOOP = J - NRWTP0
          DO 70 I=LOOP,NRWBLK
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*BINCRJ
   70     CONTINUE
   80   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 110 J=NRWBK1,NBKTOP
          INCRJ = INCR + J
          JRWTOP = J - NRWTOP
          X(INCRJ) = B(INCRJ)/ARRAY(JRWTOP,J,K)
          IF (J.EQ.NBKTOP) GO TO 100
          XINCRJ = -X(INCRJ)
          LOOP = J - NRWTP0
          DO 90 I=LOOP,NRWBLK
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
   90     CONTINUE
  100     CONTINUE
  110   CONTINUE
        INCR = INCR + NRWBLK
  120 CONTINUE
C
C       ********************************************************
C
C          ***  FINALLY, IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INCRTP = INCR + NRWTOP
      DO 140 J=1,NRWTOP
        INCRJ = INCR + J
        XINCRJ = -X(INCRJ)
        DO 130 I=1,NRWBOT
          INCRI = INCRTP + I
          B(INCRI) = B(INCRI) + BOTBLK(I,J)*XINCRJ
  130   CONTINUE
  140 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWBOT.EQ.1) GO TO 180
      DO 170 J=NRWTP1,NVRLP0
        INCRJ = INCR + J
        JPIVOT = PIVOT(INCRJ)
        IF (JPIVOT.EQ.INCRJ) GO TO 150
        SWAP = B(INCRJ)
        B(INCRJ) = B(JPIVOT)
        B(JPIVOT) = SWAP
  150   CONTINUE
        BINCRJ = -B(INCRJ)
        LOOP = J - NRWTP0
        DO 160 I=LOOP,NRWBOT
          INCRI = INCRTP + I
          B(INCRI) = B(INCRI) + BOTBLK(I,J)*BINCRJ
  160   CONTINUE
  170 CONTINUE
  180 CONTINUE
C
C***************************************************************
C
C               ****  BACKWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 210 LL=1,NRWBOT
        J = NVRLP1 - LL
        INCRJ = INCR + J
        NRWBTL = NRWBT1 - LL
        X(INCRJ) = B(INCRJ)/BOTBLK(NRWBTL,J)
        IF (LL.EQ.NRWBOT) GO TO 200
        XINCRJ = -X(INCRJ)
        LOOP = NRWBOT - LL
        DO 190 I=1,LOOP
          INCRI = INCRTP + I
          B(INCRI) = B(INCRI) + BOTBLK(I,J)*XINCRJ
  190   CONTINUE
  200   CONTINUE
  210 CONTINUE
C
C       ********************************************************
C
C          ***  THEN IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 300 L=1,NBLOKS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        K = NBLKS1 - L
        INCR = INCR - NRWBLK
        DO 240 L1=NRWEL1,NRWBLK
          I = NRWBLK + NRWEL1 - L1
          IPLUSN = I + NRWTOP
          LOOP = IPLUSN + 1
          INCRN = INCR + IPLUSN
          DOTPRD = X(INCRN)
          DO 220 J=LOOP,NCLBLK
            INCRJ = INCR + J
            DOTPRD = DOTPRD - ARRAY(I,J,K)*X(INCRJ)
  220     CONTINUE
          X(INCRN) = DOTPRD
          IPVTN = PIVOT(INCRN)
          IF (INCRN.EQ.IPVTN) GO TO 230
          SWAP = X(INCRN)
          X(INCRN) = X(IPVTN)
          X(IPVTN) = SWAP
  230     CONTINUE
  240   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCRTP = INCR + NRWTOP
        DO 260 J=NRWBK1,NCLBLK
          INCRJ = INCR + J
          XINCRJ = -X(INCRJ)
          DO 250 I=1,NROWEL
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
  250     CONTINUE
  260   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 290 LL=1,NROWEL
          J = NRWBK1 - LL
          INCRJ = INCR + J
          NRWELL = NRWEL1 - LL
          X(INCRJ) = B(INCRJ)/ARRAY(NRWELL,J,K)
          IF (LL.EQ.NROWEL) GO TO 280
          XINCRJ = -X(INCRJ)
          LOOP = NROWEL - LL
          DO 270 I=1,LOOP
            INCRI = INCRTP + I
            B(INCRI) = B(INCRI) + ARRAY(I,J,K)*XINCRJ
  270     CONTINUE
  280     CONTINUE
  290   CONTINUE
  300 CONTINUE
C
C       ********************************************************
C
C          ***  IN TOPBLK FINISH WITH....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 330 L=1,NRWTOP
        I = NRWTP1 - L
        LOOP = I + 1
        DOTPRD = X(I)
        DO 310 J=LOOP,NOVRLP
          DOTPRD = DOTPRD - TOPBLK(I,J)*X(J)
  310   CONTINUE
        X(I) = DOTPRD
        IPVTI = PIVOT(I)
        IF (I.EQ.IPVTI) GO TO 320
        SWAP = X(I)
        X(I) = X(IPVTI)
        X(IPVTI) = SWAP
  320   CONTINUE
  330 CONTINUE
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN   10
C                                                                       MAN   20
C       THIS DRIVER RUNS A TEST PROBLEM USING THE ROUTINE               MAN   30
C       A R C E C O. THE ARRAYS MTR (FOR THE STRUCTURE OF THE           MAN   40
C       LINEAR SYSTEM), AR (THE COEFFICIENT MATRIX), AND                MAN   50
C       B (THE RIGHT HAND SIDE) ARE SET UP IN DATA STATE-               MAN   60
C       MENTS.                                                          MAN   70
C                                                                       MAN   80
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN   90
      DOUBLE PRECISION AR, X, B, ERROR, ERR                             MAN  100
      DIMENSION AR(114), MTR(3,5), B(18), X(18)                         MAN  110
      INTEGER PIVOT(18)                                                 MAN  120
      DATA N /18/, NMBLKS /5/                                           MAN  130
      DATA MTR(1,1), MTR(2,1), MTR(3,1), MTR(1,2), MTR(2,2), MTR(3,2),  MAN  140
     * MTR(1,3), MTR(2,3), MTR(3,3), MTR(1,4), MTR(2,4), MTR(3,4),      MAN  150
     * MTR(1,5), MTR(2,5), MTR(3,5) /2,4,3,4,7,4,5,8,2,3,6,3,4,5,0/     MAN  160
      DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7), AR(8),      MAN  170
     * AR(9), AR(10), AR(11), AR(12), AR(13), AR(14), AR(15), AR(16),   MAN  180
     * AR(17), AR(18), AR(19), AR(20), AR(21), AR(22), AR(23), AR(24)   MAN  190
     * /-1.00D0,-1.00D0,-0.98D0,0.25D0,-0.79D0,-0.87D0,-0.15D0,0.35D0,  MAN  200
     * 0.78D0,-0.82D0,-0.83D0,-0.21D0,0.31D0,0.12D0,-0.98D0,-0.93D0,    MAN  210
     * -0.85D0,-0.01D0,-0.58D0,-0.84D0,0.89D0,0.75D0,0.04D0,0.37D0/     MAN  220
      DATA AR(25), AR(26), AR(27), AR(28), AR(29), AR(30), AR(31),      MAN  230
     * AR(32), AR(33), AR(34), AR(35), AR(36), AR(37), AR(38), AR(39),  MAN  240
     * AR(40), AR(41), AR(42), AR(43), AR(44), AR(45), AR(46), AR(47),  MAN  250
     * AR(48) /-0.69D0,0.32D0,0.87D0,-0.94D0,-0.98D0,-1.00D0,0.38D0,    MAN  260
     * -0.96D0,-0.76D0,-0.53D0,-1.00D0,-1.00D0,-0.99D0,-0.87D0,-0.93D0, MAN  270
     * 0.85D0,0.17D0,-0.91D0,-0.14D0,-0.91D0,-0.39D0,-1.37D0,-0.28D0,   MAN  280
     * -1.00D0/                                                         MAN  290
      DATA AR(49), AR(50), AR(51), AR(52), AR(53), AR(54), AR(55),      MAN  300
     * AR(56), AR(57), AR(58), AR(59), AR(60), AR(61), AR(62), AR(63),  MAN  310
     * AR(64), AR(65), AR(66), AR(67), AR(68), AR(69), AR(70), AR(71),  MAN  320
     * AR(72) /0.10D0,0.79D0,1.29D0,0.90D0,-0.59D0,-0.89D0,-0.71D0,     MAN  330
     * -1.59D0,0.78D0,-0.99D0,-0.68D0,0.39D0,1.10D0,-0.93D0,0.21D0,     MAN  340
     * -0.09D0,-0.99D0,-1.63D0,-0.76D0,-0.73D0,-0.58D0,-0.12D0,-1.01D0, MAN  350
     * 0.48D0/                                                          MAN  360
      DATA AR(73), AR(74), AR(75), AR(76), AR(77), AR(78), AR(79),      MAN  370
     * AR(80), AR(81), AR(82), AR(83), AR(84), AR(85), AR(86), AR(87),  MAN  380
     * AR(88), AR(89), AR(90), AR(91), AR(92), AR(93), AR(94), AR(95),  MAN  390
     * AR(96) /-0.48D0,-0.21D0,-0.75D0,-0.27D0,0.08D0,-0.67D0,-0.24D0,  MAN  400
     * 0.61D0,0.56D0,-0.41D0,0.54D0,-0.99D0,0.40D0,-0.41D0,0.16D0,      MAN  410
     * -0.93D0,0.16D0,-0.16D0,0.70D0,-0.46D0,0.98D0,0.43D0,0.71D0,      MAN  420
     * -0.47D0/                                                         MAN  430
      DATA AR(97), AR(98), AR(99), AR(100), AR(101), AR(102), AR(103),  MAN  440
     * AR(104), AR(105), AR(106), AR(107), AR(108), AR(109), AR(110),   MAN  450
     * AR(111), AR(112), AR(113), AR(114) /-0.25D0,0.89D0,-0.97D0,      MAN  460
     * -0.98D0,-0.92D0,-0.94D0,-0.60D0,-0.73D0,-0.52D0,-0.54D0,-0.30D0, MAN  470
     * 0.07D0,-0.46D0,-1.00D0,0.18D0,0.04D0,-0.58D0,-0.36D0/            MAN  480
      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), MAN  490
     * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18)           MAN  500
     * /-2.92D0,-1.27D0,-1.30D0,-1.17D0,-2.10D0,-4.51D0,-1.71D0,-4.59D0,MAN  510
     * -4.19D0,-0.93D0,-3.31D0,0.52D0,-0.12D0,-0.05D0,-0.98D0,-2.07D0,  MAN  520
     * -2.73D0,-1.95D0/                                                 MAN  530
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN  540
C                                                                       MAN  550
C       THE INPUT MATRIX HAS 5 BLOCKS, AND HAS THE FOLLOWING            MAN  560
C       STRUCTURE, WHERE THE INTEGER ENTRIES ARE THE INDICES            MAN  570
C       OF THE ELEMENTS IN THE ARRAY AR, GIVEN MODULO 100.              MAN  580
C                                                                       MAN  590
C                                                                       MAN  600
C                                                                       MAN  610
C       1  3  5  7                                                      MAN  620
C       2  4  6  8                                                      MAN  630
C          9 13 17 21 25 29 33                                          MAN  640
C         10 14 18 22 26 30 34                                          MAN  650
C         11 15 19 23 27 31 35                                          MAN  660
C         12 16 20 24 28 32 36                                          MAN  670
C                  37 42 47 52 57 62 67 72                              MAN  680
C                  38 43 48 53 58 63 68 73                              MAN  690
C                  39 44 49 54 59 64 69 74                              MAN  700
C                  40 45 50 55 60 65 70 75                              MAN  710
C                  41 46 51 56 61 66 71 76                              MAN  720
C                                    77 80 83 86 89 92                  MAN  730
C                                    78 81 84 87 90 93                  MAN  740
C                                    79 82 85 88 91 94                  MAN  750
C                                             95 99  3  7 11            MAN  760
C                                             96  0  4  8 12            MAN  770
C                                             97  1  5  9 13            MAN  780
C                                             98  2  6 10 14            MAN  790
C                                                                       MAN  800
C                                                                       MAN  810
C       THE RIGHT HAND SIDE IS GIVEN BY :                               MAN  820
C                                                                       MAN  830
C       B =                                                             MAN  840
C                                                                       MAN  850
C       ( -2.92,-1.27,-1.30,-1.17,-2.10,-4.51,-1.71,-4.59,-4.19,        MAN  860
C       -0.93,-3.31, 0.52,-0.12,-0.05,-0.98,-2.07,-2.73,-1.95)          MAN  870
C                                                                       MAN  880
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY :                       MAN  890
C                                                                       MAN  900
C       X = (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)                       MAN  910
C                                                                       MAN  920
C                                                                       MAN  930
C       THE STRUCTURE MATRIX, MTR(1:3,1:5), IS GIVEN BY :               MAN  940
C                                                                       MAN  950
C       MTR =                                                           MAN  960
C                                                                       MAN  970
C                 2   4   5   3   4                                     MAN  980
C                 4   7   8   6   5                                     MAN  990
C                 3   4   2   3   0                                     MAN 1000
C                                                                       MAN 1010
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAN 1020
      CALL ARCECO(N, AR, MTR, NMBLKS, PIVOT, B, X, IFLAG)               MAN 1030
      IF (IFLAG.NE.0) GO TO 20                                          MAN 1040
      ERROR = 0.D0                                                      MAN 1050
      DO 10 I=1,18                                                      MAN 1060
        ERR = 1.D0 - X(I)                                               MAN 1070
        ERROR = DMAX1(ERROR,DABS(ERR))                                  MAN 1080
        WRITE (6,99998) X(I), ERR                                       MAN 1090
   10 CONTINUE                                                          MAN 1100
      WRITE (6,99999) ERROR                                             MAN 1110
99999 FORMAT (12H MAX ERROR =, D15.7)                                   MAN 1120
99998 FORMAT (1H , F15.7, D15.7)                                        MAN 1130
      RETURN                                                            MAN 1140
   20 CONTINUE                                                          MAN 1150
      WRITE (6,99997) IFLAG                                             MAN 1160
99997 FORMAT (9H IFLAG = , I3)                                          MAN 1170
      RETURN                                                            MAN 1180
      END                                                               MAN 1190
      SUBROUTINE ARCECO(N, ARRAY, MTRSTR, NMBLKS, PIVOT, B, X, IFLAG)   ARC   10
C
C***************************************************************
C
C  THIS PROGRAM SOLVES THE LINEAR SYSTEM  A*X = B  WHERE  A IS
C  AN ALMOST BLOCK DIAGONAL MATRIX.  THE METHOD IMPLEMENTED IS
C  BASED ON GAUSS ELIMINATION WITH ALTERNATE ROW AND COLUMN
C  ELIMINATION WITH PARTIAL PIVOTING, WHICH PRODUCES A STABLE
C  DECOMPOSITION OF THE MATRIX  A  WITHOUT INTRODUCING FILL-IN.
C
C***************************************************************
C
C    TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE,
C    REMOVE ALL DOUBLE PRECISION STATEMENTS. THERE IS ONE
C    SUCH STATEMENT IN EACH OF :-  ARCECO,ARCEDC,ARCEPR,
C    ARCESL,ARCEFE,ARCEFS,ARCEFM,ARCEBM,ARCEBS,ARCEBE; AND
C    TWO IN ARCEPC.
C    IN ADDITION, REFERENCES TO THE BUILT-IN FUNCTIONS DABS
C    AND DMAX1 MUST BE REPLACED BY ABS, AND AMAX1, RESPECTIVELY,
C        DABS OCCURS TWICE IN ARCEPR (JUST BEFORE THE LOOP
C    DO 20, AND INSIDE THAT LOOP), AND TWICE IN ARCEPC (AGAIN
C    JUST BEFORE THE LOOP DO 20, AND INSIDE THAT LOOP).
C        DMAX1 OCCURS TWICE IN THE PACKAGE, ONCE IN ARCEPR AFTER
C    THE TEST FOR SINGULARTY) AND ONCE IN ARCEPC, ALSO AFTER
C    THE TEST FOR SINGULARITY.
C    FINALLY, ZERO IS INITIALISED TO 0.D0 IN A DATA STATEMENT
C    IN ARCEDC. THIS MUST BE REPLACED BY :
C        DATA ZERO/0.0/
C
C***************************************************************
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C                N     - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM, WHERE
C                           N = SUM(MTRSTR(1,K),K=1,NMBLKS)
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         WHERE
C                          NUMELS = SUM(MTRSTR(1,K)*MTRSTR(2,K);
C                               K=1,NMBLKS).
C                         CONTAINS THE ENTRIES OF THE ALMOST
C                         BLOCK DIAGONAL MATRIX  A  WHOSE BLOCK
C                         STRUCTURE IS GIVEN BY THE INTEGER ARRAY
C                         MTRSTR. THE ELEMENTS OF A ARE STORED BY
C                         COLUMNS, IN BLOCKS CORRESPONDING TO THE
C                         GIVEN STRUCTURE.
C
C               MTRSTR - INTEGER(3,NMBLKS)
C                         DESCRIBES THE BLOCK STRUCTURE OF  A:
C                           MTRSTR(1,K)  = NUMBER OF ROWS IN
C                                            BLOCK K.
C                           MTRSTR(2,K)  = NUMBER OF COLUMNS IN
C                                            BLOCK K.
C                           MTRSTR(3,K)  = NUMBER OF COLUMNS
C                                          OVERLAPPED BY BLOCK K
C                                          AND BLOCK (K+1).
C                        MTRSTR MUST SATISFY SOME RESTRICTIONS.
C                         IN ORDER THAT A BE SQUARE, WE NEED
C                         SUM(MTRSTR(1,K(,K=1,NMBLKS) = N =
C                         SUM((MTRSTR(2,K)-MTRSTR(3,K)),K=1,NMBLKS).
C                        IN ADDITION, TO ENSURE THAT THREE SUCCESS-
C                         IVE BLOCKS DO NOT HAVE COLUMNS IN COMMON,
C                         MTRSTR MUST SATISFY
C                         MTRSTR(3,K-1)+MTRSTR(3,K) <= MTRSTR(2,K),
C                         FOR K = 2,NMBLKS.
C                        FINALLY,  A R C E C O,  SETS
C                         MTRSTR(3,NMBLKS) = 0, IN ARCECD.
C
C               NMBLKS - INTEGER
C                         TOTAL NUMBER OF BLOCKS IN  A
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         CONTAINS THE MODIFIED ALTERNATE ROW
C                         AND COLUMN DECOMPOSITION OF  A (IF
C                         IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR (IF IFLAG = 0)
C
C                IFLAG - INTEGER
C                         =  1,IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
C               *****  AUXILIARY PROGRAMS  *****
C
C            ARCEDC(ARRAY,MTRSTR,NMBLKS,PIVOT,IFLAG)
C               - DECOMPOSES THE MATRIX  A  USING MODIFIED
C                 ALTERNATE ROW AND COLUMN ELIMINATION
C                 WITH PARTIAL PIVOTING, AND IS USED FOR
C                 THIS PURPOSE IN  A R C E C O.
C                 THE ARGUMENTS ARE ALL AS IN A R C E C O.
C
C            ARCESL(ARRAY,MTRSTR,NMBLKS,PIVOT,B,X)
C               - SOLVES THE SYSTEM  A*X = B  ONCE  A  IS
C                 DECOMPOSED.
C                 THE ARGUMENTS ARE ALL AS IN A R C E C O .
C
C***************************************************************
C
C               *****  BLOCK STRUCTURE OF  A  *****
C
C  THE NMBLKS BLOCKS OF A ARE STORED CONSECUTIVELY IN THE ONE
C  DIMENSIONAL MATRIX  ARRAY, THE ENTRIES OF  A  BEING STORED
C  AS FOLLOWS:
C
C       IN ARRAY(1)     THE (1,1) ENTRY OF THE TOP BLOCK,
C
C       IN ARRAY(INDEX) THE (1,1) ENTRY OF THE ITH BLOCK WHERE
C                       INDEX = 1 + SUM(MTRSTR(1,J)*MTRSTR(2,J);
C                               J=1,I-1), I=2,NMBLKS.
C
C***************************************************************
C
C       THE SUBROUTINE  A R C E C O  AUTOMATICALLY SOLVES THE
C  INPUT SYSTEM WHEN IFLAG=0.  A R C E C O  IS CALLED ONLY ONCE
C  FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT
C  HAND SIDES CAN BE OBTAINED BY ONE CALL TO  A R C E C O  AND
C  P-1 CALLS TO ARCESL ONLY. SINCE THE ARRAYS ARRAY AND
C  PIVOT CONTAIN, RESPECTIVELY, THE DECOMPOSITION OF THE GIVEN
C  COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM
C  A R C E C O , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE
C  CALLS TO ARCESL WITH THE SAME RIGHT HAND SIDES. FOR THE
C  SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT
C  MATRIX, THE ARRAY ARRAY MUST BE COPIED BEFORE A CALL
C  TO  A R C E C O .
C
C**********************************************************************
C
C               *****  SAMPLE CALLING PROGRAM  *****
C
C       THE FOLLOWING PROGRAM WILL EXERCISE  A R C E C O , IN
C       THE CASE WHEN THE COEFFICIENT MATRIX IS NONSINGULAR.
C
C       DOUBLE PRECISION AR,X,B,ERROR,ERR
C       DIMENSION AR(114),MTR(3,5),B(18),X(18)
C       INTEGER PIVOT(18)
C       DATA MTR(1,1),MTR(2,1),MTR(3,1),MTR(1,2),MTR(2,2),MTR(3,2),
C    *       MTR(1,3),MTR(2,3),MTR(3,3),MTR(1,4),MTR(2,4),MTR(3,4),
C    *       MTR(1,5),MTR(2,5),MTR(3,5)/
C    *  2,4,3,4,7,4,5,8,2,3,6,3,4,5,0/
C       DATA AR(1),AR(2),AR(3),AR(4),AR(5),AR(6),AR(7),AR(8),AR(9),
C    *       AR(10),AR(11),AR(12),AR(13),AR(14),AR(15),AR(16),AR(17),
C    *       AR(18),AR(19),AR(20),AR(21),AR(22),AR(23),AR(24)/
C    *-1.00D0,-1.00D0,-0.98D0, 0.25D0,-0.79D0,-0.87D0,-0.15D0, 0.35D0,
C    * 0.78D0,-0.82D0,-0.83D0,-0.21D0, 0.31D0, 0.12D0,-0.98D0,-0.93D0,
C    *-0.85D0,-0.01D0,-0.58D0,-0.84D0, 0.89D0, 0.75D0, 0.04D0, 0.37D0/
C       DATA AR(25),AR(26),AR(27),AR(28),AR(29),AR(30),AR(31),AR(32),
C    *       AR(33),AR(34),AR(35),AR(36),AR(37),AR(38),AR(39),AR(40),
C    *       AR(41),AR(42),AR(43),AR(44),AR(45),AR(46),AR(47),AR(48)/
C    *-0.69D0, 0.32D0, 0.87D0,-0.94D0,-0.98D0,-1.00D0, 0.38D0,-0.96D0,
C    *-0.76D0,-0.53D0,-1.00D0,-1.00D0,-0.99D0,-0.87D0,-0.93D0, 0.85D0,
C    * 0.17D0,-0.91D0,-0.14D0,-0.91D0,-0.39D0,-1.37D0,-0.28D0,-1.00D0/
C       DATA AR(49),AR(50),AR(51),AR(52),AR(53),AR(54),AR(55),AR(56),
C    *       AR(57),AR(58),AR(59),AR(60),AR(61),AR(62),AR(63),AR(64),
C    *       AR(65),AR(66),AR(67),AR(68),AR(69),AR(70),AR(71),AR(72)/
C    * 0.10D0, 0.79D0, 1.29D0, 0.90D0,-0.59D0,-0.89D0,-0.71D0,-1.59D0,
C    * 0.78D0,-0.99D0,-0.68D0, 0.39D0, 1.10D0,-0.93D0, 0.21D0,-0.09D0,
C    *-0.99D0,-1.63D0,-0.76D0,-0.73D0,-0.58D0,-0.12D0,-1.01D0, 0.48D0/
C       DATA AR(73),AR(74),AR(75),AR(76),AR(77),AR(78),AR(79),AR(80),
C    *       AR(81),AR(82),AR(83),AR(84),AR(85),AR(86),AR(87),AR(88),
C    *       AR(89),AR(90),AR(91),AR(92),AR(93),AR(94),AR(95),AR(96)/
C    *-0.48D0,-0.21D0,-0.75D0,-0.27D0, 0.08D0,-0.67D0,-0.24D0, 0.61D0,
C    * 0.56D0,-0.41D0, 0.54D0,-0.99D0, 0.40D0,-0.41D0, 0.16D0,-0.93D0,
C    * 0.16D0,-0.16D0, 0.70D0,-0.46D0, 0.98D0, 0.43D0, 0.71D0,-0.47D0/
C       DATA AR(97),AR(98),AR(99),AR(100),AR(101),AR(102),AR(103),
C    *       AR(104),AR(105),AR(106),AR(107),AR(108),AR(109),AR(110),
C    *       AR(111),AR(112),AR(113),AR(114)/
C    *-0.25D0, 0.89D0,-0.97D0,-0.98D0,-0.92D0,-0.94D0,-0.60D0,-0.73D0,
C    *-0.52D0,-0.54D0,-0.30D0, 0.07D0,-0.46D0,-1.00D0, 0.18D0, 0.04D0,
C    *-0.58D0,-0.36D0/
C       DATA B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9),B(10),B(11),
C    *       B(12),B(13),B(14),B(15),B(16),B(17),B(18)/
C    *-2.92D0,-1.27D0,-1.30D0,-1.17D0,-2.10D0,-4.51D0,-1.71D0,-4.59D0,
C    *-4.19D0,-0.93D0,-3.31D0, 0.52D0,-0.12D0,-0.05D0,-0.98D0,-2.07D0,
C    *-2.73D0,-1.95D0/
C       DATA N/18/,NMBLKS/5/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       THE MATRIX AR HAS 5 BLOCKS, WITH THE FOLLOWING STRUCTURE, THE
C       ENTRIES BEING THE INDICES OF THE ELEMENTS OF AR, MODULO 100.
C
C
C       1  3  5  7
C       2  4  6  8
C          9 13 17 21 25 29 33
C         10 14 18 22 26 30 34
C         11 15 19 23 27 31 35
C         12 16 20 24 28 32 36
C                  37 42 47 52 57 62 67 72
C                  38 43 48 53 58 63 68 73
C                  39 44 49 54 59 64 69 74
C                  40 45 50 55 60 65 70 75
C                  41 46 51 56 61 66 71 76
C                                    77 80 83 86 89 92
C                                    78 81 84 87 90 93
C                                    79 82 85 88 91 94
C                                             95 99  3  7 11
C                                             96  0  4  8 12
C                                             97  1  5  9 13
C                                             98  2  6 10 14
C
C
C       THE RIGHT HAND SIDE IS GIVEN BY :
C
C       B = ( -2.92,-1.27,-1.30,-1.17,-2.10,-4.51,-1.71,-4.59,-4.19,
C       -0.93,-3.31, 0.52,-0.12,-0.05,-0.98,-2.07,-2.73,-1.95)
C
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY :
C
C       X = (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
C
C       THE STRUCTURE MATRIX, MTR(1:3,1:5), IS GIVEN BY :
C
C                 2   4   5   3   4
C       MTR =     4   7   8   6   5
C                 3   4   2   3   0
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       CALL ARCECO(N,AR,MTR,NMBLKS,PIVOT,B,X,IFLAG)
C       IF(IFLAG.NE.0)GO TO 1000
C       ERROR = 0.D0
C       DO 10 I=1,18
C          ERR = 1.D0 - X(I)
C          ERROR = DMAX1(ERROR,DABS(ERR))
C          WRITE(5,100)X(I),ERR
C  10   CONTINUE
C       WRITE(5,200)ERROR
C 200   FORMAT(12H MAX ERROR = ,D15.7)
C 100   FORMAT(1H ,F15.7,D15.7)
C       RETURN
C1000   CONTINUE
C       WRITE(5,300)IFLAG
C 300   FORMAT(9H IFLAG =  ,I3)
C       RETURN
C       END
      DOUBLE PRECISION ARRAY, B, X
      INTEGER MTRSTR(3,1), PIVOT(1)
      DIMENSION ARRAY(1), B(1), X(1)
      CALL ARCEDC(N, ARRAY, MTRSTR, NMBLKS, PIVOT, IFLAG)
      IF (IFLAG.NE.0) RETURN
      CALL ARCESL(ARRAY, MTRSTR, NMBLKS, PIVOT, B, X)
      RETURN
      END
      SUBROUTINE ARCEDC(N, ARRAY, MTRSTR, NMBLKS, PIVOT, IFLAG)         ARC   10
C
C***************************************************************
C
C  A R C E D C SUPERVISES THE MODIFIED ALTERNATE ROW AND COLUMN
C  DECOMPOSITION WITH PARTIAL PIVOTING OF THE ALMOST BLOCK
C  DIAGONAL MATRIX  A   STORED IN THE ARRAYS A R R A Y  AND
C  M T R S T R .
C
C***************************************************************
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C                N     - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM, WHERE
C                           N = SUM(MTRSTR(1,K),K=1,NMBLKS)
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         WHERE
C                          NUMELS = SUM(MTRSTR(1,K)*MTRSTR(2,K);
C                               K=1,NMBLKS).
C                         CONTAINS THE ENTRIES OF THE ALMOST
C                         BLOCK DIAGONAL MATRIX  A  WHOSE BLOCK
C                         STRUCTURE IS GIVEN BY THE INTEGER ARRAY
C                         MTRSTR. THE ELEMENTS OF A ARE STORED BY
C                         COLUMNS, IN BLOCKS CORRESPONDING TO THE
C                         GIVEN STRUCTURE.
C               MTRSTR - INTEGER(3,NMBLKS)
C                         DESCRIBES THE BLOCK STRUCTURE OF  A:
C                           MTRSTR(1,K)  = NUMBER OF ROWS IN
C                                            BLOCK K.
C                           MTRSTR(2,K)  = NUMBER OF COLUMNS IN
C                                            BLOCK K.
C                           MTRSTR(3,K)  = NUMBER OF COLUMNS
C                                          OVERLAPPED BY BLOCK K
C                                          AND BLOCK (K+1).
C                        MTRSTR MUST SATISFY SOME RESTRICTIONS.
C                         IN ORDER THAT A BE SQUARE, WE NEED
C                         SUM(MTRSTR(1,K(,K=1,NMBLKS) = N =
C                         SUM((MTRSTR(2,K)-MTRSTR(3,K)),K=1,NMBLKS).
C                        IN ADDITION, TO ENSURE THAT THREE SUCCESS-
C                         IVE BLOCKS DO NOT HAVE COLUMNS IN COMMON,
C                         MTRSTR MUST SATISFY
C                         MTRSTR(3,K-1)+MTRSTR(3,K) <= MTRSTR(2,K),
C                         FOR K = 2,NMBLKS.
C                        FINALLY,  A R C E C O,  SETS
C                         MTRSTR(3,NMBLKS) = 0, IN ARCECD.
C
C               NMBLKS - INTEGER
C                         TOTAL NUMBER OF BLOCKS
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         CONTAINS THE MODIFIED ALTERNATE ROW
C                         AND COLUMN DECOMPOSITION OF  A (IF
C                         IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                IFLAG - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
C               *****  AUXILIARY PROGRAMS  *****
C
C       ARCEPR(BLOCK,NRWBLK,NCLBLK,NRWPIV,PIVOT,PIVMAX,IFLAG)
C        CARRIES OUT THE ROW ELIMINATIONS
C
C       ARCEPC(TOPBLK,NRWTOP,NOVRLP,BOTBLK,NRWBOT,NCLPIV,
C            PIVOT,PIVMAX,IFLAG)
C        CARRIES OUT THE COLUMN ELIMINATIONS
C
C***************************************************************
C
      DOUBLE PRECISION ARRAY, PIVMAX, ZERO
      INTEGER PIVOT(1)
      DIMENSION ARRAY(1), MTRSTR(3,1)
      DATA ZERO /0.D0/
C
C***************************************************************
C
C       ****  CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C            IF PARAMETERS ARE INVALID THEN TERMINATE AT 7;
C                                      ELSE CONTINUE AT 8.
C
C***************************************************************
C
C
      MTRSTR(3,NMBLKS) = 0
      DO 10 K=2,NMBLKS
        IF (MTRSTR(3,K-1)+MTRSTR(3,K).GT.MTRSTR(2,K)) GO TO 30
   10 CONTINUE
      ISUM1 = 0
      ISUM2 = 0
      DO 20 K=1,NMBLKS
        ISUM1 = ISUM1 + MTRSTR(1,K)
        ISUM2 = ISUM2 + MTRSTR(2,K) - MTRSTR(3,K)
   20 CONTINUE
      IF (ISUM1.NE.ISUM2) GO TO 30
      IF (ISUM1.NE.N) GO TO 30
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE ACCEPTABLE - CONTINUE AT 8
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      GO TO 40
   30 CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C        PARAMETERS ARE INVALID.  SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = 1
      RETURN
   40 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C  INTERNAL PARAMETERS:
C       C
C       C       INDEX1: POINTER TO THE ELEMENT IN THE COLUMN
C       C               WHERE ROW PIVOTING STARTS.
C       C
C       C       INDEX2: POINTER TO THE ELEMENT IN THE COLUMN
C       C               WHERE COLUMN PIVOTING STARTS.
C       C
C       C       INDEX3: POINTER TO 1ST ELEMENT IN 1ST COLUMN
C       C               OF NEXT BLOCK.
C       C
C       C       INDPIV: POINTER TO 1ST ELEMENT OF BLOCK OF PIVOT
C       C
C       C       NRWBLK: NUMBER OF ROWS IN BLOCK.
C       C
C       C       NRWBK2: NUMBER OF ROWS IN NEXT BLOCK.
C       C
C       C       NRWPIV: NUMBER OF ROW ELIMINATIONS.
C       C
C       C       NCLBLK: NUMBER OF COLUMNS IN BLOCK TO BE
C       C               ROW PIVOTED.
C       C
C       C       NCLPIV: NUMBER OF COLUMN ELIMINATIONS.
C       C
C       C       NOVRLP: NUMBER OF COLUMNS OVERLAPPED BY THE
C       C               CURRENT BLOCK AND THE NEXT BLOCK.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      PIVMAX = ZERO
      IFLAG = 0
      INDEX1 = 1
      INDPIV = 1
      NRWBLK = MTRSTR(1,1)
      NCLBLK = MTRSTR(2,1)
      NOVRLP = MTRSTR(3,1)
      NRWPIV = NCLBLK - NOVRLP
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       CALL ARCEPR TO PERFORM NRWPIV ROW ELIMINATIONS
C       C               ON TOP BLOCK.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWPIV.GT.0) CALL ARCEPR(ARRAY(INDEX1), NRWBLK, NCLBLK,
     * NRWPIV, PIVOT(INDPIV), PIVMAX, IFLAG)
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       IF MATRIX SINGULAR RETURN.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (IFLAG.LT.0) RETURN
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       NOW DO DECOMPOSITION PROCEEDING ONE BLOCK AT A
C       C               TIME.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 70 K=2,NMBLKS
        INDPIV = INDPIV + NRWPIV
        INDEX2 = INDEX1 + NRWBLK*NRWPIV
        INDEX3 = INDEX2 + NRWBLK*NOVRLP
        NCLPIV = NRWBLK - NRWPIV
        NRWBK2 = MTRSTR(1,K)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C  CALL ARCEPC TO PERFORM NCLPIV COLUMN ELIMINATIONS.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.EQ.0) GO TO 50
        CALL ARCEPC(ARRAY(INDEX2), NRWBLK, NOVRLP, ARRAY(INDEX3),
     *   NRWBK2, NCLPIV, PIVOT(INDPIV), PIVMAX, IFLAG)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C             C
C             C    IF MATRIX IS SINGULAR RETURN.
C             C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (IFLAG.LT.0) RETURN
   50   CONTINUE
        NRWBLK = NRWBK2
        INDEX1 = INDEX3 + NRWBLK*NCLPIV
        NCLBLK = MTRSTR(2,K) - NCLPIV
        NOVRLP = MTRSTR(3,K)
        NRWPIV = NCLBLK - NOVRLP
        INDPIV = INDPIV + NCLPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEPR TO PERFORM NRWPIV ROW ELIMINATIONS.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NRWPIV.EQ.0) GO TO 60
        CALL ARCEPR(ARRAY(INDEX1), NRWBLK, NCLBLK, NRWPIV,
     *   PIVOT(INDPIV), PIVMAX, IFLAG)
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C             C
C             C    IF MATRIX IS SINGULAR RETURN.
C             C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (IFLAG.LT.0) RETURN
   60   CONTINUE
   70 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEPR(BLOCK, NRWBLK, NCLBLK, NRWPIV, PIVOT, PIVMAX,   ARC   10
     * IFLAG)
C
C***************************************************************
C
C  A R C E P R  PERFORMS  NRWPIV  ROW ELIMINATIONS ON THE MATRIX
C               BLOCK
C
C***************************************************************
C
      INTEGER PIVOT(NRWBLK)
      DOUBLE PRECISION BLOCK, ROWMAX, PIVMAX, TEMPIV, ROWPIV, SWAP
      DIMENSION BLOCK(NRWBLK,NCLBLK)
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C        PERFORM NRWPIV ROW ELIMINATIONS...
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 90 J=1,NRWPIV
        JPLUS1 = J + 1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    DETERMINE ROW PIVOT AND PIVOT INDEX
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        MAX = J
        ROWMAX = DABS(BLOCK(J,J))
        IF (J.EQ.NRWBLK) GO TO 30
        DO 20 I1=JPLUS1,NRWBLK
          TEMPIV = DABS(BLOCK(I1,J))
          IF (TEMPIV.LE.ROWMAX) GO TO 10
          ROWMAX = TEMPIV
          MAX = I1
   10     CONTINUE
   20   CONTINUE
   30   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    TEST FOR SINGULARITY:
C          C
C          C             IF SINGULAR THEN TERMINATE AT 90;
C          C                         ELSE CONTINUE.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 100
        PIVMAX = DMAX1(PIVMAX,ROWMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    IF NECESSARY INTERCHANGE ROWS
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        PIVOT(J) = MAX
        IF (J.EQ.MAX) GO TO 50
        DO 40 J1=J,NCLBLK
          SWAP = BLOCK(MAX,J1)
          BLOCK(MAX,J1) = BLOCK(J,J1)
          BLOCK(J,J1) = SWAP
   40   CONTINUE
   50   CONTINUE
        IF (J.EQ.NRWBLK) RETURN
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    COMPUTE THE MULTIPLIERS
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        ROWPIV = BLOCK(J,J)
        DO 60 I1=JPLUS1,NRWBLK
          BLOCK(I1,J) = BLOCK(I1,J)/ROWPIV
   60   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    PERFORM ROW ELIMINATIONS WITH COLUMN INDEXING
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 80 J1=JPLUS1,NCLBLK
          DO 70 L1=JPLUS1,NRWBLK
            BLOCK(L1,J1) = BLOCK(L1,J1) - BLOCK(L1,J)*BLOCK(J,J1)
   70     CONTINUE
   80   CONTINUE
   90 CONTINUE
      RETURN
  100 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       MATRIX IS SINGULAR - SET IFLAG = -1.
C       C                            TERMINATE AT 90.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = -1
      RETURN
      END
      SUBROUTINE ARCEPC(TOPBLK, NRWTOP, NOVRLP, BOTBLK, NRWBOT, NCLPIV, ARC   10
     * PIVOT, PIVMAX, IFLAG)
C
C***************************************************************
C
C     A R C E P C PERFORMS NCLPIV COLUMN ELIMINATIONS ON THE
C                 MATRICES TOPBLK AND BOTBLK
C***************************************************************
C
      DOUBLE PRECISION TOPBLK, BOTBLK, COLMAX, PIVMAX, COLMLT
      DOUBLE PRECISION TEMPIV, SWAP
      INTEGER PIVOT(NRWTOP)
      DIMENSION TOPBLK(NRWTOP,NOVRLP), BOTBLK(NRWBOT,NOVRLP)
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       PERFORM THE COLUMN ELIMINATIONS ON A LOOP.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 110 J=1,NCLPIV
        I = NRWTOP - NCLPIV + J
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    DETERMINE COLUMN PIVOT AND PIVOT INDEX
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        MAX = J
        COLMAX = DABS(TOPBLK(I,J))
        IF (J.EQ.NOVRLP) GO TO 30
        JPLUS1 = J + 1
        DO 20 J1=JPLUS1,NOVRLP
          TEMPIV = DABS(TOPBLK(I,J1))
          IF (TEMPIV.LE.COLMAX) GO TO 10
          COLMAX = TEMPIV
          MAX = J1
   10     CONTINUE
   20   CONTINUE
   30   CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    TEST FOR SINGULARITY:
C          C
C          C             IF SINGULAR THEN TERMINATE AT 110;
C          C                         ELSE CONTINUE.
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 120
        PIVMAX = DMAX1(PIVMAX,COLMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    IF NECESSARY INTERCHANGE COLUMNS
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        PIVOT(J) = MAX
        IF (J.EQ.MAX) GO TO 60
        DO 40 I1=I,NRWTOP
          SWAP = TOPBLK(I1,J)
          TOPBLK(I1,J) = TOPBLK(I1,MAX)
          TOPBLK(I1,MAX) = SWAP
   40   CONTINUE
        DO 50 I2=1,NRWBOT
          SWAP = BOTBLK(I2,J)
          BOTBLK(I2,J) = BOTBLK(I2,MAX)
          BOTBLK(I2,MAX) = SWAP
   50   CONTINUE
   60   CONTINUE
        IF (J.EQ.NOVRLP) RETURN
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    COMPUTE MULTIPLIERS AND PERFORM COLUMN
C          C            ELIMINATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 100 J1=JPLUS1,NOVRLP
          COLMLT = TOPBLK(I,J1)/TOPBLK(I,J)
          TOPBLK(I,J1) = COLMLT
          IF (I.EQ.NRWTOP) GO TO 80
          IPLUS1 = I + 1
          DO 70 L1=IPLUS1,NRWTOP
            TOPBLK(L1,J1) = TOPBLK(L1,J1) - COLMLT*TOPBLK(L1,J)
   70     CONTINUE
   80     CONTINUE
          DO 90 L1=1,NRWBOT
            BOTBLK(L1,J1) = BOTBLK(L1,J1) - COLMLT*BOTBLK(L1,J)
   90     CONTINUE
  100   CONTINUE
  110 CONTINUE
      RETURN
  120 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       MATRIX IS SINGULAR - SET IFLAG = -1.
C       C                            TERMINATE AT 110.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IFLAG = -1
      RETURN
      END
      SUBROUTINE ARCESL(ARRAY, MTRSTR, NMBLKS, PIVOT, B, X)             ARC   10
C
C***************************************************************
C
C  A R C E S L  SUPERVISES THE SOLUTION OF THE LINEAR SYSTEM
C                          A*X = B
C  USING THE DECOMPOSITION OF THE MATRIX  A  ALREADY GENERATED
C  IN  A R C E D C.  IT INVOLVES TWO LOOPS, THE FORWARD LOOP,
C  CONSISTING OF FORWARD SOLUTION, FORWARD MODIFICATION, AND
C  FORWARD ELIMINATION, AND THE BACKWARD LOOP, CONSISTING OF
C  BACKWARD SOLUTION, BACKWARD MODIFICATION, AND BACKWARD
C  ELIMINATION.
C
C***************************************************************
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C                ARRAY - DOUBLE PRECISION(NUMELS)
C                         WHERE
C                          NUMELS = SUM(MTRSTR(1,K)*MTRSTR(2,K);
C                               K=1,NMBLKS).
C                         OUTPUT FROM A R C E D C
C
C               MTRSTR - INTEGER(3,NMBLKS)
C                         DESCRIBES THE BLOCK STRUCTURE OF  A:
C                           MTRSTR(1,K)  = NUMBER OF ROWS IN
C                                            BLOCK K.
C                           MTRSTR(2,K)  = NUMBER OF COLUMNS IN
C                                            BLOCK K.
C                           MTRSTR(3,K)  = NUMBER OF COLUMNS
C                                          OVERLAPPED BY BLOCK K
C                                          AND BLOCK (K+1).
C
C                        THE LINEAR SYSTEM IS OF ORDER:
C                           N  = SUM(MTRSTR(1,K),K=1,NMBLKS)
C
C               NMBLKS - INTEGER
C                         TOTAL NUMBER OF BLOCKS IN  A
C
C                PIVOT - INTEGER(N)
C                         OUTPUT FROM A R C E D C
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR
C
C***************************************************************
C
C               *****  AUXILIARY PROGRAMS  *****
C
C
C       ARCEFS - PERFORMS FORWARD SOLUTION STEP
C
C       ARCEFM - PERFORMS FORWARD MODIFICATION STEP
C
C       ARCEFE - PERFORMS FORWARD ELIMINATION STEP
C
C       ARCEBS - PERFORMS BACKWARD SOLUTION STEP
C
C       ARCEBM - PERFORMS BACKWARD MODIFICATION STEP
C
C       ARCEBE - PERFORMS BACKWARD ELIMINATION STEP
C
C***************************************************************
C
      DOUBLE PRECISION ARRAY, B, X
      INTEGER PIVOT(1)
      DIMENSION ARRAY(1), MTRSTR(3,1), B(1), X(1)
      INDPIV = 1
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C  INTERNAL PARAMETERS:
C       C
C       C       INDEXA: POINTER TO 1ST ELEMENT OF BLOCK OF  A  .
C       C
C       C       INDEXB: POINTER TO 1ST ELEMENT OF BLOCK OF B.
C       C
C       C       INDPIV,NRWBLK,NRWPIV,NCLBLK,NCLPIV,NOVRLP
C       C       ARE AS IN ARCEDC.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INDEXA = 1
      NRWBLK = MTRSTR(1,1)
      NCLBLK = MTRSTR(2,1)
      NOVRLP = MTRSTR(3,1)
      NRWPIV = NCLBLK - NOVRLP
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       CALL ARCEFE TO PERFORM FORWARD ELIMINATION.
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWPIV.GT.0) CALL ARCEFE(ARRAY(INDEXA), NRWBLK, NRWPIV,
     * PIVOT(INDPIV), B(INDPIV))
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       FORWARD LOOP
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 10 K=2,NMBLKS
        INDEXA = INDEXA + NRWBLK*NRWPIV
        NCLPIV = NRWBLK - NRWPIV
        INDPIV = INDPIV + NRWPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEFS TO PERFORM FORWARD SOLUTION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.GT.0) CALL ARCEFS(ARRAY(INDEXA), NRWBLK, NCLPIV,
     *   NOVRLP, B(INDPIV), X(INDPIV))
        INDEXA = INDEXA + NOVRLP*NRWBLK
        NRWBLK = MTRSTR(1,K)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEFM TO PERFORM FORWARD MODIFICATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.GT.0) CALL ARCEFM(ARRAY(INDEXA), NRWBLK, NCLPIV,
     *   B(INDPIV), X(INDPIV))
        INDEXA = INDEXA + NRWBLK*NCLPIV
        NCLBLK = MTRSTR(2,K) - NCLPIV
        NOVRLP = MTRSTR(3,K)
        NRWPIV = NCLBLK - NOVRLP
        INDPIV = INDPIV + NCLPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEFE TO PERFORM FORWARD ELIMINATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NRWPIV.GT.0) CALL ARCEFE(ARRAY(INDEXA), NRWBLK, NRWPIV,
     *   PIVOT(INDPIV), B(INDPIV))
   10 CONTINUE
      INDEXB = INDPIV + NRWPIV - 1
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       BACKWARD LOOP
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DO 30 LL=2,NMBLKS
        K = NMBLKS - LL + 1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEBM TO PERFORM BACKWARD MODIFICATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NRWPIV.EQ.0) GO TO 20
        IF (NRWPIV.NE.NCLBLK) CALL ARCEBM(ARRAY(INDEXA), NRWBLK,
     *   NCLBLK, NRWPIV, B(INDPIV), X(INDPIV))
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEBS TO PERFORM BACKWARD SOLUTION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        CALL ARCEBS(ARRAY(INDEXA), NRWBLK, NCLBLK, NRWPIV, B(INDPIV),
     *   X(INDPIV))
   20   CONTINUE
        INDEXA = INDEXA - NRWBLK*NCLPIV
        NRWBLK = MTRSTR(1,K)
        NOVRLP = MTRSTR(3,K)
        INDEXA = INDEXA - NRWBLK*NOVRLP
        INDPIV = INDPIV - NCLPIV
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C          C
C          C    CALL ARCEBE TO PERFORM BACKWARD ELIMINATION
C          C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF (NCLPIV.GT.0) CALL ARCEBE(ARRAY(INDEXA), NRWBLK, NCLPIV,
     *   NOVRLP, PIVOT(INDPIV), X(INDPIV))
        NRWPIV = NRWBLK - NCLPIV
        NCLBLK = NOVRLP + NRWPIV
        INDEXA = INDEXA - NRWBLK*NRWPIV
        INDPIV = INDPIV - NRWPIV
        NCLPIV = MTRSTR(2,K) - NCLBLK
   30 CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       C
C       C       IF ROW ELIMINATIONS WERE DONE IN TOPBLOCK ,CALL
C       C       ARCEBS TO PERFORM BACKWARD SOLUTION
C       C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (NRWPIV.EQ.0) RETURN
      IF (NRWPIV.NE.NCLBLK) CALL ARCEBM(ARRAY(INDEXA), NRWBLK, NCLBLK,
     * NRWPIV, B(INDPIV), X(INDPIV))
      CALL ARCEBS(ARRAY(INDEXA), NRWBLK, NCLBLK, NRWPIV, B(INDPIV),
     * X(INDPIV))
      RETURN
      END
      SUBROUTINE ARCEFE(BLOCK, NRWBLK, NRWPIV, PIVOT, B)                ARC   10
C
C***************************************************************
C
C  A R C E F E  PERFORMS THE FORWARD ELIMINATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, BI, SWAP
      INTEGER PIVOT(NRWPIV), PIVOTI
      DIMENSION BLOCK(NRWBLK,NRWPIV), B(1)
      DO 30 I=1,NRWPIV
        PIVOTI = PIVOT(I)
        IF (PIVOTI.EQ.I) GO TO 10
        SWAP = B(I)
        B(I) = B(PIVOTI)
        B(PIVOTI) = SWAP
   10   CONTINUE
        IF (I.EQ.NRWBLK) RETURN
        BI = B(I)
        IPLUS1 = I + 1
        DO 20 L=IPLUS1,NRWBLK
          B(L) = B(L) - BLOCK(L,I)*BI
   20   CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEFS(BLOCK, NRWBLK, NCLPIV, NOVRLP, B, X)            ARC   10
C
C***************************************************************
C
C  A R C E F S PERFORMS THE FORWARD SOLUTION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NOVRLP), B(1), X(1)
      DO 20 J=1,NCLPIV
        I = NRWBLK - NCLPIV + J
        X(J) = B(J)/BLOCK(I,J)
        IF (I.EQ.NRWBLK) RETURN
        IPLUS1 = I + 1
        LONG = NRWBLK - I
        JPLUS1 = J + 1
        XJ = X(J)
        DO 10 L=1,LONG
          IPLUSL = I + L
          JPLUSL = J + L
          B(JPLUSL) = B(JPLUSL) - BLOCK(IPLUSL,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEFM(BLOCK, NRWBLK, NCLPIV, B, X)                    ARC   10
C
C***************************************************************
C
C  A R C E F M  PERFORMS THE FORWARD MODIFICATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NCLPIV), B(1), X(1)
      NCLPV1 = NCLPIV + 1
      DO 20 J=1,NCLPIV
        XJ = X(J)
        DO 10 L=1,NRWBLK
          NCLPVL = NCLPIV + L
          B(NCLPVL) = B(NCLPVL) - BLOCK(L,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEBM(BLOCK, NRWBLK, NCLBLK, NRWPIV, B, X)            ARC   10
C
C***************************************************************
C
C  A R C E B M  PERFORMS THE BACKWARD MODIFICATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NCLBLK), B(1), X(1)
      NRWPV1 = NRWPIV + 1
      DO 20 J=NRWPV1,NCLBLK
        XJ = X(J)
        DO 10 L=1,NRWPIV
          B(L) = B(L) - BLOCK(L,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEBS(BLOCK, NRWBLK, NCLBLK, NRWPIV, B, X)            ARC   10
C
C***************************************************************
C
C  A R C E B S PERFORMS THE BACKWARD SOLUTION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, B, X, XJ
      DIMENSION BLOCK(NRWBLK,NCLBLK), B(1), X(1)
      DO 20 NJ=1,NRWPIV
        J = NRWPIV - NJ + 1
        X(J) = B(J)/BLOCK(J,J)
        IF (J.EQ.1) RETURN
        JMIN1 = J - 1
        XJ = X(J)
        DO 10 L=1,JMIN1
          B(L) = B(L) - BLOCK(L,J)*XJ
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE ARCEBE(BLOCK, NRWBLK, NCLPIV, NOVRLP, PIVOT, X)        ARC   10
C
C***************************************************************
C
C  A R C E B E PERFORMS THE BACKWARD ELIMINATION STEP IN THE
C  SOLUTION PHASE OF  A R C E C O.
C
C***************************************************************
C
      DOUBLE PRECISION BLOCK, X, DOTPRD, SWAP
      INTEGER PIVOT(NRWBLK), PIVOTJ
      DIMENSION BLOCK(NRWBLK,NOVRLP), X(1)
      DO 40 NJ=1,NCLPIV
        J = NCLPIV + 1 - NJ
        I = NRWBLK + 1 - NJ
        DOTPRD = X(J)
        IF (J.EQ.NOVRLP) GO TO 20
        JPLUS1 = J + 1
        DO 10 J1=JPLUS1,NOVRLP
          DOTPRD = DOTPRD - X(J1)*BLOCK(I,J1)
   10   CONTINUE
   20   CONTINUE
        X(J) = DOTPRD
        PIVOTJ = PIVOT(J)
        IF (PIVOTJ.EQ.J) GO TO 30
        SWAP = X(PIVOTJ)
        X(PIVOTJ) = X(J)
        X(J) = SWAP
   30   CONTINUE
   40 CONTINUE
      RETURN
      END
