Notebooks · lib/minpack/mgh

Test problems for least squares, from Moré, Garbow and Hillstrom

Eleven of the eighteen problems of MINPACK's own test driver for lmder (ref/file17.f: SSQFCN, SSQJAC and INITPT), ported with the Fortran's arithmetic in its order, so the Viv port of lmder can be compared with a gfortran build of the original iterate by iterate. The problems are those of Moré, Garbow and Hillstrom, "Testing Unconstrained Optimization Software", ACM TOMS 7 (1981) 17-41. The numbers are MINPACK's:

 4  Rosenbrock                   n  2, m  2
 5  helical valley               n  3, m  3
 6  Powell singular              n  4, m  4
 7  Freudenstein and Roth        n  2, m  2
 8  Bard                         n  3, m 15
 9  Kowalik and Osborne          n  4, m 11
10  Meyer                        n  3, m 16
12  Box three-dimensional        n  3, m 10
13  Jennrich and Sampson         n  2, m 10
17  Osborne 1                    n  5, m 33
18  Osborne 2                    n 11, m 65

mgh_fcn is the function lmder calls; its data array is [nprob, m, n] (mgh_data). mgh_start(nprob, n, factor) is the standard starting point times factor, as INITPT gives it.

Used with import minpack.mgh.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 31–339

# the measurements of the data-fitting problems, from SSQFCN's data statements
MGH_V = [4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625]
MGH_Y1 = [0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37, 0.58, 0.73, 0.96, 1.34, 2.1, 4.39]
MGH_Y2 = [0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246]
MGH_Y3 = [34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0,
          8261.0, 7030.0, 6005.0, 5147.0, 4427.0, 3820.0, 3307.0, 2872.0]
MGH_Y4 = [0.844, 0.908, 0.932, 0.936, 0.925, 0.908, 0.881, 0.85, 0.818, 0.784, 0.751,
          0.718, 0.685, 0.658, 0.628, 0.603, 0.58, 0.558, 0.538, 0.522, 0.506, 0.49,
          0.478, 0.467, 0.457, 0.448, 0.438, 0.431, 0.424, 0.42, 0.414, 0.411, 0.406]
MGH_Y5 = [1.366, 1.191, 1.112, 1.013, 0.991, 0.885, 0.831, 0.847, 0.786, 0.725, 0.746,
          0.679, 0.608, 0.655, 0.616, 0.606, 0.602, 0.626, 0.651, 0.724, 0.649, 0.649,
          0.694, 0.644, 0.624, 0.661, 0.612, 0.558, 0.533, 0.495, 0.5, 0.423, 0.395,
          0.375, 0.372, 0.391, 0.396, 0.405, 0.428, 0.429, 0.523, 0.562, 0.607, 0.653,
          0.672, 0.708, 0.633, 0.668, 0.645, 0.632, 0.591, 0.559, 0.597, 0.625, 0.739,
          0.71, 0.729, 0.72, 0.636, 0.581, 0.428, 0.292, 0.162, 0.098, 0.054]

#| The data array mgh_fcn reads: which problem, and its dimensions.
def mgh_data(nprob, m, n) {
    return [nprob, m, n]
}

#| The residuals of problem nprob at x, into fvec (SSQFCN).
def mgh_residuals(nprob, m, x: Arr, fvec: Arr) {
    if nprob == 4 {
        fvec[0] = 10.0 * (x[1] - x[0] * x[0])
        fvec[1] = 1.0 - x[0]
    } else if nprob == 5 {
        tpi = 8.0 * atan(1.0)
        # dsign(0.25, x(2)); only x(1) = 0 reaches it, and never with x(2) = -0
        tmp1 = 0.25
        if x[1] < 0.0 {
            tmp1 = -0.25
        }
        if x[0] > 0.0 {
            tmp1 = atan(x[1] / x[0]) / tpi
        }
        if x[0] < 0.0 {
            tmp1 = atan(x[1] / x[0]) / tpi + 0.5
        }
        tmp2 = sqrt(x[0] * x[0] + x[1] * x[1])
        fvec[0] = 10.0 * (x[2] - 10.0 * tmp1)
        fvec[1] = 10.0 * (tmp2 - 1.0)
        fvec[2] = x[2]
    } else if nprob == 6 {
        fvec[0] = x[0] + 10.0 * x[1]
        fvec[1] = sqrt(5.0) * (x[2] - x[3])
        d = x[1] - 2.0 * x[2]
        fvec[2] = d * d
        d = x[0] - x[3]
        fvec[3] = sqrt(10.0) * (d * d)
    } else if nprob == 7 {
        fvec[0] = -13.0 + x[0] + ((5.0 - x[1]) * x[1] - 2.0) * x[1]
        fvec[1] = -29.0 + x[0] + ((1.0 + x[1]) * x[1] - 14.0) * x[1]
    } else if nprob == 8 {
        for i to len(MGH_Y1) {
            tmp1 = i + 1.0
            tmp2 = 15.0 - i
            tmp3 = tmp1
            if i + 1.0 > 8.0 {
                tmp3 = tmp2
            }
            fvec[i] = MGH_Y1[i] - (x[0] + tmp1 / (x[1] * tmp2 + x[2] * tmp3))
        }
    } else if nprob == 9 {
        for i to len(MGH_V) {
            tmp1 = MGH_V[i] * (MGH_V[i] + x[1])
            tmp2 = MGH_V[i] * (MGH_V[i] + x[2]) + x[3]
            fvec[i] = MGH_Y2[i] - x[0] * tmp1 / tmp2
        }
    } else if nprob == 10 {
        for i to len(MGH_Y3) {
            temp = 5.0 * (i + 1.0) + 45.0 + x[2]
            tmp1 = x[1] / temp
            tmp2 = exp(tmp1)
            fvec[i] = x[0] * tmp2 - MGH_Y3[i]
        }
    } else if nprob == 12 {
        for i to m {
            temp = i + 1.0
            tmp1 = temp / 10.0
            fvec[i] = exp(-(tmp1 * x[0])) - exp(-(tmp1 * x[1])) + (exp(-temp) - exp(-tmp1)) * x[2]
        }
    } else if nprob == 13 {
        for i to m {
            temp = i + 1.0
            fvec[i] = 2.0 + 2.0 * temp - exp(temp * x[0]) - exp(temp * x[1])
        }
    } else if nprob == 17 {
        for i to len(MGH_Y4) {
            temp = 10.0 * i
            tmp1 = exp(-(x[3] * temp))
            tmp2 = exp(-(x[4] * temp))
            fvec[i] = MGH_Y4[i] - (x[0] + x[1] * tmp1 + x[2] * tmp2)
        }
    } else if nprob == 18 {
        for i to len(MGH_Y5) {
            temp = i / 10.0
            tmp1 = exp(-(x[4] * temp))
            d9 = temp - x[8]
            d10 = temp - x[9]
            d11 = temp - x[10]
            tmp2 = exp(-(x[5] * (d9 * d9)))
            tmp3 = exp(-(x[6] * (d10 * d10)))
            tmp4 = exp(-(x[7] * (d11 * d11)))
            fvec[i] = MGH_Y5[i] - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4)
        }
    }
    return 0
}

#| The Jacobian of problem nprob at x, into fjac (SSQJAC).
def mgh_jacobian(nprob, m, x: Arr, fjac: Arr) {
    if nprob == 4 {
        fjac[0, 0] = -(20.0 * x[0])
        fjac[0, 1] = 10.0
        fjac[1, 0] = -1.0
        fjac[1, 1] = 0.0
    } else if nprob == 5 {
        tpi = 8.0 * atan(1.0)
        temp = x[0] * x[0] + x[1] * x[1]
        tmp1 = tpi * temp
        tmp2 = sqrt(temp)
        fjac[0, 0] = 100.0 * x[1] / tmp1
        fjac[0, 1] = -(100.0 * x[0] / tmp1)
        fjac[0, 2] = 10.0
        fjac[1, 0] = 10.0 * x[0] / tmp2
        fjac[1, 1] = 10.0 * x[1] / tmp2
        fjac[1, 2] = 0.0
        fjac[2, 0] = 0.0
        fjac[2, 1] = 0.0
        fjac[2, 2] = 1.0
    } else if nprob == 6 {
        for j to len(x) {
            for i to m {
                fjac[i, j] = 0.0
            }
        }
        fjac[0, 0] = 1.0
        fjac[0, 1] = 10.0
        fjac[1, 2] = sqrt(5.0)
        fjac[1, 3] = -fjac[1, 2]
        fjac[2, 1] = 2.0 * (x[1] - 2.0 * x[2])
        fjac[2, 2] = -(2.0 * fjac[2, 1])
        fjac[3, 0] = 2.0 * sqrt(10.0) * (x[0] - x[3])
        fjac[3, 3] = -fjac[3, 0]
    } else if nprob == 7 {
        fjac[0, 0] = 1.0
        fjac[0, 1] = x[1] * (10.0 - 3.0 * x[1]) - 2.0
        fjac[1, 0] = 1.0
        fjac[1, 1] = x[1] * (2.0 + 3.0 * x[1]) - 14.0
    } else if nprob == 8 {
        for i to len(MGH_Y1) {
            tmp1 = i + 1.0
            tmp2 = 15.0 - i
            tmp3 = tmp1
            if i + 1.0 > 8.0 {
                tmp3 = tmp2
            }
            s = x[1] * tmp2 + x[2] * tmp3
            tmp4 = s * s
            fjac[i, 0] = -1.0
            fjac[i, 1] = tmp1 * tmp2 / tmp4
            fjac[i, 2] = tmp1 * tmp3 / tmp4
        }
    } else if nprob == 9 {
        for i to len(MGH_V) {
            tmp1 = MGH_V[i] * (MGH_V[i] + x[1])
            tmp2 = MGH_V[i] * (MGH_V[i] + x[2]) + x[3]
            fjac[i, 0] = -(tmp1 / tmp2)
            fjac[i, 1] = -(MGH_V[i] * x[0] / tmp2)
            fjac[i, 2] = fjac[i, 0] * fjac[i, 1]
            fjac[i, 3] = fjac[i, 2] / MGH_V[i]
        }
    } else if nprob == 10 {
        for i to len(MGH_Y3) {
            temp = 5.0 * (i + 1.0) + 45.0 + x[2]
            tmp1 = x[1] / temp
            tmp2 = exp(tmp1)
            fjac[i, 0] = tmp2
            fjac[i, 1] = x[0] * tmp2 / temp
            fjac[i, 2] = -(tmp1 * fjac[i, 1])
        }
    } else if nprob == 12 {
        for i to m {
            temp = i + 1.0
            tmp1 = temp / 10.0
            fjac[i, 0] = -(tmp1 * exp(-(tmp1 * x[0])))
            fjac[i, 1] = tmp1 * exp(-(tmp1 * x[1]))
            fjac[i, 2] = exp(-temp) - exp(-tmp1)
        }
    } else if nprob == 13 {
        for i to m {
            temp = i + 1.0
            fjac[i, 0] = -(temp * exp(temp * x[0]))
            fjac[i, 1] = -(temp * exp(temp * x[1]))
        }
    } else if nprob == 17 {
        for i to len(MGH_Y4) {
            temp = 10.0 * i
            tmp1 = exp(-(x[3] * temp))
            tmp2 = exp(-(x[4] * temp))
            fjac[i, 0] = -1.0
            fjac[i, 1] = -tmp1
            fjac[i, 2] = -tmp2
            fjac[i, 3] = temp * x[1] * tmp1
            fjac[i, 4] = temp * x[2] * tmp2
        }
    } else if nprob == 18 {
        for i to len(MGH_Y5) {
            temp = i / 10.0
            d9 = temp - x[8]
            d10 = temp - x[9]
            d11 = temp - x[10]
            tmp1 = exp(-(x[4] * temp))
            tmp2 = exp(-(x[5] * (d9 * d9)))
            tmp3 = exp(-(x[6] * (d10 * d10)))
            tmp4 = exp(-(x[7] * (d11 * d11)))
            fjac[i, 0] = -tmp1
            fjac[i, 1] = -tmp2
            fjac[i, 2] = -tmp3
            fjac[i, 3] = -tmp4
            fjac[i, 4] = temp * x[0] * tmp1
            fjac[i, 5] = x[1] * (d9 * d9) * tmp2
            fjac[i, 6] = x[2] * (d10 * d10) * tmp3
            fjac[i, 7] = x[3] * (d11 * d11) * tmp4
            fjac[i, 8] = -(2.0 * x[1] * x[5] * d9 * tmp2)
            fjac[i, 9] = -(2.0 * x[2] * x[6] * d10 * tmp3)
            fjac[i, 10] = -(2.0 * x[3] * x[7] * d11 * tmp4)
        }
    }
    return 0
}

#| The function lmder calls (the driver's fcn): residuals for iflag 1, the
#| Jacobian for iflag 2, nothing for iflag 0.
def mgh_fcn(data: Arr, x: Arr, fvec: Arr, fjac: Arr, iflag) {
    if iflag == 1 {
        mgh_residuals(data[0], data[1], x, fvec)
    }
    if iflag == 2 {
        mgh_jacobian(data[0], data[1], x, fjac)
    }
    return iflag
}

#| The standard starting point of problem nprob, times factor (INITPT).
def mgh_start(nprob, n, factor) {
    x = zeros(n)
    if nprob == 4 {
        x[0] = -1.2
        x[1] = 1.0
    } else if nprob == 5 {
        x[0] = -1.0
    } else if nprob == 6 {
        x[0] = 3.0
        x[1] = -1.0
        x[3] = 1.0
    } else if nprob == 7 {
        x[0] = 0.5
        x[1] = -2.0
    } else if nprob == 8 {
        x[0] = 1.0
        x[1] = 1.0
        x[2] = 1.0
    } else if nprob == 9 {
        x[0] = 0.25
        x[1] = 0.39
        x[2] = 0.415
        x[3] = 0.39
    } else if nprob == 10 {
        x[0] = 0.02
        x[1] = 4000.0
        x[2] = 250.0
    } else if nprob == 12 {
        x[1] = 10.0
        x[2] = 20.0
    } else if nprob == 13 {
        x[0] = 0.3
        x[1] = 0.4
    } else if nprob == 17 {
        x[0] = 0.5
        x[1] = 1.5
        x[2] = -1.0
        x[3] = 0.01
        x[4] = 0.02
    } else if nprob == 18 {
        x[0] = 1.3
        x[1] = 0.65
        x[2] = 0.65
        x[3] = 0.7
        x[4] = 0.6
        x[5] = 3.0
        x[6] = 5.0
        x[7] = 7.0
        x[8] = 2.0
        x[9] = 4.5
        x[10] = 5.5
    }
    if factor != 1.0 {
        for j to n {
            x[j] = factor * x[j]
        }
    }
    return x
}

# SPDX-License-Identifier: LicenseRef-Minpack
# SPDX-FileCopyrightText: 1999 University of Chicago (Argonne National Laboratory, MINPACK)
# SPDX-FileCopyrightText: 2026 The Viv authors

# the measurements of the data-fitting problems, from SSQFCN's data statements
MGH_V = [4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625]
MGH_Y1 = [0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37, 0.58, 0.73, 0.96, 1.34, 2.1, 4.39]
MGH_Y2 = [0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246]
MGH_Y3 = [34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0,
          8261.0, 7030.0, 6005.0, 5147.0, 4427.0, 3820.0, 3307.0, 2872.0]
MGH_Y4 = [0.844, 0.908, 0.932, 0.936, 0.925, 0.908, 0.881, 0.85, 0.818, 0.784, 0.751,
          0.718, 0.685, 0.658, 0.628, 0.603, 0.58, 0.558, 0.538, 0.522, 0.506, 0.49,
          0.478, 0.467, 0.457, 0.448, 0.438, 0.431, 0.424, 0.42, 0.414, 0.411, 0.406]
MGH_Y5 = [1.366, 1.191, 1.112, 1.013, 0.991, 0.885, 0.831, 0.847, 0.786, 0.725, 0.746,
          0.679, 0.608, 0.655, 0.616, 0.606, 0.602, 0.626, 0.651, 0.724, 0.649, 0.649,
          0.694, 0.644, 0.624, 0.661, 0.612, 0.558, 0.533, 0.495, 0.5, 0.423, 0.395,
          0.375, 0.372, 0.391, 0.396, 0.405, 0.428, 0.429, 0.523, 0.562, 0.607, 0.653,
          0.672, 0.708, 0.633, 0.668, 0.645, 0.632, 0.591, 0.559, 0.597, 0.625, 0.739,
          0.71, 0.729, 0.72, 0.636, 0.581, 0.428, 0.292, 0.162, 0.098, 0.054]

#| The data array mgh_fcn reads: which problem, and its dimensions.
def mgh_data(nprob, m, n) {
    return [nprob, m, n]
}

#| The residuals of problem nprob at x, into fvec (SSQFCN).
def mgh_residuals(nprob, m, x: Arr, fvec: Arr) {
    if nprob == 4 {
        fvec[0] = 10.0 * (x[1] - x[0] * x[0])
        fvec[1] = 1.0 - x[0]
    } else if nprob == 5 {
        tpi = 8.0 * atan(1.0)
        # dsign(0.25, x(2)); only x(1) = 0 reaches it, and never with x(2) = -0
        tmp1 = 0.25
        if x[1] < 0.0 {
            tmp1 = -0.25
        }
        if x[0] > 0.0 {
            tmp1 = atan(x[1] / x[0]) / tpi
        }
        if x[0] < 0.0 {
            tmp1 = atan(x[1] / x[0]) / tpi + 0.5
        }
        tmp2 = sqrt(x[0] * x[0] + x[1] * x[1])
        fvec[0] = 10.0 * (x[2] - 10.0 * tmp1)
        fvec[1] = 10.0 * (tmp2 - 1.0)
        fvec[2] = x[2]
    } else if nprob == 6 {
        fvec[0] = x[0] + 10.0 * x[1]
        fvec[1] = sqrt(5.0) * (x[2] - x[3])
        d = x[1] - 2.0 * x[2]
        fvec[2] = d * d
        d = x[0] - x[3]
        fvec[3] = sqrt(10.0) * (d * d)
    } else if nprob == 7 {
        fvec[0] = -13.0 + x[0] + ((5.0 - x[1]) * x[1] - 2.0) * x[1]
        fvec[1] = -29.0 + x[0] + ((1.0 + x[1]) * x[1] - 14.0) * x[1]
    } else if nprob == 8 {
        for i to len(MGH_Y1) {
            tmp1 = i + 1.0
            tmp2 = 15.0 - i
            tmp3 = tmp1
            if i + 1.0 > 8.0 {
                tmp3 = tmp2
            }
            fvec[i] = MGH_Y1[i] - (x[0] + tmp1 / (x[1] * tmp2 + x[2] * tmp3))
        }
    } else if nprob == 9 {
        for i to len(MGH_V) {
            tmp1 = MGH_V[i] * (MGH_V[i] + x[1])
            tmp2 = MGH_V[i] * (MGH_V[i] + x[2]) + x[3]
            fvec[i] = MGH_Y2[i] - x[0] * tmp1 / tmp2
        }
    } else if nprob == 10 {
        for i to len(MGH_Y3) {
            temp = 5.0 * (i + 1.0) + 45.0 + x[2]
            tmp1 = x[1] / temp
            tmp2 = exp(tmp1)
            fvec[i] = x[0] * tmp2 - MGH_Y3[i]
        }
    } else if nprob == 12 {
        for i to m {
            temp = i + 1.0
            tmp1 = temp / 10.0
            fvec[i] = exp(-(tmp1 * x[0])) - exp(-(tmp1 * x[1])) + (exp(-temp) - exp(-tmp1)) * x[2]
        }
    } else if nprob == 13 {
        for i to m {
            temp = i + 1.0
            fvec[i] = 2.0 + 2.0 * temp - exp(temp * x[0]) - exp(temp * x[1])
        }
    } else if nprob == 17 {
        for i to len(MGH_Y4) {
            temp = 10.0 * i
            tmp1 = exp(-(x[3] * temp))
            tmp2 = exp(-(x[4] * temp))
            fvec[i] = MGH_Y4[i] - (x[0] + x[1] * tmp1 + x[2] * tmp2)
        }
    } else if nprob == 18 {
        for i to len(MGH_Y5) {
            temp = i / 10.0
            tmp1 = exp(-(x[4] * temp))
            d9 = temp - x[8]
            d10 = temp - x[9]
            d11 = temp - x[10]
            tmp2 = exp(-(x[5] * (d9 * d9)))
            tmp3 = exp(-(x[6] * (d10 * d10)))
            tmp4 = exp(-(x[7] * (d11 * d11)))
            fvec[i] = MGH_Y5[i] - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4)
        }
    }
    return 0
}

#| The Jacobian of problem nprob at x, into fjac (SSQJAC).
def mgh_jacobian(nprob, m, x: Arr, fjac: Arr) {
    if nprob == 4 {
        fjac[0, 0] = -(20.0 * x[0])
        fjac[0, 1] = 10.0
        fjac[1, 0] = -1.0
        fjac[1, 1] = 0.0
    } else if nprob == 5 {
        tpi = 8.0 * atan(1.0)
        temp = x[0] * x[0] + x[1] * x[1]
        tmp1 = tpi * temp
        tmp2 = sqrt(temp)
        fjac[0, 0] = 100.0 * x[1] / tmp1
        fjac[0, 1] = -(100.0 * x[0] / tmp1)
        fjac[0, 2] = 10.0
        fjac[1, 0] = 10.0 * x[0] / tmp2
        fjac[1, 1] = 10.0 * x[1] / tmp2
        fjac[1, 2] = 0.0
        fjac[2, 0] = 0.0
        fjac[2, 1] = 0.0
        fjac[2, 2] = 1.0
    } else if nprob == 6 {
        for j to len(x) {
            for i to m {
                fjac[i, j] = 0.0
            }
        }
        fjac[0, 0] = 1.0
        fjac[0, 1] = 10.0
        fjac[1, 2] = sqrt(5.0)
        fjac[1, 3] = -fjac[1, 2]
        fjac[2, 1] = 2.0 * (x[1] - 2.0 * x[2])
        fjac[2, 2] = -(2.0 * fjac[2, 1])
        fjac[3, 0] = 2.0 * sqrt(10.0) * (x[0] - x[3])
        fjac[3, 3] = -fjac[3, 0]
    } else if nprob == 7 {
        fjac[0, 0] = 1.0
        fjac[0, 1] = x[1] * (10.0 - 3.0 * x[1]) - 2.0
        fjac[1, 0] = 1.0
        fjac[1, 1] = x[1] * (2.0 + 3.0 * x[1]) - 14.0
    } else if nprob == 8 {
        for i to len(MGH_Y1) {
            tmp1 = i + 1.0
            tmp2 = 15.0 - i
            tmp3 = tmp1
            if i + 1.0 > 8.0 {
                tmp3 = tmp2
            }
            s = x[1] * tmp2 + x[2] * tmp3
            tmp4 = s * s
            fjac[i, 0] = -1.0
            fjac[i, 1] = tmp1 * tmp2 / tmp4
            fjac[i, 2] = tmp1 * tmp3 / tmp4
        }
    } else if nprob == 9 {
        for i to len(MGH_V) {
            tmp1 = MGH_V[i] * (MGH_V[i] + x[1])
            tmp2 = MGH_V[i] * (MGH_V[i] + x[2]) + x[3]
            fjac[i, 0] = -(tmp1 / tmp2)
            fjac[i, 1] = -(MGH_V[i] * x[0] / tmp2)
            fjac[i, 2] = fjac[i, 0] * fjac[i, 1]
            fjac[i, 3] = fjac[i, 2] / MGH_V[i]
        }
    } else if nprob == 10 {
        for i to len(MGH_Y3) {
            temp = 5.0 * (i + 1.0) + 45.0 + x[2]
            tmp1 = x[1] / temp
            tmp2 = exp(tmp1)
            fjac[i, 0] = tmp2
            fjac[i, 1] = x[0] * tmp2 / temp
            fjac[i, 2] = -(tmp1 * fjac[i, 1])
        }
    } else if nprob == 12 {
        for i to m {
            temp = i + 1.0
            tmp1 = temp / 10.0
            fjac[i, 0] = -(tmp1 * exp(-(tmp1 * x[0])))
            fjac[i, 1] = tmp1 * exp(-(tmp1 * x[1]))
            fjac[i, 2] = exp(-temp) - exp(-tmp1)
        }
    } else if nprob == 13 {
        for i to m {
            temp = i + 1.0
            fjac[i, 0] = -(temp * exp(temp * x[0]))
            fjac[i, 1] = -(temp * exp(temp * x[1]))
        }
    } else if nprob == 17 {
        for i to len(MGH_Y4) {
            temp = 10.0 * i
            tmp1 = exp(-(x[3] * temp))
            tmp2 = exp(-(x[4] * temp))
            fjac[i, 0] = -1.0
            fjac[i, 1] = -tmp1
            fjac[i, 2] = -tmp2
            fjac[i, 3] = temp * x[1] * tmp1
            fjac[i, 4] = temp * x[2] * tmp2
        }
    } else if nprob == 18 {
        for i to len(MGH_Y5) {
            temp = i / 10.0
            d9 = temp - x[8]
            d10 = temp - x[9]
            d11 = temp - x[10]
            tmp1 = exp(-(x[4] * temp))
            tmp2 = exp(-(x[5] * (d9 * d9)))
            tmp3 = exp(-(x[6] * (d10 * d10)))
            tmp4 = exp(-(x[7] * (d11 * d11)))
            fjac[i, 0] = -tmp1
            fjac[i, 1] = -tmp2
            fjac[i, 2] = -tmp3
            fjac[i, 3] = -tmp4
            fjac[i, 4] = temp * x[0] * tmp1
            fjac[i, 5] = x[1] * (d9 * d9) * tmp2
            fjac[i, 6] = x[2] * (d10 * d10) * tmp3
            fjac[i, 7] = x[3] * (d11 * d11) * tmp4
            fjac[i, 8] = -(2.0 * x[1] * x[5] * d9 * tmp2)
            fjac[i, 9] = -(2.0 * x[2] * x[6] * d10 * tmp3)
            fjac[i, 10] = -(2.0 * x[3] * x[7] * d11 * tmp4)
        }
    }
    return 0
}

#| The function lmder calls (the driver's fcn): residuals for iflag 1, the
#| Jacobian for iflag 2, nothing for iflag 0.
def mgh_fcn(data: Arr, x: Arr, fvec: Arr, fjac: Arr, iflag) {
    if iflag == 1 {
        mgh_residuals(data[0], data[1], x, fvec)
    }
    if iflag == 2 {
        mgh_jacobian(data[0], data[1], x, fjac)
    }
    return iflag
}

#| The standard starting point of problem nprob, times factor (INITPT).
def mgh_start(nprob, n, factor) {
    x = zeros(n)
    if nprob == 4 {
        x[0] = -1.2
        x[1] = 1.0
    } else if nprob == 5 {
        x[0] = -1.0
    } else if nprob == 6 {
        x[0] = 3.0
        x[1] = -1.0
        x[3] = 1.0
    } else if nprob == 7 {
        x[0] = 0.5
        x[1] = -2.0
    } else if nprob == 8 {
        x[0] = 1.0
        x[1] = 1.0
        x[2] = 1.0
    } else if nprob == 9 {
        x[0] = 0.25
        x[1] = 0.39
        x[2] = 0.415
        x[3] = 0.39
    } else if nprob == 10 {
        x[0] = 0.02
        x[1] = 4000.0
        x[2] = 250.0
    } else if nprob == 12 {
        x[1] = 10.0
        x[2] = 20.0
    } else if nprob == 13 {
        x[0] = 0.3
        x[1] = 0.4
    } else if nprob == 17 {
        x[0] = 0.5
        x[1] = 1.5
        x[2] = -1.0
        x[3] = 0.01
        x[4] = 0.02
    } else if nprob == 18 {
        x[0] = 1.3
        x[1] = 0.65
        x[2] = 0.65
        x[3] = 0.7
        x[4] = 0.6
        x[5] = 3.0
        x[6] = 5.0
        x[7] = 7.0
        x[8] = 2.0
        x[9] = 4.5
        x[10] = 5.5
    }
    if factor != 1.0 {
        for j to n {
            x[j] = factor * x[j]
        }
    }
    return x
}

# SPDX-License-Identifier: LicenseRef-Minpack
# SPDX-FileCopyrightText: 1999 University of Chicago (Argonne National Laboratory, MINPACK)
# SPDX-FileCopyrightText: 2026 The Viv authors