Open ivan-pi opened 2 years ago
The CMinpack version of enorm
decided to redefine the rdwarf
and rgiant
constants so that:
rdwarf = sqrt(dpmpar(2)*1.5) * 10
rgiant = sqrt(dpmpar(3)) * 0.1
following the MPFIT library.
They also offer the option of using the BLAS norm with a preprocessor statement:
#ifdef USE_CBLAS
return __cminpack_cblas__(nrm2)(n, x, 1);
#else /* !USE_CBLAS */
// ...
Honestly, all of these feel like magic numbers to me. The Fortran intrinsic norm2
is hopefully the most portable option moving forward, shifting the burden of detecting overflow/underflow to the compiler developers. It also works with any real type.
I agree with using the standard norm2. I experimented a bit to mimic the situation where the intrinsic function is not available and would have to be replaced by an external function. That is a trifle subtle, but with a bit of manipulation that is quite doable.
Op za 12 feb. 2022 om 17:00 schreef Ivan Pribec @.***>:
The CMinpack https://github.com/devernay/cminpack/blob/master/enorm.c version of enorm also decided to redefine the rdwarf and rgiant constants so that:
rdwarf = sqrt(dpmpar(2)1.5) 10 rgiant = sqrt(dpmpar(3)) * 0.1
following the MPFIT http://cow.physics.wisc.edu/~craigm/idl/fitting.html library.
They also offer the option of using the BLAS norm with a preprocessor statement:
ifdef USE_CBLAS
return __cminpack_cblas__(nrm2)(n, x, 1);
else / !USE_CBLAS /
// ...
— Reply to this email directly, view it on GitHub https://github.com/fortran-lang/minpack/issues/36#issuecomment-1037268621, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR5S55E2X52H2OK63Q3U2Z72XANCNFSM5OFWTWOQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.
You are receiving this because you are subscribed to this thread.Message ID: @.***>
In the case norm2
is not available, I would probably just link it as an external function:
#if NO_F2008_NORM2
external :: norm2
#else
intrinsic :: norm2
#endif
One downside is in case you need the optional dim
argument from norm2(array[, dim])
. The built-in version will give a warning about rank-mismatch. The current interface of enorm
implies it is only used for rank-1 arrays.
Edit: Since we are now using a module, it might make more sense to play with importing or over-writing the intrinsic procedure, or even some other preprocessor trickery. In any case it is not worth pursuing unless users ask for it kindly.
A quick grep shows there are 43 occurrences where enorm
appears on the right-hand side:
$ grep -n .*=.*enorm\(.*\) minpack.f90
282: qnorm = enorm(n, Wa2)
301: gnorm = enorm(n, Wa1)
321: temp = enorm(n, Wa2)
333: bnorm = enorm(n, Qtb)
701: fnorm = enorm(n, Fvec)
746: xnorm = enorm(n, Wa3)
818: pnorm = enorm(n, Wa3)
830: fnorm1 = enorm(n, Wa4)
848: temp = enorm(n, Wa3)
880: xnorm = enorm(n, Wa2)
1143: fnorm = enorm(n, Fvec)
1185: xnorm = enorm(n, Wa3)
1258: pnorm = enorm(n, Wa3)
1270: fnorm1 = enorm(n, Wa4)
1288: temp = enorm(n, Wa3)
1322: xnorm = enorm(n, Wa2)
1614: fnorm = enorm(m, Fvec)
1660: xnorm = enorm(n, Wa3)
1728: pnorm = enorm(n, Wa3)
1740: fnorm1 = enorm(m, Wa4)
1758: temp1 = enorm(n, Wa3)/fnorm
1793: xnorm = enorm(n, Wa2)
2083: fnorm = enorm(m, Fvec)
2129: xnorm = enorm(n, Wa3)
2198: pnorm = enorm(n, Wa3)
2210: fnorm1 = enorm(m, Wa4)
2228: temp1 = enorm(n, Wa3)/fnorm
2266: xnorm = enorm(n, Wa2)
2503: dxnorm = enorm(n, Wa2)
2530: temp = enorm(n, Wa1)
2544: gnorm = enorm(n, Wa1)
2570: dxnorm = enorm(n, Wa2)
2599: temp = enorm(n, Wa1)
2768: fnorm = enorm(m, Fvec)
2813: Wa2(j) = enorm(j, Fjac(1, j))
2849: xnorm = enorm(n, Wa3)
2896: pnorm = enorm(n, Wa3)
2908: fnorm1 = enorm(m, Wa4)
2926: temp1 = enorm(n, Wa3)/fnorm
2963: xnorm = enorm(n, Wa2)
3235: Acnorm(j) = enorm(m, a(1, j))
3270: ajnorm = enorm(m - j + 1, a(j, j))
3296: Rdiag(k) = enorm(m - j, a(jp1, k))
In most cases the number of elements matches the array size so no problem there. Probably only four occurrences reference rank-2 arrays:
$ grep -n .*=.*enorm\(.*\(.*\)\) minpack.f90
2813: Wa2(j) = enorm(j, Fjac(1, j))
3235: Acnorm(j) = enorm(m, a(1, j))
3270: ajnorm = enorm(m - j + 1, a(j, j))
3296: Rdiag(k) = enorm(m - j, a(jp1, k))
Here's the wider context of line 2813:
sing = .false.
do j = 1, n
if (Fjac(j, j) == zero) sing = .true.
Ipvt(j) = j
Wa2(j) = enorm(j, Fjac(1, j))
end do
The norm is taken over columns in the upper triangular part of the matrix Fjac
, so you can't use the dim
argument anyways. Here's a refactored version with the three assignments separated:
! check if Jacobian is singular
sing = any(diag(fjac) == zero)
! initialize pivot array
Ipvt = [(j, j = 1, n)]
! calculate norm of upper triangular columns
do j = 1, n
Wa2(j) = norm2(fjac(1:j,j))
end do
While it communicates the intent better, I can't decide whether it's a modification worth pursuing.
Here's the context of line 3235:
! compute the initial column norms and initialize several arrays.
do j = 1, n
Acnorm(j) = enorm(m, a(1, j))
Rdiag(j) = Acnorm(j)
Wa(j) = Rdiag(j)
if (Pivot) Ipvt(j) = j
end do
In this case we could use the dim
argument:
! compute the initial column norms and initialize several arrays
acnorm = norm2(a, dim=2)
rdiag = acnorm
wa = rdiag
if (pivot) Ipvt = [(j, j = 1, n)]
What do you guys think:
norm2
) in favor of clarity/simplicity an acceptable goal of modernization? I feel like this is part of a bigger discussion whether the "modern" MINPACK uses raw loops or Fortran 90 array operations.Do we need a set of refactoring/style guidelines or will these establish themselves naturally through submitting pull requests? Personally, I already miss the consistent style of the original F77 code.
cc @certik @zbeekman @milancurcic @awvwgk
Before I forget, the two instances:
3270: ajnorm = enorm(m - j + 1, a(j, j))
3296: Rdiag(k) = enorm(m - j, a(jp1, k))
can be rewritten as
ajnorm = norm2(a(j:m,j))
Rdiag(k) = norm2(a(jp1:m,k))
In both cases, the norm is taken over a section of an array column, so no need for the dim
option either. Just an observation, the array section syntax makes it obvious we are taking these along the first dimension. The older method of passing the address to the first element obfuscates this.
The function
enorm
provides the Euclidean norm of a vector. From the MINPACK user guide:The function
enorm
is obsolescent as far as I'm concerned. The Fortran 2008 standard and later provide anorm2
intrinsic function. The standard recommends that processors compute the result without undue overflow or underflow.Other notable norm implementations include:
enorm
: https://www.netlib.org/minpack/enorm.fdnrm2
: https://www.netlib.org/slatec/lin/dnrm2.fdnrm2
: http://www.netlib.org/lapack/explore-3.1.1-html/dnrm2.f.htmlnorm2
: gfortran docs, ifort docsHere's a quick demonstration program:
An example run:
A potential issue of removing
enorm
in favor ofnorm2
are the workarounds needed to support older compilers. One approach could be to use a preprocessor block:Personally, I'd just replace
enorm
withnorm2
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