The current tests are not actually testing anything other than running the examples. They provide no way to check the results other than by manual inspection, which makes them not reliable as regression tests.

Required steps:

• split the objective functions in separate callbacks rather than one big select case via a common block variable
• have a resource module to hold the objective functions and initializers to avoid code duplication
• actually test the outcome of the examples is correct up to a given tolerance

I think these are test problems for MINPACK: https://ftp.mcs.anl.gov/pub/MINPACK-2/tprobs/

Need to modernize them and add them to the unit tests.

Tracking issue for first release of the modernized minpack v2.0.0.

Currently the Python bindings can be built with setuptools or meson, the latter is preferred since we have more control with meson over the build environment. For building the Python wheels we therefore use meson-python a PEP517 compliant build system wrapping meson.

Some upstream issues with preferably should be resolved first:

• possibility to include subprojects in sdist: https://github.com/FFY00/meson-python/issues/59
• support for editable / development installations: https://github.com/FFY00/meson-python/issues/47

From our side we need to do some testing on MacOS, Windows and *BSD.

• successful build on MacOS/x86_64
• successful build on Windows
• successful build on *BSD

Here's a section of code from rwupdt:

        do j = 1, n
            rowj = w(j)
            jm1 = j - 1

            ! apply the previous transformations to
            ! r(i,j), i=1,2,...,j-1, and to w(j).

            if (jm1 >= 1) then
                do i = 1, jm1
                    temp = Cos(i)*r(i, j) + Sin(i)*rowj
                    rowj = -Sin(i)*r(i, j) + Cos(i)*rowj
                    r(i, j) = temp
                end do
            end if

As you may notice for j = 1, jm1 = 0, meaning the if loop surrounding the inner do loop is totally superfluous.

As @arjenmarkus kindly explained in his reply at Discourse, this is a "legacy" issue related to Fortran 66 DO semantics. Edit: See Section 7.1.2.8 in the ANSI X 3.9 Fortran 66 standard. The requirement to surround the DO loop with an IF block came from the line:

At time of execution of the DO statement, m1, m2, and m3 must be greater than zero.

The MINPACK source code in the present state remains full of this pattern. Maybe CamFort, plusFORT, or fpt have a tool which could fix this?

It would be a good idea to start a document or wiki- page collecting such legacy "gotchas" if we plan to do more "modernization".

Using GitHub actions, on Linux, macOS and Windows.

The C bindings pass a nested function as callback to allow the C API to carry along user data via an opaque pointer (void*). This mechanism is also used in the Python API to propagate exceptions storing the small object in a handle and pass it as void* to the Minpack C API.

This issue has been fixed for GCC 12, see iains/gcc-darwin-arm64#26.

Potential workaround available in #30

While minpack is an important package among the Fortran legacy codes it's age is apparent.

The Ceres Solver is a non-linear least squares package written in C++ and developed at Google. They also performed a comparison of some solvers in their release note: https://groups.google.com/g/ceres-solver/c/UcicgMPgbXw

The comparison is performed using a set of 54 test problems compiled by NIST: https://www.itl.nist.gov/div898/strd/nls/nls_info.shtml

Both Ceres and HBN (not sure if this is equal to immoptibox?) have a higher total score than minpack.

How should we aim to promote minpack when better solvers are available (at least for some problems)?

Does it make sense to "restart" minpack development and bring the new algorithms into minpack? (Ceres is BSD-licensed.)

From Recommendations To Write (Slightly More) Readable And (Thus) Robust Code 7:

Exit early to avoid waterfall code. Especially for longer code blocks, this minimizes cognitive overhead.

(originally a suggestion from @interkosmos posted at Discourse)

The current modularized minpack source code, includes sections for checking the input parameters like this one taken from lmder1:

        Info = 0

        ! check the input parameters for errors.

        if (n > 0 .and. m >= n .and. Ldfjac >= m .and. Tol >= zero .and. &
            Lwa >= 5*n + m) then
            ! call lmder.
            maxfev = 100*(n + 1)
            ftol = Tol
            xtol = Tol
            gtol = zero
            mode = 1
            nprint = 0
            call lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, ftol, xtol, gtol, maxfev,   &
                     & Wa(1), mode, factor, nprint, Info, nfev, njev, Ipvt, Wa(n + 1)&
                     & , Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1))
            if (Info == 8) Info = 4
        end if

Here in the "easy" driver there is only level of indentation, so it's not as problematic, but for the actual routines which do the work, getting rid of one level of indentation would make them easier on the eye.

The solution is to return early, per the advice above:

info = 0  ! flag for improper input parameters

! check the input parameters for errors

if (n < 1 .or. m < n) return
if (ldfjac < m) return
if (tol < zero) return
if (lwa < 5*n + m) return


! default settings

maxfev = 100*(n+1)
ftol = tol
xtol = tol
gtol = zero
mode = 1
nprint = 0

call lmder( ... )
if (info == 8) info = 4

Since Fortran does not short-circuit expressions I guess it also makes sense to separate the condition. (For correct input parameters, all the conditionals should be false, but in the case not, we might save a few evaluations and exit early.) This might also just by a personal preference of mine, that is I dislike long logical expressions and/or line continuation:

if (n < 1 .or. m < n .or. ldfjac < m .or. tol < zero .or. lwa < 5*n + m) &
    return

Currently, we are using double precision. Consider having the option of:

1. exporting multiple precisions
2. specify the precision via preprocessor flag

For 2, I usually would do something like this:

#ifdef REAL32
    integer,parameter,public :: minpack_rk = real32   !! real kind used by this module [4 bytes]
#elif REAL64
    integer,parameter,public :: minpack_rk = real64   !! real kind used by this module [8 bytes]
#elif REAL128
    integer,parameter,public :: minpack_rk = real128  !! real kind used by this module [16 bytes]
#else
    integer,parameter,public :: minpack_rk = real64   !! real kind used by this module [8 bytes]
#endif

    integer,parameter,private :: wp = minpack_rk  !! local copy of `minpack_rk` with a shorter name

So, if you do nothing, you get double precision. The module also doesn't export wp, since that is too short a variable name to be exporting (if you want it, you will get minpack_rk).

For 1: Look at what I did for quadpack. Not sure we need that for minpack though. Are there use cases for having multiple versions in the same project with different precisions?

Just found this report which tested MINPACK and NL2SOL on a problem originating from a heart model:

A New Nonlinear Equations Test Problem, J.E. Dennis, Jr., David M. Gay, and
Phuong Ahn Vu, Technical Report 83-1.P, December 1985.

A Google search gives a PDF copy at the following address: https://scholarship.rice.edu/bitstream/handle/1911/101557/TR83-16.pdf

Aside from the report, Fortran source code for the problem can also be found in opt/dqed.f on Netlib.

It would be nice if data could be passed to the objective function. This could be done with a void C pointer, type(c_ptr), or any other method. Right now data must be passed globally, which makes programs not thread-safe.

Reported by @xecej4:

The INTENT of arguments fvec and fjac of subroutine FCN should be changed from OUT to IN OUT. As of now, when the subroutine is entered with iflag = 1, fjac becomes undefined; when entered with iflag = another value, fvec becomes undefined.

This subtle difference matters when flags such as NAG's -C=undefined or FTN95's /undef are used.

Transferred from certik/minpack#7

See: https://www.netlib.org/minpack/ex/file03

(Note that column 73 is missing from that file). Good grief.

The function enorm provides the Euclidean norm of a vector. From the MINPACK user guide:

The algorithm used in ENORM is a simplified version of Blue's [1978] algorithm. An advantage of the MINPACK-1 version is that it does not require machine constants; a disadvantage is that nondestructive underflows are allowed.

The function enorm is obsolescent as far as I'm concerned. The Fortran 2008 standard and later provide a norm2 intrinsic function. The standard recommends that processors compute the result without undue overflow or underflow.

Other notable norm implementations include:

• Original MINPACK enorm: https://www.netlib.org/minpack/enorm.f
• SLATEC dnrm2: https://www.netlib.org/slatec/lin/dnrm2.f
• Reference BLAS dnrm2: http://www.netlib.org/lapack/explore-3.1.1-html/dnrm2.f.html
• Fortran intrinsic function norm2: gfortran docs, ifort docs

Here's a quick demonstration program:

program main
  implicit none  
  integer, parameter :: dp = kind(1.0d0)
  real(dp) :: x(20)
  real(dp) :: dnrm2, enorm
  call random_number(x)
  print *, "enorm ", enorm(size(x), x)
  print *, "dnrm2 ", dnrm2(size(x), x, 1)
  print *, "norm2 ", norm2(x)
end program

An example run:

$ gfortran dnrm2.f enorm.f norm_example.f90 
$ ./a.out
 enorm    2.1356777438951129     
 dnrm2    2.1356777438951133     
 norm2    2.1356777438951129     

A potential issue of removing enorm in favor of norm2 are the workarounds needed to support older compilers. One approach could be to use a preprocessor block:

pure real(wp) function enorm(n, x)
  integer, intent(in) :: n
  real(wp), intent(in) :: x(n)
#if MINPACK1_ENORM
  ! ... current MINPACK-1 algorithm ...
#else
  enorm = norm2(x)
#endif
end function

Personally, I'd just replace enorm with norm2 everywhere it occurs. If deemed necessary, the original implementation can be preserved as an external subroutine in a "compatibility" folder for older compilers.

I'd also be interested in learning what type of tests can we come up with to prove the "worthiness" of the intrinsic norm2 over the MINPACK-1 algorithm.

Literature

• A Portable Fortran Program to Find the Euclidean Norm of a Vector (1978) | https://doi.org/10.1145/355769.355771
• Efficient Calculations of Faithfully Rounded l2-Norms of n-Vectors (2015) | https://doi.org/10.1145/2699469
• Remark on Algorithm 539: A Modern Fortran Reference Implementation for Carefully Computing the Euclidean Norm (2017) | https://dl.acm.org/doi/10.1145/3134441
• Algorithm 978: Safe Scaling in the Level 1 BLAS (2017) | https://doi.org/10.1145/3061665
• A Rapid Euclidean Norm Calculation Algorithm that Reduces Overflow and Underflow (2021) | https://doi.org/10.1007/978-3-030-86653-2_7
• Accurate calculation of Euclidean Norms using Double-word arithmetic (2022) | https://hal.archives-ouvertes.fr/hal-03482567/

The documentation provided with the Python bindings should be collected and deployed together with the Fortran API docs.

We should rename the default branch from master to main. See fortran-lang/fpm#421 for details.

Only two explicit references to master as default branch are in the README, so this should be straight-forward.

❯ rg master .
./README.md
27: * The API documentation for the current ```master``` branch can be found [here](https://fortran-lang.github.io/minpack/).  This is generated by processing the source files with [FORD](https://github.com/Fortran-FOSS-Programmers/ford).
31:The Minpack source code and related files and documentation are distributed under a permissive free software [license](https://github.com/fortran-lang/minpack/blob/master/LICENSE.txt) (BSD-style).

The C bindings are currently not documented outside of the comments in the header. Maybe we can make use of doxygen to create API docs for those.

Not relevant yet, but once we provide Python bindings (see #43), we probably have to face the topic of creating wheels as well and distributing over PyPI.

From my initial attempts for creating wheels with compiled extension modules, I can say it is quite tedious and somewhat of pain, but maybe somebody made better experience with it and knows how to get it done correctly. The more robust distribution channel for Python extension models seems to be conda, especially via conda-forge.

Nvidia's Fortran compiler can compile minpack, however it segfaults in the examples

❯ nvfortran --version

nvfortran 22.1-0 64-bit target on x86-64 Linux -tp haswell 
NVIDIA Compilers and Tools
Copyright (c) 2022, NVIDIA CORPORATION & AFFILIATES.  All rights reserved.

❯ fpm run --example example_primes --compiler nvfortran --runner gdb
...
(gdb) run
Starting program: /home/awvwgk/projects/src/git/minpack/build/nvfortran_15E77B831143CB1F/example/example_primes 
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/usr/lib/libthread_db.so.1".

Program received signal SIGSEGV, Segmentation fault.
0x00000000004028f4 in fcn () at examples/example_primes.f90:60
60	fvec = data_y - expr(data_x, x)
(gdb) bt
#0  0x00000000004028f4 in fcn () at examples/example_primes.f90:60
#1  0x000000000040ba45 in minpack_module::lmdif () at ./src/minpack.f90:2080
#2  0x000000000040d951 in minpack_module::lmdif1 () at ./src/minpack.f90:2383
#3  0x000000000040279c in find_fit_module::find_fit () at examples/example_primes.f90:48
#4  0x0000000000402ceb in example_primes () at examples/example_primes.f90:81

See https://fortran-lang.discourse.group/t/moving-minpack-under-fortran-lang/2631

See also: https://github.com/certik/minpack/

I'm open to moving this under fortran-lang.

The driver programs do not assign zero initial values to the variable arrays x, Fvec and Fjac. Only the nonzero initial values are programmed. As a result, if the Fortran compiler in use does not provide a switch to do zero-initialization, or the user does not use the switch, the results will most probably be completely wrong.

Correcting this deficiency should also take care of issue #5.

I have archived the https://github.com/certik/minpack repository and moved the open issue from there to here. What remains is to port all the improvements that were implemented there over to this repository.

The tree from the other repo is available at https://github.com/fortran-lang/minpack/tree/c0dc9cd as well.

Valgrind says that some of the tests have errors. For example test_hybrd

==7124== 
==7124== HEAP SUMMARY:
==7124==     in use at exit: 0 bytes in 0 blocks
==7124==   total heap usage: 27 allocs, 27 frees, 26,772 bytes allocated
==7124== 
==7124== All heap blocks were freed -- no leaks are possible
==7124== 
==7124== Use --track-origins=yes to see where uninitialised values come from
==7124== For lists of detected and suppressed errors, rerun with: -s
==7124== ERROR SUMMARY: 88315 errors from 436 contexts (suppressed: 0 from 0)

example error:

==7124== Conditional jump or move depends on uninitialised value(s)
==7124==    at 0x116888: __minpack_module_MOD_enorm (minpack.f90:392)
==7124==    by 0x10AC6F: MAIN__ (test_hybrd.f90:66)
==7124==    by 0x10B003: main (test_hybrd.f90:3)

The iflag arguments to fdjac1 and fdjac2 should be intent(inout).

The meta data in the package manifest seems somewhat inaccurate

For comparison the version by Ondřej @certik specifies

I think using a semantic version 1.0.0 would be more appropriate, however if the API was changed it should use 2.0.0 instead. Also, the author and copyright field seem incorrect in this version.

SciPy contains minpack: https://github.com/scipy/scipy/tree/main/scipy/optimize/minpack; we should develop this fortran-lang/minpack in a way so that SciPy could use it. We should add some improvements (perhaps some new minimization algorithms), and then include them in SciPy. If we can pull this off, that would be a huge win, for both Fortran and SciPy. SciPy gets a maintainable version, we can all contribute to it in Fortran, and modern Fortran gets a "user" in a high profile Python library.

Fortran Forum discussion thread: https://fortran-lang.discourse.group/t/lets-get-fortran-lang-minpack-into-scipy/2722

Follow on to #40. Consolidate all the duplicated code from the various units tests.

In Subroutine FCN of file test_lmstr.f90, the 2-D array TEMP is filled when FCN is called with IFLAG=2. When FCN is called with IFLAG=3,4,..., the argument array FJROW is filled with a row of TEMP. This will work only if TEMP is given the SAVE attribute. Please see the comments in subroutine FCN.

Use iso_c_binding to make a C-style interface to the routines now in a module. See also #14.

A typical usage case for non-linear least squares is curve fitting.

In the archived version of @certik, there is already a demo example of a find_fit subprogram: https://github.com/certik/minpack/blob/b46766bd42ec6f08114caf5f6d811d815f78a809/examples/example_primes.f90#L23

Both SciPy and MATLAB provide convenience functions for this purpose:

• SciPy: scipy.optimize.curve_fit
• MATLAB: lsqcurvefit

These are effectively just wrappers of the general non-linear least squares routines, which take the pair (xdata, ydata) alongside the function to be fitted.

The code has about 30 occurrences of goto. We should replace it with modern Fortran's exit, cycle, return and other flow constructs.

This is very important to make people more excited to use the new code base (such as for SciPy #14), see e.g. this comment:

I’m not sure if that’s what actually modern Fortran looks like - I had expected not, because there’s still a lot of goto's in the code (which is bad).

There is a MINPACK wrapper for R called minpack.lm.

It might be interesting to check how many times it has been downloaded and how many R packages depend on it.

Having a C interface also makes the task of writing an R wrapper easier. Perhaps we can get in touch with the authors of that package at some point and ask them to use the new maintained version instead.

We could provide Python bindings directly in this project. Dependent Python projects would be able to just import minpack rather than wrapping the C API themselves.

• Numba compatibility of Python bindings possible with ctypes and cffi but not with Cython:
  https://fortran-lang.discourse.group/t/fortranlang-a-python-package/2870

Related to #14

Seems like the Intel compiler can almost compile minpack, except for one of the examples

❯ ifx -V
Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 2022.0.0 Build 20211123
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.
❯ fpm test --compiler ifx
minpack.f90                            done.
minpack_testing.f90                    done.
example_hybrd1.f90                     done.
example_lmdif1.f90                     done.
example_hybrd.f90                      done.
example_lmder1.f90                     done.
example_primes.f90                     failed.
test_chkder.f90                        done.
test_hybrd.f90                         done.
test_hybrj.f90                         done.
[  38%]Compiling...
examples/example_primes.f90: remark #5415: Feature not yet implemented: Some 'check' options temporarily disabled.
 #0 0x0000000001084382 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x1084382)
 #1 0x0000000001084356 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x1084356)
 #2 0x0000000000fdf5b6 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0xfdf5b6)
 #3 0x000000000100082a (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x100082a)
 #4 0x0000000000ffec51 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0xffec51)
 #5 0x0000000000fe1e3d (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0xfe1e3d)
 #6 0x000000000105bc37 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x105bc37)
 #7 0x000000000105e391 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x105e391)
 #8 0x000000000105f5d8 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x105f5d8)
 #9 0x00000000010d1e75 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x10d1e75)
#10 0x00000000010cfae2 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x10cfae2)
#11 0x00000000010d01b2 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x10d01b2)
#12 0x00000000010cd9a2 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x10cd9a2)
#13 0x00000000010cfae2 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x10cfae2)
#14 0x00000000010cd05a (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x10cd05a)
#15 0x00000000010cfae2 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x10cfae2)
#16 0x0000000000f9d436 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0xf9d436)
#17 0x0000000000f9cb71 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0xf9cb71)
#18 0x0000000001116541 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0x1116541)
#19 0x00007fc61d475b25 __libc_start_main (/usr/lib/libc.so.6+0x27b25)
#20 0x0000000000de27a9 (/opt/intel/oneapi/compiler/2022.0.2/linux/bin-llvm/xfortcom+0xde27a9)

examples/example_primes.f90(31): error #5533: Feature found on this line is not yet supported in ifx
    function expr(x, pars) result(y)
-------------^
compilation aborted for examples/example_primes.f90 (code 3)
<ERROR> Compilation failed for object " examples_example_primes.f90.o "
<ERROR>stopping due to failed compilation
STOP 1

But passes with a small workaround

diff --git a/examples/example_primes.f90 b/examples/example_primes.f90
index c014bd3..57b2531 100644
--- a/examples/example_primes.f90
+++ b/examples/example_primes.f90
@@ -18,6 +18,14 @@ implicit none
 private
 public find_fit
 
+interface
+    function expr_interface(x, pars) result(y)
+    use types, only: dp
+    implicit none
+    real(dp), intent(in) :: x(:), pars(:)
+    real(dp) :: y(size(x))
+    end function
+end interface
 contains
 
 subroutine find_fit(data_x, data_y, expr, pars)
@@ -27,15 +35,8 @@ subroutine find_fit(data_x, data_y, expr, pars)
 ! array 'x'. The arrays 'data_x' and 'data_y' must have the same
 ! length.
 real(dp), intent(in) :: data_x(:), data_y(:)
-interface
-    function expr(x, pars) result(y)
-    use types, only: dp
-    implicit none
-    real(dp), intent(in) :: x(:), pars(:)
-    real(dp) :: y(size(x))
-    end function
-end interface
 real(dp), intent(inout) :: pars(:)
+procedure(expr_interface) :: expr
 
 real(dp) :: tol, fvec(size(data_x))
 integer :: iwa(size(pars)), info, m, n

Just add the normal fortran-lang MIT style license on top of the other two existing ones? with: Copyright (c) 2022 fortran-lang minpack contributors?

The upload of the release artifacts should be automatic, except for the minpack-2.0.0-source.tar.xz, which can be created by meson dist all other artifacts are already created and available via CI.

In order for the 18 problems in test_lmder.f90 to work correctly, please note that the values of n (number of fit coefficients) and m (number of "data" points to fit) should not be fixed at 40 and 65, respectively, as doing so will cause some test problems to fail to run to completion. For example, with nprob = 13, m should be 10, but calculating as if m were 65 causes overflow when the argument to exp() becomes larger than permissible, since that argument is i*x(1), and i is allowed to exceed m = 55, with x(1) = 12.998.

Please change the outermost DO loop in the program unit as follows:

    do 
       read *,nprob,n,m,ntries
       if (NPROb == 0) then

where the data are read from the file provided at Netlib: https://netlib.org/minpack/ex/file22 .

One of the suggestions in #14 (comment) was to port the MINPACK tests from SciPy to our test suite.

I've tried to locate some of the relevant files:

• https://github.com/scipy/scipy/blob/820fe86bf3ad8ba1c5647632b131659ed7bd447e/scipy/optimize/tests/test_nonlin.py (nice to see @certik among the authors)
• https://github.com/scipy/scipy/blob/820fe86bf3ad8ba1c5647632b131659ed7bd447e/scipy/optimize/tests/test_minpack.py
• https://github.com/scipy/scipy/blob/820fe86bf3ad8ba1c5647632b131659ed7bd447e/scipy/optimize/tests/test_least_squares.py
• https://github.com/scipy/scipy/blob/820fe86bf3ad8ba1c5647632b131659ed7bd447e/scipy/optimize/tests/test__root.py

Note that some of the test may only be testing properties of the SciPy interface, and not have much to do with Fortran at all.

Similar to what was done for the tests.

Should we move the examples into the test folder? Or is there a reason to keep them separate?

The Python bindings should be Numba compatible.

• CFFI with Numba

Currently the Levenberg-Marquardt routine in MINPACK calls a "home-brewed" QR factorization derived from LINPACK:

https://github.com/jacobwilliams/minpack/blob/b4e671d1692b6f9eb048ba920e97f8a4eda195e3/src/minpack.f90#L2578

This is understandable since MINPACK (released in 1980) precedes LAPACK (released in 1992) by more than a decade. To bring minpack into the 21st century it would be nice to support the LAPACK factorization routines, and only use the original ones as a fallback.

I've looked into this before, but unfortunately I lack the domain knowledge needed to do the modernization. I also had the feeling, minpack does a row by row factorization to somehow reduce storage costs. This might also have to do with rank-defective Jacobians. I see LAPACK comes in two flavors exactly for this purpose:

• QR Factorization
• QR Factorization with Column Pivoting (useful when a matrix is not full rank)

An overview of the orthogonal factorization routines in LAPACK is also provided in the Intel oneAPI MKL documentation

• Orthogonal Factorizations: LAPACK Computational Routines

Given the reliability and trustworthiness of minpack routines I think a plan is needed how to perform the refactoring without introducing bugs and guaranteeing the accuracy of the results is maintained (or even improved).

Currently, MINPACK exports 3 API's:

• the original Fortran API
• the C API
• the Python API

There are some minimal tests for each of the API's and a big messy collection of Fortran test functions.

Clearly, some minimal tests for each of the different interfaces are needed in order to test the correctness of the binding, with respect to argument passing, types, side-effects, etc. For the actual algorithms in Fortran however, it doesn't matter in which languages / interface the tests are written.

Does it make sense to write the algorithm tests (#10, #35) in Python?

This would test all three interfaces simultaneously. It also has the advantage of faster development due to dynamic typing, simpler output formatting, and other high-level properties of Python. Moreover the larger Python community could support us with their CI/CD knowledge, e.g. from SciPy and other big Python libraries.

cc @awvwgk @ilayn @certik

In the newly added Python bindings, the expected callback is of form fcn(x, fvec) -> None where fvec is an output vector, of the same size as x, which is modified in-place. The full set of Python callback interfaces is defined here: python/minpack/typing.py

The SciPy least-squares on the other hand, expects a callback with the following description:

Function which computes the vector of residuals, with the signature fun(x, *args, **kwargs), i.e., the minimization proceeds with respect to its first argument. The argument x passed to this function is an ndarray of shape (n,) (never a scalar, even for n=1). It must allocate and return a 1-D array_like of shape (m,) or a scalar. If the argument x is complex or the function fun returns complex residuals, it must be wrapped in a real function of real arguments, as shown at the end of the Examples section.

Is the reason for the different interfaces since the present Fortran MINPACK interface cannot handle extra arguments?

Is the idea that the SciPy layer on top of the new Python binding, would nest a small adaptor function?

    def _callback_hy(x,fev) -> None:
        fev[:] = fun(x, *args, **kwargs)

A search for "minpack" on StackExchange gives 138 results mostly on StackOverflow and Computational Science StackExchange pages.

If anyone has some spare time, sifting through the results could offer some insights on the issues people encounter while using MINPACK.

The version is here: https://github.com/scipy/scipy/tree/a0535bee9d4da58e6d41e963e907c1b9f7684585/scipy/optimize/minpack and we need to check to see what changes they have made to the original code and port those to this version.

Towards #14.

Testing with NAG's Fortran compiler yields runtime errors due to arithmetic exceptions

❯ nagfor
NAG Fortran Compiler Release 7.1(Hanzomon) Build 7104
❯ fpm test --compiler nagfor test_lmder
...
Runtime Error: *** Arithmetic exception: Floating overflow - aborting
test/test_lmder.f90, line 904: Error occurred in TEST:SSQFCN
test/test_lmder.f90, line 118: Called by TEST:FCN
././src/minpack.f90, line 1737: Called by MINPACK_MODULE:LMDER
././src/minpack.f90, line 1921: Called by MINPACK_MODULE:LMDER1
test/test_lmder.f90, line 63: Called by TEST
 <ERROR> Execution failed for "test_lmder"
STOP 1
❯ fpm test --compiler nagfor test_lmdif
...
Runtime Error: *** Arithmetic exception: Floating overflow - aborting
test/test_lmdif.f90, line 370: Error occurred in SSQFCN
test/test_lmdif.f90, line 122: Called by FCN
././src/minpack.f90, line 2207: Called by MINPACK_MODULE:LMDIF
././src/minpack.f90, line 2383: Called by MINPACK_MODULE:LMDIF1
test/test_lmdif.f90, line 74: Called by TEST
 <ERROR> Execution failed for "test_lmdif"
STOP 1

Add an object oriented interface that is a higher-level wrapper to the functions. (for example, see how this is done in bspline-fortran. See also: #4

Many computational problems will have banded matrices. In those cases a performance increase can be expected from using specialized routines for banded storage.

This question appeared previously on scicomp.stackexchange: https://scicomp.stackexchange.com/questions/36703/nonlinear-root-solving-libraries-which-accept-a-jacobian-in-band-storage

MINPACK uses a QR factorization for the Jacobians. Unfortunately, LAPACK does not provide QR routines for banded Jacobians, however a library does exist for this purpose:

• LAPACK-Style Codes for the QR Factorization of Banded Matrices
• BAND_QR (contains sources for the routines of Remon et al.)

A similar specialization could be pursued for sparse matrices.