bchfac Subroutine

The factorization has the form

  C = L * D * L'

with L unit lower triangular and D diagonal, for a given matrix C of 
order NROW, where C is symmetric positive semidefinite and banded, 
having NBANDS diagonals at and below the main diagonal.

Gauss elimination is used, adapted to the symmetry and bandedness of C.

Near-zero pivots are handled in a special way.  The diagonal 
element C(N,N) = W(1,N) is saved initially in DIAG(N), all N.

At the N-th elimination step, the current pivot element, W(1,N), 
is compared with its original value, DIAG(N).  If, as the result 
of prior elimination steps, this element has been reduced by about 
a word length, that is, if W(1,N) + DIAG(N) <= DIAG(N), then the pivot 
is declared to be zero, and the entire N-th row is declared to
be linearly dependent on the preceding rows.  This has the effect 
of producing X(N) = 0 when solving C * X = B for X, regardless of B.

Justification for this is as follows.  In contemplated applications 
of this program, the given equations are the normal equations for 
some least-squares approximation problem, DIAG(N) = C(N,N) gives 
the norm-square of the N-th basis function, and, at this point, 
W(1,N) contains the norm-square of the error in the least-squares 
approximation to the N-th basis function by linear combinations 
of the first N-1.

Having W(1,N)+DIAG(N) <= DIAG(N) signifies that the N-th function 
is linearly dependent to machine accuracy on the first N-1 
functions, therefore can safely be left out from the basis of 
approximating functions.

The solution of a linear system C * X = B is effected by the 
succession of the following two calls:

  call bchfac ( w, nbands, nrow, diag )

  call bchslv ( w, nbands, nrow, b, x )

14 February 2007

Carl DeBoor

Carl DeBoor,
A Practical Guide to Splines,
Springer, 2001,
ISBN: 0387953663,
LC: QA1.A647.v27.

Input/output, real ( kind = 8 ) W(NBANDS,NROW).
On input, W contains the NBANDS diagonals in its rows, 
with the main diagonal in row 1.  Precisely, W(I,J) 
contains C(I+J-1,J), I=1,...,NBANDS, J=1,...,NROW.
For example, the interesting entries of a seven diagonal
symmetric matrix C of order 9 would be stored in W as

  11 22 33 44 55 66 77 88 99
  21 32 43 54 65 76 87 98  *
  31 42 53 64 75 86 97  *  *
  41 52 63 74 85 96  *  *  *

Entries of the array not associated with an
entry of C are never referenced.
On output, W contains the Cholesky factorization 
C = L*D*L', with W(1,I) containing 1/D(I,I) and W(I,J) 
containing L(I-1+J,J), I=2,...,NBANDS.

Input, integer ( kind = 4 ) NBANDS, indicates the bandwidth of the
matrix C, that is, C(I,J) = 0 for NBANDS < abs(I-J).

Input, integer ( kind = 4 ) NROW, is the order of the matrix C.

Work array, real ( kind = 8 ) DIAG(NROW).