--- mpfit.c.orig 2014-06-05 10:02:14.819618639 -0500
+++ mpfit.c 2014-06-05 10:02:21.131645649 -0500
@@ -13,7 +13,7 @@
*/
/* Main mpfit library routines (double precision)
- $Id: mpfit.c,v 1.20 2010/11/13 08:15:35 craigm Exp $
+ $Id: mpfit.c,v 1.23 2013/04/23 04:24:16 craigm Exp $
*/
#include <stdio.h>
@@ -29,7 +29,8 @@
double *wa, void *priv, int *nfev,
double *step, double *dstep, int *dside,
int *qulimited, double *ulimit,
- int *ddebug, double *ddrtol, double *ddatol);
+ int *ddebug, double *ddrtol, double *ddatol,
+ double *wa2, double **dvecptr);
static void mp_qrfac(int m, int n, double *a, int lda,
int pivot, int *ipvt, int lipvt,
double *rdiag, double *acnorm, double *wa);
@@ -274,7 +275,7 @@
{
mp_config conf;
int i, j, info, iflag, nfree, npegged, iter;
- int qanylim = 0, qanypegged = 0;
+ int qanylim = 0;
int ij,jj,l;
double actred,delta,dirder,fnorm,fnorm1,gnorm, orignorm;
@@ -297,6 +298,7 @@
double *fvec = 0, *qtf = 0;
double *x = 0, *xnew = 0, *fjac = 0, *diag = 0;
double *wa1 = 0, *wa2 = 0, *wa3 = 0, *wa4 = 0;
+ double **dvecptr = 0;
int *ipvt = 0;
int ldfjac;
@@ -323,6 +325,7 @@
if (config->nprint >= 0) conf.nprint = config->nprint;
if (config->epsfcn > 0) conf.epsfcn = config->epsfcn;
if (config->maxiter > 0) conf.maxiter = config->maxiter;
+ if (config->maxiter == MP_NO_ITER) conf.maxiter = 0;
if (config->douserscale != 0) conf.douserscale = config->douserscale;
if (config->covtol > 0) conf.covtol = config->covtol;
if (config->nofinitecheck > 0) conf.nofinitecheck = config->nofinitecheck;
@@ -334,6 +337,7 @@
nfree = 0;
npegged = 0;
+ /* Basic error checking */
if (funct == 0) {
return MP_ERR_FUNC;
}
@@ -442,6 +446,7 @@
mp_malloc(wa3, double, npar);
mp_malloc(wa4, double, m);
mp_malloc(ipvt, int, npar);
+ mp_malloc(dvecptr, double *, npar);
/* Evaluate user function with initial parameter values */
iflag = mp_call(funct, m, npar, xall, fvec, 0, private_data);
@@ -483,30 +488,34 @@
iflag = mp_fdjac2(funct, m, nfree, ifree, npar, xnew, fvec, fjac, ldfjac,
conf.epsfcn, wa4, private_data, &nfev,
step, dstep, mpside, qulim, ulim,
- ddebug, ddrtol, ddatol);
+ ddebug, ddrtol, ddatol, wa2, dvecptr);
if (iflag < 0) {
goto CLEANUP;
}
/* Determine if any of the parameters are pegged at the limits */
- qanypegged = 0;
if (qanylim) {
for (j=0; j<nfree; j++) {
int lpegged = (qllim[j] && (x[j] == llim[j]));
int upegged = (qulim[j] && (x[j] == ulim[j]));
sum = 0;
-
+
+ /* If the parameter is pegged at a limit, compute the gradient
+ direction */
if (lpegged || upegged) {
- qanypegged = 1;
ij = j*ldfjac;
for (i=0; i<m; i++, ij++) {
sum += fvec[i] * fjac[ij];
}
}
+ /* If pegged at lower limit and gradient is toward negative then
+ reset gradient to zero */
if (lpegged && (sum > 0)) {
ij = j*ldfjac;
for (i=0; i<m; i++, ij++) fjac[ij] = 0;
}
+ /* If pegged at upper limit and gradient is toward positive then
+ reset gradient to zero */
if (upegged && (sum < 0)) {
ij = j*ldfjac;
for (i=0; i<m; i++, ij++) fjac[ij] = 0;
@@ -624,7 +633,10 @@
*/
if (gnorm <= conf.gtol) info = MP_OK_DIR;
if (info != 0) goto L300;
- if (conf.maxiter == 0) goto L300;
+ if (conf.maxiter == 0) {
+ info = MP_MAXITER;
+ goto L300;
+ }
/*
* rescale if necessary.
@@ -965,7 +977,7 @@
if (qulim) free(qulim);
if (llim) free(llim);
if (ulim) free(ulim);
-
+ if (dvecptr) free(dvecptr);
return info;
}
@@ -980,7 +992,8 @@
double *wa, void *priv, int *nfev,
double *step, double *dstep, int *dside,
int *qulimited, double *ulimit,
- int *ddebug, double *ddrtol, double *ddatol)
+ int *ddebug, double *ddrtol, double *ddatol,
+ double *wa2, double **dvec)
{
/*
* **********
@@ -1063,17 +1076,15 @@
int iflag = 0;
double eps,h,temp;
static double zero = 0.0;
- double **dvec = 0;
int has_analytical_deriv = 0, has_numerical_deriv = 0;
int has_debug_deriv = 0;
temp = mp_dmax1(epsfcn,MP_MACHEP0);
eps = sqrt(temp);
ij = 0;
- ldfjac = 0; /* Prevents compiler warning */
+ ldfjac = 0; /* Prevent compiler warning */
+ if (ldfjac){} /* Prevent compiler warning */
- dvec = (double **) malloc(sizeof(double **)*npar);
- if (dvec == 0) return MP_ERR_MEMORY;
for (j=0; j<npar; j++) dvec[j] = 0;
/* Initialize the Jacobian derivative matrix */
@@ -1164,12 +1175,12 @@
(fjold == 0)?(0):((fjold-fjac[ij])/fjold));
}
}
- }
+ } /* end debugging */
- } else {
+ } else { /* dside > 2 */
/* COMPUTE THE TWO-SIDED DERIVATIVE */
- for (i=0; i<m; i++, ij++) {
- fjac[ij] = wa[i]; /* Store temp data: fjac[i+m*j] */
+ for (i=0; i<m; i++) {
+ wa2[i] = wa[i];
}
/* Evaluate at x - h */
@@ -1180,15 +1191,16 @@
x[ifree[j]] = temp;
/* Now compute derivative as (f(x+h) - f(x-h))/(2h) */
- ij -= m;
if (! debug ) {
+ /* Non-debug path for speed */
for (i=0; i<m; i++, ij++) {
fjac[ij] = (fjac[ij] - wa[i])/(2*h); /* fjac[i+m*j] */
}
} else {
+ /* Debug path for correctness */
for (i=0; i<m; i++, ij++) {
double fjold = fjac[ij];
- fjac[ij] = (fjac[ij] - wa[i])/(2*h); /* fjac[i+m*j] */
+ fjac[ij] = (wa2[i] - wa[i])/(2*h); /* fjac[i+m*j] */
if ((da == 0 && dr == 0 && (fjold != 0 || fjac[ij] != 0)) ||
((da != 0 || dr != 0) && (fabs(fjold-fjac[ij]) > da + fabs(fjold)*dr))) {
printf(" %10d %10.4g %10.4g %10.4g %10.4g %10.4g\n",
@@ -1196,17 +1208,16 @@
(fjold == 0)?(0):((fjold-fjac[ij])/fjold));
}
}
- }
+ } /* end debugging */
- }
- }
+ } /* if (dside > 2) */
+ } /* if (has_numerical_derivative) */
if (has_debug_deriv) {
printf("FJAC DEBUG END\n");
}
DONE:
- if (dvec) free(dvec);
if (iflag < 0) return iflag;
return 0;
/*
@@ -1305,8 +1316,10 @@
static double one = 1.0;
static double p05 = 0.05;
- lda = 0; /* Prevent compiler warning */
- lipvt = 0; /* Prevent compiler warning */
+ lda = 0; /* Prevent compiler warning */
+ lipvt = 0; /* Prevent compiler warning */
+ if (lda) {} /* Prevent compiler warning */
+ if (lipvt) {} /* Prevent compiler warning */
/*
* compute the initial column norms and initialize several arrays.