From 9352192411f858f7b35223fc77ecb55f829e50a1 Mon Sep 17 00:00:00 2001 From: "weslley.spereira" Date: Fri, 19 Feb 2021 14:14:19 -0300 Subject: [PATCH 01/10] Add safe scaling xLASSQ routines from https://doi.org/10.1145/3061665; Let the compiler determine the Fortran layout by the file extension --- CMakeLists.txt | 3 - SRC/CMakeLists.txt | 12 +-- SRC/classq.f90 | 170 +++++++++++++++++++++++++++++++++++++++++ SRC/dlassq.f90 | 161 ++++++++++++++++++++++++++++++++++++++ SRC/la_constants.f90 | 40 ++++++++++ SRC/la_constants32.f90 | 40 ++++++++++ SRC/la_xisnan.f90 | 59 ++++++++++++++ SRC/slassq.f90 | 161 ++++++++++++++++++++++++++++++++++++++ SRC/zlassq.f90 | 170 +++++++++++++++++++++++++++++++++++++++++ 9 files changed, 808 insertions(+), 8 deletions(-) create mode 100644 SRC/classq.f90 create mode 100644 SRC/dlassq.f90 create mode 100644 SRC/la_constants.f90 create mode 100644 SRC/la_constants32.f90 create mode 100644 SRC/la_xisnan.f90 create mode 100644 SRC/slassq.f90 create mode 100644 SRC/zlassq.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 73ab55f3d2..b36fec1540 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,9 +16,6 @@ set(CMAKE_MODULE_PATH "${LAPACK_SOURCE_DIR}/CMAKE" ${CMAKE_MODULE_PATH}) # Export all symbols on Windows when building shared libraries SET(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS TRUE) -# Tell CMake that our Fortran sources are written in fixed format. -set(CMAKE_Fortran_FORMAT FIXED) - # Set a default build type if none was specified if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) message(STATUS "Setting build type to 'Release' as none was specified.") diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index bd3be3c6e8..e15eca69bc 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -36,11 +36,12 @@ ####################################################################### set(ALLAUX ilaenv.f ilaenv2stage.f ieeeck.f lsamen.f iparmq.f iparam2stage.F - ilaprec.f ilatrans.f ilauplo.f iladiag.f chla_transtype.f + ilaprec.f ilatrans.f ilauplo.f iladiag.f chla_transtype.f la_xisnan.f90 ../INSTALL/ilaver.f ../INSTALL/lsame.f xerbla.f xerbla_array.f ../INSTALL/slamch.f) set(SCLAUX + la_constants32.f90 sbdsdc.f sbdsqr.f sdisna.f slabad.f slacpy.f sladiv.f slae2.f slaebz.f slaed0.f slaed1.f slaed2.f slaed3.f slaed4.f slaed5.f slaed6.f @@ -53,12 +54,13 @@ set(SCLAUX slasd0.f slasd1.f slasd2.f slasd3.f slasd4.f slasd5.f slasd6.f slasd7.f slasd8.f slasda.f slasdq.f slasdt.f slaset.f slasq1.f slasq2.f slasq3.f slasq4.f slasq5.f slasq6.f - slasr.f slasrt.f slassq.f slasv2.f spttrf.f sstebz.f sstedc.f + slasr.f slasrt.f slassq.f90 slasv2.f spttrf.f sstebz.f sstedc.f ssteqr.f ssterf.f slaisnan.f sisnan.f slartgp.f slartgs.f ${SECOND_SRC}) set(DZLAUX + la_constants.f90 dbdsdc.f dbdsqr.f ddisna.f dlabad.f dlacpy.f dladiv.f dlae2.f dlaebz.f dlaed0.f dlaed1.f dlaed2.f dlaed3.f dlaed4.f dlaed5.f dlaed6.f @@ -71,7 +73,7 @@ set(DZLAUX dlasd0.f dlasd1.f dlasd2.f dlasd3.f dlasd4.f dlasd5.f dlasd6.f dlasd7.f dlasd8.f dlasda.f dlasdq.f dlasdt.f dlaset.f dlasq1.f dlasq2.f dlasq3.f dlasq4.f dlasq5.f dlasq6.f - dlasr.f dlasrt.f dlassq.f dlasv2.f dpttrf.f dstebz.f dstedc.f + dlasr.f dlasrt.f dlassq.f90 dlasv2.f dpttrf.f dstebz.f dstedc.f dsteqr.f dsterf.f dlaisnan.f disnan.f dlartgp.f dlartgs.f ../INSTALL/dlamch.f ${DSECOND_SRC}) @@ -209,7 +211,7 @@ set(CLASRC claqsp.f claqsy.f clar1v.f clar2v.f ilaclr.f ilaclc.f clarf.f clarfb.f clarfb_gett.f clarfg.f clarfgp.f clarft.f clarfx.f clarfy.f clargv.f clarnv.f clarrv.f clartg.f clartv.f - clarz.f clarzb.f clarzt.f clascl.f claset.f clasr.f classq.f + clarz.f clarzb.f clarzt.f clascl.f claset.f clasr.f classq.f90 claswp.f clasyf.f clasyf_rook.f clasyf_rk.f clasyf_aa.f clatbs.f clatdf.f clatps.f clatrd.f clatrs.f clatrz.f clauu2.f clauum.f cpbcon.f cpbequ.f cpbrfs.f cpbstf.f cpbsv.f @@ -405,7 +407,7 @@ set(ZLASRC zlarfg.f zlarfgp.f zlarft.f zlarfx.f zlarfy.f zlargv.f zlarnv.f zlarrv.f zlartg.f zlartv.f zlarz.f zlarzb.f zlarzt.f zlascl.f zlaset.f zlasr.f - zlassq.f zlaswp.f zlasyf.f zlasyf_rook.f zlasyf_rk.f zlasyf_aa.f + zlassq.f90 zlaswp.f zlasyf.f zlasyf_rook.f zlasyf_rk.f zlasyf_aa.f zlatbs.f zlatdf.f zlatps.f zlatrd.f zlatrs.f zlatrz.f zlauu2.f zlauum.f zpbcon.f zpbequ.f zpbrfs.f zpbstf.f zpbsv.f zpbsvx.f zpbtf2.f zpbtrf.f zpbtrs.f zpocon.f zpoequ.f zporfs.f diff --git a/SRC/classq.f90 b/SRC/classq.f90 new file mode 100644 index 0000000000..75def12baf --- /dev/null +++ b/SRC/classq.f90 @@ -0,0 +1,170 @@ +subroutine CLASSQ( n, x, incx, scl, sumsq ) + use LA_CONSTANTS32, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_XISNAN +! +! LAPACK auxiliary routine +! Based on Blue's algorithm from ACM TOMS, March 1978 +! E. Anderson +! August 9, 2016 +! +! .. Scalar Arguments .. + integer :: incx, n + real(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + complex(wp) :: x(*) +! .. +! +! Purpose +! ======= +! +! CLASSQ returns the values scl and smsq such that +! +! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +! +! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +! assumed to be non-negative and scl returns the value +! +! scl = max( scale, abs( x( i ) ) ). +! +! scale and sumsq must be supplied in SCL and SUMSQ and +! scl and smsq are overwritten on SCL and SUMSQ respectively. +! +! Arguments +! ========= +! +! N (input) INTEGER +! The number of elements of the vector x. +! +! X (input) COMPLEX array, dimension (1+(N-1)*abs(INCX)) +! The n-element vector x. +! +! INCX (input) INTEGER +! The increment between successive values of the vector 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. +! +! SCL (input/output) REAL +! On entry, the value scale in the equation above. +! On exit, SCL is overwritten with scl , the scaling factor +! for the sum of squares. +! +! SUMSQ (input/output) REAL +! On entry, the value sumsq in the equation above. +! On exit, SUMSQ is overwritten with smsq , the basic sum of +! squares from which scl has been factored out. +! +! ===================================================================== +! +! .. Local Scalars .. + integer :: i, ix + logical :: notbig + real(wp) :: abig, amed, asml, ax, ymax, ymin +! .. +! +! Quick return if possible +! + if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return + if( sumsq == zero ) scl = one + if( scl == zero ) then + scl = one + sumsq = zero + end if + if (n <= 0) then + return + end if +! +! 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 +! +! Put the existing sum of squares into one of the accumulators +! + if( sumsq > zero ) then + ax = scl*sqrt( sumsq ) + 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 + end if +! +! 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. LA_ISNAN(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. LA_ISNAN(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 + return +end subroutine diff --git a/SRC/dlassq.f90 b/SRC/dlassq.f90 new file mode 100644 index 0000000000..a9621ea01b --- /dev/null +++ b/SRC/dlassq.f90 @@ -0,0 +1,161 @@ +subroutine DLASSQ( n, x, incx, scl, sumsq ) + use LA_CONSTANTS, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_XISNAN +! +! LAPACK auxiliary routine +! Based on Blue's algorithm from ACM TOMS, March 1978 +! E. Anderson +! August 9, 2016 +! +! .. Scalar Arguments .. + integer :: incx, n + real(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + real(wp) :: x(*) +! .. +! +! Purpose +! ======= +! +! DLASSQ returns the values scl and smsq such that +! +! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +! +! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +! assumed to be non-negative and scl returns the value +! +! scl = max( scale, abs( x( i ) ) ). +! +! scale and sumsq must be supplied in SCL and SUMSQ and +! scl and smsq are overwritten on SCL and SUMSQ respectively. +! +! Arguments +! ========= +! +! N (input) INTEGER +! The number of elements of the vector x. +! +! X (input) REAL array, dimension (1+(N-1)*abs(INCX)) +! The n-element vector x. +! +! INCX (input) INTEGER +! The increment between successive values of the vector 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. +! +! SCL (input/output) REAL +! On entry, the value scale in the equation above. +! On exit, SCL is overwritten with scl , the scaling factor +! for the sum of squares. +! +! SUMSQ (input/output) REAL +! On entry, the value sumsq in the equation above. +! On exit, SUMSQ is overwritten with smsq , the basic sum of +! squares from which scl has been factored out. +! +! ===================================================================== +! +! .. Local Scalars .. + integer :: i, ix + logical :: notbig + real(wp) :: abig, amed, asml, ax, ymax, ymin +! .. +! +! Quick return if possible +! + if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return + if( sumsq == zero ) scl = one + if( scl == zero ) then + scl = one + sumsq = zero + end if + if (n <= 0) then + return + end if +! +! 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 +! +! Put the existing sum of squares into one of the accumulators +! + if( sumsq > zero ) then + ax = scl*sqrt( sumsq ) + 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 + end if +! +! 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. LA_ISNAN(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. LA_ISNAN(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 + return +end subroutine diff --git a/SRC/la_constants.f90 b/SRC/la_constants.f90 new file mode 100644 index 0000000000..8306d1ed83 --- /dev/null +++ b/SRC/la_constants.f90 @@ -0,0 +1,40 @@ +module LA_CONSTANTS +! +! -- BLAS/LAPACK module -- +! May 06, 2016 +! +! Standard constants +! + integer, parameter :: wp = 8 + real(wp), parameter :: zero = 0.0_wp + real(wp), parameter :: half = 0.5_wp + real(wp), parameter :: one = 1.0_wp + real(wp), parameter :: two = 2.0_wp + real(wp), parameter :: three = 3.0_wp + real(wp), parameter :: four = 4.0_wp + real(wp), parameter :: eight = 8.0_wp + real(wp), parameter :: ten = 10.0_wp + complex(wp), parameter :: czero = ( 0.0_wp, 0.0_wp ) + complex(wp), parameter :: chalf = ( 0.5_wp, 0.0_wp ) + complex(wp), parameter :: cone = ( 1.0_wp, 0.0_wp ) + character*1, parameter :: sprefix = 'D' + character*1, parameter :: cprefix = 'Z' +! +! Model parameters +! + real(wp), parameter :: eps = 0.11102230246251565404E-015_wp + real(wp), parameter :: ulp = 0.22204460492503130808E-015_wp + real(wp), parameter :: safmin = 0.22250738585072013831E-307_wp + real(wp), parameter :: safmax = 0.44942328371557897693E+308_wp + real(wp), parameter :: smlnum = 0.10020841800044863890E-291_wp + real(wp), parameter :: bignum = 0.99792015476735990583E+292_wp + real(wp), parameter :: rtmin = 0.10010415475915504622E-145_wp + real(wp), parameter :: rtmax = 0.99895953610111751404E+146_wp +! +! Blue's scaling constants +! + real(wp), parameter :: tsml = 0.14916681462400413487E-153_wp + real(wp), parameter :: tbig = 0.19979190722022350281E+147_wp + real(wp), parameter :: ssml = 0.44989137945431963828E+162_wp + real(wp), parameter :: sbig = 0.11113793747425387417E-161_wp +end module LA_CONSTANTS diff --git a/SRC/la_constants32.f90 b/SRC/la_constants32.f90 new file mode 100644 index 0000000000..d495709e92 --- /dev/null +++ b/SRC/la_constants32.f90 @@ -0,0 +1,40 @@ +module LA_CONSTANTS32 +! +! -- BLAS/LAPACK module -- +! May 06, 2016 +! +! Standard constants +! + integer, parameter :: wp = 4 + real(wp), parameter :: zero = 0.0_wp + real(wp), parameter :: half = 0.5_wp + real(wp), parameter :: one = 1.0_wp + real(wp), parameter :: two = 2.0_wp + real(wp), parameter :: three = 3.0_wp + real(wp), parameter :: four = 4.0_wp + real(wp), parameter :: eight = 8.0_wp + real(wp), parameter :: ten = 10.0_wp + complex(wp), parameter :: czero = ( 0.0_wp, 0.0_wp ) + complex(wp), parameter :: chalf = ( 0.5_wp, 0.0_wp ) + complex(wp), parameter :: cone = ( 1.0_wp, 0.0_wp ) + character*1, parameter :: sprefix = 'S' + character*1, parameter :: cprefix = 'C' +! +! Model parameters +! + real(wp), parameter :: eps = 0.5960464478E-07_wp + real(wp), parameter :: ulp = 0.1192092896E-06_wp + real(wp), parameter :: safmin = 0.1175494351E-37_wp + real(wp), parameter :: safmax = 0.8507059173E+38_wp + real(wp), parameter :: smlnum = 0.9860761315E-31_wp + real(wp), parameter :: bignum = 0.1014120480E+32_wp + real(wp), parameter :: rtmin = 0.3140184864E-15_wp + real(wp), parameter :: rtmax = 0.3184525782E+16_wp +! +! Blue's scaling constants +! + real(wp), parameter :: tsml = 0.1084202172E-18_wp + real(wp), parameter :: tbig = 0.4503599627E+16_wp + real(wp), parameter :: ssml = 0.3777893186E+23_wp + real(wp), parameter :: sbig = 0.1323488980E-22_wp +end module LA_CONSTANTS32 diff --git a/SRC/la_xisnan.f90 b/SRC/la_xisnan.f90 new file mode 100644 index 0000000000..6ba91de0bd --- /dev/null +++ b/SRC/la_xisnan.f90 @@ -0,0 +1,59 @@ +module LA_XISNAN + interface LA_ISNAN + + module procedure SISNAN + module procedure DISNAN + + end interface + +contains + + logical function SISNAN( x ) + use LA_CONSTANTS32, only: wp +#ifdef USE_IEEE_INTRINSIC + use, intrinsic :: ieee_arithmetic +#elif USE_ISNAN + intrinsic :: isnan +#endif + real(wp) :: x +#ifdef USE_IEEE_INTRINSIC + sisnan = ieee_is_nan(x) +#elif USE_ISNAN + sisnan = isnan(x) +#else + sisnan = SLAISNAN(x,x) + + contains + logical function SLAISNAN( x, y ) + use LA_CONSTANTS32, only: wp + real(wp) :: x, y + SLAISNAN = ( x.ne.y ) + end function SLAISNAN +#endif + end function SISNAN + + logical function DISNAN( x ) + use LA_CONSTANTS, only: wp +#ifdef USE_IEEE_INTRINSIC + use, intrinsic :: ieee_arithmetic +#elif USE_ISNAN + intrinsic :: isnan +#endif + real(wp) :: x +#ifdef USE_IEEE_INTRINSIC + DISNAN = ieee_is_nan(x) +#elif USE_ISNAN + DISNAN = isnan(x) +#else + DISNAN = DLAISNAN(x,x) + + contains + logical function DLAISNAN( x, y ) + use LA_CONSTANTS, only: wp + real(wp) :: x, y + DLAISNAN = ( x.ne.y ) + end function DLAISNAN +#endif + end function DISNAN + +end module LA_XISNAN diff --git a/SRC/slassq.f90 b/SRC/slassq.f90 new file mode 100644 index 0000000000..f8a7f6b4a4 --- /dev/null +++ b/SRC/slassq.f90 @@ -0,0 +1,161 @@ +subroutine SLASSQ( n, x, incx, scl, sumsq ) + use LA_CONSTANTS32, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_XISNAN +! +! LAPACK auxiliary routine +! Based on Blue's algorithm from ACM TOMS, March 1978 +! E. Anderson +! August 9, 2016 +! +! .. Scalar Arguments .. + integer :: incx, n + real(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + real(wp) :: x(*) +! .. +! +! Purpose +! ======= +! +! SLASSQ returns the values scl and smsq such that +! +! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +! +! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +! assumed to be non-negative and scl returns the value +! +! scl = max( scale, abs( x( i ) ) ). +! +! scale and sumsq must be supplied in SCL and SUMSQ and +! scl and smsq are overwritten on SCL and SUMSQ respectively. +! +! Arguments +! ========= +! +! N (input) INTEGER +! The number of elements of the vector x. +! +! X (input) REAL array, dimension (1+(N-1)*abs(INCX)) +! The n-element vector x. +! +! INCX (input) INTEGER +! The increment between successive values of the vector 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. +! +! SCL (input/output) REAL +! On entry, the value scale in the equation above. +! On exit, SCL is overwritten with scl , the scaling factor +! for the sum of squares. +! +! SUMSQ (input/output) REAL +! On entry, the value sumsq in the equation above. +! On exit, SUMSQ is overwritten with smsq , the basic sum of +! squares from which scl has been factored out. +! +! ===================================================================== +! +! .. Local Scalars .. + integer :: i, ix + logical :: notbig + real(wp) :: abig, amed, asml, ax, ymax, ymin +! .. +! +! Quick return if possible +! + if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return + if( sumsq == zero ) scl = one + if( scl == zero ) then + scl = one + sumsq = zero + end if + if (n <= 0) then + return + end if +! +! 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 +! +! Put the existing sum of squares into one of the accumulators +! + if( sumsq > zero ) then + ax = scl*sqrt( sumsq ) + 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 + end if +! +! 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. LA_ISNAN(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. LA_ISNAN(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 + return +end subroutine diff --git a/SRC/zlassq.f90 b/SRC/zlassq.f90 new file mode 100644 index 0000000000..002e198eda --- /dev/null +++ b/SRC/zlassq.f90 @@ -0,0 +1,170 @@ +subroutine ZLASSQ( n, x, incx, scl, sumsq ) + use LA_CONSTANTS, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_XISNAN +! +! LAPACK auxiliary routine +! Based on Blue's algorithm from ACM TOMS, March 1978 +! E. Anderson +! August 9, 2016 +! +! .. Scalar Arguments .. + integer :: incx, n + real(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + complex(wp) :: x(*) +! .. +! +! Purpose +! ======= +! +! ZLASSQ returns the values scl and smsq such that +! +! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +! +! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +! assumed to be non-negative and scl returns the value +! +! scl = max( scale, abs( x( i ) ) ). +! +! scale and sumsq must be supplied in SCL and SUMSQ and +! scl and smsq are overwritten on SCL and SUMSQ respectively. +! +! Arguments +! ========= +! +! N (input) INTEGER +! The number of elements of the vector x. +! +! X (input) COMPLEX array, dimension (1+(N-1)*abs(INCX)) +! The n-element vector x. +! +! INCX (input) INTEGER +! The increment between successive values of the vector 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. +! +! SCL (input/output) REAL +! On entry, the value scale in the equation above. +! On exit, SCL is overwritten with scl , the scaling factor +! for the sum of squares. +! +! SUMSQ (input/output) REAL +! On entry, the value sumsq in the equation above. +! On exit, SUMSQ is overwritten with smsq , the basic sum of +! squares from which scl has been factored out. +! +! ===================================================================== +! +! .. Local Scalars .. + integer :: i, ix + logical :: notbig + real(wp) :: abig, amed, asml, ax, ymax, ymin +! .. +! +! Quick return if possible +! + if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return + if( sumsq == zero ) scl = one + if( scl == zero ) then + scl = one + sumsq = zero + end if + if (n <= 0) then + return + end if +! +! 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 +! +! Put the existing sum of squares into one of the accumulators +! + if( sumsq > zero ) then + ax = scl*sqrt( sumsq ) + 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 + end if +! +! 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. LA_ISNAN(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. LA_ISNAN(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 + return +end subroutine From a0fc2cc93ba779a34d11cb1fd4145eae1ca9af39 Mon Sep 17 00:00:00 2001 From: "weslley.spereira" Date: Fri, 19 Feb 2021 17:54:50 -0300 Subject: [PATCH 02/10] Fix Makefile and Meson builds; Rename original xlassq.f files to xlassq_lapackv390.f --- SRC/Makefile | 6 +++++- SRC/{classq.f => classq_lapackv390.f} | 0 SRC/{dlassq.f => dlassq_lapackv390.f} | 0 SRC/{slassq.f => slassq_lapackv390.f} | 0 SRC/{zlassq.f => zlassq_lapackv390.f} | 0 5 files changed, 5 insertions(+), 1 deletion(-) rename SRC/{classq.f => classq_lapackv390.f} (100%) rename SRC/{dlassq.f => dlassq_lapackv390.f} (100%) rename SRC/{slassq.f => slassq_lapackv390.f} (100%) rename SRC/{zlassq.f => zlassq_lapackv390.f} (100%) diff --git a/SRC/Makefile b/SRC/Makefile index 3895fe5963..f860394ed2 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -57,9 +57,11 @@ TOPSRCDIR = .. include $(TOPSRCDIR)/make.inc -.SUFFIXES: .F .o +.SUFFIXES: .F .f90 .o .F.o: $(FC) $(FFLAGS) -c -o $@ $< +.f90.o: + $(FC) $(FFLAGS) -c -o $@ $< ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \ iparmq.o iparam2stage.o \ @@ -67,6 +69,7 @@ ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \ ../INSTALL/ilaver.o ../INSTALL/lsame.o ../INSTALL/slamch.o SCLAUX = \ + la_constants32.o \ sbdsdc.o \ sbdsqr.o sdisna.o slabad.o slacpy.o sladiv.o slae2.o slaebz.o \ slaed0.o slaed1.o slaed2.o slaed3.o slaed4.o slaed5.o slaed6.o \ @@ -85,6 +88,7 @@ SCLAUX = \ ../INSTALL/second_$(TIMER).o DZLAUX = \ + la_constants.o \ dbdsdc.o \ dbdsqr.o ddisna.o dlabad.o dlacpy.o dladiv.o dlae2.o dlaebz.o \ dlaed0.o dlaed1.o dlaed2.o dlaed3.o dlaed4.o dlaed5.o dlaed6.o \ diff --git a/SRC/classq.f b/SRC/classq_lapackv390.f similarity index 100% rename from SRC/classq.f rename to SRC/classq_lapackv390.f diff --git a/SRC/dlassq.f b/SRC/dlassq_lapackv390.f similarity index 100% rename from SRC/dlassq.f rename to SRC/dlassq_lapackv390.f diff --git a/SRC/slassq.f b/SRC/slassq_lapackv390.f similarity index 100% rename from SRC/slassq.f rename to SRC/slassq_lapackv390.f diff --git a/SRC/zlassq.f b/SRC/zlassq_lapackv390.f similarity index 100% rename from SRC/zlassq.f rename to SRC/zlassq_lapackv390.f From ad21f0c3b96e4e22b5dca4be57a87e4e7c41f46d Mon Sep 17 00:00:00 2001 From: "weslley.spereira" Date: Fri, 19 Feb 2021 20:16:00 -0300 Subject: [PATCH 03/10] Fix makefile --- SRC/Makefile | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/SRC/Makefile b/SRC/Makefile index f860394ed2..fc601bb2a1 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -57,14 +57,22 @@ TOPSRCDIR = .. include $(TOPSRCDIR)/make.inc -.SUFFIXES: .F .f90 .o -.F.o: +ALLMOD = la_xisnan.mod la_constants.mod la_constants32.mod + +.SUFFIXES: .f .F .f90 .o .mod +%.o: %.f $(ALLMOD) $(FC) $(FFLAGS) -c -o $@ $< -.f90.o: +%.o: %.F $(ALLMOD) $(FC) $(FFLAGS) -c -o $@ $< - +%.o: %.f90 $(ALLMOD) + $(FC) $(FFLAGS) -c -o $@ $< +%.o: %.F90 $(ALLMOD) + $(FC) $(FFLAGS) -c -o $@ $< +.o.mod: + @true + ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \ - iparmq.o iparam2stage.o \ + iparmq.o iparam2stage.o la_xisnan.o \ ilaprec.o ilatrans.o ilauplo.o iladiag.o chla_transtype.o \ ../INSTALL/ilaver.o ../INSTALL/lsame.o ../INSTALL/slamch.o @@ -612,7 +620,7 @@ endif .PHONY: clean cleanobj cleanlib clean: cleanobj cleanlib cleanobj: - rm -f *.o DEPRECATED/*.o + rm -f *.o *.mod DEPRECATED/*.o DEPRECATED/*.mod cleanlib: rm -f $(LAPACKLIB) @@ -622,3 +630,11 @@ sla_wwaddw.o: sla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< dla_wwaddw.o: dla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< cla_wwaddw.o: cla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< zla_wwaddw.o: zla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< + +# Modules +la_xisnan.o: la_xisnan.f90 la_constants32.mod la_constants.mod + $(FC) $(FFLAGS) -c -o $@ $< +la_constants32.o: la_constants32.f90 + $(FC) $(FFLAGS) -c -o $@ $< +la_constants.o: la_constants.f90 + $(FC) $(FFLAGS) -c -o $@ $< \ No newline at end of file From ef4a7529f68b7fe56cd0d6463a430a008e8c278e Mon Sep 17 00:00:00 2001 From: "weslley.spereira" Date: Mon, 22 Feb 2021 11:10:08 -0300 Subject: [PATCH 04/10] Move original xLASSQ to DEPRECATED; Add EOF in SRC/Makefile; Correct la_xisnan extension --- SRC/CMakeLists.txt | 2 +- SRC/{ => DEPRECATED}/classq_lapackv390.f | 3 +++ SRC/{ => DEPRECATED}/dlassq_lapackv390.f | 3 +++ SRC/{ => DEPRECATED}/slassq_lapackv390.f | 3 +++ SRC/{ => DEPRECATED}/zlassq_lapackv390.f | 3 +++ SRC/Makefile | 7 ++++--- SRC/{la_xisnan.f90 => la_xisnan.F90} | 0 7 files changed, 17 insertions(+), 4 deletions(-) rename SRC/{ => DEPRECATED}/classq_lapackv390.f (96%) rename SRC/{ => DEPRECATED}/dlassq_lapackv390.f (96%) rename SRC/{ => DEPRECATED}/slassq_lapackv390.f (96%) rename SRC/{ => DEPRECATED}/zlassq_lapackv390.f (96%) rename SRC/{la_xisnan.f90 => la_xisnan.F90} (100%) diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index e15eca69bc..deddeb4d8a 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -36,7 +36,7 @@ ####################################################################### set(ALLAUX ilaenv.f ilaenv2stage.f ieeeck.f lsamen.f iparmq.f iparam2stage.F - ilaprec.f ilatrans.f ilauplo.f iladiag.f chla_transtype.f la_xisnan.f90 + ilaprec.f ilatrans.f ilauplo.f iladiag.f chla_transtype.f la_xisnan.F90 ../INSTALL/ilaver.f ../INSTALL/lsame.f xerbla.f xerbla_array.f ../INSTALL/slamch.f) diff --git a/SRC/classq_lapackv390.f b/SRC/DEPRECATED/classq_lapackv390.f similarity index 96% rename from SRC/classq_lapackv390.f rename to SRC/DEPRECATED/classq_lapackv390.f index a1c13309bc..90e323b64c 100644 --- a/SRC/classq_lapackv390.f +++ b/SRC/DEPRECATED/classq_lapackv390.f @@ -99,6 +99,9 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * +*> This implementation of CLASSQ has been deprecated with LAPACKv3.10. +*> A better version of CLASSQ was contributed by Ed Anderson and released in 3.10. +* *> \ingroup complexOTHERauxiliary * * ===================================================================== diff --git a/SRC/dlassq_lapackv390.f b/SRC/DEPRECATED/dlassq_lapackv390.f similarity index 96% rename from SRC/dlassq_lapackv390.f rename to SRC/DEPRECATED/dlassq_lapackv390.f index 5c24d00350..d3126951ca 100644 --- a/SRC/dlassq_lapackv390.f +++ b/SRC/DEPRECATED/dlassq_lapackv390.f @@ -96,6 +96,9 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * +*> This implementation of DLASSQ has been deprecated with LAPACKv3.10. +*> A better version of DLASSQ was contributed by Ed Anderson and released in 3.10. +* *> \ingroup OTHERauxiliary * * ===================================================================== diff --git a/SRC/slassq_lapackv390.f b/SRC/DEPRECATED/slassq_lapackv390.f similarity index 96% rename from SRC/slassq_lapackv390.f rename to SRC/DEPRECATED/slassq_lapackv390.f index 72503cb77a..73f63b5150 100644 --- a/SRC/slassq_lapackv390.f +++ b/SRC/DEPRECATED/slassq_lapackv390.f @@ -96,6 +96,9 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * +*> This implementation of SLASSQ has been deprecated with LAPACKv3.10. +*> A better version of SLASSQ was contributed by Ed Anderson and released in 3.10. +* *> \ingroup OTHERauxiliary * * ===================================================================== diff --git a/SRC/zlassq_lapackv390.f b/SRC/DEPRECATED/zlassq_lapackv390.f similarity index 96% rename from SRC/zlassq_lapackv390.f rename to SRC/DEPRECATED/zlassq_lapackv390.f index 02b0d7a14a..0133f2456f 100644 --- a/SRC/zlassq_lapackv390.f +++ b/SRC/DEPRECATED/zlassq_lapackv390.f @@ -99,6 +99,9 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * +*> This implementation of ZLASSQ has been deprecated with LAPACKv3.10. +*> A better version of ZLASSQ was contributed by Ed Anderson and released in 3.10. +* *> \ingroup complex16OTHERauxiliary * * ===================================================================== diff --git a/SRC/Makefile b/SRC/Makefile index fc601bb2a1..28d2fa9954 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -59,7 +59,7 @@ include $(TOPSRCDIR)/make.inc ALLMOD = la_xisnan.mod la_constants.mod la_constants32.mod -.SUFFIXES: .f .F .f90 .o .mod +.SUFFIXES: .f .F .f90 .F90 .o .mod %.o: %.f $(ALLMOD) $(FC) $(FFLAGS) -c -o $@ $< %.o: %.F $(ALLMOD) @@ -632,9 +632,10 @@ cla_wwaddw.o: cla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< zla_wwaddw.o: zla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< # Modules -la_xisnan.o: la_xisnan.f90 la_constants32.mod la_constants.mod +la_xisnan.o: la_xisnan.F90 la_constants32.mod la_constants.mod $(FC) $(FFLAGS) -c -o $@ $< la_constants32.o: la_constants32.f90 $(FC) $(FFLAGS) -c -o $@ $< la_constants.o: la_constants.f90 - $(FC) $(FFLAGS) -c -o $@ $< \ No newline at end of file + $(FC) $(FFLAGS) -c -o $@ $< + diff --git a/SRC/la_xisnan.f90 b/SRC/la_xisnan.F90 similarity index 100% rename from SRC/la_xisnan.f90 rename to SRC/la_xisnan.F90 From 3edef3b5e0acda99dc6bff04ab796be38e81b559 Mon Sep 17 00:00:00 2001 From: "weslley.spereira" Date: Mon, 1 Mar 2021 17:38:46 -0300 Subject: [PATCH 05/10] Add the module la_constants from @ecanesc with the suggestions of @zerothi --- SRC/CMakeLists.txt | 2 +- SRC/Makefile | 10 +-- SRC/la_constants.f90 | 171 ++++++++++++++++++++++++++++++++++--------- 3 files changed, 142 insertions(+), 41 deletions(-) diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index deddeb4d8a..070eeb8100 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -41,7 +41,7 @@ set(ALLAUX ilaenv.f ilaenv2stage.f ieeeck.f lsamen.f iparmq.f iparam2stage.F ../INSTALL/slamch.f) set(SCLAUX - la_constants32.f90 + la_constants.f90 sbdsdc.f sbdsqr.f sdisna.f slabad.f slacpy.f sladiv.f slae2.f slaebz.f slaed0.f slaed1.f slaed2.f slaed3.f slaed4.f slaed5.f slaed6.f diff --git a/SRC/Makefile b/SRC/Makefile index 28d2fa9954..20a2da5cac 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -57,7 +57,7 @@ TOPSRCDIR = .. include $(TOPSRCDIR)/make.inc -ALLMOD = la_xisnan.mod la_constants.mod la_constants32.mod +ALLMOD = la_xisnan.mod la_constants.mod .SUFFIXES: .f .F .f90 .F90 .o .mod %.o: %.f $(ALLMOD) @@ -70,14 +70,14 @@ ALLMOD = la_xisnan.mod la_constants.mod la_constants32.mod $(FC) $(FFLAGS) -c -o $@ $< .o.mod: @true - + ALLAUX = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \ iparmq.o iparam2stage.o la_xisnan.o \ ilaprec.o ilatrans.o ilauplo.o iladiag.o chla_transtype.o \ ../INSTALL/ilaver.o ../INSTALL/lsame.o ../INSTALL/slamch.o SCLAUX = \ - la_constants32.o \ + la_constants.o \ sbdsdc.o \ sbdsqr.o sdisna.o slabad.o slacpy.o sladiv.o slae2.o slaebz.o \ slaed0.o slaed1.o slaed2.o slaed3.o slaed4.o slaed5.o slaed6.o \ @@ -632,9 +632,7 @@ cla_wwaddw.o: cla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< zla_wwaddw.o: zla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< # Modules -la_xisnan.o: la_xisnan.F90 la_constants32.mod la_constants.mod - $(FC) $(FFLAGS) -c -o $@ $< -la_constants32.o: la_constants32.f90 +la_xisnan.o: la_xisnan.F90 la_constants.mod $(FC) $(FFLAGS) -c -o $@ $< la_constants.o: la_constants.f90 $(FC) $(FFLAGS) -c -o $@ $< diff --git a/SRC/la_constants.f90 b/SRC/la_constants.f90 index 8306d1ed83..0bf798195d 100644 --- a/SRC/la_constants.f90 +++ b/SRC/la_constants.f90 @@ -1,40 +1,143 @@ +!> \brief \b LA_CONSTANTS is a module for the scaling constants for the compiled Fortran single and double precisions +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \date May 2016 +! +!> \ingroup OTHERauxiliary +! +!> \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 +! module LA_CONSTANTS ! -! -- BLAS/LAPACK module -- -! May 06, 2016 -! -! Standard constants -! - integer, parameter :: wp = 8 - real(wp), parameter :: zero = 0.0_wp - real(wp), parameter :: half = 0.5_wp - real(wp), parameter :: one = 1.0_wp - real(wp), parameter :: two = 2.0_wp - real(wp), parameter :: three = 3.0_wp - real(wp), parameter :: four = 4.0_wp - real(wp), parameter :: eight = 8.0_wp - real(wp), parameter :: ten = 10.0_wp - complex(wp), parameter :: czero = ( 0.0_wp, 0.0_wp ) - complex(wp), parameter :: chalf = ( 0.5_wp, 0.0_wp ) - complex(wp), parameter :: cone = ( 1.0_wp, 0.0_wp ) - character*1, parameter :: sprefix = 'D' - character*1, parameter :: cprefix = 'Z' -! -! Model parameters -! - real(wp), parameter :: eps = 0.11102230246251565404E-015_wp - real(wp), parameter :: ulp = 0.22204460492503130808E-015_wp - real(wp), parameter :: safmin = 0.22250738585072013831E-307_wp - real(wp), parameter :: safmax = 0.44942328371557897693E+308_wp - real(wp), parameter :: smlnum = 0.10020841800044863890E-291_wp - real(wp), parameter :: bignum = 0.99792015476735990583E+292_wp - real(wp), parameter :: rtmin = 0.10010415475915504622E-145_wp - real(wp), parameter :: rtmax = 0.99895953610111751404E+146_wp +! -- LAPACK auxiliary module (version 3.9.0) -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! February 2021 +! +! Standard constants for + integer, parameter :: sp = kind(1.e0) +! + real(sp), parameter :: szero = 0.0_sp + real(sp), parameter :: shalf = 0.5_sp + real(sp), parameter :: sone = 1.0_sp + real(sp), parameter :: stwo = 2.0_sp + real(sp), parameter :: sthree = 3.0_sp + real(sp), parameter :: sfour = 4.0_sp + real(sp), parameter :: seight = 8.0_sp + real(sp), parameter :: sten = 10.0_sp + complex(sp), parameter :: czero = ( 0.0_sp, 0.0_sp ) + complex(sp), parameter :: chalf = ( 0.5_sp, 0.0_sp ) + complex(sp), parameter :: cone = ( 1.0_sp, 0.0_sp ) + character*1, parameter :: sprefix = 'S' + character*1, parameter :: cprefix = 'C' +! +! Scaling constants +! + real(sp), parameter :: sulp = epsilon(0._sp) + real(sp), parameter :: seps = sulp * 0.5_sp + real(sp), parameter :: ssafmin = real(radix(0._sp),sp)**max( & + minexponent(0._sp)-1, & + 1-maxexponent(0._sp) & + ) + real(sp), parameter :: ssafmax = sone / ssafmin + real(sp), parameter :: ssmlnum = ssafmin / sulp + real(sp), parameter :: sbignum = ssafmax * sulp + real(sp), parameter :: srtmin = sqrt(ssmlnum) + real(sp), parameter :: srtmax = sqrt(sbignum) +! +! Blue's scaling constants +! + real(sp), parameter :: stsml = real(radix(0._sp), sp)**ceiling( & + real(( minexponent(0._sp) - 1_sp ) / 2, sp) & + ) + real(sp), parameter :: stbig = real(radix(0._sp), sp)**floor( & + real(( maxexponent(0._sp) - digits(0._sp) + 1_sp) / 2, sp) & + ) +! ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771 + real(sp), parameter :: sssml = real(radix(0._sp), sp)**( - floor( & + real(( minexponent(0._sp) - 1_sp ) / 2 ), sp) & + ) +! ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771 + real(sp), parameter :: ssbig = real(radix(0._sp), sp)**( - ceiling( & + real(( maxexponent(0._sp) - digits(0._sp) + 1_sp) / 2 ), sp) & + ) +! +! +! Standard constants for + integer, parameter :: dp = kind(1.d0) +! + real(dp), parameter :: dzero = 0.0_dp + real(dp), parameter :: dhalf = 0.5_dp + real(dp), parameter :: done = 1.0_dp + real(dp), parameter :: dtwo = 2.0_dp + real(dp), parameter :: dthree = 3.0_dp + real(dp), parameter :: dfour = 4.0_dp + real(dp), parameter :: deight = 8.0_dp + real(dp), parameter :: dten = 10.0_dp + complex(dp), parameter :: zzero = ( 0.0_dp, 0.0_dp ) + complex(dp), parameter :: zhalf = ( 0.5_dp, 0.0_dp ) + complex(dp), parameter :: zone = ( 1.0_dp, 0.0_dp ) + character*1, parameter :: dprefix = 'D' + character*1, parameter :: zprefix = 'Z' +! +! Scaling constants +! + real(dp), parameter :: dulp = epsilon(0._dp) + real(dp), parameter :: deps = dulp * 0.5_dp + real(dp), parameter :: dsafmin = real(radix(0._dp),dp)**max( & + minexponent(0._dp)-1, & + 1-maxexponent(0._dp) & + ) + real(dp), parameter :: dsafmax = done / dsafmin + real(dp), parameter :: dsmlnum = dsafmin / dulp + real(dp), parameter :: dbignum = dsafmax * dulp + real(dp), parameter :: drtmin = sqrt(dsmlnum) + real(dp), parameter :: drtmax = sqrt(dbignum) ! ! Blue's scaling constants ! - real(wp), parameter :: tsml = 0.14916681462400413487E-153_wp - real(wp), parameter :: tbig = 0.19979190722022350281E+147_wp - real(wp), parameter :: ssml = 0.44989137945431963828E+162_wp - real(wp), parameter :: sbig = 0.11113793747425387417E-161_wp + real(dp), parameter :: dtsml = real(radix(0._dp), dp)**ceiling( & + real(( minexponent(0._dp) - 1_sp ) / 2, dp) & + ) + real(dp), parameter :: dtbig = real(radix(0._dp), dp)**floor( & + real(( maxexponent(0._dp) - digits(0._dp) + 1_sp) / 2, dp) & + ) +! ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771 + real(dp), parameter :: dssml = real(radix(0._dp) ,dp)**( - floor( & + real(( minexponent(0._dp) - 1_sp ) / 2 ), dp) & + ) +! ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771 + real(dp), parameter :: dsbig = real(radix(0._dp), dp)**( - ceiling( & + real(( maxexponent(0._dp) - digits(0._dp) + 1_sp) / 2 ), dp) & + ) +! end module LA_CONSTANTS From 44da37564ece33c1865394e31c52b37caecd2db1 Mon Sep 17 00:00:00 2001 From: "weslley.spereira" Date: Mon, 1 Mar 2021 17:41:15 -0300 Subject: [PATCH 06/10] Adds .f90 extension to Doxyfile --- DOCS/Doxyfile | 1 + 1 file changed, 1 insertion(+) diff --git a/DOCS/Doxyfile b/DOCS/Doxyfile index f1d9c3d400..e6053c43bc 100644 --- a/DOCS/Doxyfile +++ b/DOCS/Doxyfile @@ -886,6 +886,7 @@ INPUT_ENCODING = UTF-8 FILE_PATTERNS = *.c \ *.f \ + *.f90 \ *.h # The RECURSIVE tag can be used to specify whether or not subdirectories should From 4746efca1ea568ae313daa4b6a75a98fc4da3c98 Mon Sep 17 00:00:00 2001 From: "Weslley S. Pereira" Date: Mon, 1 Mar 2021 17:52:45 -0300 Subject: [PATCH 07/10] Includes Nick Papior as a contributor --- SRC/la_constants.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/SRC/la_constants.f90 b/SRC/la_constants.f90 index 0bf798195d..07b37bc109 100644 --- a/SRC/la_constants.f90 +++ b/SRC/la_constants.f90 @@ -18,6 +18,7 @@ ! ================== !> !> Weslley Pereira, University of Colorado Denver, USA +!> Nick Papior, Technical University of Denmark, DNK ! !> \par Further Details: ! ===================== From ff8ab2694430dfa87b6d485cf62e8eeeb1b5ebfb Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Tue, 2 Mar 2021 09:47:39 +0100 Subject: [PATCH 08/10] mnt: cleanup and removal of erroneous kind specifications --- SRC/la_constants.f90 | 52 ++++++++++++++++---------------------------- 1 file changed, 19 insertions(+), 33 deletions(-) diff --git a/SRC/la_constants.f90 b/SRC/la_constants.f90 index 07b37bc109..0b9af3739d 100644 --- a/SRC/la_constants.f90 +++ b/SRC/la_constants.f90 @@ -18,7 +18,7 @@ ! ================== !> !> Weslley Pereira, University of Colorado Denver, USA -!> Nick Papior, Technical University of Denmark, DNK +!> Nick Papior, Technical University of Denmark, DK ! !> \par Further Details: ! ===================== @@ -38,15 +38,14 @@ !> \endverbatim ! module LA_CONSTANTS -! ! -- LAPACK auxiliary module (version 3.9.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! February 2021 -! + ! Standard constants for integer, parameter :: sp = kind(1.e0) -! + real(sp), parameter :: szero = 0.0_sp real(sp), parameter :: shalf = 0.5_sp real(sp), parameter :: sone = 1.0_sp @@ -60,9 +59,8 @@ module LA_CONSTANTS complex(sp), parameter :: cone = ( 1.0_sp, 0.0_sp ) character*1, parameter :: sprefix = 'S' character*1, parameter :: cprefix = 'C' -! + ! Scaling constants -! real(sp), parameter :: sulp = epsilon(0._sp) real(sp), parameter :: seps = sulp * 0.5_sp real(sp), parameter :: ssafmin = real(radix(0._sp),sp)**max( & @@ -74,28 +72,22 @@ module LA_CONSTANTS real(sp), parameter :: sbignum = ssafmax * sulp real(sp), parameter :: srtmin = sqrt(ssmlnum) real(sp), parameter :: srtmax = sqrt(sbignum) -! + ! Blue's scaling constants -! real(sp), parameter :: stsml = real(radix(0._sp), sp)**ceiling( & - real(( minexponent(0._sp) - 1_sp ) / 2, sp) & - ) + (minexponent(0._sp) - 1) * 0.5_sp) real(sp), parameter :: stbig = real(radix(0._sp), sp)**floor( & - real(( maxexponent(0._sp) - digits(0._sp) + 1_sp) / 2, sp) & - ) + (maxexponent(0._sp) - digits(0._sp) + 1) * 0.5_sp) ! ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771 real(sp), parameter :: sssml = real(radix(0._sp), sp)**( - floor( & - real(( minexponent(0._sp) - 1_sp ) / 2 ), sp) & - ) + (minexponent(0._sp) - 1) * 0.5_sp)) ! ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771 real(sp), parameter :: ssbig = real(radix(0._sp), sp)**( - ceiling( & - real(( maxexponent(0._sp) - digits(0._sp) + 1_sp) / 2 ), sp) & - ) -! -! + (maxexponent(0._sp) - digits(0._sp) + 1) * 0.5_sp)) + ! Standard constants for integer, parameter :: dp = kind(1.d0) -! + real(dp), parameter :: dzero = 0.0_dp real(dp), parameter :: dhalf = 0.5_dp real(dp), parameter :: done = 1.0_dp @@ -109,9 +101,8 @@ module LA_CONSTANTS complex(dp), parameter :: zone = ( 1.0_dp, 0.0_dp ) character*1, parameter :: dprefix = 'D' character*1, parameter :: zprefix = 'Z' -! + ! Scaling constants -! real(dp), parameter :: dulp = epsilon(0._dp) real(dp), parameter :: deps = dulp * 0.5_dp real(dp), parameter :: dsafmin = real(radix(0._dp),dp)**max( & @@ -123,22 +114,17 @@ module LA_CONSTANTS real(dp), parameter :: dbignum = dsafmax * dulp real(dp), parameter :: drtmin = sqrt(dsmlnum) real(dp), parameter :: drtmax = sqrt(dbignum) -! + ! Blue's scaling constants -! real(dp), parameter :: dtsml = real(radix(0._dp), dp)**ceiling( & - real(( minexponent(0._dp) - 1_sp ) / 2, dp) & - ) + (minexponent(0._dp) - 1) * 0.5_dp) real(dp), parameter :: dtbig = real(radix(0._dp), dp)**floor( & - real(( maxexponent(0._dp) - digits(0._dp) + 1_sp) / 2, dp) & - ) + (maxexponent(0._dp) - digits(0._dp) + 1) * 0.5_dp) ! ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771 - real(dp), parameter :: dssml = real(radix(0._dp) ,dp)**( - floor( & - real(( minexponent(0._dp) - 1_sp ) / 2 ), dp) & - ) + real(dp), parameter :: dssml = real(radix(0._dp), dp)**( - floor( & + (minexponent(0._dp) - 1) * 0.5_dp)) ! ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771 real(dp), parameter :: dsbig = real(radix(0._dp), dp)**( - ceiling( & - real(( maxexponent(0._dp) - digits(0._dp) + 1_sp) / 2 ), dp) & - ) -! + (maxexponent(0._dp) - digits(0._dp) + 1) * 0.5_dp)) + end module LA_CONSTANTS From 36060232053c8d80c6d9a63804035cab47c18309 Mon Sep 17 00:00:00 2001 From: "Weslley S. Pereira" Date: Tue, 2 Mar 2021 12:22:54 -0300 Subject: [PATCH 09/10] Update la_constants.f90 --- SRC/la_constants.f90 | 3 +-- SRC/la_constants32.f90 | 40 ---------------------------------------- 2 files changed, 1 insertion(+), 42 deletions(-) delete mode 100644 SRC/la_constants32.f90 diff --git a/SRC/la_constants.f90 b/SRC/la_constants.f90 index 0b9af3739d..3b1ab5e8c3 100644 --- a/SRC/la_constants.f90 +++ b/SRC/la_constants.f90 @@ -38,10 +38,9 @@ !> \endverbatim ! module LA_CONSTANTS -! -- LAPACK auxiliary module (version 3.9.0) -- +! -- LAPACK auxiliary module -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -! February 2021 ! Standard constants for integer, parameter :: sp = kind(1.e0) diff --git a/SRC/la_constants32.f90 b/SRC/la_constants32.f90 deleted file mode 100644 index d495709e92..0000000000 --- a/SRC/la_constants32.f90 +++ /dev/null @@ -1,40 +0,0 @@ -module LA_CONSTANTS32 -! -! -- BLAS/LAPACK module -- -! May 06, 2016 -! -! Standard constants -! - integer, parameter :: wp = 4 - real(wp), parameter :: zero = 0.0_wp - real(wp), parameter :: half = 0.5_wp - real(wp), parameter :: one = 1.0_wp - real(wp), parameter :: two = 2.0_wp - real(wp), parameter :: three = 3.0_wp - real(wp), parameter :: four = 4.0_wp - real(wp), parameter :: eight = 8.0_wp - real(wp), parameter :: ten = 10.0_wp - complex(wp), parameter :: czero = ( 0.0_wp, 0.0_wp ) - complex(wp), parameter :: chalf = ( 0.5_wp, 0.0_wp ) - complex(wp), parameter :: cone = ( 1.0_wp, 0.0_wp ) - character*1, parameter :: sprefix = 'S' - character*1, parameter :: cprefix = 'C' -! -! Model parameters -! - real(wp), parameter :: eps = 0.5960464478E-07_wp - real(wp), parameter :: ulp = 0.1192092896E-06_wp - real(wp), parameter :: safmin = 0.1175494351E-37_wp - real(wp), parameter :: safmax = 0.8507059173E+38_wp - real(wp), parameter :: smlnum = 0.9860761315E-31_wp - real(wp), parameter :: bignum = 0.1014120480E+32_wp - real(wp), parameter :: rtmin = 0.3140184864E-15_wp - real(wp), parameter :: rtmax = 0.3184525782E+16_wp -! -! Blue's scaling constants -! - real(wp), parameter :: tsml = 0.1084202172E-18_wp - real(wp), parameter :: tbig = 0.4503599627E+16_wp - real(wp), parameter :: ssml = 0.3777893186E+23_wp - real(wp), parameter :: sbig = 0.1323488980E-22_wp -end module LA_CONSTANTS32 From d71dd6e76c40db681c484bcf41da18dd76bdca7a Mon Sep 17 00:00:00 2001 From: "weslley.spereira" Date: Thu, 1 Apr 2021 09:48:24 -0300 Subject: [PATCH 10/10] Remove old lassq methods Update Doxygen headers Fix typos --- SRC/DEPRECATED/classq_lapackv390.f | 168 --------------------------- SRC/DEPRECATED/dlassq_lapackv390.f | 155 ------------------------- SRC/DEPRECATED/slassq_lapackv390.f | 155 ------------------------- SRC/DEPRECATED/zlassq_lapackv390.f | 168 --------------------------- SRC/Makefile | 4 +- SRC/classq.f90 | 180 +++++++++++++++++++++-------- SRC/dlassq.f90 | 180 +++++++++++++++++++++-------- SRC/la_xisnan.F90 | 8 +- SRC/slassq.f90 | 180 +++++++++++++++++++++-------- SRC/zlassq.f90 | 180 +++++++++++++++++++++-------- 10 files changed, 525 insertions(+), 853 deletions(-) delete mode 100644 SRC/DEPRECATED/classq_lapackv390.f delete mode 100644 SRC/DEPRECATED/dlassq_lapackv390.f delete mode 100644 SRC/DEPRECATED/slassq_lapackv390.f delete mode 100644 SRC/DEPRECATED/zlassq_lapackv390.f diff --git a/SRC/DEPRECATED/classq_lapackv390.f b/SRC/DEPRECATED/classq_lapackv390.f deleted file mode 100644 index 90e323b64c..0000000000 --- a/SRC/DEPRECATED/classq_lapackv390.f +++ /dev/null @@ -1,168 +0,0 @@ -*> \brief \b CLASSQ updates a sum of squares represented in scaled form. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download CLASSQ + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE CLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* .. Scalar Arguments .. -* INTEGER INCX, N -* REAL SCALE, SUMSQ -* .. -* .. Array Arguments .. -* COMPLEX X( * ) -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> CLASSQ returns the values scl and ssq such that -*> -*> ( scl**2 )*ssq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -*> -*> where x( i ) = abs( X( 1 + ( i - 1 )*INCX ) ). The value of sumsq is -*> assumed to be at least unity and the value of ssq will then satisfy -*> -*> 1.0 <= ssq <= ( sumsq + 2*n ). -*> -*> scale is assumed to be non-negative and scl returns the value -*> -*> scl = max( scale, abs( real( x( i ) ) ), abs( aimag( x( i ) ) ) ), -*> i -*> -*> scale and sumsq must be supplied in SCALE and SUMSQ respectively. -*> SCALE and SUMSQ are overwritten by scl and ssq respectively. -*> -*> The routine makes only one pass through the vector X. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] N -*> \verbatim -*> N is INTEGER -*> The number of elements to be used from the vector X. -*> \endverbatim -*> -*> \param[in] X -*> \verbatim -*> X is COMPLEX array, dimension (1+(N-1)*INCX) -*> The vector x as described above. -*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. -*> \endverbatim -*> -*> \param[in] INCX -*> \verbatim -*> INCX is INTEGER -*> The increment between successive values of the vector X. -*> INCX > 0. -*> \endverbatim -*> -*> \param[in,out] SCALE -*> \verbatim -*> SCALE is REAL -*> On entry, the value scale in the equation above. -*> On exit, SCALE is overwritten with the value scl . -*> \endverbatim -*> -*> \param[in,out] SUMSQ -*> \verbatim -*> SUMSQ is REAL -*> On entry, the value sumsq in the equation above. -*> On exit, SUMSQ is overwritten with the value ssq . -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> This implementation of CLASSQ has been deprecated with LAPACKv3.10. -*> A better version of CLASSQ was contributed by Ed Anderson and released in 3.10. -* -*> \ingroup complexOTHERauxiliary -* -* ===================================================================== - SUBROUTINE CLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* -- LAPACK auxiliary routine -- -* -- LAPACK is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX, N - REAL SCALE, SUMSQ -* .. -* .. Array Arguments .. - COMPLEX X( * ) -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER ( ZERO = 0.0E+0 ) -* .. -* .. Local Scalars .. - INTEGER IX - REAL TEMP1 -* .. -* .. External Functions .. - LOGICAL SISNAN - EXTERNAL SISNAN -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, AIMAG, REAL -* .. -* .. Executable Statements .. -* - IF( N.GT.0 ) THEN - DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX - TEMP1 = ABS( REAL( X( IX ) ) ) - IF( TEMP1.GT.ZERO.OR.SISNAN( TEMP1 ) ) THEN - IF( SCALE.LT.TEMP1 ) THEN - SUMSQ = 1 + SUMSQ*( SCALE / TEMP1 )**2 - SCALE = TEMP1 - ELSE - SUMSQ = SUMSQ + ( TEMP1 / SCALE )**2 - END IF - END IF - TEMP1 = ABS( AIMAG( X( IX ) ) ) - IF( TEMP1.GT.ZERO.OR.SISNAN( TEMP1 ) ) THEN - IF( SCALE.LT.TEMP1 .OR. SISNAN( TEMP1 ) ) THEN - SUMSQ = 1 + SUMSQ*( SCALE / TEMP1 )**2 - SCALE = TEMP1 - ELSE - SUMSQ = SUMSQ + ( TEMP1 / SCALE )**2 - END IF - END IF - 10 CONTINUE - END IF -* - RETURN -* -* End of CLASSQ -* - END diff --git a/SRC/DEPRECATED/dlassq_lapackv390.f b/SRC/DEPRECATED/dlassq_lapackv390.f deleted file mode 100644 index d3126951ca..0000000000 --- a/SRC/DEPRECATED/dlassq_lapackv390.f +++ /dev/null @@ -1,155 +0,0 @@ -*> \brief \b DLASSQ updates a sum of squares represented in scaled form. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download DLASSQ + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* .. Scalar Arguments .. -* INTEGER INCX, N -* DOUBLE PRECISION SCALE, SUMSQ -* .. -* .. Array Arguments .. -* DOUBLE PRECISION X( * ) -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> DLASSQ returns the values scl and smsq such that -*> -*> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -*> -*> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is -*> assumed to be non-negative and scl returns the value -*> -*> scl = max( scale, abs( x( i ) ) ). -*> -*> scale and sumsq must be supplied in SCALE and SUMSQ and -*> scl and smsq are overwritten on SCALE and SUMSQ respectively. -*> -*> The routine makes only one pass through the vector x. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] N -*> \verbatim -*> N is INTEGER -*> The number of elements to be used from the vector X. -*> \endverbatim -*> -*> \param[in] X -*> \verbatim -*> X is DOUBLE PRECISION array, dimension (1+(N-1)*INCX) -*> The vector for which a scaled sum of squares is computed. -*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. -*> \endverbatim -*> -*> \param[in] INCX -*> \verbatim -*> INCX is INTEGER -*> The increment between successive values of the vector X. -*> INCX > 0. -*> \endverbatim -*> -*> \param[in,out] SCALE -*> \verbatim -*> SCALE is DOUBLE PRECISION -*> On entry, the value scale in the equation above. -*> On exit, SCALE is overwritten with scl , the scaling factor -*> for the sum of squares. -*> \endverbatim -*> -*> \param[in,out] SUMSQ -*> \verbatim -*> SUMSQ is DOUBLE PRECISION -*> On entry, the value sumsq in the equation above. -*> On exit, SUMSQ is overwritten with smsq , the basic sum of -*> squares from which scl has been factored out. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> This implementation of DLASSQ has been deprecated with LAPACKv3.10. -*> A better version of DLASSQ was contributed by Ed Anderson and released in 3.10. -* -*> \ingroup OTHERauxiliary -* -* ===================================================================== - SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* -- LAPACK auxiliary routine -- -* -- LAPACK is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX, N - DOUBLE PRECISION SCALE, SUMSQ -* .. -* .. Array Arguments .. - DOUBLE PRECISION X( * ) -* .. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D+0 ) -* .. -* .. Local Scalars .. - INTEGER IX - DOUBLE PRECISION ABSXI -* .. -* .. External Functions .. - LOGICAL DISNAN - EXTERNAL DISNAN -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS -* .. -* .. Executable Statements .. -* - IF( N.GT.0 ) THEN - DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX - ABSXI = ABS( X( IX ) ) - IF( ABSXI.GT.ZERO.OR.DISNAN( ABSXI ) ) THEN - IF( SCALE.LT.ABSXI ) THEN - SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2 - SCALE = ABSXI - ELSE - SUMSQ = SUMSQ + ( ABSXI / SCALE )**2 - END IF - END IF - 10 CONTINUE - END IF - RETURN -* -* End of DLASSQ -* - END diff --git a/SRC/DEPRECATED/slassq_lapackv390.f b/SRC/DEPRECATED/slassq_lapackv390.f deleted file mode 100644 index 73f63b5150..0000000000 --- a/SRC/DEPRECATED/slassq_lapackv390.f +++ /dev/null @@ -1,155 +0,0 @@ -*> \brief \b SLASSQ updates a sum of squares represented in scaled form. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download SLASSQ + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE SLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* .. Scalar Arguments .. -* INTEGER INCX, N -* REAL SCALE, SUMSQ -* .. -* .. Array Arguments .. -* REAL X( * ) -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> SLASSQ returns the values scl and smsq such that -*> -*> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -*> -*> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is -*> assumed to be non-negative and scl returns the value -*> -*> scl = max( scale, abs( x( i ) ) ). -*> -*> scale and sumsq must be supplied in SCALE and SUMSQ and -*> scl and smsq are overwritten on SCALE and SUMSQ respectively. -*> -*> The routine makes only one pass through the vector x. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] N -*> \verbatim -*> N is INTEGER -*> The number of elements to be used from the vector X. -*> \endverbatim -*> -*> \param[in] X -*> \verbatim -*> X is REAL array, dimension (1+(N-1)*INCX) -*> The vector for which a scaled sum of squares is computed. -*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. -*> \endverbatim -*> -*> \param[in] INCX -*> \verbatim -*> INCX is INTEGER -*> The increment between successive values of the vector X. -*> INCX > 0. -*> \endverbatim -*> -*> \param[in,out] SCALE -*> \verbatim -*> SCALE is REAL -*> On entry, the value scale in the equation above. -*> On exit, SCALE is overwritten with scl , the scaling factor -*> for the sum of squares. -*> \endverbatim -*> -*> \param[in,out] SUMSQ -*> \verbatim -*> SUMSQ is REAL -*> On entry, the value sumsq in the equation above. -*> On exit, SUMSQ is overwritten with smsq , the basic sum of -*> squares from which scl has been factored out. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> This implementation of SLASSQ has been deprecated with LAPACKv3.10. -*> A better version of SLASSQ was contributed by Ed Anderson and released in 3.10. -* -*> \ingroup OTHERauxiliary -* -* ===================================================================== - SUBROUTINE SLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* -- LAPACK auxiliary routine -- -* -- LAPACK is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX, N - REAL SCALE, SUMSQ -* .. -* .. Array Arguments .. - REAL X( * ) -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER ( ZERO = 0.0E+0 ) -* .. -* .. Local Scalars .. - INTEGER IX - REAL ABSXI -* .. -* .. External Functions .. - LOGICAL SISNAN - EXTERNAL SISNAN -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS -* .. -* .. Executable Statements .. -* - IF( N.GT.0 ) THEN - DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX - ABSXI = ABS( X( IX ) ) - IF( ABSXI.GT.ZERO.OR.SISNAN( ABSXI ) ) THEN - IF( SCALE.LT.ABSXI ) THEN - SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2 - SCALE = ABSXI - ELSE - SUMSQ = SUMSQ + ( ABSXI / SCALE )**2 - END IF - END IF - 10 CONTINUE - END IF - RETURN -* -* End of SLASSQ -* - END diff --git a/SRC/DEPRECATED/zlassq_lapackv390.f b/SRC/DEPRECATED/zlassq_lapackv390.f deleted file mode 100644 index 0133f2456f..0000000000 --- a/SRC/DEPRECATED/zlassq_lapackv390.f +++ /dev/null @@ -1,168 +0,0 @@ -*> \brief \b ZLASSQ updates a sum of squares represented in scaled form. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download ZLASSQ + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE ZLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* .. Scalar Arguments .. -* INTEGER INCX, N -* DOUBLE PRECISION SCALE, SUMSQ -* .. -* .. Array Arguments .. -* COMPLEX*16 X( * ) -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> ZLASSQ returns the values scl and ssq such that -*> -*> ( scl**2 )*ssq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -*> -*> where x( i ) = abs( X( 1 + ( i - 1 )*INCX ) ). The value of sumsq is -*> assumed to be at least unity and the value of ssq will then satisfy -*> -*> 1.0 <= ssq <= ( sumsq + 2*n ). -*> -*> scale is assumed to be non-negative and scl returns the value -*> -*> scl = max( scale, abs( real( x( i ) ) ), abs( aimag( x( i ) ) ) ), -*> i -*> -*> scale and sumsq must be supplied in SCALE and SUMSQ respectively. -*> SCALE and SUMSQ are overwritten by scl and ssq respectively. -*> -*> The routine makes only one pass through the vector X. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] N -*> \verbatim -*> N is INTEGER -*> The number of elements to be used from the vector X. -*> \endverbatim -*> -*> \param[in] X -*> \verbatim -*> X is COMPLEX*16 array, dimension (1+(N-1)*INCX) -*> The vector x as described above. -*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. -*> \endverbatim -*> -*> \param[in] INCX -*> \verbatim -*> INCX is INTEGER -*> The increment between successive values of the vector X. -*> INCX > 0. -*> \endverbatim -*> -*> \param[in,out] SCALE -*> \verbatim -*> SCALE is DOUBLE PRECISION -*> On entry, the value scale in the equation above. -*> On exit, SCALE is overwritten with the value scl . -*> \endverbatim -*> -*> \param[in,out] SUMSQ -*> \verbatim -*> SUMSQ is DOUBLE PRECISION -*> On entry, the value sumsq in the equation above. -*> On exit, SUMSQ is overwritten with the value ssq . -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> This implementation of ZLASSQ has been deprecated with LAPACKv3.10. -*> A better version of ZLASSQ was contributed by Ed Anderson and released in 3.10. -* -*> \ingroup complex16OTHERauxiliary -* -* ===================================================================== - SUBROUTINE ZLASSQ( N, X, INCX, SCALE, SUMSQ ) -* -* -- LAPACK auxiliary routine -- -* -- LAPACK is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - INTEGER INCX, N - DOUBLE PRECISION SCALE, SUMSQ -* .. -* .. Array Arguments .. - COMPLEX*16 X( * ) -* .. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D+0 ) -* .. -* .. Local Scalars .. - INTEGER IX - DOUBLE PRECISION TEMP1 -* .. -* .. External Functions .. - LOGICAL DISNAN - EXTERNAL DISNAN -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, DIMAG -* .. -* .. Executable Statements .. -* - IF( N.GT.0 ) THEN - DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX - TEMP1 = ABS( DBLE( X( IX ) ) ) - IF( TEMP1.GT.ZERO.OR.DISNAN( TEMP1 ) ) THEN - IF( SCALE.LT.TEMP1 ) THEN - SUMSQ = 1 + SUMSQ*( SCALE / TEMP1 )**2 - SCALE = TEMP1 - ELSE - SUMSQ = SUMSQ + ( TEMP1 / SCALE )**2 - END IF - END IF - TEMP1 = ABS( DIMAG( X( IX ) ) ) - IF( TEMP1.GT.ZERO.OR.DISNAN( TEMP1 ) ) THEN - IF( SCALE.LT.TEMP1 ) THEN - SUMSQ = 1 + SUMSQ*( SCALE / TEMP1 )**2 - SCALE = TEMP1 - ELSE - SUMSQ = SUMSQ + ( TEMP1 / SCALE )**2 - END IF - END IF - 10 CONTINUE - END IF -* - RETURN -* -* End of ZLASSQ -* - END diff --git a/SRC/Makefile b/SRC/Makefile index 20a2da5cac..a54420afc7 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -59,9 +59,7 @@ include $(TOPSRCDIR)/make.inc ALLMOD = la_xisnan.mod la_constants.mod -.SUFFIXES: .f .F .f90 .F90 .o .mod -%.o: %.f $(ALLMOD) - $(FC) $(FFLAGS) -c -o $@ $< +.SUFFIXES: .F .f90 .F90 .o .mod %.o: %.F $(ALLMOD) $(FC) $(FFLAGS) -c -o $@ $< %.o: %.f90 $(ALLMOD) diff --git a/SRC/classq.f90 b/SRC/classq.f90 index 75def12baf..331f6179ae 100644 --- a/SRC/classq.f90 +++ b/SRC/classq.f90 @@ -1,11 +1,136 @@ +!> \brief \b CLASSQ updates a sum of squares represented in scaled form. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +!> \htmlonly +!> Download CLASSQ + dependencies +!> +!> [TGZ] +!> +!> [ZIP] +!> +!> [TXT] +!> \endhtmlonly +! +! Definition: +! =========== +! +! SUBROUTINE CLASSQ( N, X, INCX, SCALE, SUMSQ ) +! +! .. Scalar Arguments .. +! INTEGER INCX, N +! REAL SCALE, SUMSQ +! .. +! .. Array Arguments .. +! COMPLEX X( * ) +! .. +! +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> CLASSQ returns the values scl and smsq such that +!> +!> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +!> +!> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +!> assumed to be non-negative and scl returns the value +!> +!> scl = max( scale, abs( x( i ) ) ). +!> +!> scale and sumsq must be supplied in SCALE and SUMSQ and +!> scl and smsq are overwritten on SCALE and SUMSQ respectively. +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] N +!> \verbatim +!> N is INTEGER +!> The number of elements to be used from the vector x. +!> \endverbatim +!> +!> \param[in] X +!> \verbatim +!> X is COMPLEX array, dimension (1+(N-1)*abs(INCX)) +!> The vector for which a scaled sum of squares is computed. +!> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. +!> \endverbatim +!> +!> \param[in] INCX +!> \verbatim +!> INCX is INTEGER +!> The increment between successive values of the vector 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 +!> +!> \param[in,out] SCALE +!> \verbatim +!> SCALE is REAL +!> On entry, the value scale in the equation above. +!> On exit, SCALE is overwritten with scl , the scaling factor +!> for the sum of squares. +!> \endverbatim +!> +!> \param[in,out] SUMSQ +!> \verbatim +!> SUMSQ is REAL +!> On entry, the value sumsq in the equation above. +!> On exit, SUMSQ is overwritten with smsq , the basic sum of +!> squares from which scl has been factored out. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +!> Nick Papior, Technical University of Denmark, DK +! +!> \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 +! +!> \ingroup OTHERauxiliary +! +! ===================================================================== subroutine CLASSQ( n, x, incx, scl, sumsq ) - use LA_CONSTANTS32, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_CONSTANTS, & + only: wp=>sp, zero=>szero, one=>sone, & + sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml use LA_XISNAN ! -! LAPACK auxiliary routine -! Based on Blue's algorithm from ACM TOMS, March 1978 -! E. Anderson -! August 9, 2016 +! -- LAPACK auxiliary routine -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! ! .. Scalar Arguments .. integer :: incx, n @@ -14,51 +139,6 @@ subroutine CLASSQ( n, x, incx, scl, sumsq ) ! .. Array Arguments .. complex(wp) :: x(*) ! .. -! -! Purpose -! ======= -! -! CLASSQ returns the values scl and smsq such that -! -! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -! -! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is -! assumed to be non-negative and scl returns the value -! -! scl = max( scale, abs( x( i ) ) ). -! -! scale and sumsq must be supplied in SCL and SUMSQ and -! scl and smsq are overwritten on SCL and SUMSQ respectively. -! -! Arguments -! ========= -! -! N (input) INTEGER -! The number of elements of the vector x. -! -! X (input) COMPLEX array, dimension (1+(N-1)*abs(INCX)) -! The n-element vector x. -! -! INCX (input) INTEGER -! The increment between successive values of the vector 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. -! -! SCL (input/output) REAL -! On entry, the value scale in the equation above. -! On exit, SCL is overwritten with scl , the scaling factor -! for the sum of squares. -! -! SUMSQ (input/output) REAL -! On entry, the value sumsq in the equation above. -! On exit, SUMSQ is overwritten with smsq , the basic sum of -! squares from which scl has been factored out. -! -! ===================================================================== -! ! .. Local Scalars .. integer :: i, ix logical :: notbig diff --git a/SRC/dlassq.f90 b/SRC/dlassq.f90 index a9621ea01b..47265e199d 100644 --- a/SRC/dlassq.f90 +++ b/SRC/dlassq.f90 @@ -1,11 +1,136 @@ +!> \brief \b DLASSQ updates a sum of squares represented in scaled form. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +!> \htmlonly +!> Download DLASSQ + dependencies +!> +!> [TGZ] +!> +!> [ZIP] +!> +!> [TXT] +!> \endhtmlonly +! +! Definition: +! =========== +! +! SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) +! +! .. Scalar Arguments .. +! INTEGER INCX, N +! DOUBLE PRECISION SCALE, SUMSQ +! .. +! .. Array Arguments .. +! DOUBLE PRECISION X( * ) +! .. +! +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> DLASSQ returns the values scl and smsq such that +!> +!> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +!> +!> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +!> assumed to be non-negative and scl returns the value +!> +!> scl = max( scale, abs( x( i ) ) ). +!> +!> scale and sumsq must be supplied in SCALE and SUMSQ and +!> scl and smsq are overwritten on SCALE and SUMSQ respectively. +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] N +!> \verbatim +!> N is INTEGER +!> The number of elements to be used from the vector x. +!> \endverbatim +!> +!> \param[in] X +!> \verbatim +!> X is DOUBLE PRECISION array, dimension (1+(N-1)*abs(INCX)) +!> The vector for which a scaled sum of squares is computed. +!> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. +!> \endverbatim +!> +!> \param[in] INCX +!> \verbatim +!> INCX is INTEGER +!> The increment between successive values of the vector 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 +!> +!> \param[in,out] SCALE +!> \verbatim +!> SCALE is DOUBLE PRECISION +!> On entry, the value scale in the equation above. +!> On exit, SCALE is overwritten with scl , the scaling factor +!> for the sum of squares. +!> \endverbatim +!> +!> \param[in,out] SUMSQ +!> \verbatim +!> SUMSQ is DOUBLE PRECISION +!> On entry, the value sumsq in the equation above. +!> On exit, SUMSQ is overwritten with smsq , the basic sum of +!> squares from which scl has been factored out. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +!> Nick Papior, Technical University of Denmark, DK +! +!> \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 +! +!> \ingroup OTHERauxiliary +! +! ===================================================================== subroutine DLASSQ( n, x, incx, scl, sumsq ) - use LA_CONSTANTS, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_CONSTANTS, & + only: wp=>dp, zero=>dzero, one=>done, & + sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml use LA_XISNAN ! -! LAPACK auxiliary routine -! Based on Blue's algorithm from ACM TOMS, March 1978 -! E. Anderson -! August 9, 2016 +! -- LAPACK auxiliary routine -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! ! .. Scalar Arguments .. integer :: incx, n @@ -14,51 +139,6 @@ subroutine DLASSQ( n, x, incx, scl, sumsq ) ! .. Array Arguments .. real(wp) :: x(*) ! .. -! -! Purpose -! ======= -! -! DLASSQ returns the values scl and smsq such that -! -! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -! -! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is -! assumed to be non-negative and scl returns the value -! -! scl = max( scale, abs( x( i ) ) ). -! -! scale and sumsq must be supplied in SCL and SUMSQ and -! scl and smsq are overwritten on SCL and SUMSQ respectively. -! -! Arguments -! ========= -! -! N (input) INTEGER -! The number of elements of the vector x. -! -! X (input) REAL array, dimension (1+(N-1)*abs(INCX)) -! The n-element vector x. -! -! INCX (input) INTEGER -! The increment between successive values of the vector 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. -! -! SCL (input/output) REAL -! On entry, the value scale in the equation above. -! On exit, SCL is overwritten with scl , the scaling factor -! for the sum of squares. -! -! SUMSQ (input/output) REAL -! On entry, the value sumsq in the equation above. -! On exit, SUMSQ is overwritten with smsq , the basic sum of -! squares from which scl has been factored out. -! -! ===================================================================== -! ! .. Local Scalars .. integer :: i, ix logical :: notbig diff --git a/SRC/la_xisnan.F90 b/SRC/la_xisnan.F90 index 6ba91de0bd..50966a5c18 100644 --- a/SRC/la_xisnan.F90 +++ b/SRC/la_xisnan.F90 @@ -9,7 +9,7 @@ module LA_XISNAN contains logical function SISNAN( x ) - use LA_CONSTANTS32, only: wp + use LA_CONSTANTS, only: wp=>sp #ifdef USE_IEEE_INTRINSIC use, intrinsic :: ieee_arithmetic #elif USE_ISNAN @@ -25,7 +25,7 @@ logical function SISNAN( x ) contains logical function SLAISNAN( x, y ) - use LA_CONSTANTS32, only: wp + use LA_CONSTANTS, only: wp=>sp real(wp) :: x, y SLAISNAN = ( x.ne.y ) end function SLAISNAN @@ -33,7 +33,7 @@ end function SLAISNAN end function SISNAN logical function DISNAN( x ) - use LA_CONSTANTS, only: wp + use LA_CONSTANTS, only: wp=>dp #ifdef USE_IEEE_INTRINSIC use, intrinsic :: ieee_arithmetic #elif USE_ISNAN @@ -49,7 +49,7 @@ logical function DISNAN( x ) contains logical function DLAISNAN( x, y ) - use LA_CONSTANTS, only: wp + use LA_CONSTANTS, only: wp=>dp real(wp) :: x, y DLAISNAN = ( x.ne.y ) end function DLAISNAN diff --git a/SRC/slassq.f90 b/SRC/slassq.f90 index f8a7f6b4a4..c1231eb02b 100644 --- a/SRC/slassq.f90 +++ b/SRC/slassq.f90 @@ -1,11 +1,136 @@ +!> \brief \b SLASSQ updates a sum of squares represented in scaled form. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +!> \htmlonly +!> Download SLASSQ + dependencies +!> +!> [TGZ] +!> +!> [ZIP] +!> +!> [TXT] +!> \endhtmlonly +! +! Definition: +! =========== +! +! SUBROUTINE SLASSQ( N, X, INCX, SCALE, SUMSQ ) +! +! .. Scalar Arguments .. +! INTEGER INCX, N +! REAL SCALE, SUMSQ +! .. +! .. Array Arguments .. +! REAL X( * ) +! .. +! +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> SLASSQ returns the values scl and smsq such that +!> +!> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +!> +!> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +!> assumed to be non-negative and scl returns the value +!> +!> scl = max( scale, abs( x( i ) ) ). +!> +!> scale and sumsq must be supplied in SCALE and SUMSQ and +!> scl and smsq are overwritten on SCALE and SUMSQ respectively. +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] N +!> \verbatim +!> N is INTEGER +!> The number of elements to be used from the vector x. +!> \endverbatim +!> +!> \param[in] X +!> \verbatim +!> X is REAL array, dimension (1+(N-1)*abs(INCX)) +!> The vector for which a scaled sum of squares is computed. +!> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. +!> \endverbatim +!> +!> \param[in] INCX +!> \verbatim +!> INCX is INTEGER +!> The increment between successive values of the vector 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 +!> +!> \param[in,out] SCALE +!> \verbatim +!> SCALE is REAL +!> On entry, the value scale in the equation above. +!> On exit, SCALE is overwritten with scl , the scaling factor +!> for the sum of squares. +!> \endverbatim +!> +!> \param[in,out] SUMSQ +!> \verbatim +!> SUMSQ is REAL +!> On entry, the value sumsq in the equation above. +!> On exit, SUMSQ is overwritten with smsq , the basic sum of +!> squares from which scl has been factored out. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +!> Nick Papior, Technical University of Denmark, DK +! +!> \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 +! +!> \ingroup OTHERauxiliary +! +! ===================================================================== subroutine SLASSQ( n, x, incx, scl, sumsq ) - use LA_CONSTANTS32, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_CONSTANTS, & + only: wp=>sp, zero=>szero, one=>sone, & + sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml use LA_XISNAN ! -! LAPACK auxiliary routine -! Based on Blue's algorithm from ACM TOMS, March 1978 -! E. Anderson -! August 9, 2016 +! -- LAPACK auxiliary routine -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! ! .. Scalar Arguments .. integer :: incx, n @@ -14,51 +139,6 @@ subroutine SLASSQ( n, x, incx, scl, sumsq ) ! .. Array Arguments .. real(wp) :: x(*) ! .. -! -! Purpose -! ======= -! -! SLASSQ returns the values scl and smsq such that -! -! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -! -! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is -! assumed to be non-negative and scl returns the value -! -! scl = max( scale, abs( x( i ) ) ). -! -! scale and sumsq must be supplied in SCL and SUMSQ and -! scl and smsq are overwritten on SCL and SUMSQ respectively. -! -! Arguments -! ========= -! -! N (input) INTEGER -! The number of elements of the vector x. -! -! X (input) REAL array, dimension (1+(N-1)*abs(INCX)) -! The n-element vector x. -! -! INCX (input) INTEGER -! The increment between successive values of the vector 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. -! -! SCL (input/output) REAL -! On entry, the value scale in the equation above. -! On exit, SCL is overwritten with scl , the scaling factor -! for the sum of squares. -! -! SUMSQ (input/output) REAL -! On entry, the value sumsq in the equation above. -! On exit, SUMSQ is overwritten with smsq , the basic sum of -! squares from which scl has been factored out. -! -! ===================================================================== -! ! .. Local Scalars .. integer :: i, ix logical :: notbig diff --git a/SRC/zlassq.f90 b/SRC/zlassq.f90 index 002e198eda..3b41c1b2b5 100644 --- a/SRC/zlassq.f90 +++ b/SRC/zlassq.f90 @@ -1,11 +1,136 @@ +!> \brief \b ZLASSQ updates a sum of squares represented in scaled form. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +!> \htmlonly +!> Download ZLASSQ + dependencies +!> +!> [TGZ] +!> +!> [ZIP] +!> +!> [TXT] +!> \endhtmlonly +! +! Definition: +! =========== +! +! SUBROUTINE ZLASSQ( N, X, INCX, SCALE, SUMSQ ) +! +! .. Scalar Arguments .. +! INTEGER INCX, N +! DOUBLE PRECISION SCALE, SUMSQ +! .. +! .. Array Arguments .. +! DOUBLE COMPLEX X( * ) +! .. +! +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> ZLASSQ returns the values scl and smsq such that +!> +!> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, +!> +!> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is +!> assumed to be non-negative and scl returns the value +!> +!> scl = max( scale, abs( x( i ) ) ). +!> +!> scale and sumsq must be supplied in SCALE and SUMSQ and +!> scl and smsq are overwritten on SCALE and SUMSQ respectively. +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] N +!> \verbatim +!> N is INTEGER +!> The number of elements to be used from the vector x. +!> \endverbatim +!> +!> \param[in] X +!> \verbatim +!> X is DOUBLE COMPLEX array, dimension (1+(N-1)*abs(INCX)) +!> The vector for which a scaled sum of squares is computed. +!> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. +!> \endverbatim +!> +!> \param[in] INCX +!> \verbatim +!> INCX is INTEGER +!> The increment between successive values of the vector 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 +!> +!> \param[in,out] SCALE +!> \verbatim +!> SCALE is DOUBLE PRECISION +!> On entry, the value scale in the equation above. +!> On exit, SCALE is overwritten with scl , the scaling factor +!> for the sum of squares. +!> \endverbatim +!> +!> \param[in,out] SUMSQ +!> \verbatim +!> SUMSQ is DOUBLE PRECISION +!> On entry, the value sumsq in the equation above. +!> On exit, SUMSQ is overwritten with smsq , the basic sum of +!> squares from which scl has been factored out. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +!> Nick Papior, Technical University of Denmark, DK +! +!> \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 +! +!> \ingroup OTHERauxiliary +! +! ===================================================================== subroutine ZLASSQ( n, x, incx, scl, sumsq ) - use LA_CONSTANTS, only: wp, zero, one, sbig, ssml, tbig, tsml + use LA_CONSTANTS, & + only: wp=>dp, zero=>dzero, one=>done, & + sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml use LA_XISNAN ! -! LAPACK auxiliary routine -! Based on Blue's algorithm from ACM TOMS, March 1978 -! E. Anderson -! August 9, 2016 +! -- LAPACK auxiliary routine -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! ! .. Scalar Arguments .. integer :: incx, n @@ -14,51 +139,6 @@ subroutine ZLASSQ( n, x, incx, scl, sumsq ) ! .. Array Arguments .. complex(wp) :: x(*) ! .. -! -! Purpose -! ======= -! -! ZLASSQ returns the values scl and smsq such that -! -! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, -! -! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is -! assumed to be non-negative and scl returns the value -! -! scl = max( scale, abs( x( i ) ) ). -! -! scale and sumsq must be supplied in SCL and SUMSQ and -! scl and smsq are overwritten on SCL and SUMSQ respectively. -! -! Arguments -! ========= -! -! N (input) INTEGER -! The number of elements of the vector x. -! -! X (input) COMPLEX array, dimension (1+(N-1)*abs(INCX)) -! The n-element vector x. -! -! INCX (input) INTEGER -! The increment between successive values of the vector 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. -! -! SCL (input/output) REAL -! On entry, the value scale in the equation above. -! On exit, SCL is overwritten with scl , the scaling factor -! for the sum of squares. -! -! SUMSQ (input/output) REAL -! On entry, the value sumsq in the equation above. -! On exit, SUMSQ is overwritten with smsq , the basic sum of -! squares from which scl has been factored out. -! -! ===================================================================== -! ! .. Local Scalars .. integer :: i, ix logical :: notbig