minpack
minpack copied to clipboard
fortran-lang
Replace enorm with intrinsic norm2
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.
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.
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.
Do you want to pull the latest main into this and have it run the tests? (looks like there are some conflicts)
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.
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 (https://github.com/ivan-pi/minpack-1/commit/3d78b0195ee7b725efb599a782d87c62edfb5c1d).
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: https://github.com/fortran-lang/minpack/blob/47171898646519c02ffae7d21a776d21ad41e222/test/test_hybrj.f90#L220
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
https://github.com/fortran-lang/minpack/blob/47171898646519c02ffae7d21a776d21ad41e222/test/test_hybrj.f90#L116
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.
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.