1 subroutine chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err)
94 integer(kind=4) ldfjac
103 real(kind=8) fjac(ldfjac, n)
105 real(kind=8) fvecp(m)
113 epsmch = epsilon(epsmch)
122 if (temp == 0.0d+00)
then 130 else if (mode == 2)
then 132 epsf = 100.0d+00*epsmch
139 if (temp == 0.0d+00)
then 142 err(1:m) = err(1:m) + temp*fjac(1:m, j)
149 if (fvec(i) /= 0.0d+00 .and. fvecp(i) /= 0.0d+00 .and. &
150 abs(fvecp(i) - fvec(i)) >= epsf*abs(fvec(i)))
then 151 temp = eps*abs((fvecp(i) - fvec(i))/eps - err(i)) &
152 /(abs(fvec(i)) + abs(fvecp(i)))
157 if (epsmch < temp .and. temp < eps)
then 158 err(i) = (log10(temp) - epslog)/epslog
161 if (eps <= temp)
then 171 subroutine dogleg(n, r, lr, diag, qtb, delta, x)
260 epsmch = epsilon(epsmch)
264 jj = (n*(n + 1))/2 + 1
274 sum2 = sum2 + r(l)*x(i)
280 if (temp == 0.0d+00)
then 284 temp = max(temp, abs(r(l)))
288 if (temp == 0.0d+00)
then 296 x(j) = (qtb(j) - sum2)/temp
303 wa2(1:n) = diag(1:n)*x(1:n)
304 qnorm = enorm(n, wa2)
306 if (qnorm <= delta)
then 317 wa1(i) = wa1(i) + r(l)*temp
320 wa1(j) = wa1(j)/diag(j)
326 gnorm = enorm(n, wa1)
330 if (gnorm /= 0.0d+00)
then 334 wa1(1:n) = (wa1(1:n)/gnorm)/diag(1:n)
340 sum2 = sum2 + r(l)*wa1(i)
347 sgnorm = (gnorm/temp)/temp
356 if (sgnorm < delta)
then 358 bnorm = enorm(n, qtb)
359 temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
360 temp = temp - (delta/qnorm)*(sgnorm/delta)**2 &
361 + sqrt((temp - (delta/qnorm))**2 &
362 + (1.0d+00 - (delta/qnorm)**2) &
363 *(1.0d+00 - (sgnorm/delta)**2))
365 alpha = ((delta/qnorm)*(1.0d+00 - (sgnorm/delta)**2)) &
375 temp = (1.0d+00 - alpha)*min(sgnorm, delta)
377 x(1:n) = temp*wa1(1:n) + alpha*x(1:n)
426 enorm = sqrt(sum(x(1:n)**2))
498 rdwarf = sqrt(tiny(rdwarf))
499 rgiant = sqrt(huge(rgiant))
506 agiant = rgiant/
real(n, kind=8)
512 if (xabs <= rdwarf)
then 514 if (x3max < xabs)
then 515 s3 = 1.0d+00 + s3*(x3max/xabs)**2
517 else if (xabs /= 0.0d+00)
then 518 s3 = s3 + (xabs/x3max)**2
521 else if (agiant <= xabs)
then 523 if (x1max < xabs)
then 524 s1 = 1.0d+00 + s1*(x1max/xabs)**2
527 s1 = s1 + (xabs/x1max)**2
540 if (s1 /= 0.0d+00)
then 542 enorm2 = x1max*sqrt(s1 + (s2/x1max)/x1max)
544 else if (s2 /= 0.0d+00)
then 546 if (x3max <= s2)
then 547 enorm2 = sqrt(s2*(1.0d+00 + (x3max/s2)*(x3max*s3)))
549 enorm2 = sqrt(x3max*((s2/x3max) + (x3max*s3)))
554 enorm2 = x3max*sqrt(s3)
560 subroutine fdjac1(fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, m, prms)
639 integer(kind=4) ldfjac
647 real(kind=8) fjac(ldfjac, n)
651 integer(kind=4) iflag
663 epsmch = epsilon(epsmch)
665 eps = sqrt(max(epsfcn, epsmch))
676 if (h == 0.0d+00)
then 682 call fcn(n, x, wa1, iflag, m, prms)
689 fjac(1:n, j) = (wa1(1:n) - fvec(1:n))/h
702 if (h == 0.0d+00)
then 709 call fcn(n, x, wa1, iflag, m, prms)
720 if (h == 0.0d+00)
then 724 fjac(1:n, j) = 0.0d+00
727 if (j - mu <= i .and. i <= j + ml)
then 728 fjac(i, j) = (wa1(i) - fvec(i))/h
740 subroutine fdjac2(fcn, m, n, x, xdat, ydat, fvec, fjac, ldfjac, iflag, epsfcn)
816 integer(kind=4) ldfjac
824 real(kind=8) fjac(ldfjac, n)
825 real(kind=8) fvec(m), xdat(m), ydat(m)
828 integer(kind=4) iflag
834 epsmch = epsilon(epsmch)
836 eps = sqrt(max(epsfcn, epsmch))
842 if (h == 0.0d+00)
then 848 call fcn(m, n, x, xdat, ydat, wa, iflag)
855 fjac(1:m, j) = (wa(1:m) - fvec(1:m))/h
862 subroutine hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
863 factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, m, prms)
991 integer(kind=4) ldfjac
1004 real(kind=8) fjac(ldfjac, n)
1007 real(kind=8) fvec(n)
1009 integer(kind=4) iflag
1010 integer(kind=4) info
1011 integer(kind=4) iter
1012 integer(kind=4) iwa(1)
1016 integer(kind=4) maxfev
1018 integer(kind=4) mode
1019 integer(kind=4) msum
1021 integer(kind=4) ncfail
1022 integer(kind=4) nslow1
1023 integer(kind=4) nslow2
1024 integer(kind=4) ncsuc
1025 integer(kind=4) nfev
1026 integer(kind=4) nprint
1043 real(kind=8) prms(m)
1045 epsmch = epsilon(epsmch)
1055 else if (xtol < 0.0d+00)
then 1057 else if (maxfev <= 0)
then 1059 else if (ml < 0)
then 1061 else if (mu < 0)
then 1063 else if (factor <= 0.0d+00)
then 1065 else if (ldfjac < n)
then 1067 else if (lr < (n*(n + 1))/2)
then 1074 if (diag(j) <= 0.0d+00)
then 1085 call fcn(n, x, fvec, iflag, m, prms)
1092 fnorm = enorm(n, fvec)
1096 msum = min(ml + mu + 1, n)
1115 call fdjac1(fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, m, prms)
1126 call qrfac(n, n, fjac, ldfjac, pivot, iwa, 1, wa1, wa2)
1135 diag(1:n) = wa2(1:n)
1137 if (wa2(j) == 0.0d+00)
then 1147 wa3(1:n) = diag(1:n)*x(1:n)
1148 xnorm = enorm(n, wa3)
1149 delta = factor*xnorm
1150 if (delta == 0.0d+00)
then 1158 qtf(1:n) = fvec(1:n)
1162 if (fjac(j, j) /= 0.0d+00)
then 1163 temp = -dot_product(qtf(j:n), fjac(j:n, j))/fjac(j, j)
1164 qtf(j:n) = qtf(j:n) + fjac(j:n, j)*temp
1180 if (wa1(j) == 0.0d+00)
then 1187 call qform(n, n, fjac, ldfjac)
1193 diag(j) = max(diag(j), wa2(j))
1203 if (0 < nprint)
then 1205 if (mod(iter - 1, nprint) == 0)
then 1206 call fcn(n, x, fvec, iflag, m, prms)
1215 call dogleg(n, r, lr, diag, qtf, delta, wa1)
1220 wa1(1:n) = -wa1(1:n)
1221 wa2(1:n) = x(1:n) + wa1(1:n)
1222 wa3(1:n) = diag(1:n)*wa1(1:n)
1224 pnorm = enorm(n, wa3)
1229 delta = min(delta, pnorm)
1235 call fcn(n, wa2, wa4, iflag, m, prms)
1242 fnorm1 = enorm(n, wa4)
1247 if (fnorm1 < fnorm)
then 1248 actred = 1.0d+00 - (fnorm1/fnorm)**2
1257 sum2 = sum2 + r(l)*wa1(j)
1260 wa3(i) = qtf(i) + sum2
1263 temp = enorm(n, wa3)
1265 if (temp < fnorm)
then 1266 prered = 1.0d+00 - (temp/fnorm)**2
1272 if (0.0d+00 < prered)
then 1273 ratio = actred/prered
1278 if (ratio < 0.1d+00)
then 1282 delta = 0.5d+00*delta
1289 if (0.5d+00 <= ratio .or. 1 < ncsuc)
then 1290 delta = max(delta, pnorm/0.5d+00)
1293 if (abs(ratio - 1.0d+00) <= 0.1d+00)
then 1294 delta = pnorm/0.5d+00
1304 if (0.0001d+00 <= ratio)
then 1306 wa2(1:n) = diag(1:n)*x(1:n)
1307 fvec(1:n) = wa4(1:n)
1308 xnorm = enorm(n, wa2)
1316 if (0.001d+00 <= actred)
then 1324 if (0.1d+00 <= actred)
then 1330 if (delta <= xtol*xnorm .or. fnorm == 0.0d+00)
then 1340 if (maxfev <= nfev)
then 1344 if (0.1d+00*max(0.1d+00*delta, pnorm) <= epsmch*xnorm)
then 1348 if (nslow2 == 5)
then 1352 if (nslow1 == 10)
then 1363 if (ncfail == 2)
then 1371 sum2 = dot_product(wa4(1:n), fjac(1:n, j))
1372 wa2(j) = (sum2 - wa3(j))/pnorm
1373 wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
1374 if (0.0001d+00 <= ratio)
then 1381 call r1updt(n, n, r, lr, wa1, wa2, wa3, sing)
1382 call r1mpyq(n, n, fjac, ldfjac, wa2, wa3)
1383 call r1mpyq(1, n, qtf, 1, wa2, wa3)
1406 if (0 < nprint)
then 1407 call fcn(n, x, fvec, iflag, m, prms)
1411 end subroutine hybrd 1413 subroutine hybrd1(fcn, n, x, fvec, tol, info, m, prms)
1499 real(kind=8) diag(n)
1503 real(kind=8) fjac(n, n)
1504 real(kind=8) fvec(n)
1505 integer(kind=4) info
1507 integer(kind=4) ldfjac
1509 integer(kind=4) maxfev
1511 integer(kind=4) mode
1513 integer(kind=4) nfev
1514 integer(kind=4) nprint
1516 real(kind=8) r((n*(n + 1))/2)
1520 real(kind=8) prms(m)
1527 if (tol < 0.0d+00)
then 1533 maxfev = 200*(n + 1)
1543 fjac(1:n, 1:n) = 0.0d+00
1545 r(1:(n*(n + 1))/2) = 0.0d+00
1549 call hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
1550 factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, m, prms)
1558 subroutine hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, &
1559 factor, nprint, info, nfev, njev, r, lr, qtf)
1683 integer(kind=4) ldfjac
1689 real(kind=8) diag(n)
1694 real(kind=8) fjac(ldfjac, n)
1697 real(kind=8) fvec(n)
1699 integer(kind=4) iflag
1700 integer(kind=4) info
1701 integer(kind=4) iter
1702 integer(kind=4) iwa(1)
1706 integer(kind=4) maxfev
1707 integer(kind=4) mode
1708 integer(kind=4) ncfail
1709 integer(kind=4) nslow1
1710 integer(kind=4) nslow2
1711 integer(kind=4) ncsuc
1712 integer(kind=4) nfev
1713 integer(kind=4) njev
1714 integer(kind=4) nprint
1732 epsmch = epsilon(epsmch)
1746 if (0 < nprint)
then 1747 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1752 if (ldfjac < n .or. &
1753 xtol < 0.0d+00 .or. &
1755 factor <= 0.0d+00 .or. &
1756 lr < (n*(n + 1))/2)
then 1761 if (0 < nprint)
then 1762 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1769 if (diag(j) <= 0.0d+00)
then 1774 if (0 < nprint)
then 1775 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1785 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1793 if (0 < nprint)
then 1794 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1799 fnorm = enorm(n, fvec)
1818 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1824 if (0 < nprint)
then 1825 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1833 call qrfac(n, n, fjac, ldfjac, pivot, iwa, 1, wa1, wa2)
1841 diag(1:n) = wa2(1:n)
1843 if (wa2(j) == 0.0d+00)
then 1852 wa3(1:n) = diag(1:n)*x(1:n)
1853 xnorm = enorm(n, wa3)
1854 delta = factor*xnorm
1855 if (delta == 0.0d+00)
then 1863 qtf(1:n) = fvec(1:n)
1866 if (fjac(j, j) /= 0.0d+00)
then 1869 sum2 = sum2 + fjac(i, j)*qtf(i)
1871 temp = -sum2/fjac(j, j)
1873 qtf(i) = qtf(i) + fjac(i, j)*temp
1889 if (wa1(j) == 0.0d+00)
then 1896 call qform(n, n, fjac, ldfjac)
1902 diag(j) = max(diag(j), wa2(j))
1912 if (0 < nprint)
then 1915 if (mod(iter - 1, nprint) == 0)
then 1916 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1922 if (0 < nprint)
then 1923 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1932 call dogleg(n, r, lr, diag, qtf, delta, wa1)
1937 wa1(1:n) = -wa1(1:n)
1938 wa2(1:n) = x(1:n) + wa1(1:n)
1939 wa3(1:n) = diag(1:n)*wa1(1:n)
1940 pnorm = enorm(n, wa3)
1945 delta = min(delta, pnorm)
1951 call fcn(n, wa2, wa4, fjac, ldfjac, iflag)
1957 if (0 < nprint)
then 1958 call fcn(n, x, fvec, fjac, ldfjac, iflag)
1963 fnorm1 = enorm(n, wa4)
1968 if (fnorm1 < fnorm)
then 1969 actred = 1.0d+00 - (fnorm1/fnorm)**2
1978 sum2 = sum2 + r(l)*wa1(j)
1981 wa3(i) = qtf(i) + sum2
1984 temp = enorm(n, wa3)
1986 if (temp < fnorm)
then 1987 prered = 1.0d+00 - (temp/fnorm)**2
1992 if (0.0d+00 < prered)
then 1993 ratio = actred/prered
2000 if (ratio < 0.1d+00)
then 2004 delta = 0.5d+00*delta
2011 if (0.5d+00 <= ratio .or. 1 < ncsuc)
then 2012 delta = max(delta, pnorm/0.5d+00)
2015 if (abs(ratio - 1.0d+00) <= 0.1d+00)
then 2016 delta = pnorm/0.5d+00
2028 if (0.0001d+00 <= ratio)
then 2030 wa2(1:n) = diag(1:n)*x(1:n)
2031 fvec(1:n) = wa4(1:n)
2032 xnorm = enorm(n, wa2)
2040 if (0.001d+00 <= actred)
then 2048 if (0.1d+00 <= actred)
then 2054 if (delta <= xtol*xnorm .or. fnorm == 0.0d+00)
then 2060 if (0 < nprint)
then 2061 call fcn(n, x, fvec, fjac, ldfjac, iflag)
2068 if (maxfev <= nfev)
then 2072 if (0.1d+00*max(0.1d+00*delta, pnorm) <= epsmch*xnorm)
then 2076 if (nslow2 == 5)
then 2080 if (nslow1 == 10)
then 2086 if (0 < nprint)
then 2087 call fcn(n, x, fvec, fjac, ldfjac, iflag)
2094 if (ncfail == 2)
then 2102 sum2 = dot_product(wa4(1:n), fjac(1:n, j))
2103 wa2(j) = (sum2 - wa3(j))/pnorm
2104 wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
2105 if (0.0001d+00 <= ratio)
then 2112 call r1updt(n, n, r, lr, wa1, wa2, wa3, sing)
2113 call r1mpyq(n, n, fjac, ldfjac, wa2, wa3)
2114 call r1mpyq(1, n, qtf, 1, wa2, wa3)
2126 end subroutine hybrj 2127 subroutine hybrj1(fcn, n, x, fvec, fjac, ldfjac, tol, info)
2213 integer(kind=4) ldfjac
2216 real(kind=8) diag(n)
2219 real(kind=8) fjac(ldfjac, n)
2220 real(kind=8) fvec(n)
2221 integer(kind=4) info
2224 integer(kind=4) maxfev
2225 integer(kind=4) mode
2226 integer(kind=4) nfev
2227 integer(kind=4) njev
2228 integer(kind=4) nprint
2230 real(kind=8) r((n*(n + 1))/2)
2239 else if (ldfjac < n)
then 2241 else if (tol < 0.0d+00)
then 2245 maxfev = 100*(n + 1)
2253 call hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, &
2254 factor, nprint, info, nfev, njev, r, lr, qtf)
2262 subroutine lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
2263 diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
2411 integer(kind=4) ldfjac
2417 real(kind=8) diag(n)
2423 real(kind=8) fjac(ldfjac, n)
2427 real(kind=8) fvec(m)
2431 integer(kind=4) iflag
2432 integer(kind=4) info
2433 integer(kind=4) ipvt(n)
2434 integer(kind=4) iter
2437 integer(kind=4) maxfev
2438 integer(kind=4) mode
2439 integer(kind=4) nfev
2440 integer(kind=4) njev
2441 integer(kind=4) nprint
2460 epsmch = epsilon(epsmch)
2478 .or. ftol < 0.0d+00 .or. xtol < 0.0d+00 .or. gtol < 0.0d+00 &
2479 .or. maxfev <= 0 .or. factor <= 0.0d+00)
then 2485 if (diag(j) <= 0.0d+00)
then 2494 call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2500 fnorm = enorm(m, fvec)
2514 call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2524 if (0 < nprint)
then 2526 if (mod(iter - 1, nprint) == 0)
then 2527 call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2537 call qrfac(m, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2)
2545 diag(1:n) = wa2(1:n)
2547 if (wa2(j) == 0.0d+00)
then 2556 wa3(1:n) = diag(1:n)*x(1:n)
2558 xnorm = enorm(n, wa3)
2559 delta = factor*xnorm
2560 if (delta == 0.0d+00)
then 2568 wa4(1:m) = fvec(1:m)
2572 if (fjac(j, j) /= 0.0d+00)
then 2573 sum2 = dot_product(wa4(j:m), fjac(j:m, j))
2574 temp = -sum2/fjac(j, j)
2575 wa4(j:m) = wa4(j:m) + fjac(j:m, j)*temp
2587 if (fnorm /= 0.0d+00)
then 2591 if (wa2(l) /= 0.0d+00)
then 2592 sum2 = dot_product(qtf(1:j), fjac(1:j, j))/fnorm
2593 gnorm = max(gnorm, abs(sum2/wa2(l)))
2601 if (gnorm <= gtol)
then 2610 diag(j) = max(diag(j), wa2(j))
2620 call lmpar(n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2)
2624 wa1(1:n) = -wa1(1:n)
2625 wa2(1:n) = x(1:n) + wa1(1:n)
2626 wa3(1:n) = diag(1:n)*wa1(1:n)
2628 pnorm = enorm(n, wa3)
2633 delta = min(delta, pnorm)
2639 call fcn(m, n, wa2, wa4, fjac, ldfjac, iflag)
2647 fnorm1 = enorm(m, wa4)
2652 if (0.1d+00*fnorm1 < fnorm)
then 2653 actred = 1.0d+00 - (fnorm1/fnorm)**2
2663 wa3(1:j) = wa3(1:j) + fjac(1:j, j)*temp
2666 temp1 = enorm(n, wa3)/fnorm
2667 temp2 = (sqrt(par)*pnorm)/fnorm
2668 prered = temp1**2 + temp2**2/0.5d+00
2669 dirder = -(temp1**2 + temp2**2)
2673 if (prered /= 0.0d+00)
then 2674 ratio = actred/prered
2681 if (ratio <= 0.25d+00)
then 2683 if (0.0d+00 <= actred)
then 2687 if (actred < 0.0d+00)
then 2688 temp = 0.5d+00*dirder/(dirder + 0.5d+00*actred)
2691 if (0.1d+00*fnorm1 >= fnorm .or. temp < 0.1d+00)
then 2695 delta = temp*min(delta, pnorm/0.1d+00)
2700 if (par == 0.0d+00 .or. ratio >= 0.75d+00)
then 2701 delta = 2.0d+00*pnorm
2711 if (0.0001d+00 <= ratio)
then 2713 wa2(1:n) = diag(1:n)*x(1:n)
2714 fvec(1:m) = wa4(1:m)
2715 xnorm = enorm(n, wa2)
2722 if (abs(actred) <= ftol .and. &
2723 prered <= ftol .and. &
2724 0.5d+00*ratio <= 1.0d+00)
then 2728 if (delta <= xtol*xnorm)
then 2732 if (abs(actred) <= ftol .and. prered <= ftol &
2733 .and. 0.5d+00*ratio <= 1.0d+00 .and. info == 2)
then 2743 if (nfev >= maxfev)
then 2747 if (abs(actred) <= epsmch .and. prered <= epsmch &
2748 .and. 0.5d+00*ratio <= 1.0d+00)
then 2752 if (delta <= epsmch*xnorm)
then 2756 if (gnorm <= epsmch)
then 2766 if (ratio < 0.0001d+00)
then 2784 if (0 < nprint)
then 2785 call fcn(m, n, x, fvec, fjac, ldfjac, iflag)
2789 end subroutine lmder 2790 subroutine lmder1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info)
2890 integer(kind=4) ldfjac
2894 real(kind=8) diag(n)
2897 real(kind=8) fjac(ldfjac, n)
2899 real(kind=8) fvec(m)
2901 integer(kind=4) info
2902 integer(kind=4) ipvt(n)
2903 integer(kind=4) maxfev
2904 integer(kind=4) mode
2905 integer(kind=4) nfev
2906 integer(kind=4) njev
2907 integer(kind=4) nprint
2917 else if (m < n)
then 2919 else if (ldfjac < m)
then 2921 else if (tol < 0.0d+00)
then 2926 maxfev = 100*(n + 1)
2933 call lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
2934 diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
2942 subroutine lmdif(fcn, m, n, x, xdat, ydat, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
2943 diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf)
3097 integer(kind=4) ldfjac
3103 real(kind=8) diag(n)
3110 real(kind=8) fjac(ldfjac, n)
3114 real(kind=8) :: fvec(m), xdat(m), ydat(m)
3118 integer(kind=4) iflag
3119 integer(kind=4) iter
3120 integer(kind=4) info
3121 integer(kind=4) ipvt(n)
3124 integer(kind=4) maxfev
3125 integer(kind=4) mode
3126 integer(kind=4) nfev
3127 integer(kind=4) nprint
3146 epsmch = epsilon(epsmch)
3154 else if (m < n)
then 3156 else if (ldfjac < m)
then 3158 else if (ftol < 0.0d+00)
then 3160 else if (xtol < 0.0d+00)
then 3162 else if (gtol < 0.0d+00)
then 3164 else if (maxfev <= 0)
then 3166 else if (factor <= 0.0d+00)
then 3172 if (diag(j) <= 0.0d+00)
then 3181 call fcn(m, n, x, xdat, ydat, fvec, iflag)
3188 fnorm = enorm(m, fvec)
3202 call fdjac2(fcn, m, n, x, xdat, ydat, fvec, fjac, ldfjac, iflag, epsfcn)
3211 if (0 < nprint)
then 3213 if (mod(iter - 1, nprint) == 0)
then 3214 call fcn(m, n, x, xdat, ydat, fvec, iflag)
3224 call qrfac(m, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2)
3232 diag(1:n) = wa2(1:n)
3234 if (wa2(j) == 0.0d+00)
then 3243 wa3(1:n) = diag(1:n)*x(1:n)
3244 xnorm = enorm(n, wa3)
3245 delta = factor*xnorm
3246 if (delta == 0.0d+00)
then 3253 wa4(1:m) = fvec(1:m)
3257 if (fjac(j, j) /= 0.0d+00)
then 3258 sum2 = dot_product(wa4(j:m), fjac(j:m, j))
3259 temp = -sum2/fjac(j, j)
3260 wa4(j:m) = wa4(j:m) + fjac(j:m, j)*temp
3272 if (fnorm /= 0.0d+00)
then 3278 if (wa2(l) /= 0.0d+00)
then 3281 sum2 = sum2 + fjac(i, j)*(qtf(i)/fnorm)
3283 gnorm = max(gnorm, abs(sum2/wa2(l)))
3292 if (gnorm <= gtol)
then 3301 diag(j) = max(diag(j), wa2(j))
3311 call lmpar(n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2)
3316 wa1(1:n) = -wa1(1:n)
3317 wa2(1:n) = x(1:n) + wa1(1:n)
3318 wa3(1:n) = diag(1:n)*wa1(1:n)
3320 pnorm = enorm(n, wa3)
3325 delta = min(delta, pnorm)
3331 call fcn(m, n, wa2, xdat, ydat, wa4, iflag)
3336 fnorm1 = enorm(m, wa4)
3340 if (0.1d+00*fnorm1 < fnorm)
then 3341 actred = 1.0d+00 - (fnorm1/fnorm)**2
3352 wa3(1:j) = wa3(1:j) + fjac(1:j, j)*temp
3355 temp1 = enorm(n, wa3)/fnorm
3356 temp2 = (sqrt(par)*pnorm)/fnorm
3357 prered = temp1**2 + temp2**2/0.5d+00
3358 dirder = -(temp1**2 + temp2**2)
3363 if (prered /= 0.0d+00)
then 3364 ratio = actred/prered
3369 if (ratio <= 0.25d+00)
then 3371 if (actred >= 0.0d+00)
then 3375 if (actred < 0.0d+00)
then 3376 temp = 0.5d+00*dirder/(dirder + 0.5d+00*actred)
3379 if (0.1d+00*fnorm1 >= fnorm .or. temp < 0.1d+00)
then 3383 delta = temp*min(delta, pnorm/0.1d+00)
3388 if (par == 0.0d+00 .or. ratio >= 0.75d+00)
then 3389 delta = 2.0d+00*pnorm
3401 if (0.0001d+00 <= ratio)
then 3403 wa2(1:n) = diag(1:n)*x(1:n)
3404 fvec(1:m) = wa4(1:m)
3405 xnorm = enorm(n, wa2)
3412 if (abs(actred) <= ftol .and. prered <= ftol &
3413 .and. 0.5d+00*ratio <= 1.0d+00)
then 3417 if (delta <= xtol*xnorm)
then 3421 if (abs(actred) <= ftol .and. prered <= ftol &
3422 .and. 0.5d+00*ratio <= 1.0d+00 .and. info == 2) info = 3
3430 if (maxfev <= nfev)
then 3434 if (abs(actred) <= epsmch .and. prered <= epsmch &
3435 .and. 0.5d+00*ratio <= 1.0d+00)
then 3439 if (delta <= epsmch*xnorm)
then 3443 if (gnorm <= epsmch)
then 3453 if (ratio < 0.0001d+00)
then 3471 if (0 < nprint)
then 3472 call fcn(m, n, x, xdat, ydat, fvec, iflag)
3476 end subroutine lmdif 3477 subroutine lmdif1(fcn, m, n, x, xdat, ydat, fvec, tol, info)
3568 real(kind=8) diag(n)
3572 real(kind=8) fjac(m, n)
3574 real(kind=8) ::fvec(m), xdat(m), ydat(m)
3576 integer(kind=4) info
3577 integer(kind=4) ipvt(n)
3578 integer(kind=4) ldfjac
3579 integer(kind=4) maxfev
3580 integer(kind=4) mode
3581 integer(kind=4) nfev
3582 integer(kind=4) nprint
3592 else if (m < n)
then 3594 else if (tol < 0.0d+00)
then 3599 maxfev = 200*(n + 1)
3608 call lmdif(fcn, m, n, x, xdat, ydat, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
3609 diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf)
3617 subroutine lmpar(n, r, ldr, ipvt, diag, qtb, delta, par, x, sdiag)
3721 real(kind=8) diag(n)
3728 integer(kind=4) ipvt(n)
3729 integer(kind=4) iter
3733 integer(kind=4) nsing
3740 real(kind=8) r(ldr, n)
3741 real(kind=8) sdiag(n)
3760 if (r(j, j) == 0.0d+00 .and. nsing == n)
then 3770 wa1(j) = wa1(j)/r(j, j)
3772 wa1(1:j - 1) = wa1(1:j - 1) - r(1:j - 1, j)*temp
3785 wa2(1:n) = diag(1:n)*x(1:n)
3786 dxnorm = enorm(n, wa2)
3789 if (fp <= 0.1d+00*delta)
then 3804 if (n <= nsing)
then 3808 wa1(j) = diag(l)*(wa2(l)/dxnorm)
3812 sum2 = dot_product(wa1(1:j - 1), r(1:j - 1, j))
3813 wa1(j) = (wa1(j) - sum2)/r(j, j)
3816 temp = enorm(n, wa1)
3817 parl = ((fp/delta)/temp)/temp
3824 sum2 = dot_product(qtb(1:j), r(1:j, j))
3826 wa1(j) = sum2/diag(l)
3829 gnorm = enorm(n, wa1)
3832 if (paru == 0.0d+00)
then 3833 paru = dwarf/min(delta, 0.1d+00)
3839 par = max(par, parl)
3840 par = min(par, paru)
3841 if (par == 0.0d+00)
then 3853 if (par == 0.0d+00)
then 3854 par = max(dwarf, 0.001d+00*paru)
3857 wa1(1:n) = sqrt(par)*diag(1:n)
3859 call qrsolv(n, r, ldr, ipvt, wa1, qtb, x, sdiag)
3861 wa2(1:n) = diag(1:n)*x(1:n)
3862 dxnorm = enorm(n, wa2)
3868 if (abs(fp) <= 0.1d+00*delta)
then 3875 if (parl == 0.0d+00 .and. fp <= temp .and. temp < 0.0d+00)
then 3877 else if (iter == 10)
then 3885 wa1(j) = diag(l)*(wa2(l)/dxnorm)
3889 wa1(j) = wa1(j)/sdiag(j)
3891 wa1(j + 1:n) = wa1(j + 1:n) - r(j + 1:n, j)*temp
3894 temp = enorm(n, wa1)
3895 parc = ((fp/delta)/temp)/temp
3899 if (0.0d+00 < fp)
then 3900 parl = max(parl, par)
3901 else if (fp < 0.0d+00)
then 3902 paru = min(paru, par)
3907 par = max(parl, par + parc)
3920 end subroutine lmpar 3921 subroutine lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
3922 diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
4073 integer(kind=4) ldfjac
4079 real(kind=8) diag(n)
4085 real(kind=8) fjac(ldfjac, n)
4089 real(kind=8) fvec(m)
4093 integer(kind=4) iflag
4094 integer(kind=4) info
4095 integer(kind=4) ipvt(n)
4096 integer(kind=4) iter
4099 integer(kind=4) maxfev
4100 integer(kind=4) mode
4101 integer(kind=4) nfev
4102 integer(kind=4) njev
4103 integer(kind=4) nprint
4123 epsmch = epsilon(epsmch)
4135 else if (m < n)
then 4137 else if (ldfjac < n)
then 4139 else if (ftol < 0.0d+00)
then 4141 else if (xtol < 0.0d+00)
then 4143 else if (gtol < 0.0d+00)
then 4145 else if (maxfev <= 0)
then 4147 else if (factor <= 0.0d+00)
then 4153 if (diag(j) <= 0.0d+00)
then 4162 call fcn(m, n, x, fvec, wa3, iflag)
4169 fnorm = enorm(m, fvec)
4182 if (0 < nprint)
then 4184 if (mod(iter - 1, nprint) == 0)
then 4185 call fcn(m, n, x, fvec, wa3, iflag)
4197 fjac(1:n, 1:n) = 0.0d+00
4201 call fcn(m, n, x, fvec, wa3, iflag)
4206 call rwupdt(n, fjac, ldfjac, wa3, qtf, temp, wa1, wa2)
4217 if (fjac(j, j) == 0.0d+00)
then 4221 wa2(j) = enorm(j, fjac(1, j))
4227 call qrfac(n, n, fjac, ldfjac, pivot, ipvt, n, wa1, wa2)
4231 if (fjac(j, j) /= 0.0d+00)
then 4233 sum2 = dot_product(qtf(j:n), fjac(j:n, j))
4234 temp = -sum2/fjac(j, j)
4235 qtf(j:n) = qtf(j:n) + fjac(j:n, j)*temp
4255 diag(1:n) = wa2(1:n)
4257 if (wa2(j) == 0.0d+00)
then 4264 wa3(1:n) = diag(1:n)*x(1:n)
4265 xnorm = enorm(n, wa3)
4266 delta = factor*xnorm
4267 if (delta == 0.0d+00)
then 4277 if (fnorm /= 0.0d+00)
then 4281 if (wa2(l) /= 0.0d+00)
then 4282 sum2 = dot_product(qtf(1:j), fjac(1:j, j))/fnorm
4283 gnorm = max(gnorm, abs(sum2/wa2(l)))
4291 if (gnorm <= gtol)
then 4300 diag(j) = max(diag(j), wa2(j))
4310 call lmpar(n, fjac, ldfjac, ipvt, diag, qtf, delta, par, wa1, wa2)
4315 wa1(1:n) = -wa1(1:n)
4316 wa2(1:n) = x(1:n) + wa1(1:n)
4317 wa3(1:n) = diag(1:n)*wa1(1:n)
4318 pnorm = enorm(n, wa3)
4323 delta = min(delta, pnorm)
4329 call fcn(m, n, wa2, wa4, wa3, iflag)
4334 fnorm1 = enorm(m, wa4)
4338 if (0.1d+00*fnorm1 < fnorm)
then 4339 actred = 1.0d+00 - (fnorm1/fnorm)**2
4351 wa3(1:j) = wa3(1:j) + fjac(1:j, j)*temp
4354 temp1 = enorm(n, wa3)/fnorm
4355 temp2 = (sqrt(par)*pnorm)/fnorm
4356 prered = temp1**2 + temp2**2/0.5d+00
4357 dirder = -(temp1**2 + temp2**2)
4361 if (prered /= 0.0d+00)
then 4362 ratio = actred/prered
4369 if (ratio <= 0.25d+00)
then 4371 if (actred >= 0.0d+00)
then 4374 temp = 0.5d+00*dirder/(dirder + 0.5d+00*actred)
4377 if (0.1d+00*fnorm1 >= fnorm .or. temp < 0.1d+00)
then 4381 delta = temp*min(delta, pnorm/0.1d+00)
4386 if (par == 0.0d+00 .or. ratio >= 0.75d+00)
then 4387 delta = pnorm/0.5d+00
4395 if (ratio >= 0.0001d+00)
then 4397 wa2(1:n) = diag(1:n)*x(1:n)
4398 fvec(1:m) = wa4(1:m)
4399 xnorm = enorm(n, wa2)
4406 if (abs(actred) <= ftol .and. prered <= ftol &
4407 .and. 0.5d+00*ratio <= 1.0d+00)
then 4411 if (delta <= xtol*xnorm)
then 4415 if (abs(actred) <= ftol .and. prered <= ftol &
4416 .and. 0.5d+00*ratio <= 1.0d+00 .and. info == 2)
then 4424 if (nfev >= maxfev)
then 4428 if (abs(actred) <= epsmch .and. prered <= epsmch &
4429 .and. 0.5d+00*ratio <= 1.0d+00)
then 4433 if (delta <= epsmch*xnorm)
then 4437 if (gnorm <= epsmch)
then 4447 if (ratio < 0.0001d+00)
then 4465 if (0 < nprint)
then 4466 call fcn(m, n, x, fvec, wa3, iflag)
4470 end subroutine lmstr 4471 subroutine lmstr1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info)
4574 integer(kind=4) ldfjac
4578 real(kind=8) diag(n)
4581 real(kind=8) fjac(ldfjac, n)
4583 real(kind=8) fvec(m)
4585 integer(kind=4) info
4586 integer(kind=4) ipvt(n)
4587 integer(kind=4) maxfev
4588 integer(kind=4) mode
4589 integer(kind=4) nfev
4590 integer(kind=4) njev
4591 integer(kind=4) nprint
4607 if (ldfjac < n)
then 4612 if (tol < 0.0d+00)
then 4618 fjac(1:ldfjac, 1:n) = 0.0d+00
4622 maxfev = 100*(n + 1)
4633 call lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
4634 diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
4642 subroutine qform(m, n, q, ldq)
4700 integer(kind=4) minmn
4701 real(kind=8) q(ldq, m)
4708 q(1:j - 1, j) = 0.0d+00
4713 q(1:m, n + 1:m) = 0.0d+00
4730 if (wa(k) /= 0.0d+00)
then 4733 temp = dot_product(wa(k:m), q(k:m, j))/wa(k)
4734 q(k:m, j) = q(k:m, j) - temp*wa(k:m)
4742 end subroutine qform 4743 subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
4819 integer(kind=4) lipvt
4823 real(kind=8) a(lda, n)
4824 real(kind=8) acnorm(n)
4829 integer(kind=4) i4_temp
4830 integer(kind=4) ipvt(lipvt)
4833 integer(kind=4) kmax
4834 integer(kind=4) minmn
4836 real(kind=8) r8_temp(m)
4837 real(kind=8) rdiag(n)
4841 epsmch = epsilon(epsmch)
4846 acnorm(j) = enorm(m, a(1:m, j))
4849 rdiag(1:n) = acnorm(1:n)
4850 wa(1:n) = acnorm(1:n)
4871 if (rdiag(kmax) < rdiag(k))
then 4878 r8_temp(1:m) = a(1:m, j)
4879 a(1:m, j) = a(1:m, kmax)
4880 a(1:m, kmax) = r8_temp(1:m)
4882 rdiag(kmax) = rdiag(j)
4886 ipvt(j) = ipvt(kmax)
4887 ipvt(kmax) = i4_temp
4896 ajnorm = enorm(m - j + 1, a(j, j))
4898 if (ajnorm /= 0.0d+00)
then 4900 if (a(j, j) < 0.0d+00)
then 4904 a(j:m, j) = a(j:m, j)/ajnorm
4905 a(j, j) = a(j, j) + 1.0d+00
4911 temp = dot_product(a(j:m, j), a(j:m, k))/a(j, j)
4913 a(j:m, k) = a(j:m, k) - temp*a(j:m, j)
4915 if (pivot .and. rdiag(k) /= 0.0d+00)
then 4917 temp = a(j, k)/rdiag(k)
4918 rdiag(k) = rdiag(k)*sqrt(max(0.0d+00, 1.0d+00 - temp**2))
4920 if (0.05d+00*(rdiag(k)/wa(k))**2 <= epsmch)
then 4921 rdiag(k) = enorm(m - j, a(j + 1, k))
4936 end subroutine qrfac 4937 subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag)
5028 real(kind=8) diag(n)
5030 integer(kind=4) ipvt(n)
5034 integer(kind=4) nsing
5037 real(kind=8) r(ldr, n)
5039 real(kind=8) sdiag(n)
5051 r(j:n, j) = r(j, j:n)
5066 if (diag(l) /= 0.0d+00)
then 5068 sdiag(j:n) = 0.0d+00
5082 if (sdiag(k) /= 0.0d+00)
then 5084 if (abs(r(k, k)) < abs(sdiag(k)))
then 5085 cotan = r(k, k)/sdiag(k)
5086 s = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5089 t = sdiag(k)/r(k, k)
5090 c = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*t**2)
5097 r(k, k) = c*r(k, k) + s*sdiag(k)
5098 temp = c*wa(k) + s*qtbpj
5099 qtbpj = -s*wa(k) + c*qtbpj
5105 temp = c*r(i, k) + s*sdiag(i)
5106 sdiag(i) = -s*r(i, k) + c*sdiag(i)
5131 if (sdiag(j) == 0.0d+00 .and. nsing == n)
then 5142 sum2 = dot_product(wa(j + 1:nsing), r(j + 1:nsing, j))
5143 wa(j) = (wa(j) - sum2)/sdiag(j)
5155 subroutine r1mpyq(m, n, a, lda, v, w)
5215 real(kind=8) a(lda, n)
5228 if (1.0d+00 < abs(v(j)))
then 5230 s = sqrt(1.0d+00 - c**2)
5233 c = sqrt(1.0d+00 - s**2)
5237 temp = c*a(i, j) - s*a(i, n)
5238 a(i, n) = s*a(i, j) + c*a(i, n)
5248 if (abs(w(j)) > 1.0d+00)
then 5250 s = sqrt(1.0d+00 - c**2)
5253 c = sqrt(1.0d+00 - s**2)
5257 temp = c*a(i, j) + s*a(i, n)
5258 a(i, n) = -s*a(i, j) + c*a(i, n)
5266 subroutine r1updt(m, n, s, ls, u, v, w, sing)
5366 jj = (n*(2*m - n + 1))/2 - (m - n)
5381 jj = jj - (m - j + 1)
5384 if (v(j) /= 0.0d+00)
then 5389 if (abs(v(n)) < abs(v(j)))
then 5391 sin = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5394 if (abs(cos)*giant > 1.0d+00)
then 5399 cos = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*tan**2)
5407 v(n) = sin*v(j) + cos*v(n)
5414 temp = cos*s(l) - sin*w(i)
5415 w(i) = sin*s(l) + cos*w(i)
5426 w(1:m) = w(1:m) + v(n)*u(1:m)
5434 if (w(j) /= 0.0d+00)
then 5439 if (abs(s(jj)) < abs(w(j)))
then 5442 sin = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5445 if (1.0d+00 < abs(cos)*giant)
then 5454 cos = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*tan**2)
5464 temp = cos*s(l) + sin*w(i)
5465 w(i) = -sin*s(l) + cos*w(i)
5478 if (s(jj) == 0.0d+00)
then 5482 jj = jj + (m - j + 1)
5494 if (s(jj) == 0.0d+00)
then 5536 character(len=*) title
5538 write (*,
'(a)')
' ' 5539 write (*,
'(a)') trim(title)
5540 write (*,
'(a)')
' ' 5542 write (*,
'(2x,i8,2x,g16.8)') i, a(i)
5547 subroutine rwupdt(n, r, ldr, w, b, alpha, c, s)
5624 real(kind=8) r(ldr, n)
5638 temp = c(i)*r(i, j) + s(i)*rowj
5639 rowj = -s(i)*r(i, j) + c(i)*rowj
5648 if (rowj /= 0.0d+00)
then 5650 if (abs(r(j, j)) < abs(rowj))
then 5651 cotan = r(j, j)/rowj
5652 s(j) = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*cotan**2)
5656 c(j) = 0.5d+00/sqrt(0.25d+00 + 0.25d+00*tan**2)
5662 r(j, j) = c(j)*r(j, j) + s(j)*rowj
5663 temp = c(j)*b(j) + s(j)*alpha
5664 alpha = -s(j)*b(j) + c(j)*alpha
5701 character(len=8) ampm
5706 character(len=9),
parameter,
dimension(12) :: month = (/ &
5707 'January ',
'February ',
'March ',
'April ', &
5708 'May ',
'June ',
'July ',
'August ', &
5709 'September',
'October ',
'November ',
'December '/)
5712 integer(kind=4) values(8)
5715 call date_and_time(values=values)
5727 else if (h == 12)
then 5728 if (n == 0 .and. s == 0)
then 5737 else if (h == 12)
then 5738 if (n == 0 .and. s == 0)
then 5746 write (*,
'(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') &
5747 d, trim(month(m)), y, h,
':', n,
':', s,
'.', mm, trim(ampm)
subroutine lmdif1(fcn, m, n, x, xdat, ydat, fvec, tol, info)
subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag)
subroutine lmder1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info)
subroutine dogleg(n, r, lr, diag, qtb, delta, x)
subroutine lmdif(fcn, m, n, x, xdat, ydat, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf)
subroutine lmstr1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info)
subroutine hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, m, prms)
real(kind=8) function enorm2(n, x)
subroutine rwupdt(n, r, ldr, w, b, alpha, c, s)
subroutine r8vec_print(n, a, title)
subroutine r1updt(m, n, s, ls, u, v, w, sing)
subroutine hybrj1(fcn, n, x, fvec, fjac, ldfjac, tol, info)
real(kind=8) function enorm(n, x)
subroutine lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
subroutine hybrd1(fcn, n, x, fvec, tol, info, m, prms)
subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm)
subroutine r1mpyq(m, n, a, lda, v, w)
subroutine fdjac1(fcn, n, x, fvec, fjac, ldfjac, iflag, ml, mu, epsfcn, m, prms)
subroutine qform(m, n, q, ldq)
subroutine hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, r, lr, qtf)
subroutine lmpar(n, r, ldr, ipvt, diag, qtb, delta, par, x, sdiag)
subroutine lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf)
subroutine chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err)
subroutine fdjac2(fcn, m, n, x, xdat, ydat, fvec, fjac, ldfjac, iflag, epsfcn)