minpack
minpack copied to clipboard
fortran-lang
Can enorm be replaced with norm2?
The function enorm provides the Euclidean norm of a vector. From the MINPACK user guide:
The algorithm used in
ENORMis 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 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:
- Is converting to elemental operations (including array transformational functions like
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.