Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Can enorm be replaced with norm2? #36

Open
ivan-pi opened this issue Feb 11, 2022 · 5 comments · May be fixed by #45
Open

Can enorm be replaced with norm2? #36

ivan-pi opened this issue Feb 11, 2022 · 5 comments · May be fixed by #45
Labels
refactoring

Comments

@ivan-pi
Copy link
Member

ivan-pi commented Feb 11, 2022 •
edited

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:

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
@ivan-pi ivan-pi changed the title Is enorm even needed? Can enorm be replaced with norm2 Feb 11, 2022
@ivan-pi ivan-pi changed the title Can enorm be replaced with norm2 Can enorm be replaced with norm2? Feb 11, 2022
@ivan-pi
Copy link
Member Author

ivan-pi commented Feb 12, 2022 •
edited

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.

@arjenmarkus
Copy link
Member

arjenmarkus commented Feb 13, 2022 via email

@ivan-pi
Copy link
Member Author

ivan-pi commented Feb 13, 2022 •
edited

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.

@ivan-pi
Copy link
Member Author

ivan-pi commented Feb 13, 2022 •
edited

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

@ivan-pi
Copy link
Member Author

ivan-pi commented Feb 13, 2022

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.

@ivan-pi ivan-pi linked a pull request Feb 26, 2022 that will close this issue
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
refactoring
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Replace enorm with intrinsic norm2 ivan-pi/minpack-1
2 participants
@arjenmarkus @ivan-pi