C This file contains: C C subroutine ABDCND: the condition number estimator. C It calls ABDNRM, ABDDEC, DONEST and C ABDSOL. C C function ABDNRM: to calculate the L1 norm of A. C C subroutine DONEST: a double precision version of the routine C SONEST by N.J. Higham. C C subroutine ABDLEQ: the factor/solve routine from AABDPACK. C It calls ABDDEC and ABDSOL. C C subroutine ABDDEC: the factor routine. C C subroutine ABDSOL: the solve routine. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C END OF SUBROUTINE ABDSOL. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE ABDCND(N,TOPBLK,NRWTOP,NOVRLP,ARRAY1,NRWBLK,NCOLS, * ARRAY2, NBLOKS,BOTBLK,NRWBOT,EST,V,ISIGN, * WORK,PIVOT) C C*************************************************************** C C THIS PROGRAM COMPUTES THE CONDITION NUMBER OF AN ALMOST BLOCK C DIAGONAL MATRIX OF THE FORM: C C TOPBLK C ARRAY1(1) C ARRAY2(1) C ARRAY1(2) C ARRAY2(2) C . C . C . C C ARRAY1(NBLOKS-1) C ARRAY2(NBLOKS-1) C ARRAY1(NBLOKS) C BOTBLK C C WHERE C TOPBLK IS NRWTOP BY NOVRLP C ARRAY1(K), K=1,NBLOKS, ARE NRWBLK BY NRWBLK+NOVRLP C ARRAY2(K), K=1,NBLOKS-1, ARE NOVRLP BY C 2*NOVRLP+NRWBLK, AND IS OF THE FORM C A(NOVRLPxNOVRLP)B(NOVRLPxNRWBLK)I(NOVRLPxNOVRLP) C WHERE I IS THE IDENTITY, C BOTBLK IS NRWBOT BY NRWBLK+NOVRLP, C AND C NOVRLP = NRWTOP + NRWBOT C WITH C NOVRLP.LE.NRWBLK . C C THE LINEAR SYSTEM IS OF ORDER N = NBLOKS*(NRWBLK + NOVRLP). C C THE METHOD IMPLEMENTED IS BASED ON GAUSS ELIMINATION WITH C ALTERNATE ROW AND COLUMN ELIMINATION WITH PARTIAL PIVOTING, C WHICH PRODUCES A STABLE DECOMPOSITION OF THE MATRIX A C WITHOUT INTRODUCING FILL-IN. SEE ABDPACK DOCUMENTATION FOR C SAMPLE DRIVING PROGRAM ETC. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM, C GIVEN BY NBLOKS*NRWBLK + NOVRLP C C TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP) C THE FIRST BLOCK OF THE ALMOST BLOCK C DIAGONAL MATRIX A C C NRWTOP - INTEGER C NUMBER OF ROWS IN THE BLOCK TOPBLK C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN WHICH SUCC- C ESSIVE BLOCKS OVERLAP, WHERE C NOVRLP = NRWTOP + NRWBOT C ALSO THE NUMBER OF ROWS IN ARRAY2 C C ARRAY1 - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS) C ARRAY1(,,K) CONTAINS THE K-TH NRWBLK C BY NCOLS BLOCK OF THE MATRIX A C C NRWBLK - INTEGER C NUMBER OF ROWS IN K-TH ARRAY1 C C NCOLS - INTEGER C NUMBER OF COLUMNS IN K-TH ARRAY1 C C ARRAY2 - DOUBLE PRECISION(NOVRLP,NCLBLK+NOVRLP, C NBLOKS-1). ARRAY2(,,K) CONTAINS THE C K-TH NOVRLP BY NCOLS+NOVRLP BLOCK OF C THE MATRIX A C C NBLOKS - INTEGER C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN C THE MATRIX A C C BOTBLK - DOUBLE PRECISION(NRWBOT,NCOLS) C THE LAST BLOCK OF THE MATRIX A C C NRWBOT - INTEGER C NUMBER OF ROWS IN THE BLOCK BOTBLK C C ISIGN - INTEGER(N) C WORK SPACE C C V - DOUBLE PRECISION (N), = INV(A)*W, WHERE C NORM(V)/NORM(W) ESTIMATES NORM OF INV(A). C C WORK - DOUBLE PRECISION(N), WORK SPACE. C C *** ON RETURN ... C C EST - ESTIMATED CONDITION NUMBER. C V - APPROXIMATE NULL VECTOR IF A IS SINGULAR. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** AUXILIARY PROGRAMS ***** C C ABDDEC(N,TOPBLK,NRWTOP,NOVRLP,ARRAY1,NRWBLK,NBLOKS,ARRAY2, C * BOTBLK,PIVOT,IFLAG) C - DECOMPOSES THE MATRIX A USING MODIFIED C ALTERNATE ROW AND COLUMN ELIMINATON WITH C PARTIAL PIVOTING, AND IS USED FOR THIS C PURPOSE IN A B D P A C K. C THE ARGUMENTS ARE AS IN A B D P A C K. C C C ABDSOL(N,TOPBLK,NRWTOP,NOVRLP,ARRAY1,NRWBLK,NBLOKS,ARRAY2, C * BOTBLK,PIVOT,B,X,JOB) C - SOLVES THE SYSTEM A*X = B ONCE A IS DECOMPOSED. C THE ARGUMENTS ARE ALL AS IN A B D P A C K. C C ABDNRM(NBLOKS,NRWTOP,NRWBOT,NOVRLP,NRWBLK,NCOLS,TOPBLK, C * ARRAY1,ARRAY2,BOTBLK) C COMPUTES THE NORM OF THE ABD MATRIX TOP, A, BOT. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK INTEGER N,NRWTOP,NOVRLP,NRWBLK,NCOLS,NBLOKS,NRWBOT,PIVOT(*), * IFLAG,ISIGN(*) DOUBLE PRECISION TOPBLK(NRWTOP,*),ARRAY1(NRWBLK,NCOLS,*), * ARRAY2(NOVRLP,NCOLS+NOVRLP,*),BOTBLK(NRWBOT,*), * EST,V(*),WORK(*) DOUBLE PRECISION ABDNRM,NORM INTEGER KASE,NOUT INTEGER COUNT COMMON /CNDCNT/ COUNT COUNT = 0 NOUT = 6 C C FIRST, COMPUTE THE NORM OF THE MATRIX ITSELF: C NORM = ABDNRM(NBLOKS,NRWTOP,NRWBOT,NOVRLP,NRWBLK,NCOLS, * TOPBLK,ARRAY1,ARRAY2,BOTBLK) C C THEN, DO THE FACTORIZATION USING CRDCMP: C CALL ABDDEC(N,TOPBLK,NRWTOP,NOVRLP,ARRAY1,NRWBLK,NBLOKS, * ARRAY2,BOTBLK,PIVOT,IFLAG) C C AND THEN DO THE CONDITION NUMBER ESTIMATION: C KASE = 0 55 CALL DONEST(N,V,WORK,ISIGN,EST,KASE) IF (KASE .NE. 0) THEN COUNT = COUNT+1 CALL ABDSOL(N,TOPBLK,NRWTOP,NOVRLP,ARRAY1,NRWBLK,NBLOKS, * ARRAY2,BOTBLK,PIVOT,WORK,WORK,KASE-1) GOTO 55 END IF C WRITE (NOUT,*)EST,NORM,EST*NORM EST = EST*NORM RETURN END C C================================================================== C DOUBLE PRECISION FUNCTION ABDNRM(NBLOKS,NTOP,NBOT,NOVRLP, * NRWBLK,NCOLS,TOP,A1,A2,BOT) C C ABDNRM IS USED IN CONJUNCTION WITH DONEST TO COMPUTE THE C CONDITION NUMBER OF AN ALMOST BLOCK DIAGONAL MATRIX LIKE C THE ONES HANDLED BY ABDPACK. [SEE COMMENTS IN ABDPACK, ABDDEC] INTEGER NBLOKS,NTOP,NBOT,NOVRLP,NRWBLK,NCOLS DOUBLE PRECISION TOP(NTOP,*),A1(NRWBLK,NCOLS,*), * A2(NOVRLP,NCOLS+NOVRLP,*),BOT(NBOT,*) INTEGER J,K DOUBLE PRECISION NCOLS2 DOUBLE PRECISION MAX,DASUM,TEMP TEMP = 0.0D0 NCOLS2 = NCOLS+NOVRLP C C FIRST, GO OVER THE COLUMNS OF TOP AND THE FIRST BLOCK: C DO 10 J = 1,NOVRLP TEMP = MAX(TEMP, DASUM(NTOP,TOP(1,J),1) + * DASUM(NRWBLK,A1(1,J,1),1) + * DASUM(NOVRLP,A2(1,J,1),1)) 10 CONTINUE DO 40 K = 1,NBLOKS-1 C C IN EACH BLOCK: C C FIRST, THE COLUMNS FROM THE KTH BLOCK ALONE: C DO 20 J = NOVRLP+1,NCOLS TEMP = MAX(TEMP, DASUM(NRWBLK,A1(1,J,K),1) + * DASUM(NOVRLP,A2(1,J,K),1)) 20 CONTINUE C C NOW, TH COLUMNS WHICH INTERSECT BOTH THE KTH AND C (K+1)ST BLOCKS. C DO 30 J = NCOLS+1,NCOLS2 TEMP = MAX(TEMP, DASUM(NOVRLP,A2(1,J,K),1) + * DASUM(NRWBLK,A1(1,J-NCOLS,K+1),1) + * DASUM(NOVRLP,A2(1,J-NCOLS,K+1),1)) 30 CONTINUE 40 CONTINUE C C FINALLY, THE COLUMNS OF THE 2ND LAST BLOCK WHICH DO NOT OVERLAP C WITH ANYTHING. C DO 50 J = NOVRLP+1,NCOLS TEMP = MAX(TEMP, DASUM(NRWBLK,A1(1,J,NBLOKS),1) + * DASUM(NOVRLP,A2(1,J,NBLOKS),1)) 50 CONTINUE C C AND THOSE COLUMNS OVERLAPPING WITH BOTH THE LAST BLOCK OF A2 C AND BOT. C DO 60 J = NCOLS+1,NCOLS2 TEMP = MAX(TEMP, DASUM(NOVRLP,A2(1,J,NBLOKS),1) + * DASUM(NBOT,BOT(1,J-NCOLS),1)) 60 CONTINUE ABDNRM = TEMP RETURN END SUBROUTINE DONEST (N, V, X, ISGN, EST, KASE) INTEGER N, ISGN(N), KASE DOUBLE PRECISION V(N), X(N), EST C C DONEST ESTIMATES THE 1-NORM OF A SQUARE, DOUBLE PRECISION MATRIX A. C REVERSE COMMUNICATION IS USED FOR EVALUATING C MATRIX-VECTOR PRODUCTS. C C ON ENTRY C C N INTEGER C THE ORDER OF THE MATRIX. N .GE. 1. C C ISGN INTEGER(N) C USED AS WORKSPACE. C C KASE INTEGER C = 0. C C ON INTERMEDIATE RETURNS C C KASE = 1 OR 2. C C X DOUBLE PRECISION(N) C MUST BE OVERWRITTEN BY C C A*X, IF KASE=1, C TRANSPOSE(A)*X, IF KASE=2, C C AND DONEST MUST BE RE-CALLED, WITH ALL THE OTHER C PARAMETERS UNCHANGED. C C ON FINAL RETURN C C KASE = 0. C C EST DOUBLE PRECISION C CONTAINS AN ESTIMATE (A LOWER BOUND) FOR NORM(A). C C V DOUBLE PRECISION(N) C = A*W, WHERE EST = NORM(V)/NORM(W) C (W IS NOT RETURNED). C C THIS VERSION DATED MARCH 16, 1988. C NICK HIGHAM, UNIVERSITY OF MANCHESTER. C C MODIFIED FOR DOUBLE PRECISION ON JUNE 11, 1996. C C REFERENCE C N.J. HIGHAM (1987) FORTRAN CODES FOR ESTIMATING C THE ONE-NORM OF A REAL OR COMPLEX MATRIX, WITH APPLICATIONS C TO CONDITION ESTIMATION, NUMERICAL ANALYSIS REPORT NO. 135, C UNIVERSITY OF MANCHESTER, MANCHESTER M13 9PL, ENGLAND. C C SUBROUTINES AND FUNCTIONS C BLAS IDAMAX, DASUM, DCOPY C GENERIC ABS, NINT, FLOAT, SIGN C INTRINSIC FLOAT DOUBLE PRECISION FLOAT INTRINSIC ABS DOUBLE PRECISION ABS INTRINSIC SIGN DOUBLE PRECISION SIGN DOUBLE PRECISION DASUM INTEGER IDAMAX INTEGER ITMAX PARAMETER (ITMAX = 5) DOUBLE PRECISION ZERO,ONE,TWO PARAMETER (ZERO = 0.0D0) PARAMETER (ONE = 1.0D0) PARAMETER (TWO = 2.0D0) C C INTERNAL VARIABLES INTEGER I, ITER, J, JLAST, JUMP DOUBLE PRECISION ALTSGN, ESTOLD, TEMP C SAVE C IF (KASE .EQ. 0) THEN DO 10,I = 1,N X(I) = ONE/FLOAT(N) 10 CONTINUE KASE = 1 JUMP = 1 RETURN ENDIF C GOTO (100, 200, 300, 400, 500) JUMP C C ................ ENTRY (JUMP = 1) C FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. C 100 CONTINUE IF (N .EQ. 1) THEN V(1) = X(1) EST = ABS(V(1)) C ... QUIT GOTO 510 ENDIF EST = DASUM(N,X,1) C DO 110,I = 1,N X(I) = SIGN(ONE,X(I)) ISGN(I) = NINT(X(I)) 110 CONTINUE KASE = 2 JUMP = 2 RETURN C C ................ ENTRY (JUMP = 2) C FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. C 200 CONTINUE J = IDAMAX(N,X,1) ITER = 2 C C MAIN LOOP - ITERATIONS 2,3,...,ITMAX. C 220 CONTINUE DO 230,I = 1,N X(I) = ZERO 230 CONTINUE X(J) = ONE KASE = 1 JUMP = 3 RETURN C C ................ ENTRY (JUMP = 3) C X HAS BEEN OVERWRITTEN BY A*X. C 300 CONTINUE CALL DCOPY(N,X,1,V,1) ESTOLD = EST EST = DASUM(N,V,1) DO 310,I = 1,N IF ( NINT( SIGN(ONE,X(I)) ) .NE. ISGN(I) ) GOTO 320 310 CONTINUE C REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. GOTO 410 C 320 CONTINUE C TEST FOR CYCLING. IF (EST .LE. ESTOLD) GOTO 410 C DO 330,I = 1,N X(I) = SIGN(ONE,X(I)) ISGN(I) = NINT(X(I)) 330 CONTINUE KASE = 2 JUMP = 4 RETURN C C ................ ENTRY (JUMP = 4) C X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. C 400 CONTINUE JLAST = J J = IDAMAX(N,X,1) IF ( ( X(JLAST) .NE. ABS(X(J)) ) .AND. + (ITER .LT. ITMAX) ) THEN ITER = ITER + 1 GOTO 220 ENDIF C C ITERATION COMPLETE. FINAL STAGE. C 410 CONTINUE ALTSGN = ONE DO 420,I = 1,N X(I) = ALTSGN * (ONE + FLOAT(I-1)/FLOAT(N-1)) ALTSGN = -ALTSGN 420 CONTINUE KASE = 1 JUMP = 5 RETURN C C ................ ENTRY (JUMP = 5) C X HAS BEEN OVERWRITTEN BY A*X. C 500 CONTINUE TEMP = TWO*DASUM(N,X,1)/FLOAT(3*N) IF (TEMP. GT. EST) THEN CALL DCOPY(N,X,1,V,1) EST = TEMP ENDIF C 510 KASE = 0 RETURN C END 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,JOB) 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 C IF(NRWTOP.EQ.0) GO TO 930 C 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, JOB) * 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, FOR JOB = 0, AND TRANSPOSE(A)X = B FOR JOB = 1. 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 JOB - INTEGER C IF JOB = 0, THEN SOLV REGULAR MATRIX. C OTHERWISE, SOLVE THE TRANSPOSE OF THE MATRIX 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),JOB * 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 DETERMINE THE JOB TO BE DONE. C IF (JOB .NE. 0) GO TO 530 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 C IF JOB IS NON-ZERO, SOLVE TRANSPOSE(A)*X = B: C 530 CONTINUE C FORWARD ELIMINATION IN TOP BLOCK TOPBLK USING TRANSPOSE DO 540 I = 1,NRWTOP IPVTI = PIVOT(I) IF(I.NE.IPVTI) THEN SWAP = X(I) X(I) = X(IPVTI) X(IPVTI) = SWAP ENDIF LOOP = I+1 CALL DAXPY(NOVRLP-I,-X(I),TOPBLK(I,LOOP),NRWTOP,X(LOOP),1) 540 CONTINUE INCR = NRWTOP C IN EACH BLOCK, K = 1,..,NBLOCKS DO 620 K = 1,NBLOKS C FORWARD ELIMINATION USING TRANSPOSE(A) DO 560 I = 1,NRWBLK IPLUSN = NRWTOP+I INCRN = INCR+I IPVTN = PIVOT(INCRN) IF(INCRN.EQ.IPVTN) GO TO 545 SWAP = X(INCRN) X(INCRN) = X(IPVTN) X(IPVTN) = SWAP 545 CONTINUE CALL DAXPY(NCOLS-IPLUSN,-X(INCRN),ARRAY1(I,IPLUSN+1,K), * NRWBLK,X(INCRN+1),1) 560 CONTINUE INCR = INCR+NRWBLK C C FORWARD ADJUSTMENT (ORDERING) OF X'S, SOLUTION, AND C MODIFICATION USING TRANSPOSE(A) C DO 580 I = 1,NRWBOT IPLUSN = NTPA1+I LOOP = IPLUSN+1 INCRI = INCR+I C C COLUMN PIVOTING ADJUSTMENT ON THE X'S USING TRANSPOSE(A) C IPVTN = PIVOT(INCRI) IF(INCRI.EQ.IPVTN) GO TO 570 INCRI = INCRI+NRWBOT IPVTN = IPVTN+NRWBOT SWAP = X(INCRI) X(INCRI) = X(IPVTN) X(IPVTN) = SWAP INCRI = INCR+I 570 CONTINUE C C FORWARD SOLUTION C X(INCRI) = X(INCRI)/ARRAY2(I,IPLUSN,K) C C FORWARD MODIFICATION USING TRANSPOSE(A) C XINCRI = -X(INCRI) X(INCRI+NRWBOT) = X(INCRI+NRWBOT)+XINCRI CALL DAXPY(NRWBOT-1,XINCRI,ARRAY2(I,LOOP,K),NOVRLP, * X(INCRI+1),1) 580 CONTINUE INCR = INCR+NRWBOT C C FORWARD ELIMINATION IN THE IDENTITY IN ARRAY2 C DO 600 I = 1,NRWTOP ITEMP = NRWBOT+I IPLUSN = NCOLS+I INCRN = INCR+I IPVTN = PIVOT(INCRN) IF(INCRN.EQ.IPVTN) GO TO 585 SWAP = X(INCRN) X(INCRN) = X(IPVTN) X(IPVTN) = SWAP 585 CONTINUE CALL DAXPY(NRWBOT,-X(INCRN),ARRAY2(ITEMP,IPLUSN+1,K),NOVRLP, * X(INCRN+1),1) 600 CONTINUE INCR = INCR+NRWTOP 620 CONTINUE C C FORWARD SOLUTION IN BOTTOM BLOCK USING TRANSPOSE(A) C DO 650 J = 1,NRWBOT INCRJ = INCR+J LL = NOVRLP-NRWBOT+J LOOP = J-1 X(INCRJ) = X(INCRJ)-DDOT(LOOP,BOTBLK(1,LL),1,X(INCR+1),1) X(INCRJ) = X(INCRJ)/BOTBLK(J,LL) 650 CONTINUE C C C BACKWARD RECURSION FOR TRANSPOSE MATRIX C C C FIRST IN BOTTOM BLOCK BOTBLK C INCRTP = INCR C****** DOUBLE CHECK THIS: ************************************ C INCR = INCR-NRWBLK IF(NRWBOT.EQ.1)GO TO 675 DO 670 LL = (NRWTOP+1),(NOVRLP-1) J = NOVRLP-LL+NRWTOP INCRJ = INCR-NRWTOP+J LOOP = J-(NRWTOP-1) X(INCRJ) = X(INCRJ)-DDOT(NRWBOT-LOOP+1,BOTBLK(LOOP,J),1, * X(INCRTP+LOOP),1) JPIVOT = PIVOT(INCRJ) IF(JPIVOT.EQ.INCRJ) GO TO 665 SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP 665 CONTINUE 670 CONTINUE 675 CONTINUE C C BACKWARD MODIFICATION C DO 690 LL = 1,NRWTOP J = NRWTOP-LL+1 INCRJ = INCR-NRWTOP+J X(INCRJ) = X(INCRJ)-DDOT(NRWBOT,BOTBLK(1,J),1, * X(INCRTP+1),1) 690 CONTINUE C C IN EACH BLOCK GOING IN REVERSE ORDER C C DO 800 L = 1,NBLOKS K = NBLOKS+1-L C C BACKWARD SOLUTION OF ARRAY2 TRANSPOSE C INCR = INCR-NCOLS INCRTP = INCR+NRWBLK IF(NRWTOP.LT.1) GO TO 710 DO 700 LL = (NCOLS+1),NTPCL J = NTPCL+NCOLS-LL+1 INCRJ = INCR-NRWTOP+J ITEMP = J-NTPA1 LOOP = ITEMP+1 X(INCRJ) = X(INCRJ)-DDOT(NOVRLP-ITEMP,ARRAY2(LOOP,J,K), * 1,X(INCRTP+LOOP),1) X(INCRJ) = X(INCRJ)/ARRAY2(ITEMP,J,K) 700 CONTINUE 710 CONTINUE C C BACKWARD ELIMINATION C IF(NTP1.GT.NOVRLP) GO TO 740 DO 730 LL = (NTPA1+1),NCOLS J = NCOLS+NTPA1+1-LL INCRJ = INCR-NRWTOP+J LOOP = J-(NTPA1-1) X(INCRJ) = X(INCRJ)-DDOT(NOVRLP-LOOP+1,ARRAY2(LOOP,J,K), * 1,X(INCRTP+LOOP),1) JPIVOT = PIVOT(INCRJ) IF(JPIVOT.EQ.INCRJ) GO TO 720 SWAP = X(INCRJ) X(INCRJ) = X(JPIVOT) X(JPIVOT) = SWAP 720 CONTINUE 730 CONTINUE 740 CONTINUE C C BACKWARD MODIFICATION C DO 750 LL = 1,NTPA1 J = NTPA1-LL+1 INCRJ = INCR-NRWTOP+J X(INCRJ) = X(INCRJ)-DDOT(NOVRLP,ARRAY2(1,J,K),1, * X(INCRTP+1),1) 750 CONTINUE C BACKWARD SOLUTION IN ARRAY1 DO 760 LL = 1,NRWBLK J = NRWBLK-LL+1 INCRJ = INCR+J NTPJ = NRWTOP+J X(INCRJ) = X(INCRJ)-DDOT(NRWBLK-J,ARRAY1(J+1,NTPJ,K), * 1,X(INCRJ+1),1) X(INCRJ) = X(INCRJ)/ARRAY1(J,NTPJ,K) 760 CONTINUE C BACKWARD MODIFICATION IN ARRAY1 DO 770 LL = 1,NRWTOP J = NRWTOP-LL+1 INCRJ = INCR-NRWTOP+J X(INCRJ) = X(INCRJ)-DDOT(NRWBLK,ARRAY1(1,J,K),1, * X(INCR+1),1) 770 CONTINUE 800 CONTINUE C C BACKWARD SOLUTION IN IN TOP BLOCK C DO 900 LL = 1,NRWTOP J = NRWTOP-LL+1 LOOP = J+1 X(J) = X(J)-DDOT(NRWTOP-J,TOPBLK(LOOP,J),1,X(LOOP),1) X(J) = X(J)/TOPBLK(J,J) 900 CONTINUE RETURN END