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