wa1[i] -= r[j * ldr + i] * wa1[j];	}	temp = lm_enorm(n, wa1);	parc = fp / delta / temp / temp;        /** depending on the sign of the function, update parl or paru. **/	if (fp > 0)	    parl = MAX(parl, *par);	else if (fp < 0)	    paru = MIN(paru, *par);	/* the case fp==0 is precluded by the break condition  */                /** compute an improved estimate for par. **/        	*par = MAX(parl, *par + parc);            }} /*** lm_lmpar. ***/void lm_qrfac(int m, int n, double *a, int pivot, int *ipvt,	      double *rdiag, double *acnorm, double *wa){/* *     This subroutine uses householder transformations with column *     pivoting (optional) to compute a qr factorization of the *     m by n matrix a. That is, qrfac determines an orthogonal *     matrix q, a permutation matrix p, and an upper trapezoidal *     matrix r with diagonal elements of nonincreasing magnitude, *     such that a*p = q*r. The householder transformation for *     column k, k = 1,2,...,min(m,n), is of the form * *			    t *	    i - (1/u(k))*u*u * *     where u has zeroes in the first k-1 positions. The form of *     this transformation and the method of pivoting first *     appeared in the corresponding linpack subroutine. * *     Parameters: * *	m is a positive integer input variable set to the number *	  of rows of a. * *	n is a positive integer input variable set to the number *	  of columns of a. * *	a is an m by n array. On input a contains the matrix for *	  which the qr factorization is to be computed. On output *	  the strict upper trapezoidal part of a contains the strict *	  upper trapezoidal part of r, and the lower trapezoidal *	  part of a contains a factored form of q (the non-trivial *	  elements of the u vectors described above). * *	pivot is a logical input variable. If pivot is set true, *	  then column pivoting is enforced. If pivot is set false, *	  then no column pivoting is done. * *	ipvt is an integer output array of length lipvt. This array *	  defines the permutation matrix p such that a*p = q*r. *	  Column j of p is column ipvt(j) of the identity matrix. *	  If pivot is false, ipvt is not referenced. * *	rdiag is an output array of length n which contains the *	  diagonal elements of r. * *	acnorm is an output array of length n which contains the *	  norms of the corresponding columns of the input matrix a. *	  If this information is not needed, then acnorm can coincide *	  with rdiag. * *	wa is a work array of length n. If pivot is false, then wa *	  can coincide with rdiag. * */    int i, j, k, kmax, minmn;    double ajnorm, sum, temp;    static double p05 = 0.05;/*** qrfac: compute initial column norms and initialize several arrays. ***/    for (j = 0; j < n; j++) {	acnorm[j] = lm_enorm(m, &a[j * m]);	rdiag[j] = acnorm[j];	wa[j] = rdiag[j];	if (pivot)	    ipvt[j] = j;    }#if BUG    printf("qrfac\n");#endif/*** qrfac: reduce a to r with householder transformations. ***/    minmn = MIN(m, n);    for (j = 0; j < minmn; j++) {	if (!pivot)	    goto pivot_ok;        /** bring the column of largest norm into the pivot position. **/	kmax = j;	for (k = j + 1; k < n; k++)	    if (rdiag[k] > rdiag[kmax])		kmax = k;	if (kmax == j)	    goto pivot_ok;	for (i = 0; i < m; i++) {	    temp = a[j * m + i];	    a[j * m + i] = a[kmax * m + i];	    a[kmax * m + i] = temp;	}	rdiag[kmax] = rdiag[j];	wa[kmax] = wa[j];	k = ipvt[j];	ipvt[j] = ipvt[kmax];	ipvt[kmax] = k;      pivot_ok:        /** compute the Householder transformation to reduce the            j-th column of a to a multiple of the j-th unit vector. **/	ajnorm = lm_enorm(m - j, &a[j * m + j]);	if (ajnorm == 0.) {	    rdiag[j] = 0;	    continue;	}	if (a[j * m + j] < 0.)	    ajnorm = -ajnorm;	for (i = j; i < m; i++)	    a[j * m + i] /= ajnorm;	a[j * m + j] += 1;        /** apply the transformation to the remaining columns            and update the norms. **/	for (k = j + 1; k < n; k++) {	    sum = 0;	    for (i = j; i < m; i++)		sum += a[j * m + i] * a[k * m + i];	    temp = sum / a[j + m * j];	    for (i = j; i < m; i++)		a[k * m + i] -= temp * a[j * m + i];	    if (pivot && rdiag[k] != 0.) {		temp = a[m * k + j] / rdiag[k];		temp = MAX(0., 1 - temp * temp);		rdiag[k] *= sqrt(temp);		temp = rdiag[k] / wa[k];		if (p05 * SQR(temp) <= LM_MACHEP) {		    rdiag[k] = lm_enorm(m - j - 1, &a[m * k + j + 1]);		    wa[k] = rdiag[k];		}	    }	}	rdiag[j] = -ajnorm;    }}void lm_qrsolv(int n, double *r, int ldr, int *ipvt, double *diag,	       double *qtb, double *x, double *sdiag, double *wa){/* *     Given an m by n matrix a, an n by n diagonal matrix d, *     and an m-vector b, the problem is to determine an x which *     solves the system * *	    a*x = b  and  d*x = 0 * *     in the least squares sense. * *     This subroutine completes the solution of the problem *     if it is provided with the necessary information from the *     qr factorization, with column pivoting, of a. That is, if *     a*p = q*r, where p is a permutation matrix, q has orthogonal *     columns, and r is an upper triangular matrix with diagonal *     elements of nonincreasing magnitude, then qrsolv expects *     the full upper triangle of r, the permutation matrix p, *     and the first n components of (q transpose)*b. The system *     a*x = b, d*x = 0, is then equivalent to * *		   t	  t *	    r*z = q *b,  p *d*p*z = 0, * *     where x = p*z. If this system does not have full rank, *     then a least squares solution is obtained. On output qrsolv *     also provides an upper triangular matrix s such that * *	     t	 t		 t *	    p *(a *a + d*d)*p = s *s. * *     s is computed within qrsolv and may be of separate interest. * *     Parameters * *	n is a positive integer input variable set to the order of r. * *	r is an n by n array. On input the full upper triangle *	  must contain the full upper triangle of the matrix r. *	  On output the full upper triangle is unaltered, and the *	  strict lower triangle contains the strict upper triangle *	  (transposed) of the upper triangular matrix s. * *	ldr is a positive integer input variable not less than n *	  which specifies the leading dimension of the array r. * *	ipvt is an integer input array of length n which defines the *	  permutation matrix p such that a*p = q*r. Column j of p *	  is column ipvt(j) of the identity matrix. * *	diag is an input array of length n which must contain the *	  diagonal elements of the matrix d. * *	qtb is an input array of length n which must contain the first *	  n elements of the vector (q transpose)*b. * *	x is an output array of length n which contains the least *	  squares solution of the system a*x = b, d*x = 0. * *	sdiag is an output array of length n which contains the *	  diagonal elements of the upper triangular matrix s. * *	wa is a work array of length n. * */    int i, kk, j, k, nsing;    double qtbpj, sum, temp;    double _sin, _cos, _tan, _cot; /* local variables, not functions */    static double p25 = 0.25;    static double p5 = 0.5;/*** qrsolv: copy r and (q transpose)*b to preserve input and initialize s.     in particular, save the diagonal elements of r in x. ***/    for (j = 0; j < n; j++) {	for (i = j; i < n; i++)	    r[j * ldr + i] = r[i * ldr + j];	x[j] = r[j * ldr + j];	wa[j] = qtb[j];    }#if BUG    printf("qrsolv\n");#endif/*** qrsolv: eliminate the diagonal matrix d using a givens rotation. ***/    for (j = 0; j < n; j++) {/*** qrsolv: prepare the row of d to be eliminated, locating the     diagonal element using p from the qr factorization. ***/	if (diag[ipvt[j]] == 0.)	    goto L90;	for (k = j; k < n; k++)	    sdiag[k] = 0.;	sdiag[j] = diag[ipvt[j]];/*** qrsolv: the transformations to eliminate the row of d modify only      a single element of (q transpose)*b beyond the first n, which is     initially 0.. ***/	qtbpj = 0.;	for (k = j; k < n; k++) {            /** determine a givens rotation which eliminates the                appropriate element in the current row of d. **/	    if (sdiag[k] == 0.)		continue;	    kk = k + ldr * k;	    if (fabs(r[kk]) < fabs(sdiag[k])) {		_cot = r[kk] / sdiag[k];		_sin = p5 / sqrt(p25 + p25 * SQR(_cot));		_cos = _sin * _cot;	    } else {		_tan = sdiag[k] / r[kk];		_cos = p5 / sqrt(p25 + p25 * SQR(_tan));		_sin = _cos * _tan;	    }            /** compute the modified diagonal element of r and                the modified element of ((q transpose)*b,0). **/	    r[kk] = _cos * r[kk] + _sin * sdiag[k];	    temp = _cos * wa[k] + _sin * qtbpj;	    qtbpj = -_sin * wa[k] + _cos * qtbpj;	    wa[k] = temp;            /** accumulate the tranformation in the row of s. **/	    for (i = k + 1; i < n; i++) {		temp = _cos * r[k * ldr + i] + _sin * sdiag[i];		sdiag[i] = -_sin * r[k * ldr + i] + _cos * sdiag[i];		r[k * ldr + i] = temp;	    }	}      L90:        /** store the diagonal element of s and restore            the corresponding diagonal element of r. **/	sdiag[j] = r[j * ldr + j];	r[j * ldr + j] = x[j];    }/*** qrsolv: solve the triangular system for z. if the system is     singular, then obtain a least squares solution. ***/    nsing = n;    for (j = 0; j < n; j++) {	if (sdiag[j] == 0. && nsing == n)	    nsing = j;	if (nsing < n)	    wa[j] = 0;    }    for (j = nsing - 1; j >= 0; j--) {	sum = 0;	for (i = j + 1; i < nsing; i++)	    sum += r[j * ldr + i] * wa[i];	wa[j] = (wa[j] - sum) / sdiag[j];    }/*** qrsolv: permute the components of z back to components of x. ***/    for (j = 0; j < n; j++)	x[ipvt[j]] = wa[j];} /*** lm_qrsolv. ***/double lm_enorm(int n, double *x){/*     Given an n-vector x, this function calculates the *     euclidean norm of x. * *     The euclidean norm is computed by accumulating the sum of *     squares in three different sums. The sums of squares for the *     small and large components are scaled so that no overflows *     occur. Non-destructive underflows are permitted. Underflows *     and overflows do not occur in the computation of the unscaled *     sum of squares for the intermediate components. *     The definitions of small, intermediate and large components *     depend on two constants, LM_SQRT_DWARF and LM_SQRT_GIANT. The main *     restrictions on these constants are that LM_SQRT_DWARF**2 not *     underflow and LM_SQRT_GIANT**2 not overflow. * *     Parameters * *	n is a positive integer input variable. * *	x is an input array of length n. */    int i;    double agiant, s1, s2, s3, xabs, x1max, x3max, temp;    s1 = 0;    s2 = 0;    s3 = 0;    x1max = 0;    x3max = 0;    agiant = LM_SQRT_GIANT / ((double) n);    /** sum squares. **/    for (i = 0; i < n; i++) {	xabs = fabs(x[i]);	if (xabs > LM_SQRT_DWARF && xabs < agiant) {            /*  sum for intermediate components. */	    s2 += xabs * xabs;	    continue;	}	if (xabs > LM_SQRT_DWARF) {            /*  sum for large components. */	    if (xabs > x1max) {		temp = x1max / xabs;		s1 = 1 + s1 * SQR(temp);		x1max = xabs;	    } else {		temp = xabs / x1max;		s1 += SQR(temp);	    }	    continue;	}        /*  sum for small components. */	if (xabs > x3max) {	    temp = x3max / xabs;	    s3 = 1 + s3 * SQR(temp);	    x3max = xabs;	} else {	    if (xabs != 0.) {		temp = xabs / x3max;		s3 += SQR(temp);	    }	}    }    /** calculation of norm. **/    if (s1 != 0)	return x1max * sqrt(s1 + (s2 / x1max) / x1max);    if (s2 != 0) {	if (s2 >= x3max)	    return sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));	else	    return sqrt(x3max * ((s2 / x3max) + (x3max * s3)));    }    return x3max * sqrt(s3);} /*** lm_enorm. ***/