*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