C ALGORITHM 704, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 2, PP. 205-210. Tape to accompany the paper: The solution of almost block diagonal linear systems arising in spline collocation at Gaussian points with monomial bases. F. Majaess, P. Keast, G. Fairweather. The tape contains the following files, repeated three times: ABBABD-DOC This file. ABBDRDP.f The driving program for ABBPACK in double precision. ABBDRSP.f The driving program for ABBPACK in single precision. ABBLEQDP.f The package ABBPACK in double presision. ABBLEQSP.f The package ABBPACK in single precision. ABDDRDP.f The driving program for ABDPACK in double precision. ABDDRSP.f The driving program for ABDPACK in single precision. ABDLEQDP.f The package ABDPACK in double presision. ABDLEQSP.f The package ABDPACK in single precision. The files are separated by 3 lines of the following type: C C####################################################################### C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS DRIVING PROGRAM RUNS THE DOUBLE PRECISION PACKAGE ABBPACK C ON A SAMPLE PROBLEM. C C THE SAMPLE SYSTEM GIVEN IN THE DATA STATEMENTS CONSISTS OF C A SYSTEM OF 12 EQUATIONS. THERE ARE TWO MAIN BLOCKS C (CORRESPONDING TO TWO SUBINTERVALS IN THE COLLOCATION). C C THE COEFFICIENT MATRIX OF THE SYSTEM IS (WITH ONLY THE C ELEMENTS USED IN ABBLEQ SHOWN): C C 1 1 2 3 C 1 2 6 3 TOPBLK(:,:) C----------------------------------------------------------------------- C 1 2 3 4 1 0 0 0 C 2 3 4 5 0 1 0 0 C 3 4 5 1 0 0 1 0 C 4 5 1 2 0 0 0 1 GAMMA(:,:,1), WORK(:,:,1) C----------------------------------------------------------------------- C 1 2 3 4 1 0 0 0 C 2 3 4 5 0 1 0 0 C 3 4 5 1 0 0 1 0 C 4 5 1 2 0 0 0 1 GAMMA(:,:,2), WORK(:,:,2) C----------------------------------------------------------------------- C 1 2 1 1 C 2 3 40 10 BOTBLK(:,:) C----------------------------------------------------------------------- C C THE PROGRAM ABBLEQ CALLS THE FACTORIZATION PROCEDURE C ABBDEC AND THEN THE SOLVE PROCEDURE ABBSOL, PROVIDED IFLAG C IS ZERO. THE CODE FOR THE SEPARATE CALLS IS INCLUDED AS C COMMENTS AFTER THE CALL TO ABBLEQ. C C NOTE THAT THE ELEMENTS OF THE ARRAY WORK ARE NOT INITIALIZED C IN THE DRIVING PROGRAM. C C THE OUTPUT FROM THIS DRIVER ON A SUN 4/280 UNDER UNIX, USING F77 C IN DOUBLE PRECISION WAS : C C IFLAG = 0 C PIVOT VECTOR : C 4, 3, 6, 6, 7, 8, 8, 10, 11, 12, 12 C SOLUTION : C 0.10000000000000013D+01 0.99999999999999878D+00 C 0.10000000000000002D+01 0.10000000000000000D+01 C 0.10000000000000002D+01 0.99999999999999900D+00 C 0.10000000000000007D+01 0.10000000000000000D+01 C 0.99999999999999922D+00 0.10000000000000000D+01 C 0.99999999999999967D+00 0.10000000000000018D+01 C C IT SHOULD BE NOTED THAT PIVOTING TOOK PLACE IN THIS EXAMPLE C AT EVERY STAGE OF THE ALGORITHM, AS SHOWN BY THE PIVOT VECTOR. C C COLUMNS 1,4 WERE INTERCHANGED IN TOPBLK, FOLLOWED BY COLUMNS C 2, 3. C THEN ROWS 3,6 WERE INTERCHANGED IN THE FIRST BLOCK OF GAMMA, C FOLLOWED BY ROWS 4,6. C COLUMNS 5,7 WERE INTERCHANGED IN THE FIRST BLOCK OF WORK, C FOLLOWED BY COLUMNS 6,8. C THEN THE ROW PAIRS 7,8 AND 8,10 WERE INTERCHANGED IN THE C SECOND BLOCK OF GAMMA, FOLLOWED BY COLUMN PAIRS 9,11 AND C 10,12 IN THE SECOND BLOCK OF WORK. C FINALLY ROWS 11,12 IN BOTBLK WERE INTERCHANGED. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DOUBLE PRECISION TOPBLK(2,4),GAMMA(4,4,2),WORK(4,4,2), * BOTBLK(2,4),B(12),X(12) INTEGER N,NRWTOP,NOVRLP,NBLOKS,PIVOT(12),IFLAG INTEGER I,K,NOUT DATA NOUT/6/ DATA N,NRWTOP,NOVRLP,NBLOKS/12,2,4,2/ DATA TOPBLK / * 1.0D0, 1.0D0, 1.0D0, 2.0D0, 2.0D0, 6.0D0, 3.0D0, 3.0D0 / C THE BLOCKS OF GAMMA ARE GIVEN BY COLUMNS AS FOLLOWS: DATA GAMMA / * 1.0D0, 2.0D0, 3.0D0, 4.0D0, 2.0D0, 3.0D0, 4.0D0, 5.0D0, * 3.0D0, 4.0D0, 5.0D0, 1.0D0, 4.0D0, 5.0D0, 1.0D0, 2.0D0, * 1.0D0, 2.0D0, 3.0D0, 4.0D0, 2.0D0, 3.0D0, 4.0D0, 5.0D0, * 3.0D0, 4.0D0, 5.0D0, 1.0D0, 4.0D0, 5.0D0, 1.0D0, 2.0D0 / DATA BOTBLK / * 1.0D0, 2.0D0, 2.0D0, 3.0D0, 1.0D0, 40.0D0, 1.0D0, 10.0D0 / C RIGHT HAND SIDE B : DATA B / * 7.0D0, 12.0D0, 11.0D0, 15.0D0, 14.0D0, 13.0D0, * 11.0D0, 15.0D0, 14.0D0, 13.0D0, 5.0D0, 55.0D0 / CALL ABBLEQ(N,TOPBLK,NRWTOP,NOVRLP,NBLOKS,GAMMA,WORK,BOTBLK, * PIVOT,B,X,IFLAG) C C THE SAME CALCULATION COULD BE DONE BY THE FOLLOWING SEPARATE C CALLS TO ABBDEC AND ABBSOL. C C FIRST, CALL THE DECOMPOSITION ROUTINE ABBDEC. C C CALL ABBDEC(N,TOPBLK,NRWTOP,NOVRLP,NBLOKS,GAMMA,WORK, C * BOTBLK,PIVOT,IFLAG) C UNLESS A IS NUMERICALLY SINGULAR OR THE PARAMETERS ARE C IMPROPERLY DEFINED, CALL ABBSOL. C C IF ( IFLAG .EQ. 0 ) THEN C CALL ABBSOL(N,TOPBLK,NRWTOP,NOVRLP, C * NBLOKS,GAMMA,WORK,BOTBLK,PIVOT,B,X) C ENDIF WRITE(NOUT,1000)IFLAG IF ( IFLAG .EQ. 0 ) THEN WRITE(NOUT,2000) WRITE(NOUT,3000)(PIVOT(I),I=1,N-1) WRITE(NOUT,4000) WRITE(NOUT,5000)(X(K),K=1,N) ELSE IF ( IFLAG .EQ. 1 ) WRITE(NOUT,7000) IF ( IFLAG .EQ. -1) WRITE(NOUT,6000) ENDIF STOP 1000 FORMAT(8HIFLAG = ,I4) 2000 FORMAT(15HPIVOT VECTOR : ) 3000 FORMAT(15I4) 4000 FORMAT(11H SOLUTION :) 5000 FORMAT(2D25.17) 6000 FORMAT(29HMATRIX APPEARS TO BE SINGULAR) 7000 FORMAT(20H IMPROPER PARAMETERS) END C C####################################################################### C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS DRIVING PROGRAM RUNS THE SINGLE PRECISION PACKAGE ABBPACK C ON A SAMPLE PROBLEM. C C THE SAMPLE SYSTEM GIVEN IN THE DATA STATEMENTS CONSISTS OF C A SYSTEM OF 12 EQUATIONS. THERE ARE TWO MAIN BLOCKS C (CORRESPONDING TO TWO SUBINTERVALS IN THE COLLOCATION). C C THE COEFFICIENT MATRIX OF THE SYSTEM IS (WITH ONLY THE C ELEMENTS USED IN ABBLEQ SHOWN): C C 1 1 2 3 C 1 2 6 3 TOPBLK(:,:) C----------------------------------------------------------------------- C 1 2 3 4 1 0 0 0 C 2 3 4 5 0 1 0 0 C 3 4 5 1 0 0 1 0 C 4 5 1 2 0 0 0 1 GAMMA(:,:,1), WORK(:,:,1) C----------------------------------------------------------------------- C 1 2 3 4 1 0 0 0 C 2 3 4 5 0 1 0 0 C 3 4 5 1 0 0 1 0 C 4 5 1 2 0 0 0 1 GAMMA(:,:,2), WORK(:,:,2) C----------------------------------------------------------------------- C 1 2 1 1 C 2 3 40 10 BOTBLK(:,:) C----------------------------------------------------------------------- C C THE PROGRAM ABBLEQ CALLS THE FACTORIZATION PROCEDURE C ABBDEC AND THEN THE SOLVE PROCEDURE ABBSOL, PROVIDED IFLAG C IS ZERO. THE CODE FOR THE SEPARATE CALLS IS INCLUDED AS C COMMENTS AFTER THE CALL TO ABBLEQ. C C NOTE THAT THE ELEMENTS OF THE ARRAY WORK ARE NOT INITIALIZED C IN THE DRIVING PROGRAM. C C THE OUTPUT FROM THIS DRIVER ON A SUN 4/280 UNDER UNIX, USING F77. C IN SINGLE PRECISION WAS : C C IFLAG = 0 C PIVOT VECTOR : C 4, 3, 6, 6, 7, 8, 8, 10, 11, 12, 12 C SOLUTION : C 0.99999756E+00 0.10000019E+01 C 0.99999952E+00 0.10000004E+01 C 0.99999857E+00 0.99999958E+00 C 0.10000012E+01 0.10000001E+01 C 0.99999774E+00 0.99999869E+00 C 0.99999875E+00 0.10000057E+01 C C IT SHOULD BE NOTED THAT PIVOTING TOOK PLACE IN THIS EXAMPLE C AT EVERY STAGE OF THE ALGORITHM, AS SHOWN BY THE PIVOT VECTOR. C C COLUMNS 1,4 WERE INTERCHANGED IN TOPBLK, FOLLOWED BY COLUMNS C 2, 3. C THEN ROWS 3,6 WERE INTERCHANGED IN THE FIRST BLOCK OF GAMMA, C FOLLOWED BY ROWS 4,6. C COLUMNS 5,7 WERE INTERCHANGED IN THE FIRST BLOCK OF WORK, C FOLLOWED BY COLUMNS 6,8. C THEN THE ROW PAIRS 7,8 AND 8,10 WERE INTERCHANGED IN THE C SECOND BLOCK OF GAMMA, FOLLOWED BY COLUMN PAIRS 9,11 AND C 10,12 IN THE SECOND BLOCK OF WORK. C FINALLY ROWS 11,12 IN BOTBLK WERE INTERCHANGED. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REAL TOPBLK(2,4),GAMMA(4,4,2),WORK(4,4,2), * BOTBLK(2,4),B(12),X(12) INTEGER N,NRWTOP,NOVRLP,NBLOKS,PIVOT(12),IFLAG INTEGER I,K,NOUT DATA NOUT/6/ DATA N,NRWTOP,NOVRLP,NBLOKS/12,2,4,2/ DATA TOPBLK / * 1.0E0, 1.0E0, 1.0E0, 2.0E0, 2.0E0, 6.0E0, 3.0E0, 3.0E0 / C THE BLOCKS OF GAMMA ARE GIVEN BY COLUMNS AS FOLLOWS: DATA GAMMA / * 1.0E0, 2.0E0, 3.0E0, 4.0E0, 2.0E0, 3.0E0, 4.0E0, 5.0E0, * 3.0E0, 4.0E0, 5.0E0, 1.0E0, 4.0E0, 5.0E0, 1.0E0, 2.0E0, * 1.0E0, 2.0E0, 3.0E0, 4.0E0, 2.0E0, 3.0E0, 4.0E0, 5.0E0, * 3.0E0, 4.0E0, 5.0E0, 1.0E0, 4.0E0, 5.0E0, 1.0E0, 2.0E0 / DATA BOTBLK / * 1.0E0, 2.0E0, 2.0E0, 3.0E0, 1.0E0, 40.0E0, 1.0E0, 10.0E0 / C RIGHT HAND SIDE B : DATA B / * 7.0E0, 12.0E0, 11.0E0, 15.0E0, 14.0E0, 13.0E0, * 11.0E0, 15.0E0, 14.0E0, 13.0E0, 5.0E0, 55.0E0 / CALL ABBLEQ(N,TOPBLK,NRWTOP,NOVRLP,NBLOKS,GAMMA,WORK,BOTBLK, * PIVOT,B,X,IFLAG) C C THE SAME CALCULATION COULD BE DONE BY THE FOLLOWING SEPARATE C CALLS TO ABBDEC AND ABBSOL. C C FIRST, CALL THE DECOMPOSITION ROUTINE ABBDEC. C C CALL ABBDEC(N,TOPBLK,NRWTOP,NOVRLP,NBLOKS,GAMMA,WORK, C * BOTBLK,PIVOT,IFLAG) C UNLESS A IS NUMERICALLY SINGULAR OR THE PARAMETERS ARE C IMPROPERLY DEFINED, CALL ABBSOL. C C IF ( IFLAG .EQ. 0 ) THEN C CALL ABBSOL(N,TOPBLK,NRWTOP,NOVRLP, C * NBLOKS,GAMMA,WORK,BOTBLK,PIVOT,B,X) C ENDIF WRITE(NOUT,1000)IFLAG IF ( IFLAG .EQ. 0 ) THEN WRITE(NOUT,2000) WRITE(NOUT,3000)(PIVOT(I),I=1,N-1) WRITE(NOUT,4000) WRITE(NOUT,5000)(X(K),K=1,N) ELSE IF ( IFLAG .EQ. 1 ) WRITE(NOUT,7000) IF ( IFLAG .EQ. -1) WRITE(NOUT,6000) ENDIF STOP 1000 FORMAT(8HIFLAG = ,I4) 2000 FORMAT(15HPIVOT VECTOR : ) 3000 FORMAT(15I4) 4000 FORMAT(11H SOLUTION :) 5000 FORMAT(2E16.8) 6000 FORMAT(29HMATRIX APPEARS TO BE SINGULAR) 7000 FORMAT(20H IMPROPER PARAMETERS) END C C####################################################################### C SUBROUTINE ABBLEQ(N, TOPBLK, NRWTOP, NOVRLP, NBLOKS, GAMMA, * WORK, BOTBLK, PIVOT, B, X, IFLAG) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A (WITH THE ALMOST BLOCK C BI-DIAGONAL STRUCTURE ARISING WHEN CONDENSATION IS APPLIED TO THE C MATRIX OBTAINED WHEN SPLINE COLLOCATION AT GAUSSIAN POINTS WITH C MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT BOUNDARY C VALUE PROBLEMS WITH SEPARATED BOUNDARY CONDITIONS), TOGETHER WITH C THE CORRESPONDING RIGHT HAND SIDE. C C THIS PARTICULAR MATRIX STRUCTURE ALSO ARISES IN CERTAIN C IMPLICIT RUNGE-KUTTA APPLICATIONS, AND IN MULTIPLE SHOOTING. C C ABBLEQ CALLS A DECOMPOSITION ROUTINE, ABBDEC, FOLLOWED BY THE C CORRESPONDING SOLVE ROUTINE, ABBSOL, TO OBTAIN THE SOLUTION C OF THE SYSTEM. THE ALGORITHM IMPLEMENTS AN ALTERNATE COLUMN C AND ROW PIVOTING TECHNIQUE. C C THIS IS THE DOUBLE PRECISION VERSION OF ABBLEQ. THE SINGLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL DOUBLE PRECISION C DECLARATIONS BY REAL, AND BY REPLACING ALL REFERENCES TO C DOUBLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C THEIR SINGLE PRECISION EQUIVALENTS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C (NBLOKS + 1) * NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - DOUBLE PRECISION (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSITION ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE ORDER C OF GAMMA(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN GAMMA. C UNCHANGED ON EXIT. C C GAMMA - DOUBLE PRECISION (NOVRLP, NOVRLP, NBLOKS) C GAMMA(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP BLOCK OF GAMMA. C CONTAINS THE DECOMPOSITION ON EXIT. C C WORK - DOUBLE PRECISION (NOVRLP, NOVRLP, NBLOKS). C WORKSPACE USED TO STORE THE IDENTITY BLOCKS. C NOT INITIALIZED BY THE USER. C CONTAINS DECOMPOSITION ON EXIT. C C BOTBLK - DOUBLE PRECISION(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSITION ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C RECORDS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION ON EXIT. C C B - DOUBLE PRECISION (N). C THE RIGHT HAND SIDE. C UNCHANGED ON EXIT. C C X - DOUBLE PRECISION (N). C WORK SPACE FOR THE SOLUTION VECTOR. C IF IFLAG = 0, CONTAINS THE SOLUTION VECTOR ON EXIT. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C ROUTINES CALLED : C C ABBDEC, ABBSOL. C C CALLING SEQUENCES : C C CALL ABBDEC(N, TOPBLK, NRWTOP, NOVRLP, C * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, IFLAG) C IF ( IFLAG .NE. 0 ) RETURN C CALL ABBSOL(N, TOPBLK, NRWTOP, NOVRLP, C * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, B, X) C C IN ADDITION, THE ROUTINES ABBDEC AND ABBSOL CALL THE FOLLOWING C ROUTINES FROM THE BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) : C C IDAMAX, DSWAP, DSCAL, DAXPY, DDOT. C C TO CONVERT TO SINGLE PRECISION, THESE MUST BE REPLACED BY C CALLS TO : C C ISAMAX, SSWAP, SSCAL, SAXPY, SDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NRWTOP,NOVRLP,NBLOKS,PIVOT(N),IFLAG * DOUBLE PRECISION TOPBLK,GAMMA,WORK,BOTBLK,B,X DIMENSION TOPBLK(NRWTOP,*),GAMMA(NOVRLP,NOVRLP,*), * WORK(NOVRLP,NOVRLP,*),BOTBLK(NOVRLP-NRWTOP,*), * B(N),X(N) CALL ABBDEC(N,TOPBLK,NRWTOP,NOVRLP, * NBLOKS,GAMMA,WORK,BOTBLK,PIVOT,IFLAG) C C IF SINGULARITY WAS NOT DETECTED, CALL THE SOLVE ROUTINE ABBSOL. C IF ( IFLAG .NE. 0 ) RETURN CALL ABBSOL(N,TOPBLK,NRWTOP,NOVRLP, * NBLOKS,GAMMA,WORK,BOTBLK,PIVOT,B,X) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABBLEQ. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC END * * * * * SUBROUTINE ABBDEC(N, TOPBLK, NRWTOP, NOVRLP, * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, IFLAG) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A (WITH THE ALMOST BLOCK C BI-DIAGONAL STRUCTURE ARISING WHEN CONDENSATION IS APPLIED TO THE C MATRIX OBTAINED WHEN SPLINE COLLOCATION AT GAUSSIAN POINTS WITH C MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT BOUNDARY C VALUE PROBLEMS WITH SEPARATED BOUNDARY CONDITIONS). C C THIS PARTICULAR MATRIX STRUCTURE ALSO ARISES IN CERTAIN C IMPLICIT RUNGE-KUTTA APPLICATIONS, AND IN MULTIPLE SHOOTING. C C ABBDEC PERFORMS THE DECOMPOSITION, USING AN ALTERNATE COLUMN C AND ROW PIVOTING TECHNIQUE. C C THIS IS THE DOUBLE PRECISION VERSION OF ABBDEC. THE SINGLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL DOUBLE PRECISION C DECLARATIONS BY REAL, AND BY REPLACING ALL REFERENCES TO C DOUBLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C THEIR SINGLE PRECISION EQUIVALENTS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C (NBLOKS + 1) * NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - DOUBLE PRECISION (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSITION ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE ORDER C OF GAMMA(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN GAMMA. C UNCHANGED ON EXIT. C C GAMMA - DOUBLE PRECISION (NOVRLP, NOVRLP, NBLOKS) C GAMMA(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP BLOCK OF GAMMA. C CONTAINS THE DECOMPOSITION ON EXIT. C C WORK - DOUBLE PRECISION (NOVRLP, MRWBLK, NBLOKS). C WORKSPACE USED TO STORE THE IDENTITY BLOCKS. C NOT INITIALIZED BY THE USER. C CONTAINS THE DECOMPOSITION ON EXIT. C C BOTBLK - DOUBLE PRECISION(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSITION ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C RECORDS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION ON EXIT. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C ROUTINES CALLED : C C NONE C C ABBDEC USES THE FOLLOWING ROUTINES FROM THE BASIC LINEAR C ALGEBRA LIBRARY (BLAS) : C C IDAMAX, DSWAP, DSCAL, DAXPY. C C TO CONVERT TO SINGLE PRECISION, THESE MUST BE REPLACED BY C CALLS TO : C C ISAMAX, SSWAP, SSCAL, SAXPY. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NRWTOP,NOVRLP,NBLOKS,PIVOT(N),IFLAG * DOUBLE PRECISION TOPBLK,GAMMA,WORK,BOTBLK DIMENSION TOPBLK(NRWTOP,*),GAMMA(NOVRLP,NOVRLP,*), * WORK(NOVRLP,NOVRLP,*),BOTBLK(NOVRLP-NRWTOP,*) C C LOCAL VARIABLES : C INTEGER I,J,K,IPLUS1,JPLUS1,L,INCR,KPLUS1, * INCRJ,INCRN,IPVT,JPVT,IMINN,JMINN,LOOP, * JPVT2,NRWTP1,NRWTPI,NRWBOT,NRWBT1,IDAMAX * DOUBLE PRECISION COLMAX,ROWMAX,PIVMAX, * COLPIV,ROWPIV,COLMLT,ROWMLT,ZERO,ONE DATA ZERO /0.D0/, ONE /1.D0/ C C DEFINE THE CONSTANTS USED THROUGHOUT. C NRWTP1 = NRWTOP+1 NRWBOT = NOVRLP-NRWTOP NRWBT1 = NRWBOT + 1 IFLAG = 0 PIVMAX = ZERO C C CHECK THE VALIDITY OF THE INPUT PARAMETERS. C IF(N.NE.(NBLOKS+1)*NOVRLP) GO TO 930 IF(NRWTOP.EQ.0) GO TO 930 IF(NRWTOP.EQ.NOVRLP) GO TO 930 DO 30 K=1,NBLOKS DO 20 J = 1,NOVRLP DO 10 I = 1,NOVRLP WORK(I,J,K) = ZERO 10 CONTINUE WORK(J,J,K) = ONE 20 CONTINUE 30 CONTINUE C C C FIRST, IN TOPBLK C C APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN PIVOTING. C C DO 195 I = 1,NRWTOP IPLUS1 = I+1 C C DETERMINE COLUMN PIVOT AND PIVOT INDEX. C IPVT = IDAMAX(NOVRLP-I+1,TOPBLK(I,I),NRWTOP) + I - 1 COLMAX = ABS(TOPBLK(I,IPVT)) C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX .EQ. PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(I) = IPVT IF(IPVT.NE.I) THEN C C INTERCHANGE COLUMNS IN TOPBLK. C CALL DSWAP(NRWTOP-I+1,TOPBLK(I,I),1,TOPBLK(I,IPVT),1) C C INTERCHANGE COLUMNS IN GAMMA(:,:,1). C CALL DSWAP(NOVRLP,GAMMA(1,I,1),1,GAMMA(1,IPVT,1),1) ENDIF C C COMPUTE MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C COLPIV = ONE/TOPBLK(I,I) CALL DSCAL(NOVRLP-I,COLPIV,TOPBLK(I,IPLUS1),NRWTOP) DO 190 J = IPLUS1,NOVRLP COLMLT = TOPBLK(I,J) CALL DAXPY(NRWTOP-I,-COLMLT,TOPBLK(IPLUS1,I),1, * TOPBLK(IPLUS1,J),1) C C ADJUST COLUMNS IN GAMMA(:,:,1). C CALL DAXPY(NOVRLP,-COLMLT,GAMMA(1,I,1),1,GAMMA(1,J,1),1) 190 CONTINUE 195 CONTINUE INCR = NRWTOP C C C C IN EACH BLOCK... C C IF(K.LT.NBLOKS) C APPLY NNRWBOT ROW ELIMINATIONS WITH ROW PIVOTING ON GAMMA(:,:,K) C NRWTOP COLUMN ELIMINATIONS WITH COLUMN PIVOTING ON WORK(:,:,K) C DO 745 K = 1,NBLOKS KPLUS1 = K+1 C C C PERFORM NRWBOT ROW ELIMINATION WITH ROW PIVOTING C ON GAMMA(:,:,K). C C DO 530 J = NRWTP1,NOVRLP JPLUS1 = J+1 JMINN = J-NRWTOP IPVT = JMINN INCRJ = INCR + JMINN LOOP = JMINN + 1 IPVT = IDAMAX(NOVRLP-JMINN+1,GAMMA(JMINN,J,K),1)+JMINN-1 ROWMAX = ABS(GAMMA(IPVT,J,K)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C C IF NECESSARY, INTERCHANGE ROWS IN GAMMA(:,:,K) C AND WORK(:,:,K), C CORRESPONDING COLUMNS IN WORK(:,:,K), C COLUMNS IN GAMMA(:,:,K+1), IF K<=NBLOKS-1, C AND COLUMNS IN BOTBLK IF K=NBLOKS. C C PIVOT(INCRJ) = INCR+IPVT IF(IPVT.EQ.JMINN) GO TO 500 C C INTERCHANGE ROWS. C CALL DSWAP(NOVRLP-J+1,GAMMA(IPVT,J,K),NOVRLP, * GAMMA(JMINN,J,K),NOVRLP) CALL DSWAP(NOVRLP,WORK(IPVT,1,K),NOVRLP, * WORK(JMINN,1,K),NOVRLP) C C INTERCHANGE COLUMNS. C CALL DSWAP(NOVRLP-JMINN+1,WORK(JMINN,IPVT,K),1, * WORK(JMINN,JMINN,K),1) IF(K.LT.NBLOKS)THEN C C INTERCHANGE COLUMNS IN GAMMA(:,:,K+1) IF K < NBLOKS. C CALL DSWAP(NOVRLP,GAMMA(1,IPVT,KPLUS1),1, * GAMMA(1,JMINN,KPLUS1),1) * ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K=NBLOKS. C CALL DSWAP(NRWBOT,BOTBLK(1,IPVT),1,BOTBLK(1,JMINN),1) ENDIF 500 CONTINUE C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/GAMMA(JMINN,J,K) CALL DSCAL(NOVRLP-JMINN,ROWPIV,GAMMA(LOOP,J,K),1) C C PERFORM ROW ELIMINATION. C DO 520 L = JPLUS1,NOVRLP ROWMLT = GAMMA(JMINN,L,K) CALL DAXPY(NOVRLP-JMINN,-ROWMLT,GAMMA(LOOP,J,K),1, * GAMMA(LOOP,L,K),1) * 520 CONTINUE DO 525 L = 1,JMINN ROWMLT = WORK(JMINN,L,K) CALL DAXPY(NOVRLP-JMINN,-ROWMLT,GAMMA(LOOP,J,K),1, * WORK(LOOP,L,K),1) 525 CONTINUE 530 CONTINUE INCR = INCR + NRWBOT C C C IN ARRAY WORK(:,:,K) PERFORM NRWTOP COLUMN ELIMINATIONS WITH C COLUMN PIVOTING. C C DO 740 I = NRWBT1,NOVRLP NRWTPI = NRWTOP + I IPLUS1 = I+1 IMINN = I-NRWBOT C C DETERMINE COLUMN PIVOT WITH PIVOT INDEX. C INCRN = INCR + IMINN JPVT = IDAMAX(NRWBOT+1,WORK(I,IMINN,K),NOVRLP) + IMINN-1 COLMAX = ABS(WORK(I,JPVT,K)) JPVT = JPVT + NOVRLP C C CHECK FOR SINGULARITY. C PIVMAX = MAX(COLMAX,PIVMAX) IF(PIVMAX+COLMAX.EQ.PIVMAX) GO TO 920 C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(INCRN) = INCRN+JPVT-NRWTPI IF(JPVT.EQ.NRWTPI) GO TO 660 C C INTERCHANGE COLUMNS IN WORK(:,:,K). C JPVT2 = JPVT - NOVRLP CALL DSWAP(NOVRLP-I+1,WORK(I,JPVT2,K),1, * WORK(I,IMINN,K),1) IF(K.LT.NBLOKS)THEN C C INTERCHANGE COLUMNS IN GAMMA(:,:,K+1). C CALL DSWAP(NOVRLP,GAMMA(1,JPVT2,KPLUS1),1, * GAMMA(1,IMINN,KPLUS1),1) ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K=NBLOKS. C CALL DSWAP(NRWBOT,BOTBLK(1,JPVT2),1,BOTBLK(1,IMINN),1) ENDIF * 660 CONTINUE C C COMPUTE COLUMN MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C COLPIV = ONE/WORK(I,IMINN,K) CALL DSCAL(NRWBOT,COLPIV,WORK(I,IMINN+1,K),NOVRLP) DO 730 J=(IMINN+1),I COLMLT = WORK(I,J,K) CALL DAXPY(NOVRLP-I,-COLMLT,WORK(IPLUS1,IMINN,K),1, * WORK(IPLUS1,J,K),1) IF(K.LT.NBLOKS)THEN C C ADJUST COLUMNS IN GAMMA(:,:,K+1) C CALL DAXPY(NOVRLP,-COLMLT,GAMMA(1,IMINN,KPLUS1),1, * GAMMA(1,J,KPLUS1),1) ELSE C C ADJUST COLUMNS IN BOTBLK IF K = NBLOKS. C CALL DAXPY(NRWBOT,-COLMLT,BOTBLK(1,IMINN),1, * BOTBLK(1,J),1) ENDIF 730 CONTINUE 740 CONTINUE INCR = INCR + NRWTOP 745 CONTINUE C C C C PERFORM (NOVRLP-NRWTOP-1)=NRWBOT-1 ROW ELIMINATION(S) WITH ROW C PIVOTING ON BOTBLK. C C C IF(NRWBOT.EQ.1) GO TO 910 DO 900 J = NRWTP1,(NOVRLP-1) JPLUS1 = J+1 JMINN = J-NRWTOP IPVT = JMINN INCRJ = INCR + JMINN LOOP = JMINN+1 IPVT = IDAMAX(NRWBOT-JMINN+1,BOTBLK(JMINN,J),1) + JMINN - 1 ROWMAX = ABS(BOTBLK(IPVT,J)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C IF NECESSARY, INTERCHANGE ROWS IN BOTBLK. C PIVOT(INCRJ) = INCR + IPVT IF(IPVT.NE.JMINN) THEN CALL DSWAP(NOVRLP-J+1,BOTBLK(IPVT,J),NRWBOT, * BOTBLK(JMINN,J),NRWBOT) ENDIF C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/BOTBLK(JMINN,J) CALL DSCAL(NRWBOT-JMINN,ROWPIV,BOTBLK(LOOP,J),1) C C PERFORM ROW ELIMINATIONS. C DO 800 L = JPLUS1,NOVRLP ROWMLT = BOTBLK(JMINN,L) CALL DAXPY(NRWBOT-JMINN,-ROWMLT,BOTBLK(LOOP,J),1, * BOTBLK(LOOP,L),1) 800 CONTINUE 900 CONTINUE 910 CONTINUE C C DONE PROVIDED LAST ELEMENT IS NOT ZERO. C IF(PIVMAX+ABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN C C C MATRIX IS SINGULAR -SET IFLAG=-1 AND TERMINATE. C C 920 CONTINUE IFLAG = -1 RETURN C C C ERROR IN INPUT PARAMETERS. C C 930 CONTINUE IFLAG = 1 * RETURN * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABBDEC. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC END * SUBROUTINE ABBSOL(N, TOPBLK, NRWTOP, NOVRLP, * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, B, X) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT THE DECOMPOSED MATRIX A PRODUCED BY C ABBDEC, TOGETHER WITH THE CORRESPONDING RIGHT HAND SIDE, AND C RETURNS THE SOLUTION OF THE LINEAR SYSTEM. C C THIS IS THE DOUBLE PRECISION VERSION OF ABBSOL. THE SINGLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL DOUBLE PRECISION C DECLARATIONS BY REAL, AND BY REPLACING ALL REFERENCES TO C DOUBLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C THEIR SINGLE PRECISION EQUIVALENTS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C (NBLOKS + 1) * NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - DOUBLE PRECISION (NRWTOP, NOVRLP) C THE DECOMPOSED FIRST BLOCK. C UNCHANGED ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE ORDER C OF GAMMA(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY C UNCHANGED ON EXIT. C C GAMMA - DOUBLE PRECISION (NOVRLP, NOVRLP, NBLOKS) C GAMMA(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP DECOMPOSED BLOCK OF GAMMA. C UNCHANGED ON EXIT. C C WORK - DOUBLE PRECISION (NOVRLP, MRWBLK, NBLOKS). C WORKSPACE USED TO STORE THE DECOMPOSITION C OF WHAT WERE ORIGINALLY THE IDENTITY BLOCKS. C UNCHANGED ON EXIT. C C BOTBLK - DOUBLE PRECISION(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSITION ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C RECORDS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION BY ABBDEC. C UNCHANGED ON EXIT. C C B - DOUBLE PRECISION (N). C THE RIGHT HAND SIDE. C UNCHANGED ON EXIT. C C X - DOUBLE PRECISION (N). C WORK SPACE FOR THE SOLUTION VECTOR. C CONTAINS THE SOLUTION VECTOR ON EXIT. C C ROUTINES CALLED : C C NONE C C ABBSOL USES THE FOLLOWING ROUTINES FROM THE BASIC LINEAR C ALGEBRA LIBRARY (BLAS) : C C DAXPY, DDOT. C C TO CONVERT TO SINGLE PRECISION, THESE MUST BE REPLACED BY C CALLS TO : C C SAXPY, SDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NBLOKS,NRWTOP,NOVRLP,PIVOT(N) DOUBLE PRECISION TOPBLK,GAMMA,WORK,BOTBLK,B,X DIMENSION TOPBLK(NRWTOP,*), * GAMMA(NOVRLP,NOVRLP,*), * WORK(NOVRLP,NOVRLP,*), * BOTBLK(NOVRLP-NRWTOP,*),B(N),X(N) C C LOCAL VARIABLES : C INTEGER I,J,K,LOOP,INCR,INCRJ, * INCRI,INCRTP,JPIVOT,IPVTI,NRWBOT, * JMINN,NRWTPI,IPIVOT,NRWTP1,NRWBTJ,NRWBTI,NRWTP0,NRWTP2 * DOUBLE PRECISION SWAP,DDOT C C DEFINE THE CONSTANTS USED THROUGHOUT. C NRWTP0 = NRWTOP - 1 NRWTP1 = NRWTOP+1 NRWTP2 = NRWTOP+2 NRWBOT = NOVRLP - NRWTOP DO 100 J = 1,N X(J) = B(J) 100 CONTINUE C C C C FORWARD RECURSION ... C C C C FIRST IN TOPBLK: C C FORWARD SOLUTION. C C C DO 130 J = 1,NRWTP0 X(J) = X(J)/TOPBLK(J,J) LOOP = J+1 CALL DAXPY(NRWTOP-J,-X(J),TOPBLK(LOOP,J),1,X(LOOP),1) 130 CONTINUE X(NRWTOP) = X(NRWTOP)/TOPBLK(NRWTOP,NRWTOP) INCR = NRWTOP C C C IN EACH BLOCK K, K = 1,NBLOKS... C C C DO 360 K = 1,NBLOKS C C FORWARD MODIFICATION C INCRTP = INCR DO 240 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL DAXPY(NOVRLP,-X(INCRJ),GAMMA(1,J,K),1,X(INCRTP+1),1) 240 CONTINUE C C FORWARD ELIMINATION C IF(NRWTP1.LE.NOVRLP)THEN DO 260 J = NRWTP1,NOVRLP INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.NE.INCRJ) THEN SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP ENDIF LOOP = J-NRWTP0 CALL DAXPY(NOVRLP-LOOP+1,-X(INCRJ),GAMMA(LOOP,J,K),1, * X(INCRTP+LOOP),1) 260 CONTINUE ENDIF C C FORWARD SOLUTION C DO 280 J = 1,NRWTP0 NRWBTJ = NRWBOT+J INCRJ = INCR+NRWBTJ X(INCRJ) = X(INCRJ)/WORK(NRWBTJ,J,K) LOOP = NRWBTJ+1 CALL DAXPY(NOVRLP-NRWBTJ,-X(INCRJ),WORK(LOOP,J,K),1, * X(LOOP+INCRTP),1) 280 CONTINUE INCR = INCR+NOVRLP X(INCR) = X(INCR)/WORK(NOVRLP,NRWTOP,K) 360 CONTINUE C C IN BOTBLK... C C FORWARD MODIFICATION C INCRTP = INCR * DO 375 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL DAXPY(NRWBOT,-X(INCRJ),BOTBLK(1,J),1,X(INCRTP+1),1) 375 CONTINUE C C FORWARD MODIFICATION C * DO 390 J = NRWTP1,(NOVRLP-1) INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.NE.INCRJ) THEN SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP ENDIF LOOP = J - NRWTP0 CALL DAXPY(NRWBOT-LOOP+1,-X(INCRJ),BOTBLK(LOOP,J),1, * X(INCRTP+LOOP),1) 390 CONTINUE C C C C BACKWARD RECURSION C C C C FIRST IN BOTBLK... C C DO 430 J = NOVRLP,NRWTP1,-1 JMINN = J-NRWTOP INCRJ = INCR+JMINN X(INCRJ) = X(INCRJ)/BOTBLK(JMINN,J) LOOP = JMINN -1 CALL DAXPY(JMINN-1,-X(INCRJ),BOTBLK(1,J),1,X(INCR+1),1) 430 CONTINUE C C IN EACH BLOCK GOING IN REVERSE ORDER... C DO 520 K = NBLOKS,1,-1 INCR = INCR-NRWTOP C C BACKWARD ELIMINATION IN WORK(:,:,K). C DO 470 I = NRWTOP,1,-1 NRWBTI = NRWBOT+I LOOP = I+1 INCRI = INCR+I X(INCRI) = X(INCRI)-DDOT(NRWBOT,WORK(NRWBTI,LOOP,K),NOVRLP, * X(INCR+LOOP),1) IPIVOT = PIVOT(INCRI) IF(INCRI.NE.IPIVOT) THEN SWAP = X(INCRI) X(INCRI) = X(IPIVOT) X(IPIVOT) = SWAP ENDIF 470 CONTINUE INCR = INCR-NRWBOT C C BACKWARD MODIFICATION, SOLUTION AND ADJUSTMENT (ORDERING) OF C THE X'S. C DO 495 I = NRWBOT,1,-1 NRWTPI = NRWTOP+I LOOP = NRWTPI + 1 INCRI = INCR + I C C BACKWARD MODIFICATION C X(INCRI) = X(INCRI) - * DDOT(NOVRLP-NRWTPI,GAMMA(I,LOOP,K),NOVRLP,X(INCRI+1),1)- * DDOT(I-1,WORK(I,1,K),NOVRLP,X(INCRI+NOVRLP-NRWTPI+1),1) X(INCRI) = X(INCRI)-X(INCRI+NRWBOT) C C BACKWARD SOLUTION C X(INCRI) = X(INCRI)/GAMMA(I,NRWTPI,K) C C COLUMN PIVOTING ADJUSTMENT ON THE X'S. C IPIVOT = PIVOT(INCRI) IF(INCRI.NE.IPIVOT) THEN INCRI = INCRI+NRWBOT IPIVOT = IPIVOT+NRWBOT SWAP = X(INCRI) X(INCRI) = X(IPIVOT) X(IPIVOT) = SWAP ENDIF 495 CONTINUE 520 CONTINUE C C BACKWARD ELIMINATION IN TOPBLK C DO 510 I = NRWTOP,1,-1 LOOP = I + 1 X(I) = X(I) - DDOT(NOVRLP-I,TOPBLK(I,LOOP),NRWTOP,X(LOOP),1) IPVTI = PIVOT(I) IF(I.NE.IPVTI) THEN SWAP = X(I) X(I) = X(IPVTI) X(IPVTI) = SWAP ENDIF 510 CONTINUE RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABBSOL. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END C C####################################################################### C SUBROUTINE ABBLEQ(N, TOPBLK, NRWTOP, NOVRLP, NBLOKS, GAMMA, * WORK, BOTBLK, PIVOT, B, X, IFLAG) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A (WITH THE ALMOST BLOCK C BI-DIAGONAL STRUCTURE ARISING WHEN CONDENSATION IS APPLIED TO THE C MATRIX OBTAINED WHEN SPLINE COLLOCATION AT GAUSSIAN POINTS WITH C MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT BOUNDARY C VALUE PROBLEMS WITH SEPARATED BOUNDARY CONDITIONS), TOGETHER WITH C THE CORRESPONDING RIGHT HAND SIDE. C C THIS PARTICULAR MATRIX STRUCTURE ALSO ARISES IN CERTAIN C IMPLICIT RUNGE-KUTTA APPLICATIONS, AND IN MULTIPLE SHOOTING. C C ABBLEQ CALLS A DECOMPOSITION ROUTINE, ABBDEC, FOLLOWED BY THE C CORRESPONDING SOLVE ROUTINE, ABBSOL, TO OBTAIN THE SOLUTION C OF THE SYSTEM. THE ALGORITHM IMPLEMENTS AN ALTERNATE COLUMN C AND ROW PIVOTING TECHNIQUE. C C THIS IS THE SINGLE PRECISION VERSION OF ABBLEQ. THE DOUBLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL REAL C DECLARATIONS BY DOUBLE PRECISION, AND BY REPLACING ALL REFERENCES C TO SINGLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C THEIR DOUBLE PRECISION EQUIVALENTS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C (NBLOKS + 1) * NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - REAL (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSITION ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE ORDER C OF GAMMA(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN GAMMA. C UNCHANGED ON EXIT. C C GAMMA - REAL (NOVRLP, NOVRLP, NBLOKS) C GAMMA(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP BLOCK OF GAMMA. C CONTAINS THE DECOMPOSITION ON EXIT. C C WORK - REAL (NOVRLP, NOVRLP, NBLOKS). C WORKSPACE USED TO STORE THE IDENTITY BLOCKS. C NOT INITIALIZED BY THE USER. C CONTAINS DECOMPOSITION ON EXIT. C C BOTBLK - REAL(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSITION ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C RECORDS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION ON EXIT. C C B - REAL (N). C THE RIGHT HAND SIDE. C UNCHANGED ON EXIT. C C X - REAL (N). C WORK SPACE FOR THE SOLUTION VECTOR. C IF IFLAG = 0, CONTAINS THE SOLUTION VECTOR ON EXIT. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C ROUTINES CALLED : C C ABBDEC, ABBSOL. C C CALLING SEQUENCES : C C CALL ABBDEC(N, TOPBLK, NRWTOP, NOVRLP, C * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, IFLAG) C IF ( IFLAG .NE. 0 ) RETURN C CALL ABBSOL(N, TOPBLK, NRWTOP, NOVRLP, C * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, B, X) C C IN ADDITION, THE ROUTINES ABBDEC AND ABBSOL CALL THE FOLLOWING C ROUTINES FROM THE BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) : C C ISAMAX, SSWAP, SSCAL, SAXPY, SDOT. C C TO CONVERT TO DOUBLE PRECISION, THESE MUST BE REPLACED BY C CALLS TO : C C ISAMAX, SSWAP, SSCAL, SAXPY, SDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NRWTOP,NOVRLP,NBLOKS,PIVOT(N),IFLAG * REAL TOPBLK,GAMMA,WORK,BOTBLK,B,X DIMENSION TOPBLK(NRWTOP,*),GAMMA(NOVRLP,NOVRLP,*), * WORK(NOVRLP,NOVRLP,*),BOTBLK(NOVRLP-NRWTOP,*), * B(N),X(N) CALL ABBDEC(N,TOPBLK,NRWTOP,NOVRLP, * NBLOKS,GAMMA,WORK,BOTBLK,PIVOT,IFLAG) C C IF SINGULARITY WAS NOT DETECTED, CALL THE SOLVE ROUTINE ABBSOL. C IF ( IFLAG .NE. 0 ) RETURN CALL ABBSOL(N,TOPBLK,NRWTOP,NOVRLP, * NBLOKS,GAMMA,WORK,BOTBLK,PIVOT,B,X) RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABBLEQ. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC END * SUBROUTINE ABBDEC(N, TOPBLK, NRWTOP, NOVRLP, * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, IFLAG) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A (WITH THE ALMOST BLOCK C BI-DIAGONAL STRUCTURE ARISING WHEN CONDENSATION IS APPLIED TO THE C MATRIX OBTAINED WHEN SPLINE COLLOCATION AT GAUSSIAN POINTS WITH C MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT BOUNDARY C VALUE PROBLEMS WITH SEPARATED BOUNDARY CONDITIONS). C C THIS PARTICULAR MATRIX STRUCTURE ALSO ARISES IN CERTAIN C IMPLICIT RUNGE-KUTTA APPLICATIONS, AND IN MULTIPLE SHOOTING. C C ABBDEC PERFORMS THE DECOMPOSITION, USING AN ALTERNATE COLUMN C AND ROW PIVOTING TECHNIQUE. C C C THIS IS THE SINGLE PRECISION VERSION OF ABBDEC. THE DOUBLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL REAL C DECLARATIONS BY DOUBLE PRECISION, AND BY REPLACING ALL REFERENCES C TO SINGLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C THEIR DOUBLE PRECISION EQUIVALENTS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C (NBLOKS + 1) * NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - REAL (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSITION ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE ORDER C OF GAMMA(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN GAMMA. C UNCHANGED ON EXIT. C C GAMMA - REAL (NOVRLP, NOVRLP, NBLOKS) C GAMMA(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP BLOCK OF GAMMA. C CONTAINS THE DECOMPOSITION ON EXIT. C C WORK - REAL (NOVRLP, MRWBLK, NBLOKS). C WORKSPACE USED TO STORE THE IDENTITY BLOCKS. C NOT INITIALIZED BY THE USER. C CONTAINS THE DECOMPOSITION ON EXIT. C C BOTBLK - REAL(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSITION ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C RECORDS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION ON EXIT. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C ROUTINES CALLED : C C NONE C C ABBDEC USES THE FOLLOWING ROUTINES FROM THE BASIC LINEAR C ALGEBRA LIBRARY (BLAS) : C C ISAMAX, SSWAP, SSCAL, SAXPY. C C TO CONVERT TO DOUBLE PRECISION, THESE MUST BE REPLACED BY C CALLS TO : C C ISAMAX, SSWAP, SSCAL, SAXPY. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NRWTOP,NOVRLP,NBLOKS,PIVOT(N),IFLAG * REAL TOPBLK,GAMMA,WORK,BOTBLK DIMENSION TOPBLK(NRWTOP,*),GAMMA(NOVRLP,NOVRLP,*), * WORK(NOVRLP,NOVRLP,*),BOTBLK(NOVRLP-NRWTOP,*) C C LOCAL VARIABLES : C INTEGER I,J,K,IPLUS1,JPLUS1,L,INCR,KPLUS1, * INCRJ,INCRN,IPVT,JPVT,IMINN,JMINN,LOOP, * JPVT2,NRWTP1,NRWTPI,NRWBOT,NRWBT1,ISAMAX * REAL COLMAX,ROWMAX,PIVMAX, * COLPIV,ROWPIV,COLMLT,ROWMLT,ZERO,ONE DATA ZERO /0.0/, ONE /1.0/ C C DEFINE THE CONSTANTS USED THROUGHOUT. C NRWTP1 = NRWTOP+1 NRWBOT = NOVRLP-NRWTOP NRWBT1 = NRWBOT + 1 IFLAG = 0 PIVMAX = ZERO C C CHECK THE VALIDITY OF THE INPUT PARAMETERS. C IF(N.NE.(NBLOKS+1)*NOVRLP) GO TO 930 IF(NRWTOP.EQ.0) GO TO 930 IF(NRWTOP.EQ.NOVRLP) GO TO 930 DO 30 K=1,NBLOKS DO 20 J = 1,NOVRLP DO 10 I = 1,NOVRLP WORK(I,J,K) = ZERO 10 CONTINUE WORK(J,J,K) = ONE 20 CONTINUE 30 CONTINUE C C C FIRST, IN TOPBLK C C APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN PIVOTING. C C DO 195 I = 1,NRWTOP IPLUS1 = I+1 C C DETERMINE COLUMN PIVOT AND PIVOT INDEX. C IPVT = ISAMAX(NOVRLP-I+1,TOPBLK(I,I),NRWTOP) + I - 1 COLMAX = ABS(TOPBLK(I,IPVT)) C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX .EQ. PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(I) = IPVT IF(IPVT.NE.I) THEN C C INTERCHANGE COLUMNS IN TOPBLK. C CALL SSWAP(NRWTOP-I+1,TOPBLK(I,I),1,TOPBLK(I,IPVT),1) C C INTERCHANGE COLUMNS IN GAMMA(:,:,1). C CALL SSWAP(NOVRLP,GAMMA(1,I,1),1,GAMMA(1,IPVT,1),1) ENDIF C C COMPUTE MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C COLPIV = ONE/TOPBLK(I,I) CALL SSCAL(NOVRLP-I,COLPIV,TOPBLK(I,IPLUS1),NRWTOP) DO 190 J = IPLUS1,NOVRLP COLMLT = TOPBLK(I,J) CALL SAXPY(NRWTOP-I,-COLMLT,TOPBLK(IPLUS1,I),1, * TOPBLK(IPLUS1,J),1) C C ADJUST COLUMNS IN GAMMA(:,:,1). C CALL SAXPY(NOVRLP,-COLMLT,GAMMA(1,I,1),1,GAMMA(1,J,1),1) 190 CONTINUE 195 CONTINUE INCR = NRWTOP C C C C IN EACH BLOCK... C C IF(K.LT.NBLOKS) C APPLY NNRWBOT ROW ELIMINATIONS WITH ROW PIVOTING ON GAMMA(:,:,K) C NRWTOP COLUMN ELIMINATIONS WITH COLUMN PIVOTING ON WORK(:,:,K) C DO 745 K = 1,NBLOKS KPLUS1 = K+1 C C C PERFORM NRWBOT ROW ELIMINATION WITH ROW PIVOTING C ON GAMMA(:,:,K). C C DO 530 J = NRWTP1,NOVRLP JPLUS1 = J+1 JMINN = J-NRWTOP IPVT = JMINN INCRJ = INCR + JMINN LOOP = JMINN + 1 IPVT = ISAMAX(NOVRLP-JMINN+1,GAMMA(JMINN,J,K),1)+JMINN-1 ROWMAX = ABS(GAMMA(IPVT,J,K)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C C IF NECESSARY, INTERCHANGE ROWS IN GAMMA(:,:,K) C AND WORK(:,:,K), C CORRESPONDING COLUMNS IN WORK(:,:,K), C COLUMNS IN GAMMA(:,:,K+1), IF K<=NBLOKS-1, C AND COLUMNS IN BOTBLK IF K=NBLOKS. C C PIVOT(INCRJ) = INCR+IPVT IF(IPVT.EQ.JMINN) GO TO 500 C C INTERCHANGE ROWS. C CALL SSWAP(NOVRLP-J+1,GAMMA(IPVT,J,K),NOVRLP, * GAMMA(JMINN,J,K),NOVRLP) CALL SSWAP(NOVRLP,WORK(IPVT,1,K),NOVRLP, * WORK(JMINN,1,K),NOVRLP) C C INTERCHANGE COLUMNS. C CALL SSWAP(NOVRLP-JMINN+1,WORK(JMINN,IPVT,K),1, * WORK(JMINN,JMINN,K),1) IF(K.LT.NBLOKS)THEN C C INTERCHANGE COLUMNS IN GAMMA(:,:,K+1) IF K < NBLOKS. C CALL SSWAP(NOVRLP,GAMMA(1,IPVT,KPLUS1),1, * GAMMA(1,JMINN,KPLUS1),1) * ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K=NBLOKS. C CALL SSWAP(NRWBOT,BOTBLK(1,IPVT),1,BOTBLK(1,JMINN),1) ENDIF 500 CONTINUE C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/GAMMA(JMINN,J,K) CALL SSCAL(NOVRLP-JMINN,ROWPIV,GAMMA(LOOP,J,K),1) C C PERFORM ROW ELIMINATION. C DO 520 L = JPLUS1,NOVRLP ROWMLT = GAMMA(JMINN,L,K) CALL SAXPY(NOVRLP-JMINN,-ROWMLT,GAMMA(LOOP,J,K),1, * GAMMA(LOOP,L,K),1) * 520 CONTINUE DO 525 L = 1,JMINN ROWMLT = WORK(JMINN,L,K) CALL SAXPY(NOVRLP-JMINN,-ROWMLT,GAMMA(LOOP,J,K),1, * WORK(LOOP,L,K),1) 525 CONTINUE 530 CONTINUE INCR = INCR + NRWBOT C C C IN ARRAY WORK(:,:,K) PERFORM NRWTOP COLUMN ELIMINATIONS WITH C COLUMN PIVOTING. C C DO 740 I = NRWBT1,NOVRLP NRWTPI = NRWTOP + I IPLUS1 = I+1 IMINN = I-NRWBOT C C DETERMINE COLUMN PIVOT WITH PIVOT INDEX. C INCRN = INCR + IMINN JPVT = ISAMAX(NRWBOT+1,WORK(I,IMINN,K),NOVRLP) + IMINN-1 COLMAX = ABS(WORK(I,JPVT,K)) JPVT = JPVT + NOVRLP C C CHECK FOR SINGULARITY. C PIVMAX = MAX(COLMAX,PIVMAX) IF(PIVMAX+COLMAX.EQ.PIVMAX) GO TO 920 C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(INCRN) = INCRN+JPVT-NRWTPI IF(JPVT.EQ.NRWTPI) GO TO 660 C C INTERCHANGE COLUMNS IN WORK(:,:,K). C JPVT2 = JPVT - NOVRLP CALL SSWAP(NOVRLP-I+1,WORK(I,JPVT2,K),1, * WORK(I,IMINN,K),1) IF(K.LT.NBLOKS)THEN C C INTERCHANGE COLUMNS IN GAMMA(:,:,K+1). C CALL SSWAP(NOVRLP,GAMMA(1,JPVT2,KPLUS1),1, * GAMMA(1,IMINN,KPLUS1),1) ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K=NBLOKS. C CALL SSWAP(NRWBOT,BOTBLK(1,JPVT2),1,BOTBLK(1,IMINN),1) ENDIF * 660 CONTINUE C C COMPUTE COLUMN MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C COLPIV = ONE/WORK(I,IMINN,K) CALL SSCAL(NRWBOT,COLPIV,WORK(I,IMINN+1,K),NOVRLP) DO 730 J=(IMINN+1),I COLMLT = WORK(I,J,K) CALL SAXPY(NOVRLP-I,-COLMLT,WORK(IPLUS1,IMINN,K),1, * WORK(IPLUS1,J,K),1) IF(K.LT.NBLOKS)THEN C C ADJUST COLUMNS IN GAMMA(:,:,K+1) C CALL SAXPY(NOVRLP,-COLMLT,GAMMA(1,IMINN,KPLUS1),1, * GAMMA(1,J,KPLUS1),1) ELSE C C ADJUST COLUMNS IN BOTBLK IF K = NBLOKS. C CALL SAXPY(NRWBOT,-COLMLT,BOTBLK(1,IMINN),1, * BOTBLK(1,J),1) ENDIF 730 CONTINUE 740 CONTINUE INCR = INCR + NRWTOP 745 CONTINUE C C C C PERFORM (NOVRLP-NRWTOP-1)=NRWBOT-1 ROW ELIMINATION(S) WITH ROW C PIVOTING ON BOTBLK. C C C IF(NRWBOT.EQ.1) GO TO 910 DO 900 J = NRWTP1,(NOVRLP-1) JPLUS1 = J+1 JMINN = J-NRWTOP IPVT = JMINN INCRJ = INCR + JMINN LOOP = JMINN+1 IPVT = ISAMAX(NRWBOT-JMINN+1,BOTBLK(JMINN,J),1) + JMINN - 1 ROWMAX = ABS(BOTBLK(IPVT,J)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C IF NECESSARY, INTERCHANGE ROWS IN BOTBLK. C PIVOT(INCRJ) = INCR + IPVT IF(IPVT.NE.JMINN) THEN CALL SSWAP(NOVRLP-J+1,BOTBLK(IPVT,J),NRWBOT, * BOTBLK(JMINN,J),NRWBOT) ENDIF C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/BOTBLK(JMINN,J) CALL SSCAL(NRWBOT-JMINN,ROWPIV,BOTBLK(LOOP,J),1) C C PERFORM ROW ELIMINATIONS. C DO 800 L = JPLUS1,NOVRLP ROWMLT = BOTBLK(JMINN,L) CALL SAXPY(NRWBOT-JMINN,-ROWMLT,BOTBLK(LOOP,J),1, * BOTBLK(LOOP,L),1) 800 CONTINUE 900 CONTINUE 910 CONTINUE C C DONE PROVIDED LAST ELEMENT IS NOT ZERO. C IF(PIVMAX+ABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN C C C MATRIX IS SINGULAR -SET IFLAG=-1 AND TERMINATE. C C 920 CONTINUE IFLAG = -1 RETURN C C C ERROR IN INPUT PARAMETERS. C C 930 CONTINUE IFLAG = 1 * RETURN * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABBDEC. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC END * SUBROUTINE ABBSOL(N, TOPBLK, NRWTOP, NOVRLP, * NBLOKS, GAMMA, WORK, BOTBLK, PIVOT, B, X) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT THE DECOMPOSED MATRIX A PRODUCED BY C ABBDEC, TOGETHER WITH THE CORRESPONDING RIGHT HAND SIDE, AND C RETURNS THE SOLUTION OF THE LINEAR SYSTEM. C C THIS IS THE SINGLE PRECISION VERSION OF ABBSOL. THE DOUBLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL REAL C DECLARATIONS BY SINGLE PRECISION, AND BY REPLACING ALL REFERENCES C TO BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C THEIR DOUBLE PRECISION EQUIVALENTS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C (NBLOKS + 1) * NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - REAL (NRWTOP, NOVRLP) C THE DECOMPOSED FIRST BLOCK. C UNCHANGED ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE ORDER C OF GAMMA(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY C UNCHANGED ON EXIT. C C GAMMA - REAL (NOVRLP, NOVRLP, NBLOKS) C GAMMA(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP DECOMPOSED BLOCK OF GAMMA. C UNCHANGED ON EXIT. C C WORK - REAL (NOVRLP, MRWBLK, NBLOKS). C WORKSPACE USED TO STORE THE DECOMPOSITION C OF WHAT WERE ORIGINALLY THE IDENTITY BLOCKS. C UNCHANGED ON EXIT. C C BOTBLK - REAL(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSITION ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C RECORDS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION BY ABBDEC. C UNCHANGED ON EXIT. C C B - REAL (N). C THE RIGHT HAND SIDE. C UNCHANGED ON EXIT. C C X - REAL (N). C WORK SPACE FOR THE SOLUTION VECTOR. C CONTAINS THE SOLUTION VECTOR ON EXIT. C C ROUTINES CALLED : C C NONE C C ABBSOL USES THE FOLLOWING ROUTINES FROM THE BASIC LINEAR C ALGEBRA LIBRARY (BLAS) : C C SAXPY, SDOT. C C TO CONVERT TO DOUBLE PRECISION, THESE MUST BE REPLACED BY C CALLS TO : C C SAXPY, SDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NBLOKS,NRWTOP,NOVRLP,PIVOT(N) REAL TOPBLK,GAMMA,WORK,BOTBLK,B,X DIMENSION TOPBLK(NRWTOP,*), * GAMMA(NOVRLP,NOVRLP,*), * WORK(NOVRLP,NOVRLP,*), * BOTBLK(NOVRLP-NRWTOP,*),B(N),X(N) C C LOCAL VARIABLES : C INTEGER I,J,K,LOOP,INCR,INCRJ, * INCRI,INCRTP,JPIVOT,IPVTI,NRWBOT, * JMINN,NRWTPI,IPIVOT,NRWTP1,NRWBTJ,NRWBTI,NRWTP0,NRWTP2 * REAL SWAP,SDOT C C DEFINE THE CONSTANTS USED THROUGHOUT. C NRWTP0 = NRWTOP - 1 NRWTP1 = NRWTOP+1 NRWTP2 = NRWTOP+2 NRWBOT = NOVRLP - NRWTOP DO 100 J = 1,N X(J) = B(J) 100 CONTINUE C C C C FORWARD RECURSION ... C C C C FIRST IN TOPBLK: C C FORWARD SOLUTION. C C C DO 130 J = 1,NRWTP0 X(J) = X(J)/TOPBLK(J,J) LOOP = J+1 CALL SAXPY(NRWTOP-J,-X(J),TOPBLK(LOOP,J),1,X(LOOP),1) 130 CONTINUE X(NRWTOP) = X(NRWTOP)/TOPBLK(NRWTOP,NRWTOP) INCR = NRWTOP C C C IN EACH BLOCK K, K = 1,NBLOKS... C C C DO 360 K = 1,NBLOKS C C FORWARD MODIFICATION C INCRTP = INCR DO 240 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL SAXPY(NOVRLP,-X(INCRJ),GAMMA(1,J,K),1,X(INCRTP+1),1) 240 CONTINUE C C FORWARD ELIMINATION C IF(NRWTP1.LE.NOVRLP)THEN DO 260 J = NRWTP1,NOVRLP INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.NE.INCRJ) THEN SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP ENDIF LOOP = J-NRWTP0 CALL SAXPY(NOVRLP-LOOP+1,-X(INCRJ),GAMMA(LOOP,J,K),1, * X(INCRTP+LOOP),1) 260 CONTINUE ENDIF C C FORWARD SOLUTION C DO 280 J = 1,NRWTP0 NRWBTJ = NRWBOT+J INCRJ = INCR+NRWBTJ X(INCRJ) = X(INCRJ)/WORK(NRWBTJ,J,K) LOOP = NRWBTJ+1 CALL SAXPY(NOVRLP-NRWBTJ,-X(INCRJ),WORK(LOOP,J,K),1, * X(LOOP+INCRTP),1) 280 CONTINUE INCR = INCR+NOVRLP X(INCR) = X(INCR)/WORK(NOVRLP,NRWTOP,K) 360 CONTINUE C C IN BOTBLK... C C FORWARD MODIFICATION C INCRTP = INCR * DO 375 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL SAXPY(NRWBOT,-X(INCRJ),BOTBLK(1,J),1,X(INCRTP+1),1) 375 CONTINUE C C FORWARD MODIFICATION C * DO 390 J = NRWTP1,(NOVRLP-1) INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.NE.INCRJ) THEN SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP ENDIF LOOP = J - NRWTP0 CALL SAXPY(NRWBOT-LOOP+1,-X(INCRJ),BOTBLK(LOOP,J),1, * X(INCRTP+LOOP),1) 390 CONTINUE C C C C BACKWARD RECURSION C C C C FIRST IN BOTBLK... C C DO 430 J = NOVRLP,NRWTP1,-1 JMINN = J-NRWTOP INCRJ = INCR+JMINN X(INCRJ) = X(INCRJ)/BOTBLK(JMINN,J) LOOP = JMINN -1 CALL SAXPY(JMINN-1,-X(INCRJ),BOTBLK(1,J),1,X(INCR+1),1) 430 CONTINUE C C IN EACH BLOCK GOING IN REVERSE ORDER... C DO 520 K = NBLOKS,1,-1 INCR = INCR-NRWTOP C C BACKWARD ELIMINATION IN WORK(:,:,K). C DO 470 I = NRWTOP,1,-1 NRWBTI = NRWBOT+I LOOP = I+1 INCRI = INCR+I X(INCRI) = X(INCRI)-SDOT(NRWBOT,WORK(NRWBTI,LOOP,K),NOVRLP, * X(INCR+LOOP),1) IPIVOT = PIVOT(INCRI) IF(INCRI.NE.IPIVOT) THEN SWAP = X(INCRI) X(INCRI) = X(IPIVOT) X(IPIVOT) = SWAP ENDIF 470 CONTINUE INCR = INCR-NRWBOT C C BACKWARD MODIFICATION, SOLUTION AND ADJUSTMENT (ORDERING) OF C THE X'S. C DO 495 I = NRWBOT,1,-1 NRWTPI = NRWTOP+I LOOP = NRWTPI + 1 INCRI = INCR + I C C BACKWARD MODIFICATION C X(INCRI) = X(INCRI) - * SDOT(NOVRLP-NRWTPI,GAMMA(I,LOOP,K),NOVRLP,X(INCRI+1),1)- * SDOT(I-1,WORK(I,1,K),NOVRLP,X(INCRI+NOVRLP-NRWTPI+1),1) X(INCRI) = X(INCRI)-X(INCRI+NRWBOT) C C BACKWARD SOLUTION C X(INCRI) = X(INCRI)/GAMMA(I,NRWTPI,K) C C COLUMN PIVOTING ADJUSTMENT ON THE X'S. C IPIVOT = PIVOT(INCRI) IF(INCRI.NE.IPIVOT) THEN INCRI = INCRI+NRWBOT IPIVOT = IPIVOT+NRWBOT SWAP = X(INCRI) X(INCRI) = X(IPIVOT) X(IPIVOT) = SWAP ENDIF 495 CONTINUE 520 CONTINUE C C BACKWARD ELIMINATION IN TOPBLK C DO 510 I = NRWTOP,1,-1 LOOP = I + 1 X(I) = X(I) - SDOT(NOVRLP-I,TOPBLK(I,LOOP),NRWTOP,X(LOOP),1) IPVTI = PIVOT(I) IF(I.NE.IPVTI) THEN SWAP = X(I) X(I) = X(IPVTI) X(IPVTI) = SWAP ENDIF 510 CONTINUE RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABBSOL. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END C C####################################################################### C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS DRIVING PROGRAM RUNS THE DOUBLE PRECISION ABDPACK PACKAGE C ON A SAMPLE PROBLEM. C C THE SAMPLE SYSTEM GIVEN IN THE DATA STATEMENTS CONSISTS OF C A SYSTEM OF 18 EQUATIONS. THERE ARE TWO BLOCKS C (CORRESPONDING TO TWO SUBINTERVALS IN THE COLLOCATION). C C THUS, N (ORDER) = 18; NBLOKS (NO. OF BLOCKS) = 2; C NOVRLP = 4; NRWBLK (ROWS IN ARRAY1) = 3; C NRWTOP= NRWBOT = 2. C C THE COEFFICIENT MATRIX OF THE SYSTEM IS (WITH ONLY THE C ELEMENTS PASSED TO ABDLEQ SHOWN): C C 1 1 2 3 C 1 2 6 3 TOPBLK(1:2,1:4) C----------------------------------------------------------------------- C 1 2 2 1 3 4 20 C 2 3 4 1 5 50 1 C 3 4 50 1 1 2 3 ARRAY1(1:3,1:7,1) C----------------------------------------------------------------------- C 1 2 3 4 5 40 1 1 0 0 0 C 2 3 4 5 1 2 1 0 1 0 0 C 3 4 5 1 2 3 20 0 0 1 0 C 4 5 1 2 3 4 5 0 0 0 1 ARRAY2(1:4,1:11,1) C----------------------------------------------------------------------- C 1 2 2 1 3 4 20 C 2 3 4 1 5 50 1 C ARRAY1(1:3,1:7,2) 3 4 50 1 1 2 3 C----------------------------------------------------------------------- C 1 2 3 4 5 40 1 1 0 0 0 C 2 3 4 5 1 2 30 0 1 0 0 C 3 4 5 1 2 3 20 0 0 1 0 C ARRAY(1:4,1:11,2) 4 5 1 2 3 4 5 0 0 0 1 C----------------------------------------------------------------------- C 1 2 1 1 C BOTBLK(1:2,1:4) 2 3 40 10 C----------------------------------------------------------------------- C C THE RIGHT HAND SIDE HAS BEEN CHOSEN TO MAKE THE EXACT SOLUTION C EQUAL TO (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1), AND IS: C C B(1:18) = 7,12,33,66,64,57,19,39,25,33,66,64,57,48,39,25,5,55 . C----------------------------------------------------------------------- C C THE PROGRAM ABDLEQ FIRST CALLS THE DECOMPOSITION PROCEDURE C ABDDEC AND THEN THE SOLVE PROCEDURE ABDSOL PROVIDED IFLAG C IS ZERO. THE CODE FOR THE SEPARATE CALLS IS INCLUDED AS C COMMENTS AFTER THE CALL TO ABDLEQ. C C THE OUTPUT FROM THIS DRIVER ON A SUN 4/280 UNDER UNIX, USING F77 C IN DOUBLE PRECISION WAS : C C IFLAG = 0 C PIVOT VECTOR (FIRST 17 COMPONENTS): * C 4,3,7,6,7,9,8,10,11,14,13,13,14,16,17,18,18. C C SOLUTION : C C 0.1000000000000002D+01 0.1000000000000002D+01 C 0.9999999999999996D+00 0.9999999999999993D+00 C 0.1000000000000008D+01 0.9999999999999991D+00 C 0.9999999999999987D+00 0.9999999999999963D+00 C 0.9999999999999901D+00 0.1000000000000002D+01 C 0.9999999999999710D+00 0.9999999999998995D+00 C 0.1000000000000011D+01 0.1000000000000016D+01 C 0.1000000000000185D+01 0.9999999999997909D+00 C 0.9999999999999312D+00 0.1000000000000300D+01 C C IT SHOULD BE NOTED THAT PIVOTING WAS CARRIED OUT AT EVERY STEP IN C THIS EXAMPLE, AS INDICATED BY THE RESULTING PIVOT VECTOR. C C THE FOLLOWING SEQUENCE OF PIVOTING WAS CARRIED OUT: C C INTERCHANGE OF COLUMNS 1 AND 4, AND OF COLUMNS 2 AND 3 IN TOPBLK. C INTERCHANGE OF COLUMNS 3 AND 7, 4 AND 6, AND OF 5 AND 7 IN C ARRAY1(:,:,1). C INTERCHANGE OF ROWS 6 AND 9, AND OF ROWS 7 AND 8 IN C ARRAY2(:,:,1). C INTERCHANGE OF COLUMNS 8 AND 10, AND OF 9 AND 11 IN ARRAY2(:,:,1) C INTERCHANGE OF COLUMNS 10 AND 14, 11 AND 13, AND OF 12 AND 13 IN C ARRAY1(:,:,2). C INTERCHANGE OF ROWS 13 AND 14, AND OF ROWS 14 AND 16 IN C ARRAY2(:,:,2). C INTERCHANGE OF COLUMNS 15 AND 17, AND OF 16 AND 18 IN C ARRAY2(:,:,2). C INTERCHANGE OF ROWS 17 AND 18 IN BOTBLK. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC * DOUBLE PRECISION TOPBLK(2,4),ARRAY1(3,7,2),ARRAY2(4,11,2), * BOTBLK(2,4),B(18),X(18) INTEGER N,NTOP,NOVRLP,NRWBLK,NBLOKS,PIVOT(18),IFLAG INTEGER I,NOUT DATA NOUT /6 / DATA N,NTOP,NOVRLP,NRWBLK,NBLOKS /18,2,4,3,2 / DATA TOPBLK / * 1.0D0, 1.0D0, 1.0D0, 2.0D0, 2.0D0, 6.0D0, 3.0D0, 3.0D0 / C C THE ARRAYS ARRAY1, ARRAY2 ARE GIVEN IN THE FOLLOWING DATA C STATEMENTS, WITH EACH BLOCK BEING GIVEN BY COLUMNS: C DATA ARRAY1 / * 1.0D0, 2.0D0, 3.0D0, 2.0D0, 3.0D0, 4.0D0, 2.0D0, 4.0D0, * 50.0D0, 1.0D0, 1.0D0, 1.0D0, 3.0D0, 5.0D0, 1.0D0, 4.0D0, * 50.0D0, 2.0D0, 20.0D0, 1.0D0, 3.0D0, * 1.0D0, 2.0D0, 3.0D0, 2.0D0, 3.0D0, 4.0D0, 2.0D0, 4.0D0, * 50.0D0, 1.0D0, 1.0D0, 1.0D0, 3.0D0, 5.0D0, 1.0D0, 4.0D0, * 50.0D0, 2.0D0, 20.0D0, 1.0D0, 3.0D0 / DATA ARRAY2 / * 1.0D0, 2.0D0, 3.0D0, 4.0D0, 2.0D0, 3.0D0, 4.0D0, 5.0D0, * 3.0D0, 4.0D0, 5.0D0, 1.0D0, 4.0D0, 5.0D0, 1.0D0, 2.0D0, * 5.0D0, 1.0D0, 2.0D0, 3.0D0, 40.0D0, 2.0D0, 3.0D0, 4.0D0, * 1.0D0, 1.0D0, 20.0D0, 5.0D0, 1.0D0, 0.0D0, 0.0D0, 0.0D0, * 0.0D0, 1.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 1.0D0, 0.0D0, * 0.0D0, 0.0D0, 0.0D0, 1.0D0, * 1.0D0, 2.0D0, 3.0D0, 4.0D0, 2.0D0, 3.0D0, 4.0D0, 5.0D0, * 3.0D0, 4.0D0, 5.0D0, 1.0D0, 4.0D0, 5.0D0, 1.0D0, 2.0D0, * 5.0D0, 1.0D0, 2.0D0, 3.0D0, 40.0D0, 2.0D0, 3.0D0, 4.0D0, * 1.0D0, 30.0D0, 20.0D0, 5.0D0, 1.0D0, 0.0D0, 0.0D0, 0.0D0, * 0.0D0, 1.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 1.0D0, 0.0D0, * 0.0D0, 0.0D0, 0.0D0, 1.0D0 / * DATA BOTBLK / * 1.0D0, 2.0D0, 2.0D0, 3.0D0, 1.0D0, 40.0D0, 1.0D0, 10.0D0 / * C RIGHT HAND SIDE B : * DATA B / * 7.0D0, 12.0D0, 33.0D0, 66.0D0, 64.0D0, 57.0D0, 19.0D0, * 39.0D0, 25.0D0, 33.0D0, 66.0D0, 64.0D0, 57.0D0, 48.0D0, * 39.0D0, 25.0D0, 5.0D0, 55.0D0/ * CALL ABDLEQ(N,TOPBLK,NTOP,NOVRLP,ARRAY1,NRWBLK,NBLOKS,ARRAY2, * BOTBLK,PIVOT,B,X,IFLAG) * C THE SAME CALCULATION COULD BE DONE BY THE FOLLOWING SEPARATE C CALLS TO ABDDEC AND ABDSOL. C C FIRST, CALL THE DECOMPOSITION ROUTINE ABDDEC. C C CALL ABDDEC(N, TOPBLK, NTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) * C UNLESS A IS NUMERICALLY SINGULAR OR THE PARAMETERS ARE C IMPROPERLY DEFINED, CALL ABDSOL. * C IF ( IFLAG .EQ. 0 ) THEN * C CALL ABDSOL(N, TOPBLK, NTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, B,X) C ENDIF * WRITE(NOUT,1000)IFLAG IF ( IFLAG .EQ. 0 ) THEN WRITE(NOUT,2000) WRITE(NOUT,3000)(PIVOT(I),I=1,N-1) WRITE(NOUT,4000) WRITE(NOUT,5000)(X(I),I=1,N) ELSE IF ( IFLAG .EQ. 1 ) WRITE(NOUT,7000) IF ( IFLAG .EQ. -1) WRITE(NOUT,6000) ENDIF STOP 1000 FORMAT(9H IFLAG = ,I4) 2000 FORMAT(16H PIVOT VECTOR : ) 3000 FORMAT(15I4) 4000 FORMAT(12H SOLUTION : ) 5000 FORMAT(2D25.16) 6000 FORMAT(30H MATRIX APPEARS TO BE SINGULAR) 7000 FORMAT(20H IMPROPER PARAMETERS) END C C####################################################################### C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS DRIVING PROGRAM RUNS THE SINGLE PRECISION ABDPACK PACKAGE ON C A SAMPLE PROBLEM. C C THE SAMPLE SYSTEM GIVEN IN THE DATA STATEMENTS CONSISTS OF C A SYSTEM OF 18 EQUATIONS. THERE ARE TWO BLOCKS C (CORRESPONDING TO TWO SUBINTERVALS IN THE COLLOCATION). C C THUS, N (ORDER) = 18; NBLOKS (NO. OF BLOCKS) = 2; C NOVRLP = 4; NRWBLK (ROWS IN ARRAY1) = 3; C NRWTOP= NRWBOT = 2. C C THE COEFFICIENT MATRIX OF THE SYSTEM IS (WITH ONLY THE C ELEMENTS PASSED TO ABDLEQ SHOWN): C C 1 1 2 3 C 1 2 6 3 TOPBLK(1:2,1:4) C----------------------------------------------------------------------- C 1 2 2 1 3 4 20 C 2 3 4 1 5 50 1 C 3 4 50 1 1 2 3 ARRAY1(1:3,1:7,1) C----------------------------------------------------------------------- C 1 2 3 4 5 40 1 1 0 0 0 C 2 3 4 5 1 2 1 0 1 0 0 C 3 4 5 1 2 3 20 0 0 1 0 C 4 5 1 2 3 4 5 0 0 0 1 ARRAY2(1:4,1:11,1) C----------------------------------------------------------------------- C 1 2 2 1 3 4 20 C 2 3 4 1 5 50 1 C ARRAY1(1:3,1:7,2) 3 4 50 1 1 2 3 C----------------------------------------------------------------------- C 1 2 3 4 5 40 1 1 0 0 0 C 2 3 4 5 1 2 30 0 1 0 0 C 3 4 5 1 2 3 20 0 0 1 0 C ARRAY(1:4,1:11,2) 4 5 1 2 3 4 5 0 0 0 1 C----------------------------------------------------------------------- C 1 2 1 1 C BOTBLK(1:2,1:4) 2 3 40 10 C----------------------------------------------------------------------- C C THE RIGHT HAND SIDE HAS BEEN CHOSEN TO MAKE THE EXACT SOLUTION C EQUAL TO (1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1), AND IS: C C B(1:18) = 7,12,33,66,64,57,19,39,25,33,66,64,57,48,39,25,5,55 . C----------------------------------------------------------------------- C C THE PROGRAM ABDLEQ FIRST CALLS THE DECOMPOSITION PROCEDURE C ABDDEC AND THEN THE SOLVE PROCEDURE ABDSOL PROVIDED IFLAG C IS ZERO. THE CODE FOR THE SEPARATE CALLS IS INCLUDED AS C COMMENTS AFTER THE CALL TO ABDLEQ. C C THE OUTPUT FROM THIS DRIVER ON A SUN 4/280 UNDER UNIX, USING F77 C IN SINGLE PRECISION WAS : C C IFLAG = 0 C PIVOT VECTOR (FIRST 17 COMPONENTS): * C 4,3,7,6,7,9,8,10,11,14,13,13,14,16,17,18,18. C C SOLUTION : C C 0.99999660E+00 0.99999809E+00 C 0.10000005E+01 0.10000013E+01 C 0.99998510E+00 0.10000017E+01 C 0.10000021E+01 0.10000081E+01 C 0.10000142E+01 0.99999630E+00 C 0.10000473E+01 0.10001651E+01 C 0.99998212E+00 0.99997491E+00 C 0.99970317E+00 0.10003388E+01 C 0.10001128E+01 0.99950618E+00 C C IT SHOULD BE NOTED THAT PIVOTING WAS CARRIED OUT AT EVERY STEP IN C THIS EXAMPLE, AS INDICATED BY THE RESULTING PIVOT VECTOR. C C THE FOLLOWING SEQUENCE OF PIVOTING WAS CARRIED OUT: C C INTERCHANGE OF COLUMNS 1 AND 4, AND OF COLUMNS 2 AND 3 IN TOPBLK. C INTERCHANGE OF COLUMNS 3 AND 7, 4 AND 6, AND OF 5 AND 7 IN C ARRAY1(:,:,1). C INTERCHANGE OF ROWS 6 AND 9, AND OF ROWS 7 AND 8 IN C ARRAY2(:,:,1). C INTERCHANGE OF COLUMNS 8 AND 10, AND OF 9 AND 11 IN ARRAY2(:,:,1) C INTERCHANGE OF COLUMNS 10 AND 14, 11 AND 13, AND OF 12 AND 13 IN C ARRAY1(:,:,2). C INTERCHANGE OF ROWS 13 AND 14, AND OF ROWS 14 AND 16 IN C ARRAY2(:,:,2). C INTERCHANGE OF COLUMNS 15 AND 17, AND OF 16 AND 18 IN C ARRAY2(:,:,2). C INTERCHANGE OF ROWS 17 AND 18 IN BOTBLK. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC * REAL TOPBLK(2,4),ARRAY1(3,7,2),ARRAY2(4,11,2), * BOTBLK(2,4),B(18),X(18) INTEGER N,NTOP,NOVRLP,NRWBLK,NBLOKS,PIVOT(18),IFLAG INTEGER I,NOUT DATA NOUT /6 / DATA N,NTOP,NOVRLP,NRWBLK,NBLOKS /18,2,4,3,2 / DATA TOPBLK / * 1.0E0, 1.0E0, 1.0E0, 2.0E0, 2.0E0, 6.0E0, 3.0E0, 3.0E0 / C C THE ARRAYS ARRAY1, ARRAY2 ARE GIVEN IN THE FOLLOWING DATA C STATEMENTS, WITH EACH BLOCK BEING GIVEN BY COLUMNS: C DATA ARRAY1 / * 1.0E0, 2.0E0, 3.0E0, 2.0E0, 3.0E0, 4.0E0, 2.0E0, 4.0E0, * 50.0E0, 1.0E0, 1.0E0, 1.0E0, 3.0E0, 5.0E0, 1.0E0, 4.0E0, * 50.0E0, 2.0E0, 20.0E0, 1.0E0, 3.0E0, * 1.0E0, 2.0E0, 3.0E0, 2.0E0, 3.0E0, 4.0E0, 2.0E0, 4.0E0, * 50.0E0, 1.0E0, 1.0E0, 1.0E0, 3.0E0, 5.0E0, 1.0E0, 4.0E0, * 50.0E0, 2.0E0, 20.0E0, 1.0E0, 3.0E0 / DATA ARRAY2 / * 1.0E0, 2.0E0, 3.0E0, 4.0E0, 2.0E0, 3.0E0, 4.0E0, 5.0E0, * 3.0E0, 4.0E0, 5.0E0, 1.0E0, 4.0E0, 5.0E0, 1.0E0, 2.0E0, * 5.0E0, 1.0E0, 2.0E0, 3.0E0, 40.0E0, 2.0E0, 3.0E0, 4.0E0, * 1.0E0, 1.0E0, 20.0E0, 5.0E0, 1.0E0, 0.0E0, 0.0E0, 0.0E0, * 0.0E0, 1.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 1.0E0, 0.0E0, * 0.0E0, 0.0E0, 0.0E0, 1.0E0, * 1.0E0, 2.0E0, 3.0E0, 4.0E0, 2.0E0, 3.0E0, 4.0E0, 5.0E0, * 3.0E0, 4.0E0, 5.0E0, 1.0E0, 4.0E0, 5.0E0, 1.0E0, 2.0E0, * 5.0E0, 1.0E0, 2.0E0, 3.0E0, 40.0E0, 2.0E0, 3.0E0, 4.0E0, * 1.0E0, 30.0E0, 20.0E0, 5.0E0, 1.0E0, 0.0E0, 0.0E0, 0.0E0, * 0.0E0, 1.0E0, 0.0E0, 0.0E0, 0.0E0, 0.0E0, 1.0E0, 0.0E0, * 0.0E0, 0.0E0, 0.0E0, 1.0E0 / * DATA BOTBLK / * 1.0E0, 2.0E0, 2.0E0, 3.0E0, 1.0E0, 40.0E0, 1.0E0, 10.0E0 / * C RIGHT HAND SIDE B : * DATA B / * 7.0E0, 12.0E0, 33.0E0, 66.0E0, 64.0E0, 57.0E0, 19.0E0, * 39.0E0, 25.0E0, 33.0E0, 66.0E0, 64.0E0, 57.0E0, 48.0E0, * 39.0E0, 25.0E0, 5.0E0, 55.0E0/ * CALL ABDLEQ(N,TOPBLK,NTOP,NOVRLP,ARRAY1,NRWBLK,NBLOKS,ARRAY2, * BOTBLK,PIVOT,B,X,IFLAG) * C THE SAME CALCULATION COULD BE DONE BY THE FOLLOWING SEPARATE C CALLS TO ABDDEC AND ABDSOL. C C FIRST, CALL THE DECOMPOSITION ROUTINE ABDDEC. C C CALL ABDDEC(N, TOPBLK, NTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) * C UNLESS A IS NUMERICALLY SINGULAR OR THE PARAMETERS ARE C IMPROPERLY DEFINED, CALL ABDSOL. * C IF ( IFLAG .EQ. 0 ) THEN * C CALL ABDSOL(N, TOPBLK, NTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, B,X) C ENDIF * WRITE(NOUT,1000)IFLAG IF ( IFLAG .EQ. 0 ) THEN WRITE(NOUT,2000) WRITE(NOUT,3000)(PIVOT(I),I=1,N-1) WRITE(NOUT,4000) WRITE(NOUT,5000)(X(I),I=1,N) ELSE IF ( IFLAG .EQ. 1 ) WRITE(NOUT,7000) IF ( IFLAG .EQ. -1) WRITE(NOUT,6000) ENDIF STOP 1000 FORMAT(9H IFLAG = ,I4) 2000 FORMAT(16H PIVOT VECTOR : ) 3000 FORMAT(15I4) 4000 FORMAT(12H SOLUTION : ) 5000 FORMAT(2E17.8) 6000 FORMAT(30H MATRIX APPEARS TO BE SINGULAR) 7000 FORMAT(20H IMPROPER PARAMETERS) END C C####################################################################### C SUBROUTINE ABDLEQ(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, B, X, IFLAG) * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A (WITH THE ALMOST BLOCK C DIAGONAL STRUCTURE ARISING WHEN SPLINE COLLOCATION AT GAUSSIAN C POINTS WITH MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT C BOUNDARY VALUE DIFFERENTIAL EQUATIONS WITH SEPARATED BOUNDARY C CONDITIONS), TOGETHER WITH THE CORRESPONDING RIGHT HAND SIDE. C C ABDLEQ CALLS A DECOMPOSITION ROUTINE, ABDDEC, FOLLOWED BY THE C CORRESPONDING SOLVE ROUTINE, ABDSOL, TO OBTAIN THE SOLUTION C OF THE SYSTEM. THE ALGORITHM IMPLEMENTS AN ALTERNATE COLUMN C AND ROW PIVOTING TECHNIQUE. C C THIS IS THE DOUBLE PRECISION VERSION OF ABDLEQ. THE SINGLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL DOUBLE PRECISION C DECLARATIONS BY REAL, AND REPLACING THE CALLS TO DOUBLE PRECISION C BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY SINGLE PRECISION CALLS. C C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C NBLOKS * (NRWBLK + NOVRLP) + NOVRLP C UNCHANGED ON EXIT. C C TOPBLK - DOUBLE PRECISION (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE NUMBER OF C ROWS IN ARRAY2(:,:, K). C UNCHANGED ON EXIT. C C ARRAY1 - DOUBLE PRECISION (NRWBLK, NRWBLK+NOVRLP,NBLOKS) C ARRAY1(:,:, K) CONTAINS THE K-TH NRWBLK C BY NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWBLK - INTEGER C THE NUMBER OF ROWS IN ARRAY1(:,:, K) C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY1 AND ARRAY2. C UNCHANGED ON EXIT. C C ARRAY2 - DOUBLE PRECISION(NOVRLP,NOVRLP+NRWBLK+NOVRLP,NBLOKS) C ARRAY2(:,:, K) CONTAINS THE K-TH NOVRLP C BY NOVRLP+NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C IT HAS THE FORM [AR2 I], WHERE AR2 IS NOVRLP BY C NRWBLK+NOVRLP, AND I IS THE NOVRLP BY NOVRLP C IDENTITY. C NOTE: THE IDENTITY MATRIX MUST BE SPECIFIED BY C THE USER. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C BOTBLK - DOUBLE PRECISION(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C ON EXIT, CONTAINS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION. C C B - DOUBLE PRECISION(N) C THE RIGHT HAND SIDES. C UNCHANGED ON EXIT. C C X - DOUBLE PRECISION(N) C WORK SPACE FOR THE SOLUTION OF THE LINEAR SYSTEM. C IF IFLAG = 0, CONTAINS THE SOLUTION OF AX = B ON C EXIT. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C ROUTINES CALLED : C C ABDDEC : TO PERFORM THE DECOMPOSITION OF THE COEFFICIENT C MATRIX. THE CALLING SEQUENCE IS GIVEN BY : C C CALL ABDDEC(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) C C C ABDSOL : TO SOLVE THE DECOMPOSED SYSTEM. C THE CALLING SEQUENCE IS GIVEN BY : C C CALL ABDSOL(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, B, X) C C C IN ADDITION, THE DECOMPOSITION AND SOLVE ROUTINES CALL C THE BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) : C C IDAMAX, DSWAP, DSCAL, DAXPY, DDOT. C C TO CONVERT ABDLEQ TO SINGLE PRECISION, THESE MUST BE C REPLACED BY THEIR SINGLE PRECISION EQUIVALENTS : C C ISAMAX, SSWAP, SSCAL, SAXPY, SDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C * INTEGER N,NRWTOP,NOVRLP,NRWBLK,NBLOKS,PIVOT(N),IFLAG * DOUBLE PRECISION TOPBLK,ARRAY1,ARRAY2,BOTBLK,B,X * DIMENSION TOPBLK(NRWTOP,*),ARRAY1(NRWBLK, NRWBLK+NOVRLP, *), * ARRAY2(NOVRLP, NOVRLP+NRWBLK+NOVRLP, *), * BOTBLK(NOVRLP-NRWTOP, *),B(N),X(N) C C FIRST, CALL THE FACTORING ROUTINE ABDDEC TO OBTAIN THE C PLBUQ FACTORIZATION. C CALL ABDDEC(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) C C IF A IS NUMERICALLY SINGULAR OR THE PARAMETERS ARE IMPROPERLY C DEFINED, RETURN. C IF ( IFLAG .NE. 0 ) RETURN C C NOW, CALL THE SOLVE ROUTINE ABDSOL. C CALL ABDSOL(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, B,X) C RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABDLEQ. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END * SUBROUTINE ABDDEC(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A WITH THE ALMOST BLOCK C DIAGONAL STRUCTURE ARISING WHEN SPLINE COLLOCATION AT GAUSSIAN C POINTS WITH MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT C BOUNDARY VALUE PROBLEMS WITH SEPARATED BOUNDARY CONDITIONS. C C THE ALGORITHM IMPLEMENTS AN ALTERNATE COLUMN AND ROW PIVOTING C TECHNIQUE. C C THIS IS THE DOUBLE PRECISION VERSION OF ABDDEC. THE SINGLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL DOUBLE PRECISION C DECLARATIONS BY REAL, AND REPLACING THE CALLS TO DOUBLE PRECISION C BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY SINGLE PRECISION CALLS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C NBLOKS * (NRWBLK + NOVRLP) + NOVRLP C UNCHANGED ON EXIT. C C TOPBLK - DOUBLE PRECISION (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE NUMBER OF C ROWS IN ARRAY2(:,:,K). C UNCHANGED ON EXIT. C C ARRAY1 - DOUBLE PRECISION (NRWBLK, NRWBLK+NOVRLP,NBLOKS) C ARRAY1(:,:,K) CONTAINS THE K-TH NRWBLK C BY NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWBLK - INTEGER C THE NUMBER OF ROWS IN ARRAY1(:,:,K) C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY1 AND ARRAY2. C UNCHANGED ON EXIT. C C ARRAY2 - DOUBLE PRECISION(NOVRLP,NOVRLP+NRWBLK+NOVRLP,NBLOKS) C ARRAY2(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP+NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C IT HAS THE FORM [AR2 I], WHERE AR2 IS NOVRLP BY C NRWBLK+NOVRLP, AND I IS THE NOVRLP BY NOVRLP C IDENTITY. C NOTE: THE IDENTITY MATRIX MUST BE SPECIFIED BY C THE USER. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C BOTBLK - DOUBLE PRECISION(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C ON EXIT, CONTAINS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C C ROUTINES CALLED : C C THE DECOMPOSITION ROUTINE CALLS THE BASIC LINEAR ALGEBRA C SUBROUTINES (BLAS) : C C IDAMAX, DSWAP, DSCAL, DAXPY. C C TO CONVERT ABDDEC TO SINGLE PRECISION, THESE MUST BE C REPLACED BY THEIR SINGLE PRECISION EQUIVALENTS : C C ISAMAX, SSWAP, SSCAL, SAXPY. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NRWTOP,NOVRLP,NRWBLK,NBLOKS,PIVOT(N), * IFLAG * DOUBLE PRECISION TOPBLK,ARRAY1,ARRAY2,BOTBLK * DIMENSION TOPBLK(NRWTOP,*),ARRAY1(NRWBLK, NRWBLK+NOVRLP, *), * ARRAY2(NOVRLP, NOVRLP+NRWBLK+NOVRLP, *), * BOTBLK(NOVRLP-NRWTOP, *) C C LOCAL VARIABLES : C INTEGER I,J,K,IPLUS1,JPLUS1,L,INCR,NCOLS,NRWBOT,NCLBLK,KP1, * INCRJ,INCRN,IPVT,JPVT,IMINN,JMINN,LOOP,JMINN2, * JPVT2,JMINN3,NTP1,NTPI,NTPA1,IDAMAX * DOUBLE PRECISION COLMAX,ROWMAX,PIVMAX, * COLPIV,ROWPIV,COLMLT,ROWMLT,ZERO,ONE DATA ZERO/ 0.D0 /, ONE / 1.D0 / C C DEFINE THE CONSTANTS USED THROUGHOUT. C NCOLS = NOVRLP+NRWBLK NRWBOT = NOVRLP - NRWTOP NCLBLK = NCOLS+NOVRLP NTP1 = NRWTOP+1 NTPA1 = NRWTOP+NRWBLK IFLAG = 0 PIVMAX = ZERO C C CHECK THE VALIDITY OF THE INPUT PARAMETER N. C IF(N.NE.NBLOKS*(NRWBLK+NOVRLP) + NOVRLP) GO TO 930 IF(NRWTOP.EQ.0) GO TO 930 IF(NRWTOP.EQ.NOVRLP) GO TO 930 C C FIRST, IN TOPBLK... C C APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN PIVOTING. C C C DO 195 I = 1,NRWTOP IPLUS1 = I+1 C C DETERMINE COLUMN PIVOT AND PIVOT INDEX. C IPVT = IDAMAX(NOVRLP-I+1,TOPBLK(I,I),NRWTOP) + I - 1 COLMAX = ABS(TOPBLK(I,IPVT)) C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX .EQ. PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(I) = IPVT IF(IPVT.EQ.I) GO TO 150 C C INTERCHANGE COLUMNS IN TOPBLK. C CALL DSWAP(NRWTOP-I+1,TOPBLK(I,I),1,TOPBLK(I,IPVT),1) C C INTERCHANGE COLUMNS IN ARRAY1(:,:1). C CALL DSWAP(NRWBLK,ARRAY1(1,I,1),1,ARRAY1(1,IPVT,1),1) C C INTERCHANGE COLUMNS IN ARRAY2(:,:,1). C CALL DSWAP(NOVRLP,ARRAY2(1,I,1),1,ARRAY2(1,IPVT,1),1) 150 CONTINUE C C COMPUTE MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C * COLPIV = ONE/TOPBLK(I,I) CALL DSCAL(NOVRLP-I,COLPIV,TOPBLK(I,IPLUS1),NRWTOP) DO 190 J = IPLUS1,NOVRLP COLMLT = TOPBLK(I,J) CALL DAXPY(NRWTOP-I,-COLMLT,TOPBLK(IPLUS1,I),1, * TOPBLK(IPLUS1,J),1) * C C ADJUST COLUMNS IN ARRAY1(:,:,1). C CALL DAXPY(NRWBLK,-COLMLT,ARRAY1(1,I,1),1,ARRAY1(1,J,1),1) C C ADJUST COLUMNS IN ARRAY2(:,:,1). C CALL DAXPY(NOVRLP,-COLMLT,ARRAY2(1,I,1),1,ARRAY2(1,J,1),1) * 190 CONTINUE * 195 CONTINUE * INCR = NRWTOP C C C C IN EACH BLOCK... C C APPLY COLUMN ELIMINATION WITH COLUMN PIVOTING IN ARRAY1, C APPLY ROW ELIMINATION WITH ROW PIVOTING ON ARRAY2, C COLUMN ELIMINATION WITH COLUMN PIVOTING ON PART OF C IDENTITY IN ARRAY2(:,:,K) C THEN C C APPLY ROW ELIMINATION WITH ROW PIVOTING ON BOTBLK. C C C DO 745 K = 1,NBLOKS KP1 = K+1 C C C FIRST, C APPLY NRWBLK COLUMN ELIMINATIONS WITH COLUMN PIVOTING C ON ARRAY1(:,:,K). C C DO 295 I = 1,NRWBLK C C DETERMINE COLUMN PIVOT WITH PIVOT INDEX. C IPLUS1 = I+1 NTPI = NRWTOP+I IPVT = NTPI INCRN = INCR + I IPVT = IDAMAX(NCOLS-NTPI+1,ARRAY1(I,NTPI,K),NRWBLK)+NTPI-1 COLMAX = ABS(ARRAY1(I,IPVT,K)) C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(INCRN) = INCR + IPVT - NRWTOP IF(IPVT.EQ.NTPI) GO TO 245 C C START IN ARRAY1(:,:,K)... C CALL DSWAP(NRWBLK-I+1,ARRAY1(I,NTPI,K),1,ARRAY1(I,IPVT,K),1) C C INTERCHANGE COLUMNS IN ARRAY2( , , K). C CALL DSWAP(NOVRLP,ARRAY2(1,NTPI,K),1,ARRAY2(1,IPVT,K),1) * 245 CONTINUE C C COMPUTE MULTIPLIERS AND PERFORM COLUMN ELIMINATION. C COLPIV = ONE/ARRAY1(I,NTPI,K) CALL DSCAL(NCOLS-NTPI,COLPIV,ARRAY1(I,NTPI+1,K),NRWBLK) DO 290 J = (NTPI+1),NCOLS COLMLT = ARRAY1(I,J,K) CALL DAXPY(NRWBLK-I,-COLMLT,ARRAY1(IPLUS1,NTPI,K),1, * ARRAY1(IPLUS1,J,K),1) C C ADJUST COLUMNS IN ARRAY2(:,:,K). C CALL DAXPY(NOVRLP,-COLMLT,ARRAY2(1,NTPI,K),1, * ARRAY2(1,J,K),1) * 290 CONTINUE * 295 CONTINUE INCR = INCR + NRWBLK C C C PERFORM NRWBOT=(NOVRLP-NRWTOP) ROW ELIMINATIONS WITH ROW PIVOTIN C ON ARRAY2(:,:,K). C C * DO 530 J = (NTPA1+1),NCOLS JPLUS1 = J+1 JMINN = J-NTPA1 INCRJ = INCR + JMINN LOOP = JMINN + 1 IPVT = IDAMAX(NOVRLP-JMINN+1,ARRAY2(JMINN,J,K),1) + JMINN-1 ROWMAX = ABS(ARRAY2(IPVT,J,K)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C C IF NECESSARY, INTERCHANGE ROWS IN ARRAY2(:,:,K) C CORRESPONDING COLUMNS IN IDENTITY IN ARRAY(:,:,K), C COLUMNS IN ARRAY1(:,:,K+1), IF K < NBLOKS C COLUMNS IN ARRAY2(:,:,K+1), IF K < NBLOKS C AND COLUMNS IN BOTBLK IF K=NBLOKS. C C PIVOT(INCRJ) = INCR+IPVT IF(IPVT.EQ.JMINN) GO TO 500 C C INTERCHANGE ROWS. C CALL DSWAP(NCLBLK-J+1,ARRAY2(IPVT,J,K),NOVRLP, * ARRAY2(JMINN,J,K),NOVRLP) C C INTERCHANGE COLUMNS IN IDENTITY IN ARRAY(:,:,K). C JMINN2 = JMINN+NCOLS JPVT2 = IPVT+NCOLS CALL DSWAP(NOVRLP-JMINN+1,ARRAY2(JMINN,JPVT2,K),1, * ARRAY2(JMINN,JMINN2,K),1) IF(K.LT.NBLOKS)THEN C C INTERCHANGE COLUMNS IN ARRAY1(:,:,K+1). C CALL DSWAP(NRWBLK,ARRAY1(1,IPVT,KP1),1, * ARRAY1(1,JMINN,KP1),1) C C INTERCHANGE COLUMNS IN ARRAY2(:,:,K+1) IF K < NBLOKS. C CALL DSWAP(NOVRLP,ARRAY2(1,IPVT,KP1),1, * ARRAY2(1,JMINN,KP1),1) ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K = NBLOKS. C CALL DSWAP(NRWBOT,BOTBLK(1,IPVT),1,BOTBLK(1,JMINN),1) ENDIF * 500 CONTINUE C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/ARRAY2(JMINN,J,K) CALL DSCAL(NOVRLP-LOOP+1,ROWPIV,ARRAY2(LOOP,J,K),1) C C PERFORM ROW ELIMINATION. C DO 520 L = (J+1),(J+NRWBOT) ROWMLT = ARRAY2(JMINN,L,K) CALL DAXPY(NOVRLP-LOOP+1,-ROWMLT,ARRAY2(LOOP,J,K),1, * ARRAY2(LOOP,L,K),1) 520 CONTINUE * 530 CONTINUE * INCR = INCR + NRWBOT C C C IN IDENTITY IN ARRAY2(:,:,K) PERFORM NRWTOP COLUMN C ELIMINATIONS WITH COLUMN PIVOTING. C C DO 740 I = (NRWBOT+1),NOVRLP IPLUS1 = I+1 IMINN = I-NRWBOT JMINN = I+NTPA1 C C DETERMINE COLUMN PIVOT WITH PIVOT INDEX. C C JPVT = JMINN INCRN = INCR + IMINN JPVT = IDAMAX(NRWBOT+1,ARRAY2(I,JMINN,K),NOVRLP)+JMINN-1 C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(INCRN) = INCRN+JPVT-JMINN IF(JPVT.EQ.JMINN) GO TO 660 C C INTERCHANGE COLUMNS IN IDENTITY IN ARRAY2(:,:,K). C CALL DSWAP(NOVRLP-I+1,ARRAY2(I,JPVT,K),1, * ARRAY2(I,JMINN,K),1) C C INTERCHANGE COLUMNS IN ARRAY1(:,:,K+1) C JPVT2 = JPVT-JMINN+IMINN * IF(K.LT.NBLOKS)THEN * CALL DSWAP(NRWBLK,ARRAY1(1,JPVT2,KP1),1, * ARRAY1(1,IMINN,KP1),1) * C C INTERCHANGE COLUMNS IN ARRAY2(:,:,K+1). C CALL DSWAP(NOVRLP,ARRAY2(1,JPVT2,KP1),1, * ARRAY2(1,IMINN,KP1),1) ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K = NBLOKS. C CALL DSWAP(NRWBOT,BOTBLK(1,JPVT2),1,BOTBLK(1,IMINN),1) * ENDIF * 660 CONTINUE C C COMPUTE COLUMN MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C COLPIV = ONE/ARRAY2(I,JMINN,K) CALL DSCAL(NRWBOT,COLPIV,ARRAY2(I,JMINN+1,K),NOVRLP) DO 730 J=(JMINN+1),(JMINN+NRWBOT) COLMLT = ARRAY2(I,J,K) CALL DAXPY(NOVRLP-I,-COLMLT,ARRAY2(IPLUS1,JMINN,K),1, * ARRAY2(IPLUS1,J,K),1) C C ADJUST COLUMNS IN ARRAY1(:,:,K+1), ARRAY2(:,:,K+1) C OR BOTBLK. C JMINN2 = J-JMINN+IMINN JMINN3 = JMINN-NCOLS IF(K.LT.NBLOKS)THEN C C FIRST ADJUST COLUMNS IN ARRAY1(:,:,K+1) C CALL DAXPY(NRWBLK,-COLMLT,ARRAY1(1,JMINN3,KP1),1, * ARRAY1(1,JMINN2,KP1),1) * C C ADJUST COLUMNS IN ARRAY2(:,:,K+1) C CALL DAXPY(NOVRLP,-COLMLT,ARRAY2(1,JMINN3,KP1),1, * ARRAY2(1,JMINN2,KP1),1) * ELSE C C ADJUST COLUMNS IN BOTBLK. C CALL DAXPY(NRWBOT,-COLMLT,BOTBLK(1,JMINN3),1, * BOTBLK(1,JMINN2),1) * ENDIF 730 CONTINUE * 740 CONTINUE * INCR = INCR + NRWTOP * 745 CONTINUE * C C C C PERFORM (NOVRLP-NRWTOP-1)=NRWBOT-1 ROW ELIMINATION(S) WITH ROW C PIVOTING ON BOTBLK. C C C IF(NRWBOT.EQ.1) GO TO 910 * DO 900 J = (NRWTOP+1),(NOVRLP-1) JPLUS1 = J+1 JMINN = J-NRWTOP INCRJ = INCR + JMINN LOOP = JMINN+1 IPVT = IDAMAX(NRWBOT-JMINN+1,BOTBLK(JMINN,J),1) + JMINN-1 ROWMAX = ABS(BOTBLK(IPVT,J)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C IF NECESSARY, INTERCHANGE ROWS IN BOTBLK. C PIVOT(INCRJ) = INCR + IPVT IF(IPVT.EQ.JMINN) GO TO 770 CALL DSWAP(NOVRLP-J+1,BOTBLK(IPVT,J),NRWBOT, * BOTBLK(JMINN,J),NRWBOT) * 770 CONTINUE C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/BOTBLK(JMINN,J) CALL DSCAL(NRWBOT-JMINN,ROWPIV,BOTBLK(LOOP,J),1) C C PERFORM ROW ELIMINATIONS. C DO 800 L = (J+1),NOVRLP ROWMLT = BOTBLK(JMINN,L) CALL DAXPY(NRWBOT-JMINN,-ROWMLT,BOTBLK(LOOP,J),1, * BOTBLK(LOOP,L),1) 800 CONTINUE * 900 CONTINUE * 910 CONTINUE C C DONE PROVIDED LAST ELEMENT IS NOT ZERO. C IF(PIVMAX+ABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN C C C MATRIX IS SINGULAR -- SET IFLAG = -1 AND TERMINATE. C C 920 CONTINUE IFLAG = -1 RETURN C C C ERROR IN THE PARAMETER N SPECIFIED. C C 930 CONTINUE IFLAG = 1 RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABDDEC. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END * * * * * SUBROUTINE ABDSOL(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, B, X) * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C GIVEN THE DECOMPOSED FORM OF THE MATRIX A WITH THE PIVOT VECTOR C (GENERATED IN ABDDEC) AND THE VECTOR B, THIS ROUTINE SOLVES C FOR X IN AX = B. C C C THIS IS THE DOUBLE PRECISION VERSION OF ABDSOL. THE SINGLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL DOUBLE PRECISION C DECLARATIONS BY REAL, AND REPLACING THE CALLS TO DOUBLE PRECISION C BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY SINGLE PRECISION CALLS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C NBLOKS * (NRWBLK + NOVRLP) + NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - DOUBLE PRECISION (NRWTOP, NOVRLP) C THE DECOMPOSED FIRST BLOCK OF A. C UNCHANGED ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK, AND THE C NUMBER OF ROWS IN ARRAY2(:,:,K). C UNCHANGED ON EXIT. C C ARRAY1 - DOUBLE PRECISION (NRWBLK, NRWBLK+NOVRLP,NBLOKS) C ARRAY1(:,:,K) CONTAINS THE K-TH NRWBLK C BY NRWBLK+NOVRLP DECOMPOSED BLOCK OF C THE MATRIX A. C UNCHANGED ON EXIT. C C NRWBLK - INTEGER C THE NUMBER OF ROWS IN ARRAY1(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY1 AND ARRAY2. C UNCHANGED ON EXIT. C C ARRAY2 - DOUBLE PRECISION(NOVRLP,NOVRLP+NRWBLK+NOVRLP,NBLOK) C ARRAY2(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP+NRWBLK+NOVRLP DECOMPOSED BLOCK OF C THE MATRIX A. C UNCHANGED ON EXIT. C C BOTBLK - DOUBLE PRECISION(NOVRLP-NRWTOP, NOVRLP) C THE LAST DECOMPOSED BLOCK OF THE MATRIX A. C UNCHANGED ON EXIT. C C PIVOT - INTEGER(N) C CONTAINS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION ROUTINE ABDDEC. C UNCHANGED ON EXIT. C C B - DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR. C UNCHANGED ON EXIT. C C X - DOUBLE PRECISION(N) C WORK SPACE FOR THE SOLUTION OF THE LINEAR SYSTEM. C IF IFLAG = 0, CONTAINS THE SOLUTION OF AX = B ON C EXIT. C C C ROUTINES CALLED : C C THE BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) : C C DAXPY, DDOT. C C TO CONVERT ABDSOL TO SINGLE PRECISION, THESE MUST BE C REPLACED BY THEIR SINGLE PRECISION EQUIVALENTS : C C SAXPY, SDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NBLOKS,NRWTOP,NOVRLP,NRWBLK,PIVOT(N) * DOUBLE PRECISION TOPBLK,ARRAY1,ARRAY2,BOTBLK,B,X * DIMENSION TOPBLK(NRWTOP,*),ARRAY1(NRWBLK,NRWBLK+NOVRLP,*), * ARRAY2(NOVRLP,NOVRLP+NRWBLK+NOVRLP,*), * BOTBLK(NOVRLP-NRWTOP,*),B(N),X(N) C C LOCAL VARIABLES : C INTEGER NCOLS,NRWBOT,L,NCLBLK,NTPJ,I,J,K,LOOP,INCR,INCRJ, * INCRI,INCR1,INCRTP,JPIVOT,ITEMP,IPVTI,L1,INCRN,LL, * NRWBTL,IPLUSN,IPVTN,NTPA1,NTP1,NTPCL * DOUBLE PRECISION XINCRI,SWAP,DDOT C C DEFINE THE CONSTANTS USED THROUGHOUT. C * NCOLS = NOVRLP+NRWBLK NRWBOT = NOVRLP-NRWTOP NCLBLK = NCOLS+NOVRLP NTPA1 = NRWTOP+NRWBLK NTP1 = NRWTOP+1 NTPCL = NRWTOP+NCOLS DO 100 J = 1,N X(J) = B(J) 100 CONTINUE * C C C C FORWARD RECURSION... C C C C FIRST IN TOPBLK... C C FORWARD SOLUTION. C C C DO 130 J = 1,NRWTOP X(J) = X(J)/TOPBLK(J,J) LOOP = J+1 CALL DAXPY(NRWTOP-J,-X(J),TOPBLK(LOOP,J),1,X(LOOP),1) 130 CONTINUE * INCR = NRWTOP C C C C IN EACH BLOCK K, K = 1,NBLOKS... C C C DO 360 K = 1,NBLOKS C C FORWARD MODIFICATION. C INCRTP = INCR DO 140 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL DAXPY(NRWBLK,-X(INCRJ),ARRAY1(1,J,K),1,X(INCRTP+1),1) 140 CONTINUE C C FORWARD SOLUTION IN ARRAY1(:,:,K). C DO 170 J = 1,NRWBLK INCRJ = INCR+J NTPJ = NRWTOP+J X(INCRJ) = X(INCRJ)/ARRAY1(J,NTPJ,K) LOOP = J+1 CALL DAXPY(NRWBLK-J,-X(INCRJ),ARRAY1(LOOP,NTPJ,K),1, * X(INCRJ+1),1) 170 CONTINUE * C C FORWARD MODIFICATION. C INCRTP = INCR+NRWBLK * DO 240 J = 1,NTPA1 INCRJ = INCR-NRWTOP+J INCR1 = INCRTP+1 CALL DAXPY(NOVRLP,-X(INCRJ),ARRAY2(1,J,K),1,X(INCR1),1) 240 CONTINUE C C FORWARD ELIMINATION. C C DO 260 J = (NTPA1+1),NCOLS INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.EQ.INCRJ) GO TO 245 SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP 245 CONTINUE LOOP = J-(NTPA1-1) CALL DAXPY(NOVRLP-LOOP+1,-X(INCRJ),ARRAY2(LOOP,J,K),1, * X(INCRTP+LOOP),1) 260 CONTINUE C C FORWARD SOLUTION. C IF(NRWTOP.LT.1) GO TO 285 C DO 280 J = (NCOLS+1),NTPCL INCRJ = INCR-NRWTOP+J ITEMP = J-NTPA1 X(INCRJ) = X(INCRJ)/ARRAY2(ITEMP,J,K) LOOP = ITEMP+1 CALL DAXPY(NOVRLP-ITEMP,-X(INCRJ),ARRAY2(LOOP,J,K),1, * X(LOOP+INCRTP),1) 280 CONTINUE 285 CONTINUE INCR = INCR+NCOLS 360 CONTINUE * C C IN BOTBLK ... C C FORWARD MODIFICATION. C INCRTP = INCR DO 300 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL DAXPY(NRWBOT,-X(INCRJ),BOTBLK(1,J),1,X(INCRTP+1),1) 300 CONTINUE C C FORWARD ELIMINATION. C IF(NRWBOT.EQ.1) GO TO 350 DO 340 J = (NRWTOP+1),(NOVRLP-1) INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.EQ.INCRJ) GO TO 325 SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP 325 CONTINUE LOOP = J - (NRWTOP-1) CALL DAXPY(NRWBOT-LOOP+1,-X(INCRJ),BOTBLK(LOOP,J),1, * X(INCRTP+LOOP),1) 340 CONTINUE 350 CONTINUE C C C C BACKWARD RECURSION... C C C C FIRST IN BOTBLK... C C BACKWARD SOLUTION. C C DO 430 LL = 1,NRWBOT J = NOVRLP+1-LL NRWBTL = NRWBOT + 1 - LL INCRJ = INCR+NRWBTL X(INCRJ) = X(INCRJ)/BOTBLK(NRWBTL,J) LOOP = NRWBOT-LL CALL DAXPY(LOOP,-X(INCRJ),BOTBLK(1,J),1,X(INCR+1),1) 430 CONTINUE * C C IN EACH BLOCK GOING IN REVERSE ORDER... C DO 520 L = 1,NBLOKS C C BACKWARD ELIMINATION. C K = NBLOKS+1-L INCR = INCR-NRWTOP C C BACKWARD ELIMINATION IN THE IDENTITY IN ARRAY2. C DO 470 L1 = 1,NRWTOP I = NTP1-L1 ITEMP = NRWBOT+I IPLUSN = NCOLS+I INCRN = INCR+I X(INCRN) = X(INCRN) - DDOT(NRWBOT, * ARRAY2(ITEMP,IPLUSN+1,K),NOVRLP,X(INCRN+1),1) IPVTN = PIVOT(INCRN) IF(INCRN.EQ.IPVTN) GO TO 465 SWAP = X(INCRN) X(INCRN) = X(IPVTN) X(IPVTN) = SWAP 465 CONTINUE 470 CONTINUE INCR = INCR-NRWBOT C C BACKWARD MODIFICATION, SOLUTION AND ADJUSTMENT (ORDERING) OF C THE X'S. C DO 495 L1 = 1,NRWBOT I = NRWBOT+1-L1 IPLUSN = NTPA1+I LOOP = IPLUSN + 1 INCRI = INCR + I C C BACKWARD MODIFICATION. C XINCRI = X(INCRI) X(INCRI) = X(INCRI) - DDOT(NRWBOT-1,ARRAY2(I,LOOP,K), * NOVRLP,X(INCR+I+1),1) X(INCRI) = X(INCRI)-X(INCRI+NRWBOT) C C BACKWARD SOLUTION. C X(INCRI) = X(INCRI)/ARRAY2(I,IPLUSN,K) C C COLUMN PIVOTING ADJUSTMENT ON THE X'S. C IPVTN = PIVOT(INCRI) IF(INCRI.EQ.IPVTN) GO TO 490 INCRI = INCRI+NRWBOT IPVTN = IPVTN+NRWBOT SWAP = X(INCRI) X(INCRI) = X(IPVTN) X(IPVTN) = SWAP 490 CONTINUE 495 CONTINUE INCR = INCR-NRWBLK DO 450 L1 = 1,NRWBLK I = NRWBLK+1-L1 IPLUSN = NRWTOP+I INCRN = INCR+I X(INCRN) = X(INCRN)-DDOT(NCOLS-IPLUSN, * ARRAY1(I,IPLUSN+1,K),NRWBLK,X(INCR+I+1),1) IPVTN = PIVOT(INCRN) IF(INCRN.EQ.IPVTN) GO TO 445 SWAP = X(INCRN) X(INCRN) = X(IPVTN) X(IPVTN) = SWAP 445 CONTINUE 450 CONTINUE 520 CONTINUE C C BACKWARD ELIMINATION IN TOPBLK. C DO 510 L1 = 1,NRWTOP I = NTP1-L1 LOOP = I + 1 X(I) = X(I) - DDOT(NOVRLP-I,TOPBLK(I,LOOP),NRWTOP, * X(LOOP),1) IPVTI = PIVOT(I) IF(I.EQ.IPVTI) GO TO 505 SWAP = X(I) X(I) = X(IPVTI) X(IPVTI) = SWAP 505 CONTINUE 510 CONTINUE RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABDSOL. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END C C####################################################################### C SUBROUTINE ABDLEQ(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, B, X, IFLAG) * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A (WITH THE ALMOST BLOCK C DIAGONAL STRUCTURE ARISING WHEN SPLINE COLLOCATION AT GAUSSIAN C POINTS WITH MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT C BOUNDARY VALUE DIFFERENTIAL EQUATIONS WITH SEPARATED BOUNDARY C CONDITIONS), TOGETHER WITH THE CORRESPONDING RIGHT HAND SIDE. C C ABDLEQ CALLS A DECOMPOSITION ROUTINE, ABDDEC, FOLLOWED BY THE C CORRESPONDING SOLVE ROUTINE, ABDSOL, TO OBTAIN THE SOLUTION C OF THE SYSTEM. THE ALGORITHM IMPLEMENTS AN ALTERNATE COLUMN C AND ROW PIVOTING TECHNIQUE. C C THIS IS THE SINGLE PRECISION VERSION OF ABDLEQ. THE DOUBLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL REAL C DECLARATIONS BY DOUBLE PRECISION, AND REPLACING THE CALLS TO C SINGLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C DOUBLE PRECISION CALLS. C C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C NBLOKS * (NRWBLK + NOVRLP) + NOVRLP C UNCHANGED ON EXIT. C C TOPBLK - REAL (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE NUMBER OF C ROWS IN ARRAY2(:,:, K). C UNCHANGED ON EXIT. C C ARRAY1 - REAL (NRWBLK, NRWBLK+NOVRLP,NBLOKS) C ARRAY1(:,:, K) CONTAINS THE K-TH NRWBLK C BY NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWBLK - INTEGER C THE NUMBER OF ROWS IN ARRAY1(:,:, K) C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY1 AND ARRAY2. C UNCHANGED ON EXIT. C C ARRAY2 - REAL(NOVRLP,NOVRLP+NRWBLK+NOVRLP,NBLOKS) C ARRAY2(:,:, K) CONTAINS THE K-TH NOVRLP C BY NOVRLP+NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C IT HAS THE FORM [AR2 I], WHERE AR2 IS NOVRLP BY C NRWBLK+NOVRLP, AND I IS THE NOVRLP BY NOVRLP C IDENTITY. C NOTE: THE IDENTITY MATRIX MUST BE SPECIFIED BY C THE USER. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C BOTBLK - REAL(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C ON EXIT, CONTAINS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION. C C B - REAL(N) C THE RIGHT HAND SIDES. C UNCHANGED ON EXIT. C C X - REAL(N) C WORK SPACE FOR THE SOLUTION OF THE LINEAR SYSTEM. C IF IFLAG = 0, CONTAINS THE SOLUTION OF AX = B ON C EXIT. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C ROUTINES CALLED : C C ABDDEC : TO PERFORM THE DECOMPOSITION OF THE COEFFICIENT C MATRIX. THE CALLING SEQUENCE IS GIVEN BY : C C CALL ABDDEC(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) C C C ABDSOL : TO SOLVE THE DECOMPOSED SYSTEM. C THE CALLING SEQUENCE IS GIVEN BY : C C CALL ABDSOL(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, C * NBLOKS, ARRAY2, BOTBLK, PIVOT, B, X) C C C IN ADDITION, THE DECOMPOSITION AND SOLVE ROUTINES CALL C THE BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) : C C ISAMAX, SSWAP, SSCAL, SAXPY, SDOT. C C TO CONVERT ABDLEQ TO DOUBLE PRECISION, THESE MUST BE C REPLACED BY THEIR DOUBLE PRECISION EQUIVALENTS : C C IDAMAX, DSWAP, DSCAL, DAXPY, DDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C * INTEGER N,NRWTOP,NOVRLP,NRWBLK,NBLOKS,PIVOT(N),IFLAG * REAL TOPBLK,ARRAY1,ARRAY2,BOTBLK,B,X * DIMENSION TOPBLK(NRWTOP,*),ARRAY1(NRWBLK, NRWBLK+NOVRLP, *), * ARRAY2(NOVRLP, NOVRLP+NRWBLK+NOVRLP, *), * BOTBLK(NOVRLP-NRWTOP, *),B(N),X(N) C C FIRST, CALL THE FACTORING ROUTINE ABDDEC TO OBTAIN THE C PLBUQ FACTORIZATION. C CALL ABDDEC(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) C C IF A IS NUMERICALLY SINGULAR OR THE PARAMETERS ARE IMPROPERLY C DEFINED, RETURN. C IF ( IFLAG .NE. 0 ) RETURN C C NOW, CALL THE SOLVE ROUTINE ABDSOL. C CALL ABDSOL(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, B,X) C RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABDLEQ. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END * SUBROUTINE ABDDEC(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, IFLAG) * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS PROCEDURE TAKES AS INPUT A MATRIX A WITH THE ALMOST BLOCK C DIAGONAL STRUCTURE ARISING WHEN SPLINE COLLOCATION AT GAUSSIAN C POINTS WITH MONOMIAL BASIS FUNCTIONS IS APPLIED TO TWO-POINT C BOUNDARY VALUE PROBLEMS WITH SEPARATED BOUNDARY CONDITIONS. C C THE ALGORITHM IMPLEMENTS AN ALTERNATE COLUMN AND ROW PIVOTING C TECHNIQUE. C C THIS IS THE SINGLE PRECISION VERSION OF ABDDEC. THE DOUBLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL REAL C DECLARATIONS BY DOUBLE PRECISION, AND REPLACING THE CALLS TO C SINGLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) C BY DOUBLE PRECISION CALLS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C NBLOKS * (NRWBLK + NOVRLP) + NOVRLP C UNCHANGED ON EXIT. C C TOPBLK - REAL (NRWTOP, NOVRLP) C THE FIRST BLOCK OF A TO BE DECOMPOSED. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK AND THE NUMBER OF C ROWS IN ARRAY2(:,:,K). C UNCHANGED ON EXIT. C C ARRAY1 - REAL (NRWBLK, NRWBLK+NOVRLP,NBLOKS) C ARRAY1(:,:,K) CONTAINS THE K-TH NRWBLK C BY NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C NRWBLK - INTEGER C THE NUMBER OF ROWS IN ARRAY1(:,:,K) C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY1 AND ARRAY2. C UNCHANGED ON EXIT. C C ARRAY2 - REAL(NOVRLP,NOVRLP+NRWBLK+NOVRLP,NBLOKS) C ARRAY2(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP+NRWBLK+NOVRLP BLOCK OF THE MATRIX A. C IT HAS THE FORM [AR2 I], WHERE AR2 IS NOVRLP BY C NRWBLK+NOVRLP, AND I IS THE NOVRLP BY NOVRLP C IDENTITY. C NOTE: THE IDENTITY MATRIX MUST BE SPECIFIED BY C THE USER. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C BOTBLK - REAL(NOVRLP-NRWTOP, NOVRLP) C THE LAST BLOCK OF THE MATRIX A. C CONTAINS DECOMPOSED BLOCK ON EXIT. C C PIVOT - INTEGER(N) C WORK SPACE FOR RECORDING PIVOT INFORMATION. C ON EXIT, CONTAINS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION. C C IFLAG - INTEGER C ON EXIT : C = 1, INPUT PARAMETER N IS NOT VALID C OR NRWTOP = 0, OR NRWTOP = NOVRLP. C = -1, MATRIX IS SINGULAR C = 0 OTHERWISE. C C C ROUTINES CALLED : C C THE DECOMPOSITION ROUTINE CALLS THE BASIC LINEAR ALGEBRA C SUBROUTINES (BLAS) : C C ISAMAX, SSWAP, SSCAL, SAXPY. C C TO CONVERT ABDDEC TO DOUBLE PRECISION, THESE MUST BE C REPLACED BY THEIR DOUBLE PRECISION EQUIVALENTS : C C ISAMAX, DSSWAP, DSSCAL, DSAXPY. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NRWTOP,NOVRLP,NRWBLK,NBLOKS,PIVOT(N), * IFLAG * REAL TOPBLK,ARRAY1,ARRAY2,BOTBLK * DIMENSION TOPBLK(NRWTOP,*),ARRAY1(NRWBLK, NRWBLK+NOVRLP, *), * ARRAY2(NOVRLP, NOVRLP+NRWBLK+NOVRLP, *), * BOTBLK(NOVRLP-NRWTOP, *) C C LOCAL VARIABLES : C INTEGER I,J,K,IPLUS1,JPLUS1,L,INCR,NCOLS,NRWBOT,NCLBLK,KP1, * INCRJ,INCRN,IPVT,JPVT,IMINN,JMINN,LOOP,JMINN2, * JPVT2,JMINN3,NTP1,NTPI,NTPA1,ISAMAX * REAL COLMAX,ROWMAX,PIVMAX, * COLPIV,ROWPIV,COLMLT,ROWMLT,ZERO,ONE DATA ZERO/ 0.D0 /, ONE / 1.D0 / C C DEFINE THE CONSTANTS USED THROUGHOUT. C NCOLS = NOVRLP+NRWBLK NRWBOT = NOVRLP - NRWTOP NCLBLK = NCOLS+NOVRLP NTP1 = NRWTOP+1 NTPA1 = NRWTOP+NRWBLK IFLAG = 0 PIVMAX = ZERO C C CHECK THE VALIDITY OF THE INPUT PARAMETER N. C IF(N.NE.NBLOKS*(NRWBLK+NOVRLP) + NOVRLP) GO TO 930 IF(NRWTOP.EQ.0) GO TO 930 IF(NRWTOP.EQ.NOVRLP) GO TO 930 C C FIRST, IN TOPBLK... C C APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN PIVOTING. C C C DO 195 I = 1,NRWTOP IPLUS1 = I+1 C C DETERMINE COLUMN PIVOT AND PIVOT INDEX. C IPVT = ISAMAX(NOVRLP-I+1,TOPBLK(I,I),NRWTOP) + I - 1 COLMAX = ABS(TOPBLK(I,IPVT)) C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX .EQ. PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(I) = IPVT IF(IPVT.EQ.I) GO TO 150 C C INTERCHANGE COLUMNS IN TOPBLK. C CALL SSWAP(NRWTOP-I+1,TOPBLK(I,I),1,TOPBLK(I,IPVT),1) C C INTERCHANGE COLUMNS IN ARRAY1(:,:1). C CALL SSWAP(NRWBLK,ARRAY1(1,I,1),1,ARRAY1(1,IPVT,1),1) C C INTERCHANGE COLUMNS IN ARRAY2(:,:,1). C CALL SSWAP(NOVRLP,ARRAY2(1,I,1),1,ARRAY2(1,IPVT,1),1) 150 CONTINUE C C COMPUTE MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C * COLPIV = ONE/TOPBLK(I,I) CALL SSCAL(NOVRLP-I,COLPIV,TOPBLK(I,IPLUS1),NRWTOP) DO 190 J = IPLUS1,NOVRLP COLMLT = TOPBLK(I,J) CALL SAXPY(NRWTOP-I,-COLMLT,TOPBLK(IPLUS1,I),1, * TOPBLK(IPLUS1,J),1) * C C ADJUST COLUMNS IN ARRAY1(:,:,1). C CALL SAXPY(NRWBLK,-COLMLT,ARRAY1(1,I,1),1,ARRAY1(1,J,1),1) C C ADJUST COLUMNS IN ARRAY2(:,:,1). C CALL SAXPY(NOVRLP,-COLMLT,ARRAY2(1,I,1),1,ARRAY2(1,J,1),1) * 190 CONTINUE * 195 CONTINUE * INCR = NRWTOP C C C C IN EACH BLOCK... C C APPLY COLUMN ELIMINATION WITH COLUMN PIVOTING IN ARRAY1, C APPLY ROW ELIMINATION WITH ROW PIVOTING ON ARRAY2, C COLUMN ELIMINATION WITH COLUMN PIVOTING ON PART OF C IDENTITY IN ARRAY2(:,:,K) C THEN C C APPLY ROW ELIMINATION WITH ROW PIVOTING ON BOTBLK. C C C DO 745 K = 1,NBLOKS KP1 = K+1 C C C FIRST, C APPLY NRWBLK COLUMN ELIMINATIONS WITH COLUMN PIVOTING C ON ARRAY1(:,:,K). C C DO 295 I = 1,NRWBLK C C DETERMINE COLUMN PIVOT WITH PIVOT INDEX. C IPLUS1 = I+1 NTPI = NRWTOP+I IPVT = NTPI INCRN = INCR + I IPVT = ISAMAX(NCOLS-NTPI+1,ARRAY1(I,NTPI,K),NRWBLK)+NTPI-1 COLMAX = ABS(ARRAY1(I,IPVT,K)) C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(INCRN) = INCR + IPVT - NRWTOP IF(IPVT.EQ.NTPI) GO TO 245 C C START IN ARRAY1(:,:,K)... C CALL SSWAP(NRWBLK-I+1,ARRAY1(I,NTPI,K),1,ARRAY1(I,IPVT,K),1) C C INTERCHANGE COLUMNS IN ARRAY2( , , K). C CALL SSWAP(NOVRLP,ARRAY2(1,NTPI,K),1,ARRAY2(1,IPVT,K),1) * 245 CONTINUE C C COMPUTE MULTIPLIERS AND PERFORM COLUMN ELIMINATION. C COLPIV = ONE/ARRAY1(I,NTPI,K) CALL SSCAL(NCOLS-NTPI,COLPIV,ARRAY1(I,NTPI+1,K),NRWBLK) DO 290 J = (NTPI+1),NCOLS COLMLT = ARRAY1(I,J,K) CALL SAXPY(NRWBLK-I,-COLMLT,ARRAY1(IPLUS1,NTPI,K),1, * ARRAY1(IPLUS1,J,K),1) C C ADJUST COLUMNS IN ARRAY2(:,:,K). C CALL SAXPY(NOVRLP,-COLMLT,ARRAY2(1,NTPI,K),1, * ARRAY2(1,J,K),1) * 290 CONTINUE * 295 CONTINUE INCR = INCR + NRWBLK C C C PERFORM NRWBOT=(NOVRLP-NRWTOP) ROW ELIMINATIONS WITH ROW PIVOTIN C ON ARRAY2(:,:,K). C C * DO 530 J = (NTPA1+1),NCOLS JPLUS1 = J+1 JMINN = J-NTPA1 INCRJ = INCR + JMINN LOOP = JMINN + 1 IPVT = ISAMAX(NOVRLP-JMINN+1,ARRAY2(JMINN,J,K),1) + JMINN-1 ROWMAX = ABS(ARRAY2(IPVT,J,K)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C C IF NECESSARY, INTERCHANGE ROWS IN ARRAY2(:,:,K) C CORRESPONDING COLUMNS IN IDENTITY IN ARRAY(:,:,K), C COLUMNS IN ARRAY1(:,:,K+1), IF K < NBLOKS C COLUMNS IN ARRAY2(:,:,K+1), IF K < NBLOKS C AND COLUMNS IN BOTBLK IF K=NBLOKS. C C PIVOT(INCRJ) = INCR+IPVT IF(IPVT.EQ.JMINN) GO TO 500 C C INTERCHANGE ROWS. C CALL SSWAP(NCLBLK-J+1,ARRAY2(IPVT,J,K),NOVRLP, * ARRAY2(JMINN,J,K),NOVRLP) C C INTERCHANGE COLUMNS IN IDENTITY IN ARRAY(:,:,K). C JMINN2 = JMINN+NCOLS JPVT2 = IPVT+NCOLS CALL SSWAP(NOVRLP-JMINN+1,ARRAY2(JMINN,JPVT2,K),1, * ARRAY2(JMINN,JMINN2,K),1) IF(K.LT.NBLOKS)THEN C C INTERCHANGE COLUMNS IN ARRAY1(:,:,K+1). C CALL SSWAP(NRWBLK,ARRAY1(1,IPVT,KP1),1, * ARRAY1(1,JMINN,KP1),1) C C INTERCHANGE COLUMNS IN ARRAY2(:,:,K+1) IF K < NBLOKS. C CALL SSWAP(NOVRLP,ARRAY2(1,IPVT,KP1),1, * ARRAY2(1,JMINN,KP1),1) ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K = NBLOKS. C CALL SSWAP(NRWBOT,BOTBLK(1,IPVT),1,BOTBLK(1,JMINN),1) ENDIF * 500 CONTINUE C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/ARRAY2(JMINN,J,K) CALL SSCAL(NOVRLP-LOOP+1,ROWPIV,ARRAY2(LOOP,J,K),1) C C PERFORM ROW ELIMINATION. C DO 520 L = (J+1),(J+NRWBOT) ROWMLT = ARRAY2(JMINN,L,K) CALL SAXPY(NOVRLP-LOOP+1,-ROWMLT,ARRAY2(LOOP,J,K),1, * ARRAY2(LOOP,L,K),1) 520 CONTINUE * 530 CONTINUE * INCR = INCR + NRWBOT C C C IN IDENTITY IN ARRAY2(:,:,K) PERFORM NRWTOP COLUMN C ELIMINATIONS WITH COLUMN PIVOTING. C C DO 740 I = (NRWBOT+1),NOVRLP IPLUS1 = I+1 IMINN = I-NRWBOT JMINN = I+NTPA1 C C DETERMINE COLUMN PIVOT WITH PIVOT INDEX. C C JPVT = JMINN INCRN = INCR + IMINN JPVT = ISAMAX(NRWBOT+1,ARRAY2(I,JMINN,K),NOVRLP)+JMINN-1 C C TEST FOR SINGULARITY. C IF(PIVMAX+COLMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(COLMAX,PIVMAX) C C INTERCHANGE COLUMNS IF NECESSARY. C PIVOT(INCRN) = INCRN+JPVT-JMINN IF(JPVT.EQ.JMINN) GO TO 660 C C INTERCHANGE COLUMNS IN IDENTITY IN ARRAY2(:,:,K). C CALL SSWAP(NOVRLP-I+1,ARRAY2(I,JPVT,K),1, * ARRAY2(I,JMINN,K),1) C C INTERCHANGE COLUMNS IN ARRAY1(:,:,K+1) C JPVT2 = JPVT-JMINN+IMINN * IF(K.LT.NBLOKS)THEN * CALL SSWAP(NRWBLK,ARRAY1(1,JPVT2,KP1),1, * ARRAY1(1,IMINN,KP1),1) * C C INTERCHANGE COLUMNS IN ARRAY2(:,:,K+1). C CALL SSWAP(NOVRLP,ARRAY2(1,JPVT2,KP1),1, * ARRAY2(1,IMINN,KP1),1) ELSE C C INTERCHANGE COLUMNS IN BOTBLK IF K = NBLOKS. C CALL SSWAP(NRWBOT,BOTBLK(1,JPVT2),1,BOTBLK(1,IMINN),1) * ENDIF * 660 CONTINUE C C COMPUTE COLUMN MULTIPLIERS AND PERFORM COLUMN ELIMINATIONS. C COLPIV = ONE/ARRAY2(I,JMINN,K) CALL SSCAL(NRWBOT,COLPIV,ARRAY2(I,JMINN+1,K),NOVRLP) DO 730 J=(JMINN+1),(JMINN+NRWBOT) COLMLT = ARRAY2(I,J,K) CALL SAXPY(NOVRLP-I,-COLMLT,ARRAY2(IPLUS1,JMINN,K),1, * ARRAY2(IPLUS1,J,K),1) C C ADJUST COLUMNS IN ARRAY1(:,:,K+1), ARRAY2(:,:,K+1) C OR BOTBLK. C JMINN2 = J-JMINN+IMINN JMINN3 = JMINN-NCOLS IF(K.LT.NBLOKS)THEN C C FIRST ADJUST COLUMNS IN ARRAY1(:,:,K+1) C CALL SAXPY(NRWBLK,-COLMLT,ARRAY1(1,JMINN3,KP1),1, * ARRAY1(1,JMINN2,KP1),1) * C C ADJUST COLUMNS IN ARRAY2(:,:,K+1) C CALL SAXPY(NOVRLP,-COLMLT,ARRAY2(1,JMINN3,KP1),1, * ARRAY2(1,JMINN2,KP1),1) * ELSE C C ADJUST COLUMNS IN BOTBLK. C CALL SAXPY(NRWBOT,-COLMLT,BOTBLK(1,JMINN3),1, * BOTBLK(1,JMINN2),1) * ENDIF 730 CONTINUE * 740 CONTINUE * INCR = INCR + NRWTOP * 745 CONTINUE * C C C C PERFORM (NOVRLP-NRWTOP-1)=NRWBOT-1 ROW ELIMINATION(S) WITH ROW C PIVOTING ON BOTBLK. C C C IF(NRWBOT.EQ.1) GO TO 910 * DO 900 J = (NRWTOP+1),(NOVRLP-1) JPLUS1 = J+1 JMINN = J-NRWTOP INCRJ = INCR + JMINN LOOP = JMINN+1 IPVT = ISAMAX(NRWBOT-JMINN+1,BOTBLK(JMINN,J),1) + JMINN-1 ROWMAX = ABS(BOTBLK(IPVT,J)) C C TEST FOR SINGULARITY. C IF(PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 920 PIVMAX = MAX(ROWMAX,PIVMAX) C C IF NECESSARY, INTERCHANGE ROWS IN BOTBLK. C PIVOT(INCRJ) = INCR + IPVT IF(IPVT.EQ.JMINN) GO TO 770 CALL SSWAP(NOVRLP-J+1,BOTBLK(IPVT,J),NRWBOT, * BOTBLK(JMINN,J),NRWBOT) * 770 CONTINUE C C COMPUTE MULTIPLIERS. C ROWPIV = ONE/BOTBLK(JMINN,J) CALL SSCAL(NRWBOT-JMINN,ROWPIV,BOTBLK(LOOP,J),1) C C PERFORM ROW ELIMINATIONS. C DO 800 L = (J+1),NOVRLP ROWMLT = BOTBLK(JMINN,L) CALL SAXPY(NRWBOT-JMINN,-ROWMLT,BOTBLK(LOOP,J),1, * BOTBLK(LOOP,L),1) 800 CONTINUE * 900 CONTINUE * 910 CONTINUE C C DONE PROVIDED LAST ELEMENT IS NOT ZERO. C IF(PIVMAX+ABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN C C C MATRIX IS SINGULAR -- SET IFLAG = -1 AND TERMINATE. C C 920 CONTINUE IFLAG = -1 RETURN C C C ERROR IN THE PARAMETER N SPECIFIED. C C 930 CONTINUE IFLAG = 1 RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABDDEC. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END * SUBROUTINE ABDSOL(N, TOPBLK, NRWTOP, NOVRLP, ARRAY1, NRWBLK, * NBLOKS, ARRAY2, BOTBLK, PIVOT, B, X) * CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C GIVEN THE DECOMPOSED FORM OF THE MATRIX A WITH THE PIVOT VECTOR C (GENERATED IN ABDDEC) AND THE VECTOR B, THIS ROUTINE SOLVES C FOR X IN AX = B. C C C THIS IS THE SINGLE PRECISION VERSION OF ABDSOL. THE DOUBLE C PRECISION VERSION IS OBTAINED BY REPLACING ALL SINGLE PRECISION C DECLARATIONS BY DOUBLE PRECISION, AND REPLACING THE CALLS C SINGLE PRECISION BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) BY C DOUBLE PRECISION CALLS. C C *** PARAMETERS *** C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM GIVEN BY C NBLOKS * (NRWBLK + NOVRLP) + NOVRLP. C UNCHANGED ON EXIT. C C TOPBLK - REAL (NRWTOP, NOVRLP) C THE DECOMPOSED FIRST BLOCK OF A. C UNCHANGED ON EXIT. C C NRWTOP - INTEGER C THE NUMBER OF ROWS IN TOPBLK. C UNCHANGED ON EXIT. C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN TOPBLK, AND THE C NUMBER OF ROWS IN ARRAY2(:,:,K). C UNCHANGED ON EXIT. C C ARRAY1 - REAL (NRWBLK, NRWBLK+NOVRLP,NBLOKS) C ARRAY1(:,:,K) CONTAINS THE K-TH NRWBLK C BY NRWBLK+NOVRLP DECOMPOSED BLOCK OF C THE MATRIX A. C UNCHANGED ON EXIT. C C NRWBLK - INTEGER C THE NUMBER OF ROWS IN ARRAY1(:,:,K). C UNCHANGED ON EXIT. C C NBLOKS - INTEGER C THE NUMBER OF BLOCKS IN ARRAY1 AND ARRAY2. C UNCHANGED ON EXIT. C C ARRAY2 - REAL(NOVRLP,NOVRLP+NRWBLK+NOVRLP,NBLOK) C ARRAY2(:,:,K) CONTAINS THE K-TH NOVRLP C BY NOVRLP+NRWBLK+NOVRLP DECOMPOSED BLOCK OF C THE MATRIX A. C UNCHANGED ON EXIT. C C BOTBLK - REAL(NOVRLP-NRWTOP, NOVRLP) C THE LAST DECOMPOSED BLOCK OF THE MATRIX A. C UNCHANGED ON EXIT. C C PIVOT - INTEGER(N) C CONTAINS THE PIVOTING INDICES DETERMINED C IN THE DECOMPOSITION ROUTINE ABDDEC. C UNCHANGED ON EXIT. C C B - REAL(N) C THE RIGHT HAND SIDE VECTOR. C UNCHANGED ON EXIT. C C X - REAL(N) C WORK SPACE FOR THE SOLUTION OF THE LINEAR SYSTEM. C IF IFLAG = 0, CONTAINS THE SOLUTION OF AX = B ON C EXIT. C C C ROUTINES CALLED : C C THE BASIC LINEAR ALGEBRA SUBROUTINES (BLAS) : C C SAXPY, SDOT. C C TO CONVERT ABDSOL TO DOUBLE PRECISION, THESE MUST BE C REPLACED BY THEIR DOUBLE PRECISION EQUIVALENTS : C C DAXPY, DDOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DECLARATIONS : C C PARAMETERS : C INTEGER N,NBLOKS,NRWTOP,NOVRLP,NRWBLK,PIVOT(N) * REAL TOPBLK,ARRAY1,ARRAY2,BOTBLK,B,X * DIMENSION TOPBLK(NRWTOP,*),ARRAY1(NRWBLK,NRWBLK+NOVRLP,*), * ARRAY2(NOVRLP,NOVRLP+NRWBLK+NOVRLP,*), * BOTBLK(NOVRLP-NRWTOP,*),B(N),X(N) C C LOCAL VARIABLES : C INTEGER NCOLS,NRWBOT,L,NCLBLK,NTPJ,I,J,K,LOOP,INCR,INCRJ, * INCRI,INCR1,INCRTP,JPIVOT,ITEMP,IPVTI,L1,INCRN,LL, * NRWBTL,IPLUSN,IPVTN,NTPA1,NTP1,NTPCL * REAL XINCRI,SWAP,SDOT C C DEFINE THE CONSTANTS USED THROUGHOUT. C * NCOLS = NOVRLP+NRWBLK NRWBOT = NOVRLP-NRWTOP NCLBLK = NCOLS+NOVRLP NTPA1 = NRWTOP+NRWBLK NTP1 = NRWTOP+1 NTPCL = NRWTOP+NCOLS DO 100 J = 1,N X(J) = B(J) 100 CONTINUE * C C C C FORWARD RECURSION... C C C C FIRST IN TOPBLK... C C FORWARD SOLUTION. C C C DO 130 J = 1,NRWTOP X(J) = X(J)/TOPBLK(J,J) LOOP = J+1 CALL SAXPY(NRWTOP-J,-X(J),TOPBLK(LOOP,J),1,X(LOOP),1) 130 CONTINUE * INCR = NRWTOP C C C C IN EACH BLOCK K, K = 1,NBLOKS... C C C DO 360 K = 1,NBLOKS C C FORWARD MODIFICATION. C INCRTP = INCR DO 140 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL SAXPY(NRWBLK,-X(INCRJ),ARRAY1(1,J,K),1,X(INCRTP+1),1) 140 CONTINUE C C FORWARD SOLUTION IN ARRAY1(:,:,K). C DO 170 J = 1,NRWBLK INCRJ = INCR+J NTPJ = NRWTOP+J X(INCRJ) = X(INCRJ)/ARRAY1(J,NTPJ,K) LOOP = J+1 CALL SAXPY(NRWBLK-J,-X(INCRJ),ARRAY1(LOOP,NTPJ,K),1, * X(INCRJ+1),1) 170 CONTINUE * C C FORWARD MODIFICATION. C INCRTP = INCR+NRWBLK * DO 240 J = 1,NTPA1 INCRJ = INCR-NRWTOP+J INCR1 = INCRTP+1 CALL SAXPY(NOVRLP,-X(INCRJ),ARRAY2(1,J,K),1,X(INCR1),1) 240 CONTINUE C C FORWARD ELIMINATION. C C DO 260 J = (NTPA1+1),NCOLS INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.EQ.INCRJ) GO TO 245 SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP 245 CONTINUE LOOP = J-(NTPA1-1) CALL SAXPY(NOVRLP-LOOP+1,-X(INCRJ),ARRAY2(LOOP,J,K),1, * X(INCRTP+LOOP),1) 260 CONTINUE C C FORWARD SOLUTION. C IF(NRWTOP.LT.1) GO TO 285 C DO 280 J = (NCOLS+1),NTPCL INCRJ = INCR-NRWTOP+J ITEMP = J-NTPA1 X(INCRJ) = X(INCRJ)/ARRAY2(ITEMP,J,K) LOOP = ITEMP+1 CALL SAXPY(NOVRLP-ITEMP,-X(INCRJ),ARRAY2(LOOP,J,K),1, * X(LOOP+INCRTP),1) 280 CONTINUE 285 CONTINUE INCR = INCR+NCOLS 360 CONTINUE * C C IN BOTBLK ... C C FORWARD MODIFICATION. C INCRTP = INCR DO 300 J = 1,NRWTOP INCRJ = INCR-NRWTOP+J CALL SAXPY(NRWBOT,-X(INCRJ),BOTBLK(1,J),1,X(INCRTP+1),1) 300 CONTINUE C C FORWARD ELIMINATION. C IF(NRWBOT.EQ.1) GO TO 350 DO 340 J = (NRWTOP+1),(NOVRLP-1) INCRJ = INCR-NRWTOP+J JPIVOT = PIVOT(INCRJ) IF(JPIVOT.EQ.INCRJ) GO TO 325 SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP 325 CONTINUE LOOP = J - (NRWTOP-1) CALL SAXPY(NRWBOT-LOOP+1,-X(INCRJ),BOTBLK(LOOP,J),1, * X(INCRTP+LOOP),1) 340 CONTINUE 350 CONTINUE C C C C BACKWARD RECURSION... C C C C FIRST IN BOTBLK... C C BACKWARD SOLUTION. C C DO 430 LL = 1,NRWBOT J = NOVRLP+1-LL NRWBTL = NRWBOT + 1 - LL INCRJ = INCR+NRWBTL X(INCRJ) = X(INCRJ)/BOTBLK(NRWBTL,J) LOOP = NRWBOT-LL CALL SAXPY(LOOP,-X(INCRJ),BOTBLK(1,J),1,X(INCR+1),1) 430 CONTINUE * C C IN EACH BLOCK GOING IN REVERSE ORDER... C DO 520 L = 1,NBLOKS C C BACKWARD ELIMINATION. C K = NBLOKS+1-L INCR = INCR-NRWTOP C C BACKWARD ELIMINATION IN THE IDENTITY IN ARRAY2. C DO 470 L1 = 1,NRWTOP I = NTP1-L1 ITEMP = NRWBOT+I IPLUSN = NCOLS+I INCRN = INCR+I X(INCRN) = X(INCRN) - SDOT(NRWBOT, * ARRAY2(ITEMP,IPLUSN+1,K),NOVRLP,X(INCRN+1),1) IPVTN = PIVOT(INCRN) IF(INCRN.EQ.IPVTN) GO TO 465 SWAP = X(INCRN) X(INCRN) = X(IPVTN) X(IPVTN) = SWAP 465 CONTINUE 470 CONTINUE INCR = INCR-NRWBOT C C BACKWARD MODIFICATION, SOLUTION AND ADJUSTMENT (ORDERING) OF C THE X'S. C DO 495 L1 = 1,NRWBOT I = NRWBOT+1-L1 IPLUSN = NTPA1+I LOOP = IPLUSN + 1 INCRI = INCR + I C C BACKWARD MODIFICATION. C XINCRI = X(INCRI) X(INCRI) = X(INCRI) - SDOT(NRWBOT-1,ARRAY2(I,LOOP,K), * NOVRLP,X(INCR+I+1),1) X(INCRI) = X(INCRI)-X(INCRI+NRWBOT) C C BACKWARD SOLUTION. C X(INCRI) = X(INCRI)/ARRAY2(I,IPLUSN,K) C C COLUMN PIVOTING ADJUSTMENT ON THE X'S. C IPVTN = PIVOT(INCRI) IF(INCRI.EQ.IPVTN) GO TO 490 INCRI = INCRI+NRWBOT IPVTN = IPVTN+NRWBOT SWAP = X(INCRI) X(INCRI) = X(IPVTN) X(IPVTN) = SWAP 490 CONTINUE 495 CONTINUE INCR = INCR-NRWBLK DO 450 L1 = 1,NRWBLK I = NRWBLK+1-L1 IPLUSN = NRWTOP+I INCRN = INCR+I X(INCRN) = X(INCRN)-SDOT(NCOLS-IPLUSN, * ARRAY1(I,IPLUSN+1,K),NRWBLK,X(INCR+I+1),1) IPVTN = PIVOT(INCRN) IF(INCRN.EQ.IPVTN) GO TO 445 SWAP = X(INCRN) X(INCRN) = X(IPVTN) X(IPVTN) = SWAP 445 CONTINUE 450 CONTINUE 520 CONTINUE C C BACKWARD ELIMINATION IN TOPBLK. C DO 510 L1 = 1,NRWTOP I = NTP1-L1 LOOP = I + 1 X(I) = X(I) - SDOT(NOVRLP-I,TOPBLK(I,LOOP),NRWTOP, * X(LOOP),1) IPVTI = PIVOT(I) IF(I.EQ.IPVTI) GO TO 505 SWAP = X(I) X(I) = X(IPVTI) X(IPVTI) = SWAP 505 CONTINUE 510 CONTINUE RETURN C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABDSOL. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C END