/*

                              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"
/*---------------------------------------------------------------------------*/
/* Table of constant values */


/*function:external*/
extern void bchslvlot(C(Pfloat *)w, C(Pint) nbands, C(Pint) nrow, C(Pint) m, 
                      C(Pint) md, C(Pfloat *)b)
PreANSI(Pfloat *w) 
PreANSI(Pint nbands) 
PreANSI(Pint nrow) 
PreANSI(Pint m) 
PreANSI(Pint md) 
PreANSI(Pfloat *b)
/*
Documentation - to be completed
*/
{
    /* System generated locals */
    Pint w_dim1, w_offset, b_dim1, b_dim2, b_offset;

    /* Local variables */
    Pint jmax, j, n, nbndm1, ii, mm;

/* +++++++++++++++++++++++++++++++++++++++++++++++ */

/*  solves the linear system     c*x = b   of order  n r o w  for  x */
/*  provided  w  contains the cholesky factorization for the banded (sym- */
/*  metric) positive definite matrix  c  as constructed in the subroutine */
/*    b c h f a c  (quo vide). */

/* ******  i n p u t  ****** */
/*  nrow.....is the order of the matrix  c . */
/*  nbands.....indicates the bandwidth of  c . */
/*  m..........the dimension of the elements */
/*  md.........the number of right hand sides of b */
/*  w.....contains the cholesky factorization for  c . as output from */
/* 	subroutine bchfac  (quo vide). */
/*  b.....the matrix  of length  n r o w  containing the M right side. */

/* ******  o u t p u t  ****** */
/*  b.....the matrix of length  n r o w  containing the M solutions. */

/* ******  m e t h o d  ****** */
/*  with the factorization  c = l*d*l-transpose  available, where  l  is */
/*  unit lower triangular and  d  is diagonal, the triangular system */
/*  l*y = b  is solved for  y (forward substitution), y is stored in  b, */
/*  the vector  d**(-1)*y is computed and stored in  b, then the triang- */
/*  ular system  l-transpose*x = d**(-1)*y is solved for  x (backsubstit- */
/*  ution). */

/* WTH August 1984 */
/* ------------------------------------------------------------------- */

    /* Parameter adjustments */
    b_dim1 = m;
    b_dim2 = md;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    w_dim1 = nbands;
    w_offset = w_dim1 + 1;
    w -= w_offset;

    /* Function Body */
    if (nrow <= 1) {
	for (mm = 1; mm <= md; ++mm) {
	    for (ii = 1; ii <= m; ++ii) {
		b[ii + (mm + b_dim2) * b_dim1] *= w[w_dim1 + 1];
	    }
	}
	return;
    }

/* 	forward substitution. solve l*y = b for y, store in b. */
    nbndm1 = nbands - 1;
    for (n = 1; n <= nrow; ++n) {
/* Computing MIN */
	jmax = P_min(nbndm1,nrow-n);
	if (jmax >= 1) {
	    for (j = 1; j <= jmax; ++j) {
		for (mm = 1; mm <= md; ++mm) {
		    for (ii = 1; ii <= m; ++ii) {
			b[ii + (mm + (j + n) * b_dim2) * b_dim1] -= w[j + 1 + 
				n * w_dim1] * b[ii + (mm + n * b_dim2) * 
				b_dim1];
		    }
		}
	    }
	}
    }

/* 	backsubstitution. solve l-transp.x = d**(-1)*y  for s, store in b. */
    for (n = nrow; n >= 1; --n) {
	for (mm = 1; mm <= md; ++mm) {
	    for (ii = 1; ii <= m; ++ii) {
		b[ii + (mm + n * b_dim2) * b_dim1] *= w[n * w_dim1 + 1];
	    }
	}
/* Computing MIN */
	jmax = P_min(nbndm1,nrow-n);
	if (jmax >= 1) {
	    for (j = 1; j <= jmax; ++j) {
		for (mm = 1; mm <= md; ++mm) {
		    for (ii = 1; ii <= m; ++ii) {
			b[ii + (mm + n * b_dim2) * b_dim1] -= w[j + 1 + n * 
				w_dim1] * b[ii + (mm + (j + n) * b_dim2) * 
				b_dim1];
		    }
		}
	    }
	}
    }

    return;
} /* bchslvlot */
/*----------------------------------------------------------------------------*/
/*function:external*/
extern void banslv(C(Pfloat *)w, C(Pint) nroww, C(Pint) nrow, C(Pint) m, 
                   C(Pint) mmax, C(Pint) nbandl, C(Pint) nbandu, C(Pfloat *)b)
PreANSI(Pfloat *w) 
PreANSI(Pint nroww) 
PreANSI(Pint nrow) 
PreANSI(Pint m)
PreANSI(Pint mmax) 
PreANSI(Pint nbandl) 
PreANSI(Pint nbandu) 
PreANSI(Pfloat *b)
/*
Documentation - to be completed
*/
{
    /* System generated locals */
    Pint w_dim1, w_offset, b_dim1, b_offset;

    /* Local variables */
    Pint jmax, i, j, ii, nrowm1, middle;

/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* 	companion routine to banfac . it returns the solution x of the */
/* 	linear system a*x = b in place of b , given the lu-factorization */
/* 	for a in the workarray w . */

/* ******  i n p u t  ****** */
/* 	w, nroww,nrow,nbandl,nbandu.....describe the lu-factorization of a */
/* 	banded matrix a of order nrow as constructed in banfac . */
/* 	for details, see banfac . */
/* 	b.....right side of the systems to be solved . */
/*       m.....dimesnions of points in b */
/*       mmax..maximum dimension of points */
/* ******  o u t p u t  ****** */
/* 	b.....contains the solution x , of order nrow . */

/* ******  m e t h o d  ****** */
/* 	(with a = l*u, as stored in w,) the unit lower triangular system */
/* 	l(u*x) = b is solved for y = u*x, and y stored in b . then the */
/* 	upper triangular system u*x = y is solved for x . the calcul- */
/* 	ations are so arranged that the innermost loops stay within columns. */
/*------------------------------------------------------------------------*/
    /* Parameter adjustments */
    b_dim1 = mmax;
    b_offset = b_dim1 + 1;
    b -= b_offset;
    w_dim1 = nroww;
    w_offset = w_dim1 + 1;
    w -= w_offset;

    /* Function Body */
    middle = nbandu + 1;
    if (nrow == 1) {
	goto L49;
    }
    nrowm1 = nrow - 1;
    if (nbandl == 0) {
	goto L30;
    }

/* 		forward pass */
/* 	for i=1,2,...,nrow-1, subtract right side(i)*(i-th column */
/* 	of l ) from right side (below i-th row) . */

    for (i = 1; i <= nrowm1; ++i) {
/* Computing MIN */
	jmax = P_min(nbandl,nrow-i);
	for (j = 1; j <= jmax; ++j) {
	    for (ii = 1; ii <= m; ++ii) {
		b[ii + (i + j) * b_dim1] -= b[ii + i * b_dim1] * w[middle + j 
			+ i * w_dim1];
	    }
	}
    }

/* 			backward pass */
/* 	for i=nrow,nrow-1,...,1, divide right side(i) by i-th diag- */
/* 	onal entry of u, then subtract right side(i)*(i-th column */
/* 	of u) from right side (above i-th row) . */

L30:
    if (nbandu > 0) {
	goto L40;
    }

/*  a is lower triangular . */

    for (i = 1; i <= nrow; ++i) {
	for (ii = 1; ii <= m; ++ii) {
	    b[ii + i * b_dim1] /= w[i * w_dim1 + 1];
	}
    }
    return;
/* ----------------- */

L40:
    for (i = nrow; i >= 2; --i) {
	for (ii = 1; ii <= m; ++ii) {
	    b[ii + i * b_dim1] /= w[middle + i * w_dim1];
	}
/* Computing MIN */
	jmax = P_min(nbandu,i-1);
	for (j = 1; j <= jmax; ++j) {
	    for (ii = 1; ii <= m; ++ii) {
		b[ii + (i - j) * b_dim1] -= b[ii + i * b_dim1] * w[middle - j 
			+ i * w_dim1];
	    }
	}
    }
L49:
    for (ii = 1; ii <= m; ++ii) {
	b[ii + b_dim1] /= w[middle + w_dim1];
    }
    return;
} /* banslv*/

/*---------------------------------------------------------------------------*/
/*function:external*/
extern void spli2d(C(Pfloat *)tau, C(Pfloat *)gtau, C(Pfloat *)t, C(Pint) n, 
	C(Pint) k, C(Pint) m, C(Pint) mmax, C(Pint) md, 
        C(Pfloat *)work, C(Pfloat *)q, C(Pfloat *)bcoef, C(Pint *)iflag)
PreANSI(Pfloat *tau) 
PreANSI(Pfloat *gtau) 
PreANSI(Pfloat *t) 
PreANSI(Pint n)
PreANSI(Pint k) 
PreANSI(Pint m) 
PreANSI(Pint mmax) 
PreANSI(Pint md) 
PreANSI(Pfloat *work) 
PreANSI(Pfloat *q) 
PreANSI(Pfloat *bcoef) 
PreANSI(Pint *iflag)
/*
Documentation - to be completed
*/

{

    /* System generated locals */
    Pint bcoef_dim1, bcoef_dim2, bcoef_offset, gtau_dim1, gtau_dim2, 
	    gtau_offset, work_dim1, work_offset;

    /* Local variables */
    Pint left, lenq;
    Pfloat taui;
    Pint kpkm2, i, j, ilp1mx;

    Pint ii, jj;
    Pint km1, np1;
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* calls bsplvb, banfac/slv */
/* This is an extended version of splint, for use in tensor product */
/* interpolation */
/*  spli2d produces the b-spline coeff.s bcoef(j,.) of the spline of order*/
/*   k with knots  t(i), i=1,..., n + k , which takes on the value */
/*   gtau(i,j) at tau(i), i=1,..., n , j=1,..., md . */
/* ******  i n p u t  ****** */
/*  tau...array of length  n , containing data point abscissae. */
/*    a s s u m p t i o n . . .  tau  is strictly increasing */
/* gtau(.,j).....corresponding array of length  n , containing data point or-*/
/* 	 dinates */
/*  t.....knot sequence, of length n+k */
/*  n.....number of data points and dimension of spline space s(k,t) */
/*  k.....order of spline */
/*  m.....the dimemsion of each point */
/* mmax..the maximum dimension of each point (e.g. 4 but set m=3 to keep Wi = 1)*/
/*  md....number of data sets */

/* ****** W o r k  A r e a ***** */

/* work...a vector of length n */

/* ******  o u t p u t  ****** */
/*  q.....array of size  (2*k-1)*n , containing the triangular factoriz- */
/* 	 ation of the coefficient matrix of the linear system for the b- */
/* 	 coefficients of the spline interpolant. */
/* 	    the b-coeffs for the interpolant of an additional data set */
/* 	 (tau(i),htau(i)), i=1,...,n with the same data abscissae can */
/* 	 be obtained without going through all the calculations in this */
/* 	 routine, simply by loading htau into bcoef and then execut- */
/* 	 ing the  call banslv ( q, 2*k-1, n, k-1, k-1, bcoef ) */
/*  bcoef.....the b-coefficients of the interpolant, of length n */
/*  iflag.....an integer indicating success (= 1) or failure (= 2) */
/* 	 the linear system to be solved is (theoretically) invertible if */
/* 	 and only if */
/* 		t(i) .lt. tau(i) .lt.  t(i+k),  all i. */
/* 	 violation of this condition is certain to lead to  iflag = 2 . */

/* ******  m e t h o d  ****** */
/* 	the i-th equation of the linear system a*bcoef = b for the b-co- */
/*  effs of the interpolant enforces interpolation at tau(i), i=1,...,n. */
/*  hence, b(i) = gtau(i), all i, and a is a band matrix with 2k-1 */
/*   bands (if it is invertible). */
/* 	the matrix a is generated row by row and stored, diagonal by di- */
/*  agonal, in the r o w s of the array q , with the main diagonal go- */
/*  ing into row k .  see comments in the program below. */
/* 	the banded system is then solved by a call to banfac (which con- */
/*  structs the triangular factorization for a and stores it again in */
/*  q ), followed by a call to banslv (which then obtains the solution */
/*  bcoef by substitution). */
/* 	BANFAC  does no pivoting, since the total positivity of the matrix */
/*  a makes this unnecessary. */
/*-------------------------------------------------------------------------*/
/* 	dimension q(2*k-1,n), t(n+k) */
/* current fortran standard makes it impossible to specify precisely the */
/*  dimension of  q  and  t  without the introduction of otherwise super- */
/*  fluous additional arguments. */

 printf("SPLI2D data\n");
 printf("n:%3ld, k:%3ld, m:%3ld, mmax:%3ld, md:%3ld\n",n,k,m,mmax,md);
 printf("tau:"); for (i=0; i<n; i++) printf("%9.3f",tau[i]); printf("\n");
 printf("t:"); for (i=0; i<n+k; i++) printf("%9.3f",t[i]); printf("\n");
 

    /* Parameter adjustments */
    bcoef_dim1 = mmax;
    bcoef_dim2 = md;
    bcoef_offset = bcoef_dim1 * (bcoef_dim2 + 1) + 1;
    bcoef -= bcoef_offset;
    --q;
    work_dim1 = mmax;
    work_offset = work_dim1 + 1;
    work -= work_offset;
    --t;
    gtau_dim1 = mmax;
    gtau_dim2 = n;
    gtau_offset = gtau_dim1 * (gtau_dim2 + 1) + 1;
    gtau -= gtau_offset;
    --tau;

    /* Function Body */
    np1 = n + 1;
    km1 = k - 1;
    kpkm2 = km1 << 1;
    left = k;
/*  zero out all entries of q */
    lenq = n * (k + km1);
    for (i = 1; i <= lenq; ++i) {
	q[i] = 0.0;
    }

/*  ***  loop over i to construct the n interpolation equations */
    for (i = 1; i <= n; ++i) {
	taui = tau[i];
/* Computing MIN */
	ilp1mx = P_min(i+k,np1);
/*  ***	find left in the closed interval (i,i+k-1) such that */
/* 		t(left) .le. tau(i) .lt. t(left+1) */
/* 	matrix is singular if this is not possible */

	left = P_max(left,i);
	if (taui < t[left]) {
	    goto L998;
	}
L15:
	if (taui < t[left + 1]) {
	    goto L16;
	}
	++left;
	if (left < ilp1mx) {
	    goto L15;
	}
	--left;
	if (taui > t[left + 1]) {
	    goto L998;
	}

/*  ***	the i-th equation enforces interpolation at taui, hence */
/* 	a(i,j) = b(j,k,t)(taui), all j. only the k entries with j = */
/* 	left-k+1,...,left actually might be nonzero. these k numbers */
/* 	are returned, in work (used for temp.storage here), by the */
/* 	following */
L16:
	bsplvb(&t[1], k, 1, taui, left, &work[work_offset]);
/* 	we therefore want work(j) = b(left-k+j)(taui) to go into */
/* 	a(i,left-k+j), i.e., into q(i-(left+j)+2*k,(left+j)-k) since */
/* 	a(i+j,j)  is to go into q(i+k,j), all i,j, if we consider q */
/* 	as a two-dim. array , with 2*k-1 rows (see comments in */
/* 	banfac). in the present program, we treat q as an equivalent */
/* 	one-dimensional array (because of fortran restrictions on */
/* 	dimension statements) . we therefore want work(j) go to into */
/* 	entry */
/* 	i -(left+j) + 2*k + ((left+j) - k-1)*(2*k-1) */
/* 		= i-left+1 + (left-k)*(2*k-1) + (2*k-2)*j */
/* 	of q . */

	jj = i - left + 1 + (left - k) * (k + km1);
	for (j = 1; j <= k; ++j) {
	    jj += kpkm2;
	    q[jj] = work[j + work_dim1];
	}
    }

/*  ***	obtain factorization of a , stored again in q. */

    banfac(&q[1], k+km1, n, km1, km1, iflag);
    switch (*iflag) {
	case 1:  goto L40;
	case 2:  goto L999;
    }

/*  ***	solve a*bcoef = gtau by backsubstitution */

L40:
    for (j = 1; j <=md; ++j) {
	for (i = 1; i <= n; ++i) {
	    for (ii = 1; ii <= m; ++ii) {
		work[ii + i * work_dim1] = gtau[ii + (i + j * gtau_dim2) * 
			gtau_dim1];
	    }
	}
	banslv(&q[1],k+km1, n, m, mmax, km1, km1, &work[work_offset]);
	for (i = 1; i <= n; ++i) {
	    for (ii = 1; ii <= m; ++ii) {
		bcoef[ii + (j + i * bcoef_dim2) * bcoef_dim1] = work[ii + i * 
			work_dim1];
	    }
	}
    }
    return;
/* --------------------- */
L998:
    *iflag = 2;
L999:
    nrb_error("NRB_SPLI2D - Linear system not invertible");
    return;
} /* spli2d*/
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void ll2appr(C(Pfloat *)tau, C(Pfloat *)gtau, C(Pfloat *)weight, 
                    C(Pint) ntau, C(Pfloat *)t, C(Pint) n, C(Pint) k, 
                    C(Pint) m, C(Pint) md, C(Pfloat *)q, C(Pfloat *)diag, 
                    C(Pfloat *)bcoef)
PreANSI(Pfloat *tau)
PreANSI(Pfloat *gtau) 
PreANSI(Pfloat *weight) 
PreANSI(Pint ntau) 
PreANSI(Pfloat *t) 
PreANSI(Pint n) 
PreANSI(Pint k) 
PreANSI(Pint m) 
PreANSI(Pint md) 
PreANSI(Pfloat *q) 
PreANSI(Pfloat *diag) 
PreANSI(Pfloat *bcoef)
/*
Documentation - to be completed
*/
{
    /* System generated locals */
    Pint gtau_dim1, gtau_dim2, gtau_offset, bcoef_dim1, bcoef_dim2, 
	    bcoef_offset, q_dim1, q_offset;

    /* Local variables */
    Pint left, i, j;
    Pfloat biatx[20];
    Pint id, jj, ll, mm;
    Pfloat dw;
    Pint leftmk;
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* constructs the (weighted discrete) l2-approximation by splines of order*/
/*  k  with knot sequence  t(1), ..., t(n+k)  to given data points */
/*  ( tau(i), gtau(i,.) ), i=1,...,ntau. the b-spline coefficients */
/*  b c o e f   of the approximating spline are determined from the */
/*  normal equations using cholesky's method. */

/* ******  i n p u t  ****** */
/*  ntau.....number of data points */
/* (tau(.),gtau(m,.,md)), i=1,...,md   are the  md data sets to be fitted.*/
/*  weight(i), i=1,...,ntau    are the corresponding weights. */
/*  m is the dimension of each point */
/*  md is the number of data sets */
/*  t(1), ..., t(n+k)  the knot sequence */
/*  n.....the dimension of the space of splines of order k with knots t. */
/*  k.....the order */

/*  w a r n i n g  . . .  the restriction   k .le. kmax (= 20)   is impo- */
/* 	sed by the arbitrary dimension statement for  biatx  below, but */
/* 	is  n o w h e r e   c h e c k e d   for. */

/* ******  w o r k  a r r a y s  ****** */
/*  q....a work array of size (at least) k*n. its first  k  rows are used */
/* 	for the  k  lower diagonals of the gramian matrix  c . */
/*  diag.....a work array of length  n  used in bchfac . */

/* ******  o u t p u t  ****** */
/*  bcoef(m,i,.) i=1,md the b-spline coeffs. of the l2-appr. */

/* ******  m e t h o d  ****** */
/*  the b-spline coefficients of the l2-appr. are determined as the sol- */
/*  ution of the normal equations */
/* 	sum ( (b(mm,i),b(mm,j))*bcoef(mm,j) : j=1,...,n)  = (b(mm,i),g), */
/* 	i = 1, ..., n . */
/*       and mm =1, ..., m. */
/*  here,  b(mm,i)  denotes the i-th b-spline,  g  denotes the function to */
/*  be approximated, and the  i n n e r   p r o d u c t  of two funct- */
/*  ions  f  and  g  is given by */
/* 	(f,g)  :=  sum ( f(tau(i))*g(tau(i))*weight(i) : i=1,...,ntau)  . */
/*  the arrays  t a u  and  w e i g h t  are given in common block */
/*   d a t a , as is the array  g t a u  containing the sequence */
/*  g(tau(i)), i=1,...,ntau. */
/*  the relevant function values of the b-splines  b(i), i=1,...,n, are */

/*  supplied by the subprogram  b s p l v b . */
/* 	the coeff.matrix  c , with */
/* 	c(i,j)  :=  (b(i),b(j)), i,j=1,...,n, */
/*  of the normal equations is symmetric and (2*k-1)-banded, therefore */
/*  can be specified by giving its k bands at or below the diagonal. for */
/*  i=1,...,n,  we store */
/*   (b(i),b(j))  =  c(i,j)  in  q(i-j+1,j), j=i,...,min0(i+k-1,n) */
/*  and the right side */
/*   (b(i), g )  in  bcoef(i) . */
/*  since b-spline values are most efficiently generated by finding sim- */
/*  ultaneously the value of  e v e r y  nonzero b-spline at one point, */
/*  the entries of  c  (i.e., of  q ), are generated by computing, for */
/*  each ll, all the terms involving  tau(ll)  simultaneously and adding */
/*  them to all relevant entries. */
/* -------------------------------------------------------- */
/* 	dimension t(n+k) */
/* current fortran standard makes it impossible to specify the exact dimen- */
/*  sion of  t  without the introduction of additional otherwise super- */
/*  fluous arguments. */

    /* Parameter adjustments */
    bcoef_dim1 = m;
    bcoef_dim2 = md;
    bcoef_offset = bcoef_dim1 * (bcoef_dim2 + 1) + 1;
    bcoef -= bcoef_offset;
    --diag;
    q_dim1 = k;
    q_offset = q_dim1 + 1;
    q -= q_offset;
    --t;
    --weight;
    gtau_dim1 = m;
    gtau_dim2 = ntau;
    gtau_offset = gtau_dim1 * (gtau_dim2 + 1) + 1;
    gtau -= gtau_offset;
    --tau;

    /* Function Body */
    for (j = 1; j <= n; ++j) {
	for (id = 1; id <= md; ++id) {
	    for (i = 1; i <= m; ++i) {
		bcoef[i + (id + j * bcoef_dim2) * bcoef_dim1] = 0.0;
	    }
	}
	for (i = 1; i <= k; ++i) {
	    q[i + j * q_dim1] = 0.0;
	}
    }

    left = k;
    leftmk = 0;
    for (ll = 1; ll <= ntau; ++ll) {
/* 	locate  l e f t  s.t. tau(ll) in (t(left),t(left+1)) */
L10:
	if (left == n) {
	    goto L15;
	}
	if (tau[ll] < t[left + 1]) {
	    goto L15;
	}
	++left;
	++leftmk;
	goto L10;
L15:
	bsplvb(&t[1], k, 1, tau[ll], left, biatx);

/* 	biatx(mm) contains the value of b(left-k+mm) at tau(ll). */
/* 	hence, with  dw := biatx(mm)*weight(ll), the number dw*gtau(ll) */

/* 	is a summand in the inner product */
/* 	   (b(left=k+mm), g)  which goes into  bcoef(left-k+mm) */
/* 	and the number biatx(jj)*dw is a cummand in the inner product */
/* 	   (b(left-k+jj), b(left-k+mm)), into  q(jj-mm+1,left-k+mm) */
/* 	since  (left-k+jj) - (left-k+mm) + 1  =  jj - mm + 1 . */

	for (mm = 1; mm <= k; ++mm) {
	    dw = biatx[mm - 1] * weight[ll];
	    j = leftmk + mm;
	    for (id = 1; id <= md; ++id) {
		for (i = 1; i <= m; ++i) {
		    bcoef[i + (id + j * bcoef_dim2) * bcoef_dim1] = dw * gtau[
			    i + (ll + id * gtau_dim2) * gtau_dim1] + bcoef[i 
			    + (id + j * bcoef_dim2) * bcoef_dim1];
		}
	    }
	    i = 1;
	    for (jj = mm; jj <= k; ++jj) {
		q[i + j * q_dim1] = biatx[jj - 1] * dw + q[i + j * q_dim1];
/* L20: */
		++i;
	    }
	}
    }

/* 	construct cholesky factorization for  c  in  q , then use */
/* 	it to solve the normal equations */
/* 	c*x  =  bcoef */
/* 	for  x , and store  x  in bcoef . */
    bchfac(&q[q_offset], k, n, &diag[1]);
    bchslvlot(&q[q_offset], k, n, m, md, &bcoef[bcoef_offset]);
    return;
} /* ll2appr*/

