***SEE ALSO: HYPOTENUSE***

What is the best implementation for NORM in common/linalg.f90?
Here are the places where it is used. It affects the performance of NEWUOA.

common/linalg.f90:    u = u / norm(u)
common/linalg.f90:    u = v / norm(v)
cobyla/trustregion.f90:    d = (delta / norm(d)) * d
cobyla/cobylb.f90:    dnorm = min(delta, norm(d))
uobyqa/trustregion.f90:    d = (delta / norm(d)) * d
newuoa/newuob.f90:    dnorm = min(delta, norm(d))
newuoa/newuob.f90:        dnorm = min(delbar, norm(d))  ! In theory, DNORM = DELBAR in this case.
newuoa/trustregion.f90:    d = (norm(s) / norm(d)) * d
newuoa/geometry.f90:    s = (norm(d) / norm(s)) * s
newuoa/geometry.f90:    s = (s / norm(s)) * norm(d)

In particular, common/linalg.f90/PROJECT1 uses NORM, which is in turn used in NEWUOA. It is probably
the place that affects NEWUOA the most. PROJECT1 is not used by any other solver yet.

Below is the implementation of DNRM2 in BLAS. If we use this one to compute 2-norms, then NEWUOA
behaves the same as when we use SQRT(SUM OF SQUARES) naively. Therefore, maybe we should simply
implement NORM in this naive way. TO be checked.


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Function DNRM2 LAPACK 3.11.0:
!https://netlib.org/lapack/explore-html/df/d28/group__single__blas__level1_gab5393665c8f0e7d5de9bd1dd2ff0d9d0.html
function dnrm2(n, x, incx)
integer, parameter :: wp = kind(1.D0)
real(wp) :: DNRM2
!
!  -- Reference BLAS level1 routine (version 3.9.1) --
!  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     March 2021
!
!  .. Constants ..
real(wp), parameter :: zero = 0.0_WP
real(wp), parameter :: one = 1.0_WP
real(wp), parameter :: maxN = huge(0.0_WP)
!  ..
!  .. Blue's scaling constants ..
real(wp), parameter :: tsml = real(radix(0._WP), wp)**ceiling( &
                       (minexponent(0._WP) - 1) * 0.5_WP)
real(wp), parameter :: tbig = real(radix(0._WP), wp)**floor( &
                       (maxexponent(0._WP) - digits(0._WP) + 1) * 0.5_WP)
real(wp), parameter :: ssml = real(radix(0._WP), wp)**(-floor( &
                                                       (minexponent(0._WP) - digits(0._WP)) * 0.5_WP))
real(wp), parameter :: sbig = real(radix(0._WP), wp)**(-ceiling( &
                                                       (maxexponent(0._WP) + digits(0._WP) - 1) * 0.5_WP))
!  ..
!  .. Scalar Arguments ..
integer :: incx, n
!  ..
!  .. Array Arguments ..
real(wp) :: x(*)
!  ..
!  .. Local Scalars ..
integer :: i, ix
logical :: notbig
real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
!
!  Quick return if possible
!
dnrm2 = zero
if (n <= 0) return
!
scl = one
sumsq = zero
!
!  Compute the sum of squares in 3 accumulators:
!     abig -- sums of squares scaled down to avoid overflow
!     asml -- sums of squares scaled up to avoid underflow
!     amed -- sums of squares that do not require scaling
!  The thresholds and multipliers are
!     tbig -- values bigger than this are scaled down by sbig
!     tsml -- values smaller than this are scaled up by ssml
!
notbig = .true.
asml = zero
amed = zero
abig = zero
ix = 1
if (incx < 0) ix = 1 - (n - 1) * incx
do i = 1, n
    ax = abs(x(ix))
    if (ax > tbig) then
        abig = abig + (ax * sbig)**2
        notbig = .false.
    else if (ax < tsml) then
        if (notbig) asml = asml + (ax * ssml)**2
    else
        amed = amed + ax**2
    end if
    ix = ix + incx
end do
!
!  Combine abig and amed or amed and asml if more than one
!  accumulator was used.
!
if (abig > zero) then
!
!     Combine abig and amed if abig > 0.
!
    if ((amed > zero) .or. (amed > maxn) .or. (amed /= amed)) then
        abig = abig + (amed * sbig) * sbig
    end if
    scl = one / sbig
    sumsq = abig
else if (asml > zero) then
!
!     Combine amed and asml if asml > 0.
!
    if ((amed > zero) .or. (amed > maxn) .or. (amed /= amed)) then
        amed = sqrt(amed)
        asml = sqrt(asml) / ssml
        if (asml > amed) then
            ymin = amed
            ymax = asml
        else
            ymin = asml
            ymax = amed
        end if
        scl = one
        sumsq = ymax**2 * (one + (ymin / ymax)**2)
    else
        scl = one / ssml
        sumsq = asml
    end if
else
!
!     Otherwise all values are mid-range
!
    scl = one
    sumsq = amed
end if
dnrm2 = scl * sqrt(sumsq)
return
end function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
