From 831ac96880aecf0dea8ecc8c629a212295f9dc8e Mon Sep 17 00:00:00 2001 From: Weslley S Pereira Date: Wed, 24 May 2023 16:59:25 -0600 Subject: [PATCH 1/5] Adds CRSCL --- SRC/CMakeLists.txt | 2 +- SRC/cgetf2.f | 1 + SRC/crscl.f | 208 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 210 insertions(+), 1 deletion(-) create mode 100644 SRC/crscl.f diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index 5d2e072584..1556b31dc7 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -227,7 +227,7 @@ set(CLASRC cppcon.f cppequ.f cpprfs.f cppsv.f cppsvx.f cpptrf.f cpptri.f cpptrs.f cptcon.f cpteqr.f cptrfs.f cptsv.f cptsvx.f cpttrf.f cpttrs.f cptts2.f crot.f cspcon.f cspmv.f cspr.f csprfs.f cspsv.f - cspsvx.f csptrf.f csptri.f csptrs.f csrscl.f cstedc.f + cspsvx.f csptrf.f csptri.f csptrs.f csrscl.f crscl.f cstedc.f cstegr.f cstein.f csteqr.f csycon.f csymv.f csyr.f csyrfs.f csysv.f csysvx.f csytf2.f csytrf.f csytri.f csytri2.f csytri2x.f csyswapr.f diff --git a/SRC/cgetf2.f b/SRC/cgetf2.f index aac9899702..776d6deb76 100644 --- a/SRC/cgetf2.f +++ b/SRC/cgetf2.f @@ -189,6 +189,7 @@ SUBROUTINE CGETF2( M, N, A, LDA, IPIV, INFO ) A( J+I, J ) = A( J+I, J ) / A( J, J ) 20 CONTINUE END IF + ! CALL CRSCL( M-J, A( J, J ), A( J+1, J ), 1 ) END IF * ELSE IF( INFO.EQ.0 ) THEN diff --git a/SRC/crscl.f b/SRC/crscl.f new file mode 100644 index 0000000000..71449d1950 --- /dev/null +++ b/SRC/crscl.f @@ -0,0 +1,208 @@ +*> \brief \b CRSCL multiplies a vector by the reciprocal of a real scalar. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CRSCL + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CRSCL( N, A, X, INCX ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* COMPLEX A +* .. +* .. Array Arguments .. +* COMPLEX X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CRSCL multiplies an n-element complex vector x by the complex scalar +*> 1/a. This is done without overflow or underflow as long as +*> the final result x/a does not overflow or underflow. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of components of the vector x. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX +*> The scalar a which is used to divide each component of x. +*> A must not be 0, or the subroutine will divide by zero. +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is COMPLEX array, dimension +*> (1+(N-1)*abs(INCX)) +*> The n-element vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between successive values of the vector X. +*> > 0: X(1) = X(1) and X(1+(i-1)*INCX) = x(i), 1< i<= n +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup complexOTHERauxiliary +* +* ===================================================================== + SUBROUTINE CRSCL( N, A, X, INCX ) +* +* -- 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 + COMPLEX A +* .. +* .. Array Arguments .. + COMPLEX X( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + REAL BIGNUM, SMLNUM, HUGE, AR, AI, ABSR, ABSI, UR + % , UI + COMPLEX INVA +* .. +* .. External Functions .. + REAL SLAMCH + COMPLEX CLADIV + EXTERNAL SLAMCH, CLADIV +* .. +* .. External Subroutines .. + EXTERNAL CSCAL, CSSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.LE.0 ) + $ RETURN +* +* Get machine parameters +* + SMLNUM = SLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM + HUGE = SLAMCH( 'O' ) +* +* Initialize constants related to A. +* + AR = REAL( A ) + AI = AIMAG( A ) + ABSR = ABS( AR ) + ABSI = ABS( AI ) +* + IF( ABSI.EQ.ZERO ) THEN +* If alpha is real, then we can use csrscl + CALL CSRSCL( N, AR, X, INCX ) +* + ELSE IF( ABSR.EQ.ZERO ) THEN +* If alpha has a zero real part, then we follow the same rules as if +* alpha were real. + IF( ABSI.GT.BIGNUM ) THEN + INVA = CMPLX( ZERO, -BIGNUM / AI ) + CALL CSSCAL( N, SMLNUM, X, INCX ) + CALL CSCAL( N, INVA, X, INCX ) + ELSE IF( ABSI.LT.SMLNUM ) THEN + INVA = CMPLX( ZERO, -SMLNUM / AI ) + CALL CSCAL( N, INVA, X, INCX ) + CALL CSSCAL( N, BIGNUM, X, INCX ) + ELSE + INVA = CMPLX( ZERO, -ONE / AI ) + CALL CSCAL( N, INVA, X, INCX ) + END IF +* + ELSE IF( (ABSR.GE.BIGNUM).OR.(ABSI.GE.BIGNUM) ) THEN +* Either real or imaginary part is too large. + INVA = CLADIV( CMPLX( BIGNUM, ZERO ), A ) + CALL CSSCAL( N, SMLNUM, X, INCX ) + CALL CSCAL( N, INVA, X, INCX ) +* + ELSE +* The following numbers can be computed without NaNs and zeros. +* They do not overflow simultaneously. +* They are the inverse of the real and imaginary parts of 1/alpha. + UR = AR + AI * ( AI / AR ) + UI = AI + AR * ( AR / AI ) +* + IF( (ABS( UR ).LT.SMLNUM).OR.(ABS( UI ).LT.SMLNUM) ) THEN + INVA = CMPLX( SMLNUM / UR, -SMLNUM / UI ) + CALL CSCAL( N, INVA, X, INCX ) + CALL CSSCAL( N, BIGNUM, X, INCX ) + ELSE IF( ABS( UR ).GT.HUGE ) THEN + IF( ABSR.GE.ABSI ) THEN + UR = (SMLNUM * AR) + AI * (SMLNUM * (AI / AR)) + ELSE + UR = (SMLNUM * AR) + AI * ((SMLNUM * AI) / AR) + END IF + INVA = CMPLX( ONE / UR, -BIGNUM / UI ) + CALL CSSCAL( N, SMLNUM, X, INCX ) + CALL CSCAL( N, INVA, X, INCX ) + ELSE IF( ABS( UI ).GT.HUGE ) THEN + IF( ABSI.GE.ABSR ) THEN + UI = (SMLNUM * AI) + AR * (SMLNUM * (AR / AI)) + ELSE + UI = (SMLNUM * AI) + AR * ((SMLNUM * AR) / AI) + END IF + INVA = CMPLX( BIGNUM / UR, -ONE / UI ) + CALL CSSCAL( N, SMLNUM, X, INCX ) + CALL CSCAL( N, INVA, X, INCX ) + ELSE IF( (ABS( UR ).GT.BIGNUM).OR.(ABS( UI ).GT.BIGNUM) ) THEN + INVA = CMPLX( BIGNUM / UR, -BIGNUM / UI ) + CALL CSSCAL( N, SMLNUM, X, INCX ) + CALL CSCAL( N, INVA, X, INCX ) + ELSE + INVA = CMPLX( ONE / UR, -ONE / UI ) + CALL CSCAL( N, INVA, X, INCX ) + END IF + END IF +* + RETURN +* +* End of CRSCL +* + END From 59a1772336ff389b5ddc1a5850db7e1730c8d02b Mon Sep 17 00:00:00 2001 From: Weslley S Pereira Date: Wed, 31 May 2023 18:29:30 -0600 Subject: [PATCH 2/5] (Last?) version of RSCL --- SRC/Makefile | 2 +- SRC/cgetf2.f | 24 +++---------- SRC/crscl.f | 98 ++++++++++++++++++++++++---------------------------- 3 files changed, 52 insertions(+), 72 deletions(-) diff --git a/SRC/Makefile b/SRC/Makefile index 35b8c64aea..ef8ff179a7 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -260,7 +260,7 @@ CLASRC = \ cppcon.o cppequ.o cpprfs.o cppsv.o cppsvx.o cpptrf.o cpptri.o cpptrs.o \ cptcon.o cpteqr.o cptrfs.o cptsv.o cptsvx.o cpttrf.o cpttrs.o cptts2.o \ crot.o cspcon.o cspmv.o cspr.o csprfs.o cspsv.o \ - cspsvx.o csptrf.o csptri.o csptrs.o csrscl.o cstedc.o \ + cspsvx.o csptrf.o csptri.o csptrs.o csrscl.o crscl.o cstedc.o \ cstegr.o cstein.o csteqr.o \ csycon.o csymv.o \ csyr.o csyrfs.o csysv.o csysvx.o csytf2.o csytrf.o csytri.o csytri2.o csytri2x.o \ diff --git a/SRC/cgetf2.f b/SRC/cgetf2.f index 776d6deb76..011dad95f7 100644 --- a/SRC/cgetf2.f +++ b/SRC/cgetf2.f @@ -126,16 +126,14 @@ SUBROUTINE CGETF2( M, N, A, LDA, IPIV, INFO ) $ ZERO = ( 0.0E+0, 0.0E+0 ) ) * .. * .. Local Scalars .. - REAL SFMIN - INTEGER I, J, JP + INTEGER J, JP * .. * .. External Functions .. - REAL SLAMCH INTEGER ICAMAX - EXTERNAL SLAMCH, ICAMAX + EXTERNAL ICAMAX * .. * .. External Subroutines .. - EXTERNAL CGERU, CSCAL, CSWAP, XERBLA + EXTERNAL CGERU, CRSCL, CSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -161,10 +159,6 @@ SUBROUTINE CGETF2( M, N, A, LDA, IPIV, INFO ) * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN -* -* Compute machine safe minimum -* - SFMIN = SLAMCH('S') * DO 10 J = 1, MIN( M, N ) * @@ -181,16 +175,8 @@ SUBROUTINE CGETF2( M, N, A, LDA, IPIV, INFO ) * * Compute elements J+1:M of J-th column. * - IF( J.LT.M ) THEN - IF( ABS(A( J, J )) .GE. SFMIN ) THEN - CALL CSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) - ELSE - DO 20 I = 1, M-J - A( J+I, J ) = A( J+I, J ) / A( J, J ) - 20 CONTINUE - END IF - ! CALL CRSCL( M-J, A( J, J ), A( J+1, J ), 1 ) - END IF + IF( J.LT.M ) + $ CALL CRSCL( M-J, A( J, J ), A( J+1, J ), 1 ) * ELSE IF( INFO.EQ.0 ) THEN * diff --git a/SRC/crscl.f b/SRC/crscl.f index 71449d1950..22919cd62c 100644 --- a/SRC/crscl.f +++ b/SRC/crscl.f @@ -101,9 +101,8 @@ SUBROUTINE CRSCL( N, A, X, INCX ) PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. - REAL BIGNUM, SMLNUM, HUGE, AR, AI, ABSR, ABSI, UR + REAL SAFMAX, SAFMIN, OV, AR, AI, ABSR, ABSI, UR % , UI - COMPLEX INVA * .. * .. External Functions .. REAL SLAMCH @@ -111,7 +110,7 @@ SUBROUTINE CRSCL( N, A, X, INCX ) EXTERNAL SLAMCH, CLADIV * .. * .. External Subroutines .. - EXTERNAL CSCAL, CSSCAL + EXTERNAL CSCAL, CSSCAL, CSRSCL * .. * .. Intrinsic Functions .. INTRINSIC ABS @@ -125,9 +124,9 @@ SUBROUTINE CRSCL( N, A, X, INCX ) * * Get machine parameters * - SMLNUM = SLAMCH( 'S' ) - BIGNUM = ONE / SMLNUM - HUGE = SLAMCH( 'O' ) + SAFMIN = SLAMCH( 'S' ) + SAFMAX = ONE / SAFMIN + OV = SLAMCH( 'O' ) * * Initialize constants related to A. * @@ -136,68 +135,63 @@ SUBROUTINE CRSCL( N, A, X, INCX ) ABSR = ABS( AR ) ABSI = ABS( AI ) * - IF( ABSI.EQ.ZERO ) THEN + IF( AI.EQ.ZERO ) THEN * If alpha is real, then we can use csrscl CALL CSRSCL( N, AR, X, INCX ) * - ELSE IF( ABSR.EQ.ZERO ) THEN + ELSE IF( AR.EQ.ZERO ) THEN * If alpha has a zero real part, then we follow the same rules as if * alpha were real. - IF( ABSI.GT.BIGNUM ) THEN - INVA = CMPLX( ZERO, -BIGNUM / AI ) - CALL CSSCAL( N, SMLNUM, X, INCX ) - CALL CSCAL( N, INVA, X, INCX ) - ELSE IF( ABSI.LT.SMLNUM ) THEN - INVA = CMPLX( ZERO, -SMLNUM / AI ) - CALL CSCAL( N, INVA, X, INCX ) - CALL CSSCAL( N, BIGNUM, X, INCX ) + IF( ABSI.GT.SAFMAX ) THEN + CALL CSSCAL( N, SAFMIN, X, INCX ) + CALL CSCAL( N, CMPLX( ZERO, -SAFMAX / AI ), X, INCX ) + ELSE IF( ABSI.LT.SAFMIN ) THEN + CALL CSCAL( N, CMPLX( ZERO, -SAFMIN / AI ), X, INCX ) + CALL CSSCAL( N, SAFMAX, X, INCX ) ELSE - INVA = CMPLX( ZERO, -ONE / AI ) - CALL CSCAL( N, INVA, X, INCX ) + CALL CSCAL( N, CMPLX( ZERO, -ONE / AI ), X, INCX ) END IF -* - ELSE IF( (ABSR.GE.BIGNUM).OR.(ABSI.GE.BIGNUM) ) THEN -* Either real or imaginary part is too large. - INVA = CLADIV( CMPLX( BIGNUM, ZERO ), A ) - CALL CSSCAL( N, SMLNUM, X, INCX ) - CALL CSCAL( N, INVA, X, INCX ) * ELSE -* The following numbers can be computed without NaNs and zeros. -* They do not overflow simultaneously. +* The following numbers can be computed. * They are the inverse of the real and imaginary parts of 1/alpha. +* Note that a and b are always different from zero. +* NaNs are only possible if either: +* 1. alphaR or alphaI is NaN. +* 2. alphaR and alphaI are both infinite, in which case it makes sense +* to propagate a NaN. UR = AR + AI * ( AI / AR ) UI = AI + AR * ( AR / AI ) * - IF( (ABS( UR ).LT.SMLNUM).OR.(ABS( UI ).LT.SMLNUM) ) THEN - INVA = CMPLX( SMLNUM / UR, -SMLNUM / UI ) - CALL CSCAL( N, INVA, X, INCX ) - CALL CSSCAL( N, BIGNUM, X, INCX ) - ELSE IF( ABS( UR ).GT.HUGE ) THEN - IF( ABSR.GE.ABSI ) THEN - UR = (SMLNUM * AR) + AI * (SMLNUM * (AI / AR)) - ELSE - UR = (SMLNUM * AR) + AI * ((SMLNUM * AI) / AR) - END IF - INVA = CMPLX( ONE / UR, -BIGNUM / UI ) - CALL CSSCAL( N, SMLNUM, X, INCX ) - CALL CSCAL( N, INVA, X, INCX ) - ELSE IF( ABS( UI ).GT.HUGE ) THEN - IF( ABSI.GE.ABSR ) THEN - UI = (SMLNUM * AI) + AR * (SMLNUM * (AR / AI)) + IF( (ABS( UR ).LT.SAFMIN).OR.(ABS( UI ).LT.SAFMIN) ) THEN +* This means that both alphaR and alphaI are very small. + CALL CSCAL( N, CMPLX( SAFMIN / UR, -SAFMIN / UI ), X, INCX ) + CALL CSSCAL( N, SAFMAX, X, INCX ) + ELSE IF( (ABS( UR ).GT.SAFMAX).OR.(ABS( UI ).GT.SAFMAX) ) THEN + IF( (ABSR.GT.OV).OR.(ABSI.GT.OV) ) THEN +* This means that a and b are both Inf. No need for scaling. + CALL CSCAL( N, CMPLX( ONE / UR, -ONE / UI ), X, INCX ) ELSE - UI = (SMLNUM * AI) + AR * ((SMLNUM * AR) / AI) + CALL CSSCAL( N, SAFMIN, X, INCX ) + IF( (ABS( UR ).GT.OV).OR.(ABS( UI ).GT.OV) ) THEN +* Infs were generated. We do proper scaling to avoid them. + IF( ABSR.GE.ABSI ) THEN +* ABS( UR ) <= ABS( UI ) + UR = (SAFMIN * AR) + SAFMIN * (AI * ( AI / AR )) + UI = (SAFMIN * AI) + AR * ( (SAFMIN * AR) / AI ) + ELSE +* ABS( UR ) > ABS( UI ) + UR = (SAFMIN * AR) + AI * ( (SAFMIN * AI) / AR ) + UI = (SAFMIN * AI) + SAFMIN * (AR * ( AR / AI )) + END IF + CALL CSCAL( N, CMPLX( ONE / UR, -ONE / UI ), X, INCX ) + ELSE + CALL CSCAL( N, CMPLX( SAFMAX / UR, -SAFMAX / UI ), + $ X, INCX ) + END IF END IF - INVA = CMPLX( BIGNUM / UR, -ONE / UI ) - CALL CSSCAL( N, SMLNUM, X, INCX ) - CALL CSCAL( N, INVA, X, INCX ) - ELSE IF( (ABS( UR ).GT.BIGNUM).OR.(ABS( UI ).GT.BIGNUM) ) THEN - INVA = CMPLX( BIGNUM / UR, -BIGNUM / UI ) - CALL CSSCAL( N, SMLNUM, X, INCX ) - CALL CSCAL( N, INVA, X, INCX ) ELSE - INVA = CMPLX( ONE / UR, -ONE / UI ) - CALL CSCAL( N, INVA, X, INCX ) + CALL CSCAL( N, CMPLX( ONE / UR, -ONE / UI ), X, INCX ) END IF END IF * From bfcece32ec4eb864def08fcb93581b20fac9bc7c Mon Sep 17 00:00:00 2001 From: Weslley S Pereira Date: Mon, 26 Jun 2023 09:34:24 -0600 Subject: [PATCH 3/5] Adding zrscl.f --- SRC/zgetf2.f | 15 +--- SRC/zrscl.f | 201 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 205 insertions(+), 11 deletions(-) create mode 100644 SRC/zrscl.f diff --git a/SRC/zgetf2.f b/SRC/zgetf2.f index c247f8645a..4827d05987 100644 --- a/SRC/zgetf2.f +++ b/SRC/zgetf2.f @@ -127,7 +127,7 @@ SUBROUTINE ZGETF2( M, N, A, LDA, IPIV, INFO ) * .. * .. Local Scalars .. DOUBLE PRECISION SFMIN - INTEGER I, J, JP + INTEGER J, JP * .. * .. External Functions .. DOUBLE PRECISION DLAMCH @@ -135,7 +135,7 @@ SUBROUTINE ZGETF2( M, N, A, LDA, IPIV, INFO ) EXTERNAL DLAMCH, IZAMAX * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZGERU, ZSCAL, ZSWAP + EXTERNAL XERBLA, ZGERU, ZRSCL, ZSWAP * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -181,15 +181,8 @@ SUBROUTINE ZGETF2( M, N, A, LDA, IPIV, INFO ) * * Compute elements J+1:M of J-th column. * - IF( J.LT.M ) THEN - IF( ABS(A( J, J )) .GE. SFMIN ) THEN - CALL ZSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) - ELSE - DO 20 I = 1, M-J - A( J+I, J ) = A( J+I, J ) / A( J, J ) - 20 CONTINUE - END IF - END IF + IF( J.LT.M ) + $ CALL ZRSCL( M-J, A( J, J ), A( J+1, J ), 1 ) * ELSE IF( INFO.EQ.0 ) THEN * diff --git a/SRC/zrscl.f b/SRC/zrscl.f new file mode 100644 index 0000000000..8f4db4053a --- /dev/null +++ b/SRC/zrscl.f @@ -0,0 +1,201 @@ +*> \brief \b ZDRSCL multiplies a vector by the reciprocal of a real scalar. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZDRSCL + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZRSCL( N, A, X, INCX ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* COMPLEX*16 A +* .. +* .. Array Arguments .. +* COMPLEX*16 X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZRSCL multiplies an n-element complex vector x by the complex scalar +*> 1/a. This is done without overflow or underflow as long as +*> the final result x/a does not overflow or underflow. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of components of the vector x. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX*16 +*> The scalar a which is used to divide each component of x. +*> A must not be 0, or the subroutine will divide by zero. +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension +*> (1+(N-1)*abs(INCX)) +*> The n-element vector x. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between successive values of the vector SX. +*> > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= n +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup complex16OTHERauxiliary +* +* ===================================================================== + SUBROUTINE ZRSCL( N, A, X, INCX ) +* +* -- 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 + COMPLEX*16 A +* .. +* .. Array Arguments .. + COMPLEX*16 X( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + DOUBLE PRECISION SAFMAX, SAFMIN, OV, AR, AI, ABSR, ABSI, UR, UI +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH + COMPLEX ZLADIV + EXTERNAL DLAMCH, ZLADIV +* .. +* .. External Subroutines .. + EXTERNAL DSCAL, ZDSCAL, ZDRSCL +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.LE.0 ) + $ RETURN +* +* Get machine parameters +* + SAFMIN = DLAMCH( 'S' ) + SAFMAX = ONE / SAFMIN + OV = DLAMCH( 'O' ) +* +* Initialize constants related to A. +* + AR = DBLE( A ) + AI = DIMAG( A ) + ABSR = ABS( AR ) + ABSI = ABS( AI ) +* + IF( AI.EQ.ZERO ) THEN +* If alpha is real, then we can use csrscl + CALL ZDRSCL( N, AR, X, INCX ) +* + ELSE IF( AR.EQ.ZERO ) THEN +* If alpha has a zero real part, then we follow the same rules as if +* alpha were real. + IF( ABSI.GT.SAFMAX ) THEN + CALL ZDSCAL( N, SAFMIN, X, INCX ) + CALL ZSCAL( N, DCMPLX( ZERO, -SAFMAX / AI ), X, INCX ) + ELSE IF( ABSI.LT.SAFMIN ) THEN + CALL ZSCAL( N, DCMPLX( ZERO, -SAFMIN / AI ), X, INCX ) + CALL ZDSCAL( N, SAFMAX, X, INCX ) + ELSE + CALL ZSCAL( N, DCMPLX( ZERO, -ONE / AI ), X, INCX ) + END IF +* + ELSE +* The following numbers can be computed. +* They are the inverse of the real and imaginary parts of 1/alpha. +* Note that a and b are always different from zero. +* NaNs are only possible if either: +* 1. alphaR or alphaI is NaN. +* 2. alphaR and alphaI are both infinite, in which case it makes sense +* to propagate a NaN. + UR = AR + AI * ( AI / AR ) + UI = AI + AR * ( AR / AI ) +* + IF( (ABS( UR ).LT.SAFMIN).OR.(ABS( UI ).LT.SAFMIN) ) THEN +* This means that both alphaR and alphaI are very small. + CALL ZSCAL( N, DCMPLX( SAFMIN / UR, -SAFMIN / UI ), X, INCX ) + CALL ZDSCAL( N, SAFMAX, X, INCX ) + ELSE IF( (ABS( UR ).GT.SAFMAX).OR.(ABS( UI ).GT.SAFMAX) ) THEN + IF( (ABSR.GT.OV).OR.(ABSI.GT.OV) ) THEN +* This means that a and b are both Inf. No need for scaling. + CALL ZSCAL( N, DCMPLX( ONE / UR, -ONE / UI ), X, INCX ) + ELSE + CALL ZDSCAL( N, SAFMIN, X, INCX ) + IF( (ABS( UR ).GT.OV).OR.(ABS( UI ).GT.OV) ) THEN +* Infs were generated. We do proper scaling to avoid them. + IF( ABSR.GE.ABSI ) THEN +* ABS( UR ) <= ABS( UI ) + UR = (SAFMIN * AR) + SAFMIN * (AI * ( AI / AR )) + UI = (SAFMIN * AI) + AR * ( (SAFMIN * AR) / AI ) + ELSE +* ABS( UR ) > ABS( UI ) + UR = (SAFMIN * AR) + AI * ( (SAFMIN * AI) / AR ) + UI = (SAFMIN * AI) + SAFMIN * (AR * ( AR / AI )) + END IF + CALL ZSCAL( N, DCMPLX( ONE / UR, -ONE / UI ), X, INCX ) + ELSE + CALL ZSCAL( N, DCMPLX( SAFMAX / UR, -SAFMAX / UI ), + $ X, INCX ) + END IF + END IF + ELSE + CALL ZSCAL( N, DCMPLX( ONE / UR, -ONE / UI ), X, INCX ) + END IF + END IF +* + RETURN +* +* End of ZRSCL +* + END From c9a49d98cf3f96eba70c4ac9bd16401c8b37215e Mon Sep 17 00:00:00 2001 From: Weslley S Pereira Date: Mon, 26 Jun 2023 09:39:40 -0600 Subject: [PATCH 4/5] Adds zrscl in the build environments --- SRC/CMakeLists.txt | 2 +- SRC/Makefile | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index 1556b31dc7..2ac56662a8 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -427,7 +427,7 @@ set(ZLASRC zppcon.f zppequ.f zpprfs.f zppsv.f zppsvx.f zpptrf.f zpptri.f zpptrs.f zptcon.f zpteqr.f zptrfs.f zptsv.f zptsvx.f zpttrf.f zpttrs.f zptts2.f zrot.f zspcon.f zspmv.f zspr.f zsprfs.f zspsv.f - zspsvx.f zsptrf.f zsptri.f zsptrs.f zdrscl.f zstedc.f + zspsvx.f zsptrf.f zsptri.f zsptrs.f zdrscl.f zrscl.f zstedc.f zstegr.f zstein.f zsteqr.f zsycon.f zsymv.f zsyr.f zsyrfs.f zsysv.f zsysvx.f zsytf2.f zsytrf.f zsytri.f zsytri2.f zsytri2x.f zsyswapr.f diff --git a/SRC/Makefile b/SRC/Makefile index ef8ff179a7..3bdf8bdf00 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -464,7 +464,7 @@ ZLASRC = \ zppcon.o zppequ.o zpprfs.o zppsv.o zppsvx.o zpptrf.o zpptri.o zpptrs.o \ zptcon.o zpteqr.o zptrfs.o zptsv.o zptsvx.o zpttrf.o zpttrs.o zptts2.o \ zrot.o zspcon.o zspmv.o zspr.o zsprfs.o zspsv.o \ - zspsvx.o zsptrf.o zsptri.o zsptrs.o zdrscl.o zstedc.o \ + zspsvx.o zsptrf.o zsptri.o zsptrs.o zdrscl.o zrscl.o zstedc.o \ zstegr.o zstein.o zsteqr.o \ zsycon.o zsymv.o \ zsyr.o zsyrfs.o zsysv.o zsysvx.o zsytf2.o zsytrf.o zsytri.o zsytri2.o zsytri2x.o \ From d212879672abb6f4a1d14711b742da2dad785e6c Mon Sep 17 00:00:00 2001 From: Weslley S Pereira Date: Mon, 26 Jun 2023 10:12:50 -0600 Subject: [PATCH 5/5] Fix bugs in ZRSCL --- SRC/zrscl.f | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/SRC/zrscl.f b/SRC/zrscl.f index 8f4db4053a..970f6de752 100644 --- a/SRC/zrscl.f +++ b/SRC/zrscl.f @@ -105,7 +105,7 @@ SUBROUTINE ZRSCL( N, A, X, INCX ) * .. * .. External Functions .. DOUBLE PRECISION DLAMCH - COMPLEX ZLADIV + COMPLEX*16 ZLADIV EXTERNAL DLAMCH, ZLADIV * .. * .. External Subroutines .. @@ -164,7 +164,8 @@ SUBROUTINE ZRSCL( N, A, X, INCX ) * IF( (ABS( UR ).LT.SAFMIN).OR.(ABS( UI ).LT.SAFMIN) ) THEN * This means that both alphaR and alphaI are very small. - CALL ZSCAL( N, DCMPLX( SAFMIN / UR, -SAFMIN / UI ), X, INCX ) + CALL ZSCAL( N, DCMPLX( SAFMIN / UR, -SAFMIN / UI ), X, + $ INCX ) CALL ZDSCAL( N, SAFMAX, X, INCX ) ELSE IF( (ABS( UR ).GT.SAFMAX).OR.(ABS( UI ).GT.SAFMAX) ) THEN IF( (ABSR.GT.OV).OR.(ABSI.GT.OV) ) THEN @@ -183,7 +184,8 @@ SUBROUTINE ZRSCL( N, A, X, INCX ) UR = (SAFMIN * AR) + AI * ( (SAFMIN * AI) / AR ) UI = (SAFMIN * AI) + SAFMIN * (AR * ( AR / AI )) END IF - CALL ZSCAL( N, DCMPLX( ONE / UR, -ONE / UI ), X, INCX ) + CALL ZSCAL( N, DCMPLX( ONE / UR, -ONE / UI ), X, + $ INCX ) ELSE CALL ZSCAL( N, DCMPLX( SAFMAX / UR, -SAFMAX / UI ), $ X, INCX )