Replace enorm with intrinsic norm2

[ivan-pi]

Closes #36.

Was not sure whether I should remove enorm entirely, so I put it in a separate file for now.

I also haven't run any tests to analyze if the output remains the same.

[awvwgk]

I also haven't run any tests to analyze if the output remains the same.

The tests run for the C API by meson should give you some minimal safety. But we really need better tests than manual comparison on the Fortran side.

[ivan-pi]

Okay, I've removed enorm entirely. We can keep this PR waiting until we have more tests. I've been working on an executable for the heart model (#35) and I've also been setting up a repo with the NIST problems (https://github.com/ivan-pi/nist_strd_nls) but it might take me awhile before I can finish them.

[awvwgk]

I'm always happy to see negative line count patches.

[jacobwilliams]

Do you want to pull the latest main into this and have it run the tests? (looks like there are some conflicts)

[ivan-pi]

I'll see to locate some for time this in the following day or two.

[ivan-pi]

I've pulled and merged the new main branch, but I have several failed test cases. In the main branch however, all tests appear to pass. I'll have to investigate a bit.

Edit: I've tried inserting norm2 directly into enorm as follows:

    pure real(wp) function enorm(n, x)

        implicit none

        integer, intent(in) :: n !! a positive integer input variable.
        real(wp), intent(in) :: x(n) !! an input array of length n.

        enorm = norm2(x(1:n))

    end function enorm

None of the tests seems to fail. I've compiled with the original implementation, and the version above, and diffed the output after redirecting it to a file:

ivan@maxwell:~/minpack-1$ fpm test test_hybrd > hybrd_norm2.txt
ivan@maxwell:~/minpack-1$ diff hybrd_enorm.txt hybrd_norm2.txt
ivan@maxwell:~/minpack-1$

It looks like it leads to no differences whatsoever. The same holds for the output of test_lmder.

Edit: I revisited this, and it seems like I didn't save my editor buffer, or some other silly error.

[ivan-pi]

TLDR: Some of the modified test drivers are broken. (Potentially, the old ones were broken too.)

After manually replacing enorm with norm2 however, I'm again facing errors. Edit: I get the same errors in the globally replaced version, or the version where I edited enorm (ivan-pi@3d78b01).

It looks like some of the cases which are failing are affected by symmetry. One of the tests that fails after switching to norm2 is Problem 7 in test_hybrd:

      PROBLEM    7      DIMENSION    5


      INITIAL L2 NORM OF THE RESIDUALS  0.2257066D+00

      FINAL L2 NORM OF THE RESIDUALS    0.3539820D-11

      NUMBER OF FUNCTION EVALUATIONS          17

      EXIT PARAMETER                           1

      FINAL APPROXIMATE SOLUTION

       0.8375126D-01  0.3127293D+00  0.5000000D+00  0.6872707D+00  0.9162487D+00





      PROBLEM    7      DIMENSION    5


      INITIAL L2 NORM OF THE RESIDUALS  0.4117243D+07

      FINAL L2 NORM OF THE RESIDUALS    0.1109510D-08

      NUMBER OF FUNCTION EVALUATIONS         258

      EXIT PARAMETER                           1

      FINAL APPROXIMATE SOLUTION

       0.8375126D-01  0.9162487D+00  0.5000000D+00  0.6872707D+00  0.3127293D+00

Failed case


     Expected x: 

       0.6872707D+00  0.9162487D+00  0.5000000D+00  0.8375126D-01  0.3127293D+00

     Computed x: 

       0.8375126D-01  0.9162487D+00  0.5000000D+00  0.6872707D+00  0.3127293D+00

     absdiff: 

       0.6035194D+00  0.5514722D-10  0.3415872D-09  0.6035194D+00  0.2318776D-09

     reldiff: 

       0.8781393D+00  0.6018804D-10  0.6831744D-09  0.7206094D+01  0.7414642D-09

As you can see in the third test invocation of problem 7, the expected and computed x contain the same values but in different order. This is the Chebyquad function which I can only assume is somehow symmetric. Also the solution of the previous iteration is just another permutation.

Another case which fails is Problem 6 in test_hybrj. The current main branch produces the output:

     PROBLEM    6      DIMENSION    9


      INITIAL L2 NORM OF THE RESIDUALS  0.1015108D+08

      FINAL L2 NORM OF THE RESIDUALS    0.3823081D-01

      NUMBER OF FUNCTION EVALUATIONS         290

      NUMBER OF JACOBIAN EVALUATIONS          12

      EXIT PARAMETER                           4

      FINAL APPROXIMATE SOLUTION

       0.3119912D+00  0.1089342D+01  0.1379293D+01 -0.1175426D+02  0.6265968D+02
      -0.1636005D+03  0.2326239D+03 -0.1697823D+03  0.5076762D+02

The expected values are "hardwired" here:

minpack/test/test_hybrj.f90

Line 220 in 4717189

case (18); x = [0.3119911862752667E+00_wp,0.1089342426782224E+01_wp,0.1379293238883866E+01_wp,&

After switching to norm2 (https://github.com/ivan-pi/minpack-1/tree/switch_enorm/) and rerunning the test, I get the output:

      PROBLEM    6      DIMENSION    9


      INITIAL L2 NORM OF THE RESIDUALS  0.1015108D+08

      FINAL L2 NORM OF THE RESIDUALS    0.3621046D-01

      NUMBER OF FUNCTION EVALUATIONS         311

      NUMBER OF JACOBIAN EVALUATIONS          13

      EXIT PARAMETER                           4

      FINAL APPROXIMATE SOLUTION

       0.3088689D+00  0.1088198D+01  0.1328601D+01 -0.1119031D+02  0.5974643D+02
      -0.1559703D+03  0.2218520D+03 -0.1619874D+03  0.4847785D+02
Failed case


     Expected x: 

       0.3119912D+00  0.1089342D+01  0.1379293D+01 -0.1175426D+02  0.6265968D+02
      -0.1636005D+03  0.2326239D+03 -0.1697823D+03  0.5076762D+02

     Computed x: 

       0.3088689D+00  0.1088198D+01  0.1328601D+01 -0.1119031D+02  0.5974643D+02
      -0.1559703D+03  0.2218520D+03 -0.1619874D+03  0.4847785D+02

     absdiff: 

       0.3122310D-02  0.1144771D-02  0.5069181D-01  0.5639537D+00  0.2913250D+01
       0.7630221D+01  0.1077192D+02  0.7794888D+01  0.2289768D+01

     reldiff: 

       0.1000769D-01  0.1050883D-02  0.3675202D-01  0.4797864D-01  0.4649322D-01
       0.4663935D-01  0.4630614D-01  0.4591107D-01  0.4510293D-01
ERROR STOP 

hybrj exits with the same error code 4, and also the L2-norm of the residuals is actually slightly smaller. This leads me to believe that the modified test_driver is wrong in stopping here.

At the following line

minpack/test/test_hybrj.f90

Line 116 in 4717189

if (info_original(ic)<5) then ! ignore any where the original minpack failed

the logical expression should probably be info_original(ic) < 4 or even info_original(ic) == 1. The error codes for hybrj are:

After fixing this, also hybrj crashes at Problem 7, due to the same symmetry issue.

[ivan-pi]

Here's another one that failed, Problem 14 in test_lmstr.

The current output (I've increased the precision of printing the L2-norm):

      PROBLEM   14      DIMENSIONS    4   20


      INITIAL L2 NORM OF THE RESIDUALS  0.6121125223D+08

      FINAL L2 NORM OF THE RESIDUALS    0.2929543061D+03

      NUMBER OF FUNCTION EVALUATIONS         237

      NUMBER OF JACOBIAN EVALUATIONS         221

      EXIT PARAMETER                           1

      FINAL APPROXIMATE SOLUTION

      -0.1159026D+02  0.1320206D+02 -0.4036921D+00  0.2366622D+00

After switching to norm2 the test fails

      PROBLEM   14      DIMENSIONS    4   20


      INITIAL L2 NORM OF THE RESIDUALS  0.6121125223D+08

      FINAL L2 NORM OF THE RESIDUALS    0.2929543000D+03

      NUMBER OF FUNCTION EVALUATIONS         240

      NUMBER OF JACOBIAN EVALUATIONS         224

      EXIT PARAMETER                           1

      FINAL APPROXIMATE SOLUTION

      -0.1159058D+02  0.1320219D+02 -0.4036522D+00  0.2366854D+00
Failed case


     Expected x: 

      -0.1159026D+02  0.1320206D+02 -0.4036926D+00  0.2366619D+00

     Computed x: 

      -0.1159058D+02  0.1320219D+02 -0.4036522D+00  0.2366854D+00

     absdiff: 

       0.3228592D-03  0.1230603D-03  0.4037176D-04  0.2349015D-04

     reldiff: 

       0.2785607D-04  0.9321291D-05  0.1000062D-03  0.9925614D-04
ERROR STOP 

but according to the L2-norm the solution is actually slightly more accurate, which is a good reason to prefer norm2.

I guess we have to amend compare_solution to also check the absolute value of the L2-norm and not only the absolute/relative difference between the solution vectors.

Branch is outdated