diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index 62e6ec72f1..070eeb8100 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 + 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) @@ -54,7 +54,7 @@ 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}) @@ -73,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}) @@ -211,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 @@ -407,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/Makefile b/SRC/Makefile index 66b78f40c8..d15a8e628b 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -57,11 +57,9 @@ TOPSRCDIR = .. include $(TOPSRCDIR)/make.inc -ALLMOD = la_constants.mod +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) @@ -72,7 +70,7 @@ ALLMOD = la_constants.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 @@ -632,5 +630,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_constants.mod + $(FC) $(FFLAGS) -c -o $@ $< la_constants.o: la_constants.f90 $(FC) $(FFLAGS) -c -o $@ $< diff --git a/SRC/classq.f b/SRC/classq.f deleted file mode 100644 index a1c13309bc..0000000000 --- a/SRC/classq.f +++ /dev/null @@ -1,165 +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. -* -*> \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/classq.f90 b/SRC/classq.f90 new file mode 100644 index 0000000000..331f6179ae --- /dev/null +++ b/SRC/classq.f90 @@ -0,0 +1,250 @@ +!> \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_CONSTANTS, & + only: wp=>sp, zero=>szero, one=>sone, & + sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml + use LA_XISNAN +! +! -- 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(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + complex(wp) :: x(*) +! .. +! .. 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.f b/SRC/dlassq.f deleted file mode 100644 index 5c24d00350..0000000000 --- a/SRC/dlassq.f +++ /dev/null @@ -1,152 +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. -* -*> \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/dlassq.f90 b/SRC/dlassq.f90 new file mode 100644 index 0000000000..47265e199d --- /dev/null +++ b/SRC/dlassq.f90 @@ -0,0 +1,241 @@ +!> \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=>dp, zero=>dzero, one=>done, & + sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml + use LA_XISNAN +! +! -- 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(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + real(wp) :: x(*) +! .. +! .. 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 index d8146cb1bb..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.10.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_xisnan.F90 b/SRC/la_xisnan.F90 new file mode 100644 index 0000000000..50966a5c18 --- /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_CONSTANTS, only: wp=>sp +#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_CONSTANTS, only: wp=>sp + real(wp) :: x, y + SLAISNAN = ( x.ne.y ) + end function SLAISNAN +#endif + end function SISNAN + + logical function DISNAN( x ) + use LA_CONSTANTS, only: wp=>dp +#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=>dp + real(wp) :: x, y + DLAISNAN = ( x.ne.y ) + end function DLAISNAN +#endif + end function DISNAN + +end module LA_XISNAN diff --git a/SRC/slassq.f b/SRC/slassq.f deleted file mode 100644 index 72503cb77a..0000000000 --- a/SRC/slassq.f +++ /dev/null @@ -1,152 +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. -* -*> \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/slassq.f90 b/SRC/slassq.f90 new file mode 100644 index 0000000000..c1231eb02b --- /dev/null +++ b/SRC/slassq.f90 @@ -0,0 +1,241 @@ +!> \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_CONSTANTS, & + only: wp=>sp, zero=>szero, one=>sone, & + sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml + use LA_XISNAN +! +! -- 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(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + real(wp) :: x(*) +! .. +! .. 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.f b/SRC/zlassq.f deleted file mode 100644 index 02b0d7a14a..0000000000 --- a/SRC/zlassq.f +++ /dev/null @@ -1,165 +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. -* -*> \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/zlassq.f90 b/SRC/zlassq.f90 new file mode 100644 index 0000000000..3b41c1b2b5 --- /dev/null +++ b/SRC/zlassq.f90 @@ -0,0 +1,250 @@ +!> \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=>dp, zero=>dzero, one=>done, & + sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml + use LA_XISNAN +! +! -- 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(wp) :: scl, sumsq +! .. +! .. Array Arguments .. + complex(wp) :: x(*) +! .. +! .. 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