kopia lustrzana https://github.com/Dsplib/libdspl-2.0
210 wiersze
5.2 KiB
Fortran
210 wiersze
5.2 KiB
Fortran
!> \brief \b SCNRM2
|
|
!
|
|
! =========== DOCUMENTATION ===========
|
|
!
|
|
! Online html documentation available at
|
|
! http://www.netlib.org/lapack/explore-html/
|
|
!
|
|
! Definition:
|
|
! ===========
|
|
!
|
|
! REAL FUNCTION SCNRM2(N,X,INCX)
|
|
!
|
|
! .. Scalar Arguments ..
|
|
! INTEGER INCX,N
|
|
! ..
|
|
! .. Array Arguments ..
|
|
! COMPLEX X(*)
|
|
! ..
|
|
!
|
|
!
|
|
!> \par Purpose:
|
|
! =============
|
|
!>
|
|
!> \verbatim
|
|
!>
|
|
!> SCNRM2 returns the euclidean norm of a vector via the function
|
|
!> name, so that
|
|
!>
|
|
!> SCNRM2 := sqrt( x**H*x )
|
|
!> \endverbatim
|
|
!
|
|
! Arguments:
|
|
! ==========
|
|
!
|
|
!> \param[in] N
|
|
!> \verbatim
|
|
!> N is INTEGER
|
|
!> number of elements in input vector(s)
|
|
!> \endverbatim
|
|
!>
|
|
!> \param[in] X
|
|
!> \verbatim
|
|
!> X is COMPLEX array, dimension (N)
|
|
!> complex vector with N elements
|
|
!> \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 single_blas_level1
|
|
!
|
|
!> \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 SCNRM2( n, x, incx )
|
|
integer, parameter :: wp = kind(1.e0)
|
|
real(wp) :: SCNRM2
|
|
!
|
|
! -- 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 ccaling 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) - 1) * 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 ..
|
|
complex(wp) :: x(*)
|
|
! ..
|
|
! .. Local Scalars ..
|
|
integer :: i, ix
|
|
logical :: notbig
|
|
real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
|
|
!
|
|
! Quick return if possible
|
|
!
|
|
SCNRM2 = 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(real(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
|
|
ax = abs(aimag(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
|
|
SCNRM2 = scl*sqrt( sumsq )
|
|
return
|
|
end function
|