1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
|
/* lmpar.f -- translated by f2c (version 20020621).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "cminpack.h"
#include <math.h>
#define real __cminpack_real__
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
__cminpack_attr__
void __cminpack_func__(lmpar)(int n, real *r, int ldr,
const int *ipvt, const real *diag, const real *qtb, real delta,
real *par, real *x, real *sdiag, real *wa1,
real *wa2)
{
/* Initialized data */
#define p1 .1
#define p001 .001
/* System generated locals */
real d1, d2;
/* Local variables */
int i, j, k, l;
real fp;
real sum, parc, parl;
int iter;
real temp, paru, dwarf;
int nsing;
real gnorm;
real dxnorm;
/* ********** */
/* subroutine lmpar */
/* given an m by n matrix a, an n by n nonsingular diagonal */
/* matrix d, an m-vector b, and a positive number delta, */
/* the problem is to determine a value for the parameter */
/* par such that if x solves the system */
/* a*x = b , sqrt(par)*d*x = 0 , */
/* in the least squares sense, and dxnorm is the euclidean */
/* norm of d*x, then either par is zero and */
/* (dxnorm-delta) .le. 0.1*delta , */
/* or par is positive and */
/* abs(dxnorm-delta) .le. 0.1*delta . */
/* 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 lmpar expects */
/* the full upper triangle of r, the permutation matrix p, */
/* and the first n components of (q transpose)*b. on output */
/* lmpar also provides an upper triangular matrix s such that */
/* t t t */
/* p *(a *a + par*d*d)*p = s *s . */
/* s is employed within lmpar and may be of separate interest. */
/* only a few iterations are generally needed for convergence */
/* of the algorithm. if, however, the limit of 10 iterations */
/* is reached, then the output par will contain the best */
/* value obtained so far. */
/* the subroutine statement is */
/* subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, */
/* wa1,wa2) */
/* where */
/* 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. */
/* delta is a positive input variable which specifies an upper */
/* bound on the euclidean norm of d*x. */
/* par is a nonnegative variable. on input par contains an */
/* initial estimate of the levenberg-marquardt parameter. */
/* on output par contains the final estimate. */
/* x is an output array of length n which contains the least */
/* squares solution of the system a*x = b, sqrt(par)*d*x = 0, */
/* for the output par. */
/* sdiag is an output array of length n which contains the */
/* diagonal elements of the upper triangular matrix s. */
/* wa1 and wa2 are work arrays of length n. */
/* subprograms called */
/* minpack-supplied ... dpmpar,enorm,qrsolv */
/* fortran-supplied ... dabs,dmax1,dmin1,dsqrt */
/* argonne national laboratory. minpack project. march 1980. */
/* burton s. garbow, kenneth e. hillstrom, jorge j. more */
/* ********** */
/* dwarf is the smallest positive magnitude. */
dwarf = __cminpack_func__(dpmpar)(2);
/* compute and store in x the gauss-newton direction. if the */
/* jacobian is rank-deficient, obtain a least squares solution. */
nsing = n;
for (j = 0; j < n; ++j) {
wa1[j] = qtb[j];
if (r[j + j * ldr] == 0. && nsing == n) {
nsing = j;
}
if (nsing < n) {
wa1[j] = 0.;
}
/* L10: */
}
if (nsing >= 1) {
for (k = 1; k <= nsing; ++k) {
j = nsing - k;
wa1[j] /= r[j + j * ldr];
temp = wa1[j];
if (j >= 1) {
for (i = 0; i < j; ++i) {
wa1[i] -= r[i + j * ldr] * temp;
}
}
}
}
for (j = 0; j < n; ++j) {
l = ipvt[j]-1;
x[l] = wa1[j];
}
/* initialize the iteration counter. */
/* evaluate the function at the origin, and test */
/* for acceptance of the gauss-newton direction. */
iter = 0;
for (j = 0; j < n; ++j) {
wa2[j] = diag[j] * x[j];
}
dxnorm = __cminpack_func__(enorm)(n, wa2);
fp = dxnorm - delta;
if (fp <= p1 * delta) {
goto TERMINATE;
}
/* if the jacobian is not rank deficient, the newton */
/* step provides a lower bound, parl, for the zero of */
/* the function. otherwise set this bound to zero. */
parl = 0.;
if (nsing >= n) {
for (j = 0; j < n; ++j) {
l = ipvt[j]-1;
wa1[j] = diag[l] * (wa2[l] / dxnorm);
}
for (j = 0; j < n; ++j) {
sum = 0.;
if (j >= 1) {
for (i = 0; i < j; ++i) {
sum += r[i + j * ldr] * wa1[i];
}
}
wa1[j] = (wa1[j] - sum) / r[j + j * ldr];
}
temp = __cminpack_func__(enorm)(n, wa1);
parl = fp / delta / temp / temp;
}
/* calculate an upper bound, paru, for the zero of the function. */
for (j = 0; j < n; ++j) {
sum = 0.;
for (i = 0; i <= j; ++i) {
sum += r[i + j * ldr] * qtb[i];
}
l = ipvt[j]-1;
wa1[j] = sum / diag[l];
}
gnorm = __cminpack_func__(enorm)(n, wa1);
paru = gnorm / delta;
if (paru == 0.) {
paru = dwarf / min(delta,(real)p1);
}
/* if the input par lies outside of the interval (parl,paru), */
/* set par to the closer endpoint. */
*par = max(*par,parl);
*par = min(*par,paru);
if (*par == 0.) {
*par = gnorm / dxnorm;
}
/* beginning of an iteration. */
for (;;) {
++iter;
/* evaluate the function at the current value of par. */
if (*par == 0.) {
/* Computing MAX */
d1 = dwarf, d2 = p001 * paru;
*par = max(d1,d2);
}
temp = sqrt(*par);
for (j = 0; j < n; ++j) {
wa1[j] = temp * diag[j];
}
__cminpack_func__(qrsolv)(n, r, ldr, ipvt, wa1, qtb, x, sdiag, wa2);
for (j = 0; j < n; ++j) {
wa2[j] = diag[j] * x[j];
}
dxnorm = __cminpack_func__(enorm)(n, wa2);
temp = fp;
fp = dxnorm - delta;
/* if the function is small enough, accept the current value */
/* of par. also test for the exceptional cases where parl */
/* is zero or the number of iterations has reached 10. */
if (fabs(fp) <= p1 * delta || (parl == 0. && fp <= temp && temp < 0.) || iter == 10) {
goto TERMINATE;
}
/* compute the newton correction. */
for (j = 0; j < n; ++j) {
l = ipvt[j]-1;
wa1[j] = diag[l] * (wa2[l] / dxnorm);
}
for (j = 0; j < n; ++j) {
wa1[j] /= sdiag[j];
temp = wa1[j];
if (n > j+1) {
for (i = j+1; i < n; ++i) {
wa1[i] -= r[i + j * ldr] * temp;
}
}
}
temp = __cminpack_func__(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);
}
if (fp < 0.) {
paru = min(paru,*par);
}
/* compute an improved estimate for par. */
/* Computing MAX */
d1 = parl, d2 = *par + parc;
*par = max(d1,d2);
/* end of an iteration. */
}
TERMINATE:
/* termination. */
if (iter == 0) {
*par = 0.;
}
/* last card of subroutine lmpar. */
} /* lmpar_ */
|