diff options
| -rw-r--r-- | unsupported/test/NonLinear.cpp | 848 |
1 files changed, 820 insertions, 28 deletions
diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 74ea2dedb..972c57c2d 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -7,11 +7,55 @@ #include "main.h" #include <stdio.h> -#include <math.h> #include <cminpack.h> -int fcn(int m, int n, const double *x, double *fvec, double *fjac, int ldfjac, int iflag); +int fcn_chkder(int /*m*/, int /*n*/, const double *x, double *fvec, double *fjac, int ldfjac, int iflag) +{ + /* subroutine fcn for chkder example. */ + + int i; + double tmp1, tmp2, tmp3, tmp4; + double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, + 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; + + + if (iflag == 0) + { + /* insert print statements here when nprint is positive. */ + return 0; + } + + if (iflag != 2) + + for (i=1; i<=15; i++) + { + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3)); + } + else + { + for (i = 1; i <= 15; i++) + { + tmp1 = i; + tmp2 = 16 - i; + + /* error introduced into next statement for illustration. */ + /* corrected statement should read tmp3 = tmp1 . */ + + tmp3 = tmp2; + if (i > 8) tmp3 = tmp2; + tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4=tmp4*tmp4; + fjac[i-1+ ldfjac*(1-1)] = -1.; + fjac[i-1+ ldfjac*(2-1)] = tmp1*tmp2/tmp4; + fjac[i-1+ ldfjac*(3-1)] = tmp1*tmp3/tmp4; + } + } + return 0; +} void testChkder() { @@ -32,9 +76,9 @@ void testChkder() ldfjac = 15; chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, 1, err); - fcn(m, n, x, fvec, fjac, ldfjac, 1); - fcn(m, n, x, fvec, fjac, ldfjac, 2); - fcn(m, n, xp, fvecp, fjac, ldfjac, 1); + fcn_chkder(m, n, x, fvec, fjac, ldfjac, 1); + fcn_chkder(m, n, x, fvec, fjac, ldfjac, 2); + fcn_chkder(m, n, xp, fvecp, fjac, ldfjac, 1); chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, 2, err); for (i=1; i<=m; i++) @@ -68,20 +112,95 @@ void testChkder() for (i=1; i<=m; i++) VERIFY_IS_APPROX(fvec[i-1], fvec_ref[i-1]); for (i=1; i<=m; i++) VERIFY_IS_APPROX(fvecp[i-1], fvecp_ref[i-1]); for (i=1; i<=m; i++) VERIFY_IS_APPROX(err[i-1], err_ref[i-1]); - return; } -int fcn(int /*m*/, int /*n*/, const double *x, double *fvec, - double *fjac, int ldfjac, int iflag) + +int fcn_lmder1(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjac, + int ldfjac, int iflag) { - /* subroutine fcn for chkder example. */ + +/* subroutine fcn for lmder1 example. */ + + int i; + double tmp1, tmp2, tmp3, tmp4; + double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, + 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; + + if (iflag != 2) + { + for (i = 1; i <= 15; i++) + { + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3)); + } + } + else + { + for ( i = 1; i <= 15; i++) + { + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4 = tmp4*tmp4; + fjac[i-1 + ldfjac*(1-1)] = -1.; + fjac[i-1 + ldfjac*(2-1)] = tmp1*tmp2/tmp4; + fjac[i-1 + ldfjac*(3-1)] = tmp1*tmp3/tmp4; + } + } + return 0; +} + +void testLmder1() +{ + int j, m, n, ldfjac, info, lwa; + int ipvt[3]; + double tol, fnorm; + double x[3], fvec[15], fjac[15*3], wa[30]; + + m = 15; + n = 3; + +/* the following starting values provide a rough fit. */ + + x[1-1] = 1.; + x[2-1] = 1.; + x[3-1] = 1.; + + ldfjac = 15; + lwa = 30; + +/* set tol to the square root of the machine precision. */ +/* unless high solutions are required, */ +/* this is the recommended setting. */ + + tol = sqrt(dpmpar(1)); + + info = lmder1(fcn_lmder1, 0, m, n, x, fvec, fjac, ldfjac, tol, + ipvt, wa, lwa); + fnorm = enorm(m, fvec); + + VERIFY_IS_APPROX(fnorm, 0.09063596); + VERIFY(info == 1); + double x_ref[] = {0.08241058, 1.133037, 2.343695 }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); +} + + +int fcn_lmder(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjac, + int ldfjac, int iflag) +{ + + /* subroutine fcn for lmder example. */ int i; double tmp1, tmp2, tmp3, tmp4; double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; - if (iflag == 0) { /* insert print statements here when nprint is positive. */ @@ -89,38 +208,711 @@ int fcn(int /*m*/, int /*n*/, const double *x, double *fvec, } if (iflag != 2) + { + for (i=1; i <= 15; i++) + { + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3)); + } + } + else + { + for (i=1; i<=15; i++) + { + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4 = tmp4*tmp4; + fjac[i-1 + ldfjac*(1-1)] = -1.; + fjac[i-1 + ldfjac*(2-1)] = tmp1*tmp2/tmp4; + fjac[i-1 + ldfjac*(3-1)] = tmp1*tmp3/tmp4; + }; + } + return 0; +} - for (i=1; i<=15; i++) - { - tmp1 = i; - tmp2 = 16 - i; - tmp3 = tmp1; - if (i > 8) tmp3 = tmp2; - fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3)); - } +void testLmder() +{ + int i, j, m, n, ldfjac, maxfev, mode, nprint, info, nfev, njev; + int ipvt[3]; + double ftol, xtol, gtol, factor, fnorm; + double x[3], fvec[15], fjac[15*3], diag[3], qtf[3], + wa1[3], wa2[3], wa3[3], wa4[15]; + double covfac; + + m = 15; + n = 3; + +/* the following starting values provide a rough fit. */ + + x[1-1] = 1.; + x[2-1] = 1.; + x[3-1] = 1.; + + ldfjac = 15; + + /* set ftol and xtol to the square root of the machine */ + /* and gtol to zero. unless high solutions are */ + /* required, these are the recommended settings. */ + + ftol = sqrt(dpmpar(1)); + xtol = sqrt(dpmpar(1)); + gtol = 0.; + + maxfev = 400; + mode = 1; + factor = 1.e2; + nprint = 0; + + info = lmder(fcn_lmder, 0, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, + maxfev, diag, mode, factor, nprint, &nfev, &njev, + ipvt, qtf, wa1, wa2, wa3, wa4); + fnorm = enorm(m, fvec); + + VERIFY_IS_APPROX(fnorm, 0.09063596); + + VERIFY(nfev==6); + VERIFY(njev==5); + VERIFY(info==1); + double x_ref[] = {0.08241058, 1.133037, 2.343695 }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); + ftol = dpmpar(1); + covfac = fnorm*fnorm/(m-n); + covar(n, fjac, ldfjac, ipvt, ftol, wa1); + + double cov_ref[] = { + 0.0001531202, 0.002869941, -0.002656662, + 0.002869941, 0.09480935, -0.09098995, + -0.002656662, -0.09098995, 0.08778727 + }; + + for (i=1; i<=n; i++) + for (j=1; j<=n; j++) + VERIFY_IS_APPROX(fjac[(i-1)*ldfjac+j-1]*covfac, cov_ref[(i-1)*3+(j-1)]); +} + + +int fcn_hybrj1(void * /*p*/, int n, const double *x, double *fvec, double *fjac, int ldfjac, + int iflag) +{ + /* subroutine fcn for hybrj1 example. */ + + int j, k; + double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4; + + if (iflag != 2) + { + for (k = 1; k <= n; k++) + { + temp = (three - two*x[k-1])*x[k-1]; + temp1 = zero; + if (k != 1) temp1 = x[k-1-1]; + temp2 = zero; + if (k != n) temp2 = x[k+1-1]; + fvec[k-1] = temp - temp1 - two*temp2 + one; + } + } + else + { + for (k = 1; k <= n; k++) + { + for (j = 1; j <= n; j++) + { + fjac[k-1 + ldfjac*(j-1)] = zero; + } + fjac[k-1 + ldfjac*(k-1)] = three - four*x[k-1]; + if (k != 1) fjac[k-1 + ldfjac*(k-1-1)] = -one; + if (k != n) fjac[k-1 + ldfjac*(k+1-1)] = -two; + } + } + return 0; +} + + +void testHybrj1() +{ + int j, n, ldfjac, info, lwa; + double tol, fnorm; + double x[9], fvec[9], fjac[9*9], wa[99]; + + n = 9; + +/* the following starting values provide a rough solution. */ + + for (j=1; j<=9; j++) + { + x[j-1] = -1.; + } + + ldfjac = 9; + lwa = 99; + +/* set tol to the square root of the machine precision. */ +/* unless high solutions are required, */ +/* this is the recommended setting. */ + + tol = sqrt(dpmpar(1)); + + info = hybrj1(fcn_hybrj1, 0, n, x, fvec, fjac, ldfjac, tol, wa, lwa); + + fnorm = enorm(n, fvec); + + VERIFY_IS_APPROX(fnorm, 1.192636e-08); + VERIFY(info==1); + double x_ref[] = { + -0.5706545, -0.6816283, -0.7017325, + -0.7042129, -0.701369, -0.6918656, + -0.665792, -0.5960342, -0.4164121 + }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); +} + +int fcn_hybrj(void * /*p*/, int n, const double *x, double *fvec, double *fjac, int ldfjac, + int iflag) +{ + + /* subroutine fcn for hybrj example. */ + + int j, k; + double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4; + + if (iflag == 0) + { + /* insert print statements here when nprint is positive. */ + return 0; + } + + if (iflag != 2) + { + for (k=1; k <= n; k++) + { + temp = (three - two*x[k-1])*x[k-1]; + temp1 = zero; + if (k != 1) temp1 = x[k-1-1]; + temp2 = zero; + if (k != n) temp2 = x[k+1-1]; + fvec[k-1] = temp - temp1 - two*temp2 + one; + } + } + else + { + for (k = 1; k <= n; k++) + { + for (j=1; j <= n; j++) + { + fjac[k-1 + ldfjac*(j-1)] = zero; + } + fjac[k-1 + ldfjac*(k-1)] = three - four*x[k-1]; + if (k != 1) fjac[k-1 + ldfjac*(k-1-1)] = -one; + if (k != n) fjac[k-1 + ldfjac*(k+1-1)] = -two; + } + } + return 0; +} + +void testHybrj() +{ + + int j, n, ldfjac, maxfev, mode, nprint, info, nfev, njev, lr; + double xtol, factor, fnorm; + double x[9], fvec[9], fjac[9*9], diag[9], r[45], qtf[9], + wa1[9], wa2[9], wa3[9], wa4[9]; + + n = 9; + +/* the following starting values provide a rough solution. */ + + for (j=1; j<=9; j++) + { + x[j-1] = -1.; + } + + ldfjac = 9; + lr = 45; + +/* set xtol to the square root of the machine precision. */ +/* unless high solutions are required, */ +/* this is the recommended setting. */ + + xtol = sqrt(dpmpar(1)); + + maxfev = 1000; + mode = 2; + for (j=1; j<=9; j++) + { + diag[j-1] = 1.; + } + factor = 1.e2; + nprint = 0; + + info = hybrj(fcn_hybrj, 0, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, + mode, factor, nprint, &nfev, &njev, r, lr, qtf, + wa1, wa2, wa3, wa4); + fnorm = enorm(n, fvec); + + VERIFY_IS_APPROX(fnorm, 1.192636e-08); + VERIFY(nfev==11); + VERIFY(njev==1); + VERIFY(info==1); + + double x_ref[] = { + -0.5706545, -0.6816283, -0.7017325, + -0.7042129, -0.701369, -0.6918656, + -0.665792, -0.5960342, -0.4164121 + }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); + printf("\n"); +} + +int fcn_hybrd1(void * /*p*/, int n, const double *x, double *fvec, int /*iflag*/) +{ +/* subroutine fcn for hybrd1 example. */ + + int k; + double one=1, temp, temp1, temp2, three=3, two=2, zero=0; + + for (k=1; k <= n; k++) + { + temp = (three - two*x[k-1])*x[k-1]; + temp1 = zero; + if (k != 1) temp1 = x[k-1-1]; + temp2 = zero; + if (k != n) temp2 = x[k+1-1]; + fvec[k-1] = temp - temp1 - two*temp2 + one; + } + return 0; +} + +void testHybrd1() +{ + int j, n, info, lwa; + double tol, fnorm; + double x[9], fvec[9], wa[180]; + + n = 9; + +/* the following starting values provide a rough solution. */ + + for (j=1; j<=9; j++) + { + x[j-1] = -1.; + } + + lwa = 180; + +/* set tol to the square root of the machine precision. */ +/* unless high solutions are required, */ +/* this is the recommended setting. */ + + tol = sqrt(dpmpar(1)); + info = hybrd1(fcn_hybrd1, 0, n, x, fvec, tol, wa, lwa); + fnorm = enorm(n, fvec); + + VERIFY_IS_APPROX(fnorm, 1.192636e-08); + VERIFY(info==1); + + double x_ref[] = { + -0.5706545, -0.6816283, -0.7017325, + -0.7042129, -0.701369, -0.6918656, + -0.665792, -0.5960342, -0.4164121 + }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); +} + +int fcn_hybrd(void * /*p*/, int n, const double *x, double *fvec, int iflag) +{ + /* subroutine fcn for hybrd example. */ + + int k; + double one=1, temp, temp1, temp2, three=3, two=2, zero=0; + + if (iflag == 0) + { + /* insert print statements here when nprint is positive. */ + return 0; + } + for (k=1; k<=n; k++) + { + + temp = (three - two*x[k-1])*x[k-1]; + temp1 = zero; + if (k != 1) temp1 = x[k-1-1]; + temp2 = zero; + if (k != n) temp2 = x[k+1-1]; + fvec[k-1] = temp - temp1 - two*temp2 + one; + } + return 0; +} + +void testHybrd() +{ + int j, n, maxfev, ml, mu, mode, nprint, info, nfev, ldfjac, lr; + double xtol, epsfcn, factor, fnorm; + double x[9], fvec[9], diag[9], fjac[9*9], r[45], qtf[9], + wa1[9], wa2[9], wa3[9], wa4[9]; + + n = 9; + +/* the following starting values provide a rough solution. */ + + for (j=1; j<=9; j++) + { + x[j-1] = -1.; + } + + ldfjac = 9; + lr = 45; + +/* set xtol to the square root of the machine precision. */ +/* unless high solutions are required, */ +/* this is the recommended setting. */ + + xtol = sqrt(dpmpar(1)); + + maxfev = 2000; + ml = 1; + mu = 1; + epsfcn = 0.; + mode = 2; + for (j=1; j<=9; j++) + { + diag[j-1] = 1.; + } + + factor = 1.e2; + nprint = 0; + + info = hybrd(fcn_hybrd, 0, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, + diag, mode, factor, nprint, &nfev, + fjac, ldfjac, r, lr, qtf, wa1, wa2, wa3, wa4); + fnorm = enorm(n, fvec); + + VERIFY_IS_APPROX(fnorm, 1.192636e-08); + VERIFY(nfev==14); + VERIFY(info==1); + double x_ref[] = { + -0.5706545, -0.6816283, -0.7017325, + -0.7042129, -0.701369, -0.6918656, + -0.665792, -0.5960342, -0.4164121 + }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); +} + +int fcn_lmstr1(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjrow, int iflag) +{ + /* subroutine fcn for lmstr1 example. */ + int i; + double tmp1, tmp2, tmp3, tmp4; + double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, + 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; + + if (iflag < 2) + { + for (i=1; i<=15; i++) + { + tmp1=i; + tmp2 = 16-i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3)); + } + } else { + i = iflag - 1; + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4=tmp4*tmp4; + fjrow[1-1] = -1; + fjrow[2-1] = tmp1*tmp2/tmp4; + fjrow[3-1] = tmp1*tmp3/tmp4; + } + return 0; +} + +void testLmstr1() +{ + int m, n, ldfjac, info, lwa, ipvt[3]; + double tol, fnorm; + double x[30], fvec[15], fjac[9], wa[30]; + + m = 15; + n = 3; + + /* the following starting values provide a rough fit. */ + + x[0] = 1.; + x[1] = 1.; + x[2] = 1.; + + ldfjac = 3; + lwa = 30; + + /* set tol to the square root of the machine precision. + unless high precision solutions are required, + this is the recommended setting. */ + + tol = sqrt(dpmpar(1)); + + info = lmstr1(fcn_lmstr1, 0, m, n, + x, fvec, fjac, ldfjac, + tol, ipvt, wa, lwa); + + fnorm = enorm(m, fvec); + + VERIFY_IS_APPROX(fnorm, 0.09063596); + VERIFY(info==1); + double x_ref[] = {0.08241058, 1.133037, 2.343695 }; + for (m=1; m<=3; m++) VERIFY_IS_APPROX(x[m-1], x_ref[m-1]); +} + +int fcn_lmstr(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjrow, int iflag) +{ + + /* subroutine fcn for lmstr example. */ + + int i; + double tmp1, tmp2, tmp3, tmp4; + double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, + 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; + + if (iflag == 0) + { + /* insert print statements here when nprint is positive. */ + return 0; + } + if (iflag < 2) + { for (i = 1; i <= 15; i++) { tmp1 = i; tmp2 = 16 - i; - - /* error introduced into next statement for illustration. */ - /* corrected statement should read tmp3 = tmp1 . */ - - tmp3 = tmp2; + tmp3 = tmp1; if (i > 8) tmp3 = tmp2; - tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4=tmp4*tmp4; - fjac[i-1+ ldfjac*(1-1)] = -1.; - fjac[i-1+ ldfjac*(2-1)] = tmp1*tmp2/tmp4; - fjac[i-1+ ldfjac*(3-1)] = tmp1*tmp3/tmp4; - } + fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3)); +} + } +else + { + i = iflag - 1; + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4 = tmp4*tmp4; + fjrow[1-1] = -1.; + fjrow[2-1] = tmp1*tmp2/tmp4; + fjrow[3-1] = tmp1*tmp3/tmp4; + } + return 0; +} + +void testLmstr() +{ + int j, m, n, ldfjac, maxfev, mode, nprint, info, nfev, njev; + int ipvt[3]; + double ftol, xtol, gtol, factor, fnorm; + double x[3], fvec[15], fjac[3*3], diag[3], qtf[3], + wa1[3], wa2[3], wa3[3], wa4[15]; + + m = 15; + n = 3; + + /* the following starting values provide a rough fit. */ + + x[1-1] = 1.; + x[2-1] = 1.; + x[3-1] = 1.; + + ldfjac = 3; + + /* set ftol and xtol to the square root of the machine */ + /* and gtol to zero. unless high solutions are */ + /* required, these are the recommended settings. */ + + ftol = sqrt(dpmpar(1)); + xtol = sqrt(dpmpar(1)); + gtol = 0.; + + maxfev = 400; + mode = 1; + factor = 1.e2; + nprint = 0; + + info = lmstr(fcn_lmstr, 0, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, + maxfev, diag, mode, factor, nprint, &nfev, &njev, + ipvt, qtf, wa1, wa2, wa3, wa4); + fnorm = enorm(m, fvec); + + VERIFY_IS_APPROX(fnorm, 0.09063596); + VERIFY(nfev==6); + VERIFY(njev==5); + VERIFY(info==1); + + double x_ref[] = {0.08241058, 1.133037, 2.343695 }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); +} + +int fcn_lmdif1(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, int /*iflag*/) +{ + /* function fcn for lmdif1 example */ + + int i; + double tmp1,tmp2,tmp3; + double y[15]={1.4e-1,1.8e-1,2.2e-1,2.5e-1,2.9e-1,3.2e-1,3.5e-1,3.9e-1, + 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; + + for (i=0; i<15; i++) + { + tmp1 = i+1; + tmp2 = 15 - i; + tmp3 = tmp1; + + if (i >= 8) tmp3 = tmp2; + fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); + } + return 0; +} + +void testLmdif1() +{ + int m, n, info, lwa, iwa[3]; + double tol, fnorm, x[3], fvec[15], wa[75]; + + m = 15; + n = 3; + + /* the following starting values provide a rough fit. */ + + x[0] = 1.e0; + x[1] = 1.e0; + x[2] = 1.e0; + + lwa = 75; + + /* set tol to the square root of the machine precision. unless high + precision solutions are required, this is the recommended + setting. */ + + tol = sqrt(dpmpar(1)); + + info = lmdif1(fcn_lmdif1, 0, m, n, x, fvec, tol, iwa, wa, lwa); + + fnorm = enorm(m, fvec); + + VERIFY_IS_APPROX(fnorm, 0.09063596); + VERIFY(info==1); + double x_ref[] = {0.0824106, 1.1330366, 2.3436947 }; + int j; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); + +} + +int fcn_lmdif(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, int iflag) +{ + +/* subroutine fcn for lmdif example. */ + + int i; + double tmp1, tmp2, tmp3; + double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, + 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; + + if (iflag == 0) + { + /* insert print statements here when nprint is positive. */ + return 0; + } + for (i = 1; i <= 15; i++) + { + tmp1 = i; + tmp2 = 16 - i; + tmp3 = tmp1; + if (i > 8) tmp3 = tmp2; + fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3)); } return 0; } +void testLmdif() +{ + int i, j, m, n, maxfev, mode, nprint, info, nfev, ldfjac; + int ipvt[3]; + double ftol, xtol, gtol, epsfcn, factor, fnorm; + double x[3], fvec[15], diag[3], fjac[15*3], qtf[3], + wa1[3], wa2[3], wa3[3], wa4[15]; + double covfac; + + m = 15; + n = 3; + +/* the following starting values provide a rough fit. */ + + x[1-1] = 1.; + x[2-1] = 1.; + x[3-1] = 1.; + + ldfjac = 15; + + /* set ftol and xtol to the square root of the machine */ + /* and gtol to zero. unless high solutions are */ + /* required, these are the recommended settings. */ + + ftol = sqrt(dpmpar(1)); + xtol = sqrt(dpmpar(1)); + gtol = 0.; + + maxfev = 800; + epsfcn = 0.; + mode = 1; + factor = 1.e2; + nprint = 0; + + info = lmdif(fcn_lmdif, 0, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, + diag, mode, factor, nprint, &nfev, fjac, ldfjac, + ipvt, qtf, wa1, wa2, wa3, wa4); + + fnorm = enorm(m, fvec); + + VERIFY_IS_APPROX(fnorm, 0.09063596); + VERIFY(nfev==21); + VERIFY(info==1); + + double x_ref[] = {0.08241058, 1.133037, 2.343695 }; + for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); + + ftol = dpmpar(1); + covfac = fnorm*fnorm/(m-n); + covar(n, fjac, ldfjac, ipvt, ftol, wa1); + + double cov_ref[] = { + 0.0001531202, 0.002869942, -0.002656662, + 0.002869942, 0.09480937, -0.09098997, + -0.002656662, -0.09098997, 0.08778729 + }; + + for (i=1; i<=n; i++) + for (j=1; j<=n; j++) + VERIFY_IS_APPROX(fjac[(i-1)*ldfjac+j-1]*covfac, cov_ref[(i-1)*3+(j-1)]); +} void test_NonLinear() { CALL_SUBTEST(testChkder()); + CALL_SUBTEST(testLmder1()); + CALL_SUBTEST(testLmder()); + CALL_SUBTEST(testHybrj1()); + CALL_SUBTEST(testHybrj()); + CALL_SUBTEST(testHybrd1()); + CALL_SUBTEST(testHybrd()); + CALL_SUBTEST(testLmstr1()); + CALL_SUBTEST(testLmstr()); + CALL_SUBTEST(testLmdif1()); + CALL_SUBTEST(testLmdif()); } |