      SUBROUTINE BCHFAC ( W, NBANDS, NROW, DIAG )
C  FROM  * A PRACTICAL GUIDE TO SPLINES *  BY C. DE BOOR
CONSTRUCTS CHOLESKY FACTORIZATION
C                     C  =  L * D * L-TRANSPOSE
C  WITH L UNIT LOWER TRIANGULAR AND D DIAGONAL, FOR GIVEN MATRIX C OF
C  ORDER  N R O W , IN CASE  C  IS (SYMMETRIC) POSITIVE SEMIDEFINITE
C  AND  B A N D E D , HAVING  N B A N D S  DIAGONALS AT AND BELOW THE
C  MAIN DIAGONAL.
C
C******  I N P U T  ******
C  NROW.....IS THE ORDER OF THE MATRIX  C .
C  NBANDS.....INDICATES ITS BANDWIDTH, I.E.,
C          C(I,J) = 0 FOR ABS(I-J) .GT. NBANDS .
C  W.....WORKARRAY OF SIZE (NBANDS,NROW)  CONTAINING THE  NBANDS  DIAGO-
C        NALS IN ITS ROWS, WITH THE MAIN DIAGONAL IN ROW  1 . PRECISELY,
C        W(I,J)  CONTAINS  C(I+J-1,J), I=1,...,NBANDS, J=1,...,NROW.
C          FOR EXAMPLE, THE INTERESTING ENTRIES OF A SEVEN DIAGONAL SYM-
C        METRIC MATRIX  C  OF ORDER  9  WOULD BE STORED IN  W  AS
C
C                       11 22 33 44 55 66 77 88 99
C                       21 32 43 54 65 76 87 98
C                       31 42 53 64 75 86 97
C                       41 52 63 74 85 96
C
C        ALL OTHER ENTRIES OF  W  NOT IDENTIFIED IN THIS WAY WITH AN EN-
C        TRY OF  C  ARE NEVER REFERENCED .
C  DIAG.....IS A WORK ARRAY OF LENGTH  NROW .
C
C******  O U T P U T  ******
C  W.....CONTAINS THE CHOLESKY FACTORIZATION  C = L*D*L-TRANSP, WITH
C        W(1,I) CONTAINING  1/D(I,I)
C        AND  W(I,J)  CONTAINING  L(I-1+J,J), I=2,...,NBANDS.
C
C******  M E T H O D  ******
C   GAUSS ELIMINATION, ADAPTED TO THE SYMMETRY AND BANDEDNESS OF  C , IS
C   USED .
C     NEAR ZERO PIVOTS ARE HANDLED IN A SPECIAL WAY. THE DIAGONAL ELE-
C  MENT C(N,N) = W(1,N) IS SAVED INITIALLY IN  DIAG(N), ALL N. AT THE N-
C  TH ELIMINATION STEP, THE CURRENT PIVOT ELEMENT, VIZ.  W(1,N), IS COM-
C  PARED WITH ITS ORIGINAL VALUE, DIAG(N). IF, AS THE RESULT OF PRIOR
C  ELIMINATION STEPS, THIS ELEMENT HAS BEEN REDUCED BY ABOUT A WORD
C  LENGTH, (I.E., IF W(1,N)+DIAG(N) .LE. DIAG(N)), THEN THE PIVOT IS DE-
C  CLARED TO BE ZERO, AND THE ENTIRE N-TH ROW IS DECLARED TO BE LINEARLY
C  DEPENDENT ON THE PRECEDING ROWS. THIS HAS THE EFFECT OF PRODUCING
C   X(N) = 0  WHEN SOLVING  C*X = B  FOR  X, REGARDLESS OF  B. JUSTIFIC-
C  ATION FOR THIS IS AS FOLLOWS. IN CONTEMPLATED APPLICATIONS OF THIS
C  PROGRAM, THE GIVEN EQUATIONS ARE THE NORMAL EQUATIONS FOR SOME LEAST-
C  SQUARES APPROXIMATION PROBLEM, DIAG(N) = C(N,N) GIVES THE NORM-SQUARE
C  OF THE N-TH BASIS FUNCTION, AND, AT THIS POINT,  W(1,N)  CONTAINS THE
C  NORM-SQUARE OF THE ERROR IN THE LEAST-SQUARES APPROXIMATION TO THE N-
C  TH BASIS FUNCTION BY LINEAR COMBINATIONS OF THE FIRST N-1 . HAVING
C  W(1,N)+DIAG(N) .LE. DIAG(N) SIGNIFIES THAT THE N-TH FUNCTION IS LIN-
C  EARLY DEPENDENT TO MACHINE ACCURACY ON THE FIRST N-1 FUNCTIONS, THERE
C  FORE CAN SAFELY BE LEFT OUT FROM THE BASIS OF APPROXIMATING FUNCTIONS
C     THE SOLUTION OF A LINEAR SYSTEM
C                       C*X = B
C   IS EFFECTED BY THE SUCCESSION OF THE FOLLOWING  T W O  CALLS ...
C     CALL BCHFAC ( W, NBANDS, NROW, DIAG )       , TO GET FACTORIZATION
C     CALL BCHSLV ( W, NBANDS, NROW, B, X )            , TO SOLVE FOR X.
C
      REAL W(NBANDS,NROW),DIAG(NROW),   RATIO
      IF (NROW .GT. 1)                  GO TO 9
      IF (W(1,1) .GT. 0.) W(1,1) = 1./W(1,1)
                                        RETURN
C                                        STORE DIAGONAL OF  C  IN  DIAG.
    9 DO 10 N=1,NROW
   10    DIAG(N) = W(1,N)
C                                                        FACTORIZATION .
      DO 20 N=1,NROW
         IF (W(1,N)+DIAG(N) .GT. DIAG(N)) GO TO 15
         DO 14 J=1,NBANDS
   14       W(J,N) = 0.
                                        GO TO 20
   15    W(1,N) = 1./W(1,N)
         IMAX = MIN0(NBANDS-1,NROW - N)
         IF (IMAX .LT. 1)               GO TO 20
         JMAX = IMAX
         DO 18 I=1,IMAX
            RATIO = W(I+1,N)*W(1,N)
            NPI = N + I
            DO 17 J=1,JMAX
               IPJ = I + J
   17          W(J,NPI) = W(J,NPI) - W(IPJ,N)*RATIO
            JMAX = JMAX - 1
   18       W(I+1,N) = RATIO
   20    CONTINUE
                                        RETURN
      END
