*DECK DBNFAC SUBROUTINE DBNFAC (W, NROWW, NROW, NBANDL, NBANDU, IFLAG) C***BEGIN PROLOGUE DBNFAC C***SUBSIDIARY C***PURPOSE Subsidiary to DBINT4 and DBINTK C***LIBRARY SLATEC C***TYPE DOUBLE PRECISION (BNFAC-S, DBNFAC-D) C***AUTHOR (UNKNOWN) C***DESCRIPTION C C DBNFAC is the BANFAC routine from C * A Practical Guide to Splines * by C. de Boor C C DBNFAC is a double precision routine C C Returns in W the LU-factorization (without pivoting) of the banded C matrix A of order NROW with (NBANDL + 1 + NBANDU) bands or diag- C onals in the work array W . C C ***** I N P U T ****** W is double precision C W.....Work array of size (NROWW,NROW) containing the interesting C part of a banded matrix A , with the diagonals or bands of A C stored in the rows of W , while columns of A correspond to C columns of W . This is the storage mode used in LINPACK and C results in efficient innermost loops. C Explicitly, A has NBANDL bands below the diagonal C + 1 (main) diagonal C + NBANDU bands above the diagonal C and thus, with MIDDLE = NBANDU + 1, C A(I+J,J) is in W(I+MIDDLE,J) for I=-NBANDU,...,NBANDL C J=1,...,NROW . C For example, the interesting entries of A (1,2)-banded matrix C of order 9 would appear in the first 1+1+2 = 4 rows of W C as follows. C 13 24 35 46 57 68 79 C 12 23 34 45 56 67 78 89 C 11 22 33 44 55 66 77 88 99 C 21 32 43 54 65 76 87 98 C C All other entries of W not identified in this way with an en- C try of A are never referenced . C NROWW.....Row dimension of the work array W . C must be .GE. NBANDL + 1 + NBANDU . C NBANDL.....Number of bands of A below the main diagonal C NBANDU.....Number of bands of A above the main diagonal . C C ***** O U T P U T ****** W is double precision C IFLAG.....Integer indicating success( = 1) or failure ( = 2) . C If IFLAG = 1, then C W.....contains the LU-factorization of A into a unit lower triangu- C lar matrix L and an upper triangular matrix U (both banded) C and stored in customary fashion over the corresponding entries C of A . This makes it possible to solve any particular linear C system A*X = B for X by a C CALL DBNSLV ( W, NROWW, NROW, NBANDL, NBANDU, B ) C with the solution X contained in B on return . C If IFLAG = 2, then C one of NROW-1, NBANDL,NBANDU failed to be nonnegative, or else C one of the potential pivots was found to be zero indicating C that A does not have an LU-factorization. This implies that C A is singular in case it is totally positive . C C ***** M E T H O D ****** C Gauss elimination W I T H O U T pivoting is used. The routine is C intended for use with matrices A which do not require row inter- C changes during factorization, especially for the T O T A L L Y C P O S I T I V E matrices which occur in spline calculations. C The routine should NOT be used for an arbitrary banded matrix. C C***SEE ALSO DBINT4, DBINTK C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 800901 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900328 Added TYPE section. (WRB) C***END PROLOGUE DBNFAC C INTEGER IFLAG, NBANDL, NBANDU, NROW, NROWW, I, IPK, J, JMAX, K, 1 KMAX, MIDDLE, MIDMK, NROWM1 DOUBLE PRECISION W(NROWW,*), FACTOR, PIVOT C C***FIRST EXECUTABLE STATEMENT DBNFAC IFLAG = 1 MIDDLE = NBANDU + 1 C W(MIDDLE,.) CONTAINS THE MAIN DIAGONAL OF A . NROWM1 = NROW - 1 IF (NROWM1) 120, 110, 10 10 IF (NBANDL.GT.0) GO TO 30 C A IS UPPER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO . DO 20 I=1,NROWM1 IF (W(MIDDLE,I).EQ.0.0D0) GO TO 120 20 CONTINUE GO TO 110 30 IF (NBANDU.GT.0) GO TO 60 C A IS LOWER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO AND C DIVIDE EACH COLUMN BY ITS DIAGONAL . DO 50 I=1,NROWM1 PIVOT = W(MIDDLE,I) IF (PIVOT.EQ.0.0D0) GO TO 120 JMAX = MIN(NBANDL,NROW-I) DO 40 J=1,JMAX W(MIDDLE+J,I) = W(MIDDLE+J,I)/PIVOT 40 CONTINUE 50 CONTINUE RETURN C C A IS NOT JUST A TRIANGULAR MATRIX. CONSTRUCT LU FACTORIZATION 60 DO 100 I=1,NROWM1 C W(MIDDLE,I) IS PIVOT FOR I-TH STEP . PIVOT = W(MIDDLE,I) IF (PIVOT.EQ.0.0D0) GO TO 120 C JMAX IS THE NUMBER OF (NONZERO) ENTRIES IN COLUMN I C BELOW THE DIAGONAL . JMAX = MIN(NBANDL,NROW-I) C DIVIDE EACH ENTRY IN COLUMN I BELOW DIAGONAL BY PIVOT . DO 70 J=1,JMAX W(MIDDLE+J,I) = W(MIDDLE+J,I)/PIVOT 70 CONTINUE C KMAX IS THE NUMBER OF (NONZERO) ENTRIES IN ROW I TO C THE RIGHT OF THE DIAGONAL . KMAX = MIN(NBANDU,NROW-I) C SUBTRACT A(I,I+K)*(I-TH COLUMN) FROM (I+K)-TH COLUMN C (BELOW ROW I ) . DO 90 K=1,KMAX IPK = I + K MIDMK = MIDDLE - K FACTOR = W(MIDMK,IPK) DO 80 J=1,JMAX W(MIDMK+J,IPK) = W(MIDMK+J,IPK) - W(MIDDLE+J,I)*FACTOR 80 CONTINUE 90 CONTINUE 100 CONTINUE C CHECK THE LAST DIAGONAL ENTRY . 110 IF (W(MIDDLE,NROW).NE.0.0D0) RETURN 120 IFLAG = 2 RETURN END