/*

                              DISCLAIMER
                              ==========

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.   

    If the software is modified by someone else and passed on, we, the authors
    want its recipients to know that what they have is not the original, so
    that any problems introduced by others will not reflect on the original
    authors' reputations.
*/                                            

#include "nurbh.h"
/*---------------------------------------------------------------------------*/
static Pint j = 1;
static Pfloat saved, deltal[100], deltar[100];
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void bsplvb(C(Pfloat *)t, C(Pint) jhigh, C(Pint) index, C(Pfloat) x, 
                   C(Pint) left, C(Pfloat *)biatx)
PreANSI(Pfloat *t)
PreANSI(Pint jhigh) 
PreANSI(Pint index) 
PreANSI(Pfloat x) 
PreANSI(Pint left) 
PreANSI(Pfloat *biatx)
/*
Documentation - to be completed
*/
{
    /* Local variables */
    Pfloat term;
    Pint i;
    Pint jp1;

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* calculates the value of all possibly nonzero b-splines at  x  of order */
/* 		jout  =  max( jhigh , (j+1)*(index-1) ) */
/*  with knot sequence  t . */
/* ******  I N P U T  ****** */
/* t.....knot sequence, of length  left + jout  , assumed to be nondecreasing.*/
/* 		A S S U M P T I O N . . . . */
/* 			t(left) .lt. t(left + 1)   . */
/*   DIVISION  BY  ZERO  WILL RESULT IF  T(LEFT)  = T(LEFT+1) */
/*  jhigh, */
/*  index.....integers which determine the order  jout = max(jhigh, */
/* 	(j+1)*(index-1))  of the b-splines whose values at  x  are to */
/* 	be returned.  index  is used to avoid recalculations when seve- */
/* 	ral columns of the triangular array of b-spline values are nee- */
/* 	ded (e.g., in  bsplpp  or in  bsplvd ). precisely, */
/* 			if  index = 1 , */
/* 	the calculation starts from scratch and the entire triangular */
/* 	array of b-spline values of orders 1,2,...,jhigh  is generated */
/* 	order by order , i.e., column by column . */
/* 			if  index = 2 , */
/* 	only the b-spline values of order  j+1, j+2, ..., jout  are ge- */
/* 	nerated, the assumption being that  biatx , j , deltal , deltar */
/* 	are, on entry, as they were on exit at the previous call. */
/* 	   in particular, if  jhigh = 0, then  jout = j+1, i.e., just */
/* 	the next column of b-spline values is generated. */
/* W A R N I N G . . .  the restriction   jout .le. jmax (= 20)  is im- */
/* 	posed arbitrarily by the dimension statement for  deltal  and */
/* 	deltar  below, but is  NOWHERE CHECKED FOR . */
/* x.....the point at which the b-splines are to be evaluated. */
/* left.....an integer chosen (usually) so that */
/* 		t(left) .le. x .le. t(left)  . */

/* ******  O U T P U T  ****** */
/*  biatx.....array of length  jout , with  biatx(i)  containing the val- */
/* 	ue at  x  of the polynomial of order  jout  which agrees with */
/* 	the b-spline  b(left-jout+i,jout,t)  on the interval (t(left), */
/* 	t(left+1)) . */

/* ******  M E T H O D  ****** */
/*  the recurrence relation */
/* 			x - t(i)		t(i+j+1) - x */
/* 	b(i,j+1)(x)  =  -----------b(i,j)(x) + ---------------b(i+1,j)(x) */
/* 			t(i+j)-t(i)	       t(i+j+1)-t(i+1) */

/*  is used (repeatedly) to generate the (j+1)-vector  b(left-j,j+1)(x), */
/*  ...,b(left,j+1)(x)  from the j-vector  b(left-j+1,j)(x),..., */
/*  b(left,j)(x), storing the new values in  biatx  over the old. the */
/*  facts that */
/* 		b(i,1) = 1  if  t(i) .le. x .lt. t(i+1) */
/*  and that */
/* 		b(i,j)(x) = 0  unless  t(i) .le. x .lt. t(i+j) */
/*  are used. */
/* The particular organization of the calculations follows algorithm  (8) */
/* in chapter x of the text. */

/*-------------------------------------------------------------------------*/
/* 	real biatx(jhigh),t(1),x,   deltal(jmax),deltar(jmax),saved,term */
/* 	dimension biatx(jout), t(left+jout) */
/* current fortran standard makes it impossible to specify the length of */
/*  t  and of  biatx  precisely without the introduction of otherwise */
/*  superfluous additional arguments. */

    /* Parameter adjustments */
    --biatx;
    --t;

    /* Function Body */

    if (index ==1){
      j = 1;
      biatx[1] = 1.0;
      if (j >= jhigh) return;
    }

    while (j< jhigh){

      jp1 = j + 1;
      deltar[j - 1] = t[left + j] - x;
      deltal[j - 1] = x - t[left + 1 - j];
      saved = 0.0;
      for (i = 1; i <= j; ++i) {
	term = biatx[i] / (deltar[i - 1] + deltal[jp1 - i - 1]);
	biatx[i] = saved + deltar[i - 1] * term;
	saved = deltal[jp1 - i - 1] * term;
      }
      biatx[jp1] = saved;
      j = jp1;
    }

} /* bsplvb*/
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void bchfac(C(Pfloat *)w, C(Pint) nbands, C(Pint) nrow, 
                   C(Pfloat *)diag)
PreANSI(Pfloat *w) 
PreANSI(Pint nbands) 
PreANSI(Pint nrow) 
PreANSI(Pfloat *diag)
/*
Documentation - to be completed
*/
{
    /* Local variables */
    Pint imax, jmax, i, j, n;
    Pfloat ratio;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++ */
/*  constructs cholesky factorization */
/* 	c  =  l * d * l-transpose */
/*  with l unit lower triangular and d diagonal, for given matrix c of */
/*  order  n r o w , in case  c  is (symmetric) positive semidefinite */
/*  and  b a n d e d , having  n b a n d s  diagonals at and below the */
/*  main diagonal. */
/* ******  i n p u t  ****** */
/*  nrow.....is the order of the matrix  c . */
/*  nbands.....indicates its bandwidth, i.e., */
/* 	c(i,j) = 0 for  i-j   .ge.  nbands . */
/*  w.....workarray of size (nbands,nrow)  containing the  nbands  diago- */
/* 	nals 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 sym- */
/* 	metric matrix  c  or 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 */

/* 	all other entries of  w  not identified in this way with an en- */
/* 	try of  c  are never referenced . */
/*  diag.....is a work array of length  nrow . */

/* ******  o u t p u t ****** */
/*  w.....contains the cholesky factorization  c = l*d*l-transp, with */
/* 	w(1,i) containing 1/d(i,i) */
/* 	and  w(i,j)  containing  l(i-1+j,j), i=2,...,nbands. */

/* ******  m e t h o d  ****** */
/*   gauss elimination, adapted to the symmetry and bandedness of  c , is */
/*   used . */
/* 	near zero pivots are handled in a special way. the diagonal ele- */
/*  ment c(n,n) = w(1,n) is saved initially in  diag(n), all n. at the n- */
/*  th elimination step, the current pivot element, viz.  w(1,n), is com- */
/*  pared with its original value, diag(n). if, as the result of prior */
/*  elimination steps, this element has been reduced by about a word */
/*  elngth, (i.e., if w(1,n)+diag(n) .le. diag(n)), then the pivot is de- */
/*  clared 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. justific- */
/*  ation 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 .le. diag(n) signifies that the n-th function is lin- */
/*  early dependent to machine accuracy on the first n-1 functions, there */
/*  fore 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  t w o  calls: */
/* 	call bchfac ( w, nbands, nrow, diag )       , to get factorization */
/* 	call bchslv ( w, nbands, nrow, b    )            , to solve for x. */
/* ----------------------------------------------------- */
/* 	real w(nbands,nrow),diag(nrow),   ratio */

    /* Parameter adjustments */
    --diag;
    --w;

    /* Function Body */
    if (nrow <= 1){
/* WTH Optimize
    if (w[ni(1,1,nbands)] > 0.0) w[ni(1,1,nbands)] = 1.0 / w[ni(1,1,nbands)];
*/
     if (w[1] > 0.0) w[1] = 1.0 / w[1];    
      return;
    }
/* ------------------------- */
/* 	store diagonal of  c  in  diag. */

    for (n = 1; n <= nrow; ++n) {
	diag[n] = w[ni(1,n,nbands)];
    }
/* 	factorization . */
    for (n = 1; n <= nrow; ++n) {
	if (w[ni(1,n,nbands)] + diag[n] > diag[n]) {
	    goto L15;
	}
	for (j = 1; j <= nbands; ++j) {
	    w[ni(j,n,nbands)] = 0.0;
	}
	goto L20;
L15:
	w[ni(1,n,nbands)] = 1.0 / w[ni(1,n,nbands)];
/* Computing MIN */
	imax = P_min(nbands,(nrow-n));
	if (imax < 1) {
	    goto L20;
	}
	jmax = imax;
	for (i = 1; i <= imax; ++i) {
	    ratio = w[ni(i+1,n,nbands)] * w[ni(1,n,nbands)];
	    for (j = 1; j <= jmax; ++j) {
		w[ni(j,n+i,nbands)] -= w[ni(j+i,n,nbands)] * ratio;
	    }
	    --jmax;
	    w[ni(i+1,n,nbands)] = ratio;
	}
L20:
	;
    }

} /* bchfac */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void banfac(C(Pfloat *)w, C(Pint) nroww, C(Pint) nrow, C(Pint) nbandl, 
                   C(Pint) nbandu, C(Pint *)iflag)
PreANSI(Pfloat *w)
PreANSI(Pint nroww) 
PreANSI(Pint nrow) 
PreANSI(Pint nbandl) 
PreANSI(Pint nbandu)
PreANSI(Pint *iflag) 
/*
Documentation - to be completed
*/
{
    /* Local variables */
    Pint jmax, kmax, i, j, k, midmk;
    Pfloat pivot;
    Pint nrowm1, middle;
    Pfloat factor;
    Pint ipk;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/*  returns in w the lu-factorization (without pivoting) of the banded */
/*  matrix a of order nrow with (nbandl + 1 + nbandu) bands or diag- */
/*  onals in the work array w . */
/* ******  i n p u t  ****** */
/*  w.....work array of size (nroww,nrow) containing the interesting */
/* 	part of a banded matrix a , with the diagonals or bands of a */
/* 	stored in the rows of w , while columns of a correspond to */
/* 	columns of w . this is the storage mode used in linpack and */
/* 	results in efficient innermost loops. */
/* 	explicitly, a has nbandl bands below the diagonal */
/* 			+ 1  (main) diagonal */
/* 			+  nbandu bands above the diagonal */
/* 	and thus, with  middle = nbandu + 1, */
/* 	a(i+j,j) is in w(i+middle,j) for i=-nbandu,...,nbandl */
/* 				j=1,...,nrow . */
/* 	for example, the interesting entries of a (1,2)-banded matrix */
/* 	of order 9 would appear in the first 1+1+2 = 4 rows of w */
/* 	as follows. */
/* 			13 24 35 46 57 68 79 */
/* 	             12 23 34 45 56 67 78 89 */
/* 	          11 22 33 44 55 66 77 88 99 */
/* 	       21 32 43 54 65 76 87 98 */
/*            . */
/* 	all other entries of w not identified in this way with an en- */
/* 	try of a are never referenced . */
/* 	nroww.....row dimension of the work array w . */
/* 	must be .ge. nbandl + 1 + nbandu . */
/* 	nbandl.....number of bands of a below the main diagonal */
/* 	nbandu.....number of bands of a above the main diagonal . */

/* ******  o u t p u t  ****** */
/* 	iflag.....integer indicating success( = 1) or failure ( = 2) . */
/* 	if iflag = 1, then */
/* 	w.....contains the lu-factorization of a into a unit lower triangu- */

/* 	lar matrix l and an upper triangular matrix u (both banded) */
/* 	and stored in customary fashion over the corresponding entries */
/* 	of a . this makes it possible to solve any particular linear */
/* 	system a*x = b for x by a */
/* 	call banslv ( w, nroww, nrow, m, mmax, nbandl, nbandu, b ) */
/* 	with the solution x contained in b on return . */
/* 	if iflag = 2, then */
/* 	one of nrow-1, nbandl,nbandu failed to be nonnegative, or else */
/* 	one of the potential pivots was found to be zero indicating */
/* 	that a does not have an lu-factorization. this implies that */
/* 	a is singular in case it is totally positive . */
/*       m is the dimesnion of the points */

/* ******  m e t h o d  ****** */
/* 	gauss elimination  w i t h o u t  pivoting is used. the routine is */
/* 	intended for use with matrices a which do not require row inter- */
/* 	changes during factorization, especially for the t o t a l l y */
/* 	p o s i t i v e  matrices which occur in spline calculations. */
/* 	the routine should not be used for an arbitrary banded matrix. */

/* -----------------------------------------------------------------------*/

    /* Parameter adjustments */
    --w;
    
    /* Function Body */
    *iflag = 1;
    middle = nbandu + 1;

/*  w(middle,.) cibtaubs tge naub duagibak if a . */

    nrowm1 =  nrow - 1;
    if (nrowm1 < 0) {
	goto L999;
    } else if (nrowm1 == 0) {
	goto L900;
    } else {
	goto L1;
    }
L1:
    if ( nbandl > 0) {
	goto L10;
    }

/*  a is upper triangular. check that diagonal is nonzero . */

    for (i = 1; i <= nrowm1; ++i) {
	if (w[ni(middle,i,nroww)] == 0.0) {
	    goto L999;
	}
    }
    goto L900;
/* ---------------- */
L10:
    if ( nbandu > 0) {
	goto L20;
    }

/* 		a is lower triangular. check that diagonal is nonzero and */
/* 		divide each column by its diagonal . */

    for (i = 1; i <= nrowm1; ++i) {
	pivot = w[ni(middle,i,nroww)];
	if (pivot == 0.0) {
	    goto L999;
	}
/* Computing MIN */

	jmax = P_min(nbandl,(nrow-i));
	for (j = 1; j <= jmax; ++j) {
	    w[ni(middle + j, i,nroww)] /= pivot;
	}
    }
    goto L900;
/* ------------------- */

/* 	a is not just a triangular matrix. construct lu factorization */

L20:
    for (i = 1; i <= nrowm1; ++i) {

/*  w(middle,i) is pivot for i-th step . */

	pivot = w[ni(middle,i,nroww)];
	if (pivot == 0.0) {
	    goto L999;
	}
/* 		jmax is the number of (nonzero) entries in column i */
/* 			below the diagonal . */
/* Computing MIN */
	jmax = P_min(nbandl,(nrow-i));
/* 		divide each entry in column i below diagonal by pivot . */

	for (j = 1; j <= jmax; ++j) {
	    w[ni(middle + j, i,nroww)] /= pivot;
	}

/* 			kmax is the number of (nonzero) entries in row i to */
/* 			    the right of the diagonal . */
/* Computing MIN */
	kmax = P_min(nbandu,(nrow-i));

/* 		subtract a(i,i+k)*(i-th column) from (i+k)-th column */
/* 		(below row i ) . */

	for (k = 1; k <= kmax; ++k) {
	    ipk = i + k;
	    midmk = middle - k;
	    factor = w[ni(midmk,ipk,nroww)];
	    for (j = 1; j <= jmax; ++j) {
		w[ni(midmk + j,ipk,nroww)] -= w[ni(middle + j,i,nroww)] * 
			factor;
	    }
	}
    }
/*  check the last diagonal entry . */
L900:
    if (w[ni(middle,nrow,nroww)] != 0.0) {
	return;
    }
L999:
    *iflag = 2;
} /* banfac */

/*--------------------------------------------------------------------------*/

