!> \brief \b DNRM2
!
! =========== DOCUMENTATION ===========
!
! Online html documentation available at
! http://www.netlib.org/lapack/explore-html/
!
! Definition:
! ===========
!
! DOUBLE PRECISION FUNCTION DNRM2(N,X,INCX)
!
! .. Scalar Arguments ..
! INTEGER INCX,N
! ..
! .. Array Arguments ..
! DOUBLE PRECISION X(*)
! ..
!
!
!> \par Purpose:
! =============
!>
!> \verbatim
!>
!> DNRM2 returns the euclidean norm of a vector via the function
!> name, so that
!>
!> DNRM2 := sqrt( x'*x )
!> \endverbatim
!
! Arguments:
! ==========
!
!> \param[in] N
!> \verbatim
!> N is INTEGER
!> number of elements in input vector(s)
!> \endverbatim
!>
!> \param[in] X
!> \verbatim
!> X is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!> \endverbatim
!>
!> \param[in] INCX
!> \verbatim
!> INCX is INTEGER, storage spacing between elements of X
!> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
!> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
!> If INCX = 0, x isn't a vector so there is no need to call
!> this subroutine. If you call it anyway, it will count x(1)
!> in the vector norm N times.
!> \endverbatim
!
! Authors:
! ========
!
!> \author Edward Anderson, Lockheed Martin
!
!> \date August 2016
!
!> \ingroup nrm2
!
!> \par Contributors:
! ==================
!>
!> Weslley Pereira, University of Colorado Denver, USA
!
!> \par Further Details:
! =====================
!>
!> \verbatim
!>
!> Anderson E. (2017)
!> Algorithm 978: Safe Scaling in the Level 1 BLAS
!> ACM Trans Math Softw 44:1--28
!> https://doi.org/10.1145/3061665
!>
!> Blue, James L. (1978)
!> A Portable Fortran Program to Find the Euclidean Norm of a Vector
!> ACM Trans Math Softw 4:15--23
!> https://doi.org/10.1145/355769.355771
!>
!> \endverbatim
!>
! =====================================================================
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
!> \brief \b DROTG
!
! =========== DOCUMENTATION ===========
!
! Online html documentation available at
! http://www.netlib.org/lapack/explore-html/
!
!> \par Purpose:
! =============
!>
!> \verbatim
!>
!> DROTG constructs a plane rotation
!> [ c s ] [ a ] = [ r ]
!> [ -s c ] [ b ] [ 0 ]
!> satisfying c**2 + s**2 = 1.
!>
!> The computation uses the formulas
!> sigma = sgn(a) if |a| > |b|
!> = sgn(b) if |b| >= |a|
!> r = sigma*sqrt( a**2 + b**2 )
!> c = 1; s = 0 if r = 0
!> c = a/r; s = b/r if r != 0
!> The subroutine also computes
!> z = s if |a| > |b|,
!> = 1/c if |b| >= |a| and c != 0
!> = 1 if c = 0
!> This allows c and s to be reconstructed from z as follows:
!> If z = 1, set c = 0, s = 1.
!> If |z| < 1, set c = sqrt(1 - z**2) and s = z.
!> If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2).
!>
!> \endverbatim
!>
!> @see lartg, @see lartgp
!
! Arguments:
! ==========
!
!> \param[in,out] A
!> \verbatim
!> A is DOUBLE PRECISION
!> On entry, the scalar a.
!> On exit, the scalar r.
!> \endverbatim
!>
!> \param[in,out] B
!> \verbatim
!> B is DOUBLE PRECISION
!> On entry, the scalar b.
!> On exit, the scalar z.
!> \endverbatim
!>
!> \param[out] C
!> \verbatim
!> C is DOUBLE PRECISION
!> The scalar c.
!> \endverbatim
!>
!> \param[out] S
!> \verbatim
!> S is DOUBLE PRECISION
!> The scalar s.
!> \endverbatim
!
! Authors:
! ========
!
!> \author Edward Anderson, Lockheed Martin
!
!> \par Contributors:
! ==================
!>
!> Weslley Pereira, University of Colorado Denver, USA
!
!> \ingroup rotg
!
!> \par Further Details:
! =====================
!>
!> \verbatim
!>
!> Anderson E. (2017)
!> Algorithm 978: Safe Scaling in the Level 1 BLAS
!> ACM Trans Math Softw 44:1--28
!> https://doi.org/10.1145/3061665
!>
!> \endverbatim
!
! =====================================================================
subroutine DROTG( a, b, c, s )
integer, parameter :: wp = kind(1.d0)
!
! -- Reference BLAS level1 routine --
! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!
! .. Constants ..
real(wp), parameter :: zero = 0.0_wp
real(wp), parameter :: one = 1.0_wp
! ..
! .. Scaling constants ..
real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( &
minexponent(0._wp)-1, &
1-maxexponent(0._wp) &
)
real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( &
1-minexponent(0._wp), &
maxexponent(0._wp)-1 &
)
! ..
! .. Scalar Arguments ..
real(wp) :: a, b, c, s
! ..
! .. Local Scalars ..
real(wp) :: anorm, bnorm, scl, sigma, r, z
! ..
anorm = abs(a)
bnorm = abs(b)
if( bnorm == zero ) then
c = one
s = zero
b = zero
else if( anorm == zero ) then
c = zero
s = one
a = b
b = one
else
scl = min( safmax, max( safmin, anorm, bnorm ) )
if( anorm > bnorm ) then
sigma = sign(one,a)
else
sigma = sign(one,b)
end if
r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) )
c = a/r
s = b/r
if( anorm > bnorm ) then
z = s
else if( c /= zero ) then
z = one/c
else
z = one
end if
a = r
b = z
end if
return
end subroutine