SUBROUTINE SLVBLK ( BLOKS, INTEGS, NBLOKS, B, IPIVOT, X, IFLAG ) BLK00100 C THIS PROGRAM SOLVES THE LINEAR SYSTEM A*X = B WHERE A IS AN C ALMOST BLOCK DIAGONAL MATRIX. SUCH ALMOST BLOCK DIAGONAL MATRICES C ARISE NATURALLY IN PIECEWISE POLYNOMIAL INTERPOLATION OR APPROX- C IMATION AND IN FINITE ELEMENT METHODS FOR TWO-POINT BOUNDARY VALUE C PROBLEMS. THE PLU FACTORIZATION METHOD IS IMPLEMENTED HERE TO TAKE C ADVANTAGE OF THE SPECIAL STRUCTURE OF SUCH SYSTEMS FOR SAVINGS IN C COMPUTING TIME AND STORAGE REQUIREMENTS. C C PARAMETERS C BLOKS A ONE-DIMENSIONAL ARRAY, OF LENGTH C SUM( INTEGS(1,I)*INTEGS(2,I) , I = 1,NBLOKS ) C ON INPUT, CONTAINS THE BLOCKS OF THE ALMOST BLOCK DIAGONAL C MATRIX A . THE ARRAY INTEGS (SEE BELOW AND THE EXAMPLE) C DESCRIBES THE BLOCK STRUCTURE. C ON OUTPUT, CONTAINS CORRESPONDINGLY THE PLU FACTORIZATION C OF A (IF IFLAG .NE. 0). CERTAIN OF THE ENTRIES INTO BLOKS C ARE ARBITRARY (WHERE THE BLOCKS OVERLAP). C INTEGS INTEGER ARRAY DESCRIPTION OF THE BLOCK STRUCTURE OF A . C INTEGS(1,I) = NO. OF ROWS OF BLOCK I = NROW C INTEGS(2,I)= NO. OF COLUMNS OF BLOCK I = NCOL C INTEGS(3,I) = NO. OF ELIM. STEPS IN BLOCK I = LAST C I = 1,2,...,NBLOKS C THE LINEAR SYSTEM IS OF ORDER C N = SUM ( INTEGS(3,I) , I=1,...,NBLOKS ), C BUT THE TOTAL NUMBER OF ROWS IN THE BLOCKS IS C NBROWS = SUM( INTEGS(1,I) , I = 1,...,NBLOKS) C NBLOKS NUMBER OF BLOCKS C B RIGHT SIDE OF THE LINEAR SYSTEM, ARRAY OF LENGTH NBROWS. C CERTAIN OF THE ENTRIES ARE ARBITRARY, CORRESPONDING TO C ROWS OF THE BLOCKS WHICH OVERLAP (SEE BLOCK STRUCTURE AND C THE EXAMPLE BELOW). C IPIVOT ON OUTPUT, INTEGER ARRAY CONTAINING THE PIVOTING SEQUENCE C USED. LENGTH IS NBROWS C X ON OUTPUT, CONTAINS THE COMPUTED SOLUTION (IF IFLAG .NE. 0) C LENGTH IS N. C IFLAG ON OUTPUT, INTEGER C = (-1)**(NO. OF INTERCHANGES DURING FACTORIZATION) C IF A IS INVERTIBLE C = 0 IF A IS SINGULAR C C AUXILIARY PROGRAMS C FCBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,SCRTCH,IFLAG) FACTORS THE MATRIX C A , AND IS USED FOR THIS PURPOSE IN SLVBLK. ITS ARGUMENTS C ARE AS IN SLVBLK, EXCEPT FOR C SCRTCH = A WORK ARRAY OF LENGTH MAX(INTEGS(1,I)). C C SBBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,B,X) SOLVES THE SYSTEM A*X = B C ONCE A IS FACTORED. THIS IS DONE AUTOMATICALLY BY SLVBLK C FOR ONE RIGHT SIDE B, BUT SUBSEQUENT SOLUTIONS MAY BE C OBTAINED FOR ADDITIONAL B-VECTORS. THE ARGUMENTS ARE ALL C AS IN SLVBLK. C C DTBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,IFLAG,DETSGN,DETLOG) COMPUTES THE C DETERMINANT OF A ONCE SLVBLK OR FCBLOK HAS DONE THE FACT- C ORIZATION.THE FIRST FIVE ARGUMENTS ARE AS IN SLVBLK. C DETSGN = SIGN OF THE DETERMINANT C DETLOG = NATURAL LOG OF THE DETERMINANT C C ------ BLOCK STRUCTURE OF A ------ C THE NBLOKS BLOCKS ARE STORED CONSECUTIVELY IN THE ARRAY BLOKS . C THE FIRST BLOCK HAS ITS (1,1)-ENTRY AT BLOKS(1), AND, IF THE I-TH C BLOCK HAS ITS (1,1)-ENTRY AT BLOKS(INDEX(I)), THEN C INDEX(I+1) = INDEX(I) + NROW(I)*NCOL(I) . C THE BLOCKS ARE PIECED TOGETHER TO GIVE THE INTERESTING PART OF A C AS FOLLOWS. FOR I = 1,2,...,NBLOKS-1, THE (1,1)-ENTRY OF THE NEXT C BLOCK (THE (I+1)ST BLOCK ) CORRESPONDS TO THE (LAST+1,LAST+1)-ENTRY C OF THE CURRENT I-TH BLOCK. RECALL LAST = INTEGS(3,I) AND NOTE THAT C THIS MEANS THAT C A. EVERY BLOCK STARTS ON THE DIAGONAL OF A . C B. THE BLOCKS OVERLAP (USUALLY). THE ROWS OF THE (I+1)ST BLOCK C WHICH ARE OVERLAPPED BY THE I-TH BLOCK MAY BE ARBITRARILY DE- C FINED INITIALLY. THEY ARE OVERWRITTEN DURING ELIMINATION. C THE RIGHT SIDE FOR THE EQUATIONS IN THE I-TH BLOCK ARE STORED COR- C RESPONDINGLY AS THE LAST ENTRIES OF A PIECE OF B OF LENGTH NROW C (= INTEGS(1,I)) AND FOLLOWING IMMEDIATELY IN B THE CORRESPONDING C PIECE FOR THE RIGHT SIDE OF THE PRECEDING BLOCK, WITH THE RIGHT SIDE C FOR THE FIRST BLOCK STARTING AT B(1) . IN THIS, THE RIGHT SIDE FOR C AN EQUATION NEED ONLY BE SPECIFIED ONCE ON INPUT, IN THE FIRST BLOCK C IN WHICH THE EQUATION APPEARS. C C ------ EXAMPLE AND TEST DRIVER ------ C THE TEST DRIVER FOR THIS PACKAGE CONTAINS AN EXAMPLE, A LINEAR C SYSTEM OF ORDER 11, WHOSE NONZERO ENTRIES ARE INDICATED IN THE FOL- C LOWING SCHEMA BY THEIR ROW AND COLUMN INDEX MODULO 10. NEXT TO IT C ARE THE CONTENTS OF THE INTEGS ARRRAY WHEN THE MATRIX IS TAKEN TO C BE ALMOST BLOCK DIAGONAL WITH NBLOKS = 5, AND BELOW IT ARE THE FIVE C BLOCKS. C C NROW1 = 3, NCOL1 = 4 C 11 12 13 14 C 21 22 23 24 NROW2 = 3, NCOL2 = 3 C 31 32 33 34 C LAST1 = 2 43 44 45 C 53 54 55 NROW3 = 3, NCOL3 = 4 C LAST2 = 3 66 67 68 69 NROW4 = 3, NCOL4 = 4 C 76 77 78 79 NROW5 = 4, NCOL5 = 4 C 86 87 88 89 C LAST3 = 1 97 98 99 90 C LAST4 = 1 08 09 00 01 C 18 19 10 11 C LAST5 = 4 C C ACTUAL INPUT TO BLOKS SHOWN BY ROWS OF BLOCKS OF A . C (THE ** ITEMS ARE ARBITRARY, THIS STORAGE IS USED BY SLVBLK) C C 11 12 13 14 / ** ** ** / 66 67 68 69 / ** ** ** ** / ** ** ** ** C 21 22 23 24 / 43 44 45 / 76 77 78 79 / ** ** ** ** / ** ** ** ** C 31 32 33 34/ 53 54 55/ 86 87 88 89/ 97 98 99 90/ 08 09 00 01 C 18 19 10 11 C C INDEX = 1 INDEX = 13 INDEX = 22 INDEX = 34 INDEX = 46 C C ACTUAL RIGHT SIDE VALUES WITH ** FOR ARBITRARY VALUES C B1 B2 B3 ** B4 B5 B6 B7 B8 ** ** B9 ** ** B10 B11 C C (IT WOULD HAVE BEEN MORE EFFICIENT TO COMBINE BLOCK 3 WITH BLOCK 4) C INTEGER INTEGS(3,NBLOKS),IPIVOT(1),IFLAG REAL BLOKS(1),B(1),X(1) C IN THE CALL TO FCBLOK, X IS USED FOR TEMPORARY STORAGE. CALL FCBLOK(BLOKS,INTEGS,NBLOKS,IPIVOT,X,IFLAG) IF (IFLAG .EQ. 0) RETURN CALL SBBLOK(BLOKS,INTEGS,NBLOKS,IPIVOT,B,X) RETURN END SUBROUTINE FCBLOK ( BLOKS, INTEGS, NBLOKS, IPIVOT, SCRTCH, IFLAG )BLK12700 CALLS SUBROUTINES F A C T R B AND S H I F T B . C C F C B L O K SUPERVISES THE PLU FACTORIZATION WITH PIVOTING OF C SCALED ROWS OF THE ALMOST BLOCK DIAGONAL MATRIX STORED IN THE ARRAYS C B L O K S AND I N T E G S . C C FACTRB = SUBPROGRAM WHICH CARRIES OUT STEPS 1,...,LAST OF GAUSS C ELIMINATION (WITH PIVOTING) FOR AN INDIVIDUAL BLOCK. C SHIFTB = SUBPROGRAM WHICH SHIFTS THE REMAINING ROWS TO THE TOP OF C THE NEXT BLOCK C C PARAMETERS C BLOKS AN ARRAY THAT INITIALLY CONTAINS THE ALMOST BLOCK DIAGONAL C MATRIX A TO BE FACTORED, AND ON RETURN CONTAINS THE COM- C PUTED FACTORIZATION OF A . C INTEGS AN INTEGER ARRAY DESCRIBING THE BLOCK STRUCTURE OF A . C NBLOKS THE NUMBER OF BLOCKS IN A . C IPIVOT AN INTEGER ARRAY OF DIMENSION SUM (INTEGS(1,I) , I=1, C ...,NBLOKS) WHICH, ON RETURN, CONTAINS THE PIVOTING STRA- C TEGY USED. C SCRTCH WORK AREA REQUIRED, OF LENGTH MAX (INTEGS(1,I) , I=1, C ...,NBLOKS). C IFLAG OUTPUT PARAMETER, C = 0 IN CASE MATRIX WAS FOUND TO BE SINGULAR. C OTHERWISE, C = (-1)**(NUMBER OF ROW INTERCHANGES DURING FACTORIZATION) C INTEGER INTEGS(3,NBLOKS),IPIVOT(1),IFLAG, I,INDEX,INDEXB,INDEXN, * LAST,NCOL,NROW REAL BLOKS(1),SCRTCH(1) IFLAG = 1 INDEXB = 1 INDEXN = 1 I = 1 C LOOP OVER THE BLOCKS. I IS LOOP INDEX 10 INDEX = INDEXN NROW = INTEGS(1,I) NCOL = INTEGS(2,I) LAST = INTEGS(3,I) C CARRY OUT ELIMINATION ON THE I-TH BLOCK UNTIL NEXT BLOCK C ENTERS, I.E., FOR COLUMNS 1,...,LAST OF I-TH BLOCK. CALL FACTRB(BLOKS(INDEX),IPIVOT(INDEXB),SCRTCH,NROW,NCOL,LAST, * IFLAG) C CHECK FOR HAVING REACHED A SINGULAR BLOCK OR THE LAST BLOCK IF (IFLAG .EQ. 0 .OR. I .EQ. NBLOKS) * RETURN I = I+1 INDEXN = NROW*NCOL + INDEX C PUT THE REST OF THE I-TH BLOCK ONTO THE NEXT BLOCK CALL SHIFTB(BLOKS(INDEX),IPIVOT(INDEXB),NROW,NCOL,LAST, * BLOKS(INDEXN),INTEGS(1,I),INTEGS(2,I)) INDEXB = INDEXB + NROW GO TO 10 END SUBROUTINE FACTRB ( W, IPIVOT, D, NROW, NCOL, LAST, IFLAG ) BLK18200 C ADAPTED FROM P.132 OF *ELEMENTARY NUMER.ANALYSIS* BY CONTE-DE BOOR C C CONSTRUCTS A PARTIAL PLU FACTORIZATION, CORRESPONDING TO STEPS 1,..., C L A S T IN GAUSS ELIMINATION, FOR THE MATRIX W OF ORDER C ( N R O W , N C O L ), USING PIVOTING OF SCALED ROWS. C C PARAMETERS C W CONTAINS THE (NROW,NCOL) MATRIX TO BE PARTIALLY FACTORED C ON INPUT, AND THE PARTIAL FACTORIZATION ON OUTPUT. C IPIVOT AN INTEGER ARRAY OF LENGTH NROW CONTAINING A RECORD OF THE C PIVOTING STRATEGY USED. ROW IPIVOT(I) IS USED DURING THE C I-TH ELIMINATION STEP, I=1,...,LAST. C D A WORK ARRAY OF LENGTH NROW USED TO STORE ROW SIZES C TEMPORARILY. C NROW NUMBER OF ROWS OF W. C NCOL NUMBER OF COLUMNS OF W. C LAST NUMBER OF ELIMINATION STEPS TO BE CARRIED OUT. C IFLAG ON OUTPUT, EQUALS IFLAG ON INPUT TIMES (-1)**(NUMBER OF C ROW INTERCHANGES DURING THE FACTORIZATION PROCESS), IN C CASE NO ZERO PIVOT WAS ENCOUNTERED. C OTHERWISE, IFLAG = 0 ON OUTPUT. C INTEGER IPIVOT(NROW),NCOL,LAST,IFLAG, I,IPIVI,IPIVK,J,K,KP1 REAL W(NROW,NCOL),D(NROW), AWIKDI,COLMAX,RATIO,ROWMAX C INITIALIZE IPIVOT, D DO 10 I=1,NROW IPIVOT(I) = I ROWMAX = 0. DO 9 J=1,NCOL 9 ROWMAX = AMAX1(ROWMAX, ABS(W(I,J))) IF (ROWMAX .EQ. 0.) GO TO 999 10 D(I) = ROWMAX C GAUSS ELIMINATION WITH PIVOTING OF SCALED ROWS, LOOP OVER K=1,.,LAST K = 1 C AS PIVOT ROW FOR K-TH STEP, PICK AMONG THE ROWS NOT YET USED, C I.E., FROM ROWS IPIVOT(K),...,IPIVOT(NROW), THE ONE WHOSE K-TH C ENTRY (COMPARED TO THE ROW SIZE) IS LARGEST. THEN, IF THIS ROW C DOES NOT TURN OUT TO BE ROW IPIVOT(K), REDEFINE IPIVOT(K) AP- C PROPRIATELY AND RECORD THIS INTERCHANGE BY CHANGING THE SIGN C OF I F L A G . 11 IPIVK = IPIVOT(K) IF (K .EQ. NROW) GO TO 21 J = K KP1 = K+1 COLMAX = ABS(W(IPIVK,K))/D(IPIVK) C FIND THE (RELATIVELY) LARGEST PIVOT DO 15 I=KP1,NROW IPIVI = IPIVOT(I) AWIKDI = ABS(W(IPIVI,K))/D(IPIVI) IF (AWIKDI .LE. COLMAX) GO TO 15 COLMAX = AWIKDI J = I 15 CONTINUE IF (J .EQ. K) GO TO 16 IPIVK = IPIVOT(J) IPIVOT(J) = IPIVOT(K) IPIVOT(K) = IPIVK IFLAG = -IFLAG 16 CONTINUE C IF PIVOT ELEMENT IS TOO SMALL IN ABSOLUTE VALUE, DECLARE C MATRIX TO BE NONINVERTIBLE AND QUIT. IF (ABS(W(IPIVK,K))+D(IPIVK) .LE. D(IPIVK)) * GO TO 999 C OTHERWISE, SUBTRACT THE APPROPRIATE MULTIPLE OF THE PIVOT C ROW FROM REMAINING ROWS, I.E., THE ROWS IPIVOT(K+1),..., C IPIVOT(NROW), TO MAKE K-TH ENTRY ZERO. SAVE THE MULTIPLIER IN C ITS PLACE. DO 20 I=KP1,NROW IPIVI = IPIVOT(I) W(IPIVI,K) = W(IPIVI,K)/W(IPIVK,K) RATIO = -W(IPIVI,K) DO 20 J=KP1,NCOL 20 W(IPIVI,J) = RATIO*W(IPIVK,J) + W(IPIVI,J) K = KP1 C CHECK FOR HAVING REACHED THE NEXT BLOCK. IF (K .LE. LAST) GO TO 11 RETURN C IF LAST .EQ. NROW , CHECK NOW THAT PIVOT ELEMENT IN LAST ROW C IS NONZERO. 21 IF( ABS(W(IPIVK,NROW))+D(IPIVK) .GT. D(IPIVK) ) * RETURN C SINGULARITY FLAG SET 999 IFLAG = 0 RETURN END SUBROUTINE SHIFTB ( AI, IPIVOT, NROWI, NCOLI, LAST, BLK26800 * AI1, NROWI1, NCOLI1 ) C SHIFTS THE ROWS IN CURRENT BLOCK, AI, NOT USED AS PIVOT ROWS, IF C ANY, I.E., ROWS IPIVOT(LAST+1),...,IPIVOT(NROWI), ONTO THE FIRST C MMAX = NROW-LAST ROWS OF THE NEXT BLOCK, AI1, WITH COLUMN LAST+J OF C AI GOING TO COLUMN J , J=1,...,JMAX=NCOLI-LAST. THE REMAINING COL- C UMNS OF THESE ROWS OF AI1 ARE ZEROED OUT. C C PICTURE C C ORIGINAL SITUATION AFTER RESULTS IN A NEW BLOCK I+1 C LAST = 2 COLUMNS HAVE BEEN CREATED AND READY TO BE C DONE IN FACTRB (ASSUMING NO FACTORED BY NEXT FACTRB CALL. C INTERCHANGES OF ROWS) C 1 C X X 1X X X X X X X X C 1 C 0 X 1X X X 0 X X X X C BLOCK I 1 --------------- C NROWI = 4 0 0 1X X X 0 0 1X X X 0 01 C NCOLI = 5 1 1 1 C LAST = 2 0 0 1X X X 0 0 1X X X 0 01 C ------------------------------- 1 1 NEW C 1X X X X X 1X X X X X1 BLOCK C 1 1 1 I+1 C BLOCK I+1 1X X X X X 1X X X X X1 C NROWI1= 5 1 1 1 C NCOLI1= 5 1X X X X X 1X X X X X1 C ------------------------------- 1-------------1 C 1 C INTEGER IPIVOT(NROWI),LAST, IP,J,JMAX,JMAXP1,M,MMAX REAL AI(NROWI,NCOLI),AI1(NROWI1,NCOLI1) MMAX = NROWI - LAST JMAX = NCOLI - LAST IF (MMAX .LT. 1 .OR. JMAX .LT. 1) RETURN C PUT THE REMAINDER OF BLOCK I INTO AI1 DO 10 M=1,MMAX IP = IPIVOT(LAST+M) DO 10 J=1,JMAX 10 AI1(M,J) = AI(IP,LAST+J) IF (JMAX .EQ. NCOLI1) RETURN C ZERO OUT THE UPPER RIGHT CORNER OF AI1 JMAXP1 = JMAX + 1 DO 20 J=JMAXP1,NCOLI1 DO 20 M=1,MMAX 20 AI1(M,J) = 0. RETURN END SUBROUTINE SBBLOK ( BLOKS, INTEGS, NBLOKS, IPIVOT, B, X ) BLK31700 CALLS SUBROUTINES S U B F O R AND S U B B A K . C C SUPERVISES THE SOLUTION (BY FORWARD AND BACKWARD SUBSTITUTION) OF C THE LINEAR SYSTEM A*X = B FOR X, WITH THE PLU FACTORIZATION OF A C ALREADY GENERATED IN F C B L O K . INDIVIDUAL BLOCKS OF EQUATIONS C ARE SOLVED VIA S U B F O R AND S U B B A K . C C PARAMETERS C BLOKS, INTEGS, NBLOKS, IPIVOT ARE AS ON RETURN FROM FCBLOK. C B THE RIGHT SIDE, STORED CORRESPONDING TO THE STORAGE OF C THE EQUATIONS. SEE COMMENTS IN S L V B L K FOR DETAILS. C X SOLUTION VECTOR C INTEGER INTEGS(3,NBLOKS),IPIVOT(1), I,INDEX,INDEXB,INDEXX,J,LAST, * NBP1,NCOL,NROW REAL BLOKS(1),B(1),X(1) C C FORWARD SUBSTITUTION PASS C INDEX = 1 INDEXB = 1 INDEXX = 1 DO 20 I=1,NBLOKS NROW = INTEGS(1,I) LAST = INTEGS(3,I) CALL SUBFOR(BLOKS(INDEX),IPIVOT(INDEXB),NROW,LAST,B(INDEXB), * X(INDEXX)) INDEX = NROW*INTEGS(2,I) + INDEX INDEXB = INDEXB + NROW 20 INDEXX = INDEXX + LAST C C BACK SUBSTITUTION PASS C NBP1 = NBLOKS + 1 DO 30 J=1,NBLOKS I = NBP1 - J NROW = INTEGS(1,I) NCOL = INTEGS(2,I) LAST = INTEGS(3,I) INDEX = INDEX - NROW*NCOL INDEXB = INDEXB - NROW INDEXX = INDEXX - LAST 30 CALL SUBBAK(BLOKS(INDEX),IPIVOT(INDEXB),NROW,NCOL,LAST, * X(INDEXX)) RETURN END SUBROUTINE SUBFOR ( W, IPIVOT, NROW, LAST, B, X ) BLK36400 C CARRIES OUT THE FORWARD PASS OF SUBSTITUTION FOR THE CURRENT BLOCK, C I.E., THE ACTION ON THE RIGHT SIDE CORRESPONDING TO THE ELIMINATION C CARRIED OUT IN F A C T R B FOR THIS BLOCK. C AT THE END, X(J) CONTAINS THE RIGHT SIDE OF THE TRANSFORMED C IPIVOT(J)-TH EQUATION IN THIS BLOCK, J=1,...,NROW. THEN, SINCE C FOR I=1,...,NROW-LAST, B(NROW+I) IS GOING TO BE USED AS THE RIGHT C SIDE OF EQUATION I IN THE NEXT BLOCK (SHIFTED OVER THERE FROM C THIS BLOCK DURING FACTORIZATION), IT IS SET EQUAL TO X(LAST+I) HERE. C C PARAMETERS C W, IPIVOT, NROW, LAST ARE AS ON RETURN FROM FACTRB. C B(J) IS EXPECTED TO CONTAIN, ON INPUT, THE RIGHT SIDE OF J-TH C EQUATION FOR THIS BLOCK, J=1,...,NROW. C B(NROW+J) CONTAINS, ON OUTPUT, THE APPROPRIATELY MODIFIED RIGHT C SIDE FOR EQUATION J IN NEXT BLOCK, J=1,...,NROW-LAST. C X(J) CONTAINS, ON OUTPUT, THE APPROPRIATELY MODIFIED RIGHT C SIDE OF EQUATION IPIVOT(J) IN THIS BLOCK, J=1,...,LAST (AND C EVEN FOR J=LAST+1,...,NROW). C INTEGER IPIVOT(NROW), IP,JMAX,K C DIMENSION B(NROW + NROW-LAST) REAL W(NROW,LAST),B(1),X(NROW) IP = IPIVOT(1) X(1) = B(IP) IF (NROW .EQ. 1) GO TO 99 DO 15 K=2,NROW IP = IPIVOT(K) JMAX = AMIN0(K-1,LAST) SUM = 0. DO 14 J=1,JMAX 14 SUM = W(IP,J)*X(J) + SUM 15 X(K) = B(IP) - SUM C C TRANSFER MODIFIED RIGHT SIDES OF EQUATIONS IPIVOT(LAST+1),..., C IPIVOT(NROW) TO NEXT BLOCK. NROWML = NROW - LAST IF (NROWML .EQ. 0) GO TO 99 LASTP1 = LAST+1 DO 25 K=LASTP1,NROW 25 B(NROWML+K) = X(K) 99 RETURN END SUBROUTINE SUBBAK ( W, IPIVOT, NROW, NCOL, LAST, X ) BLK40700 C CARRIES OUT BACKSUBSTITUTION FOR CURRENT BLOCK. C C PARAMETERS C W, IPIVOT, NROW, NCOL, LAST ARE AS ON RETURN FROM FACTRB. C X(1),...,X(NCOL) CONTAINS, ON INPUT, THE RIGHT SIDE FOR THE C EQUATIONS IN THIS BLOCK AFTER BACKSUBSTITUTION HAS BEEN C CARRIED UP TO BUT NOT INCLUDING EQUATION IPIVOT(LAST). C MEANS THAT X(J) CONTAINS THE RIGHT SIDE OF EQUATION IPI- C VOT(J) AS MODIFIED DURING ELIMINATION, J=1,...,LAST, WHILE C FOR J .GT. LAST, X(J) IS ALREADY A COMPONENT OF THE SOLUT- C ION VECTOR. C X(1),...,X(NCOL) CONTAINS, ON OUTPUT, THE COMPONENTS OF THE SOLUT- C ION CORRESPONDING TO THE PRESENT BLOCK. C INTEGER IPIVOT(NROW),LAST, IP,J,K,KP1 REAL W(NROW,NCOL),X(NCOL), SUM K = LAST IP = IPIVOT(K) SUM = 0. IF (K .EQ. NCOL) GO TO 4 KP1 = K+1 2 DO 3 J=KP1,NCOL 3 SUM = W(IP,J)*X(J) + SUM 4 X(K) = (X(K) - SUM)/W(IP,K) IF (K .EQ. 1) RETURN KP1 = K K = K-1 IP = IPIVOT(K) SUM = 0. GO TO 2 END SUBROUTINE DTBLOK ( BLOKS, INTEGS, NBLOKS, IPIVOT, IFLAG, BLK43900 * DETSGN, DETLOG ) C COMPUTES THE DETERMINANT OF AN ALMOST BLOCK DIAGONAL MATRIX WHOSE C PLU FACTORIZATION HAS BEEN OBTAINED PREVIOUSLY IN FCBLOK. C *** THE LOGARITHM OF THE DETERMINANT IS COMPUTED INSTEAD OF THE C DETERMINANT ITSELF TO AVOID THE DANGER OF OVERFLOW OR UNDERFLOW C INHERENT IN THIS CALCULATION. C C PARAMETERS C BLOKS, INTEGS, NBLOKS, IPIVOT, IFLAG ARE AS ON RETURN FROM FCBLOK. C IN PARTICULAR, IFLAG = (-1)**(NUMBER OF INTERCHANGES DUR- C ING FACTORIZATION) IF SUCCESSFUL, OTHERWISE IFLAG = 0. C DETSGN ON OUTPUT, CONTAINS THE SIGN OF THE DETERMINANT. C DETLOG ON OUTPUT, CONTAINS THE NATURAL LOGARITHM OF THE DETERMI- C NANT IF DETERMINANT IS NOT ZERO. OTHERWISE CONTAINS 0. C INTEGER INTEGS(3,NBLOKS),IPIVOT(1),IFLAG, I,INDEXP,IP,K,LAST REAL BLOKS(1),DETSGN,DETLOG C DETSGN = IFLAG DETLOG = 0. IF (IFLAG .EQ. 0) RETURN INDEX = 0 INDEXP = 0 DO 2 I=1,NBLOKS NROW = INTEGS(1,I) LAST = INTEGS(3,I) DO 1 K=1,LAST IP = INDEX + NROW*(K-1) + IPIVOT(INDEXP+K) DETLOG = DETLOG + ALOG(ABS(BLOKS(IP))) 1 DETSGN = DETSGN*SIGN(1.,BLOKS(IP)) INDEX = NROW*INTEGS(2,I) + INDEX 2 INDEXP = INDEXP + NROW RETURN END C TEST PROGRAM FOR THE SOLVEBLOK PACKAGE SOLVES AN 11TH ORDER LINEAR BLK47400 C ALMOST BLOCK DIAGONAL SYSTEM. BLK47500 DIMENSION BLOKS(61),B(16),X(11), IPIVOT(16),INTEGS(15) BLK47600 C NROW NCOL LAST BLK47700 DATA INTEGS/ 3 , 4 , 2 BLK47800 2 , 3 , 3 , 3 BLK47900 3 , 3 , 4 , 1 BLK48000 4 , 3 , 4 , 1 BLK48100 5 , 4 , 4 , 4/ BLK48200 DATA BLOKS /.1,.2,-1., 2.,-.2,.3, -.1,-.2,-.3, -.1,4.,.3 BLK48300 2 ,0.,-.4,3., 0.,.4,.5, 0.,-5.,-.5 BLK48400 3 ,.6,.5,3., -.6,4.,.4, -.6,.5,-.4, 5.,-.5,.4 BLK48500 4 ,0.,0.,.3, 0.,0.,-.3, 0.,0.,.3, 0.,0.,7. BLK48600 5 ,0.,0.,.2,6., 0.,0.,-.2,.1, 0.,0.,-.2,-.1, BLK48700 5 0.,0.,8.,-.1/ BLK48800 DATA B /1.94,3.04,-.83, 0.,-3.54,2.75, 1.32,2.35,1.96, BLK48900 2 2*0.,1.52, 2*0.,.78,2.40 / BLK49000 NBLOKS = 5 BLK49100 N = 11 BLK49200 CALL SLVBLK ( BLOKS, INTEGS, NBLOKS, B, IPIVOT, X, IFLAG ) BLK49300 ERROR = 0. BLK49400 DO 10 I=1,N BLK49500 B(I) = FLOAT(12-I)/10. - X(I) BLK49600 10 ERROR = AMAX1(ERROR,ABS(B(I))) BLK49700 WRITE (6,610) IFLAG,ERROR, (IPIVOT(I),I=1,16), (I,X(I),B(I),I=1,N)BLK49800 610 FORMAT (28H CHECKOUT SOLVEBLOK ROUTINES/ BLK49900 * 8H IFLAG =,I2,17H, MAXIMUM ERROR =,E9.4// BLK50000 * 9H IPIVOT =,5(2X,3I2),I2// BLK50100 * 27H I X(I) ERROR(I)/(I3,F11.5,E13.5)) BLK50200 CALL DTBLOK ( BLOKS, INTEGS, NBLOKS, IPIVOT, IFLAG, BLK50300 * DETSGN, DETLOG ) BLK50400 DET = DETSGN*EXP(DETLOG) BLK50500 ERROR = DET - 2.418821899E6 BLK50600 WRITE (6,611) DET,ERROR BLK50700 611 FORMAT(//14H DETERMINANT =,E15.8,11H ERROR = ,E11.5) BLK50800 STOP BLK50900 END BLK51000