diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index 070eeb8100..82b2973825 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -50,7 +50,7 @@ set(SCLAUX slapy2.f slapy3.f slarnv.f slarra.f slarrb.f slarrc.f slarrd.f slarre.f slarrf.f slarrj.f slarrk.f slarrr.f slaneg.f - slartg.f slaruv.f slas2.f slascl.f + slartg.f90 slaruv.f slas2.f slascl.f 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 @@ -69,7 +69,7 @@ set(DZLAUX dlapy2.f dlapy3.f dlarnv.f dlarra.f dlarrb.f dlarrc.f dlarrd.f dlarre.f dlarrf.f dlarrj.f dlarrk.f dlarrr.f dlaneg.f - dlartg.f dlaruv.f dlas2.f dlascl.f + dlartg.f90 dlaruv.f dlas2.f dlascl.f 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 @@ -210,7 +210,7 @@ set(CLASRC claqr0.f claqr1.f claqr2.f claqr3.f claqr4.f claqr5.f 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 + clarfx.f clarfy.f clargv.f clarnv.f clarrv.f clartg.f90 clartv.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 @@ -405,7 +405,7 @@ set(ZLASRC zlaqsp.f zlaqsy.f zlar1v.f zlar2v.f ilazlr.f ilazlc.f zlarcm.f zlarf.f zlarfb.f zlarfb_gett.f zlarfg.f zlarfgp.f zlarft.f - zlarfx.f zlarfy.f zlargv.f zlarnv.f zlarrv.f zlartg.f zlartv.f + zlarfx.f zlarfy.f zlargv.f zlarnv.f zlarrv.f zlartg.f90 zlartv.f zlarz.f zlarzb.f zlarzt.f zlascl.f zlaset.f zlasr.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 diff --git a/SRC/clartg.f b/SRC/clartg.f deleted file mode 100644 index 308900797d..0000000000 --- a/SRC/clartg.f +++ /dev/null @@ -1,247 +0,0 @@ -*> \brief \b CLARTG generates a plane rotation with real cosine and complex sine. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download CLARTG + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE CLARTG( F, G, CS, SN, R ) -* -* .. Scalar Arguments .. -* REAL CS -* COMPLEX F, G, R, SN -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> CLARTG generates a plane rotation so that -*> -*> [ CS SN ] [ F ] [ R ] -*> [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1. -*> [ -SN CS ] [ G ] [ 0 ] -*> -*> This is a faster version of the BLAS1 routine CROTG, except for -*> the following differences: -*> F and G are unchanged on return. -*> If G=0, then CS=1 and SN=0. -*> If F=0, then CS=0 and SN is chosen so that R is real. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] F -*> \verbatim -*> F is COMPLEX -*> The first component of vector to be rotated. -*> \endverbatim -*> -*> \param[in] G -*> \verbatim -*> G is COMPLEX -*> The second component of vector to be rotated. -*> \endverbatim -*> -*> \param[out] CS -*> \verbatim -*> CS is REAL -*> The cosine of the rotation. -*> \endverbatim -*> -*> \param[out] SN -*> \verbatim -*> SN is COMPLEX -*> The sine of the rotation. -*> \endverbatim -*> -*> \param[out] R -*> \verbatim -*> R is COMPLEX -*> The nonzero component of the rotated vector. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup complexOTHERauxiliary -* -*> \par Further Details: -* ===================== -*> -*> \verbatim -*> -*> 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel -*> -*> This version has a few statements commented out for thread safety -*> (machine parameters are computed on each entry). 10 feb 03, SJH. -*> \endverbatim -*> -* ===================================================================== - SUBROUTINE CLARTG( F, G, CS, SN, R ) -* -* -- 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 .. - REAL CS - COMPLEX F, G, R, SN -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL TWO, ONE, ZERO - PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) - COMPLEX CZERO - PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) -* .. -* .. Local Scalars .. -* LOGICAL FIRST - INTEGER COUNT, I - REAL D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN, - $ SAFMN2, SAFMX2, SCALE - COMPLEX FF, FS, GS -* .. -* .. External Functions .. - REAL SLAMCH, SLAPY2 - LOGICAL SISNAN - EXTERNAL SLAMCH, SLAPY2, SISNAN -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, AIMAG, CMPLX, CONJG, INT, LOG, MAX, REAL, - $ SQRT -* .. -* .. Statement Functions .. - REAL ABS1, ABSSQ -* .. -* .. Statement Function definitions .. - ABS1( FF ) = MAX( ABS( REAL( FF ) ), ABS( AIMAG( FF ) ) ) - ABSSQ( FF ) = REAL( FF )**2 + AIMAG( FF )**2 -* .. -* .. Executable Statements .. -* - SAFMIN = SLAMCH( 'S' ) - EPS = SLAMCH( 'E' ) - SAFMN2 = SLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / - $ LOG( SLAMCH( 'B' ) ) / TWO ) - SAFMX2 = ONE / SAFMN2 - SCALE = MAX( ABS1( F ), ABS1( G ) ) - FS = F - GS = G - COUNT = 0 - IF( SCALE.GE.SAFMX2 ) THEN - 10 CONTINUE - COUNT = COUNT + 1 - FS = FS*SAFMN2 - GS = GS*SAFMN2 - SCALE = SCALE*SAFMN2 - IF( SCALE.GE.SAFMX2 .AND. COUNT .LT. 20) - $ GO TO 10 - ELSE IF( SCALE.LE.SAFMN2 ) THEN - IF( G.EQ.CZERO.OR.SISNAN( ABS( G ) ) ) THEN - CS = ONE - SN = CZERO - R = F - RETURN - END IF - 20 CONTINUE - COUNT = COUNT - 1 - FS = FS*SAFMX2 - GS = GS*SAFMX2 - SCALE = SCALE*SAFMX2 - IF( SCALE.LE.SAFMN2 ) - $ GO TO 20 - END IF - F2 = ABSSQ( FS ) - G2 = ABSSQ( GS ) - IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN -* -* This is a rare case: F is very small. -* - IF( F.EQ.CZERO ) THEN - CS = ZERO - R = SLAPY2( REAL( G ), AIMAG( G ) ) -* Do complex/real division explicitly with two real divisions - D = SLAPY2( REAL( GS ), AIMAG( GS ) ) - SN = CMPLX( REAL( GS ) / D, -AIMAG( GS ) / D ) - RETURN - END IF - F2S = SLAPY2( REAL( FS ), AIMAG( FS ) ) -* G2 and G2S are accurate -* G2 is at least SAFMIN, and G2S is at least SAFMN2 - G2S = SQRT( G2 ) -* Error in CS from underflow in F2S is at most -* UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS -* If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN, -* and so CS .lt. sqrt(SAFMIN) -* If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN -* and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS) -* Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S - CS = F2S / G2S -* Make sure abs(FF) = 1 -* Do complex/real division explicitly with 2 real divisions - IF( ABS1( F ).GT.ONE ) THEN - D = SLAPY2( REAL( F ), AIMAG( F ) ) - FF = CMPLX( REAL( F ) / D, AIMAG( F ) / D ) - ELSE - DR = SAFMX2*REAL( F ) - DI = SAFMX2*AIMAG( F ) - D = SLAPY2( DR, DI ) - FF = CMPLX( DR / D, DI / D ) - END IF - SN = FF*CMPLX( REAL( GS ) / G2S, -AIMAG( GS ) / G2S ) - R = CS*F + SN*G - ELSE -* -* This is the most common case. -* Neither F2 nor F2/G2 are less than SAFMIN -* F2S cannot overflow, and it is accurate -* - F2S = SQRT( ONE+G2 / F2 ) -* Do the F2S(real)*FS(complex) multiply with two real multiplies - R = CMPLX( F2S*REAL( FS ), F2S*AIMAG( FS ) ) - CS = ONE / F2S - D = F2 + G2 -* Do complex/real division explicitly with two real divisions - SN = CMPLX( REAL( R ) / D, AIMAG( R ) / D ) - SN = SN*CONJG( GS ) - IF( COUNT.NE.0 ) THEN - IF( COUNT.GT.0 ) THEN - DO 30 I = 1, COUNT - R = R*SAFMX2 - 30 CONTINUE - ELSE - DO 40 I = 1, -COUNT - R = R*SAFMN2 - 40 CONTINUE - END IF - END IF - END IF - RETURN -* -* End of CLARTG -* - END diff --git a/SRC/clartg.f90 b/SRC/clartg.f90 new file mode 100644 index 0000000000..f63a0f8d20 --- /dev/null +++ b/SRC/clartg.f90 @@ -0,0 +1,233 @@ +!> \brief \b CLARTG generates a plane rotation with real cosine and complex sine. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! SUBROUTINE CLARTG( F, G, C, S, R ) +! +! .. Scalar Arguments .. +! REAL(wp) C +! COMPLEX(wp) F, G, R, S +! .. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> CLARTG generates a plane rotation so that +!> +!> [ C S ] . [ F ] = [ R ] +!> [ -conjg(S) C ] [ G ] [ 0 ] +!> +!> where C is real and C**2 + |S|**2 = 1. +!> +!> The mathematical formulas used for C and S are +!> +!> sgn(x) = { x / |x|, x != 0 +!> { 1, x = 0 +!> +!> R = sgn(F) * sqrt(|F|**2 + |G|**2) +!> +!> C = |F| / sqrt(|F|**2 + |G|**2) +!> +!> S = sgn(F) * conjg(G) / sqrt(|F|**2 + |G|**2) +!> +!> When F and G are real, the formulas simplify to C = F/R and +!> S = G/R, and the returned values of C, S, and R should be +!> identical to those returned by CLARTG. +!> +!> The algorithm used to compute these quantities incorporates scaling +!> to avoid overflow or underflow in computing the square root of the +!> sum of squares. +!> +!> This is a faster version of the BLAS1 routine CROTG, except for +!> the following differences: +!> F and G are unchanged on return. +!> If G=0, then C=1 and S=0. +!> If F=0, then C=0 and S is chosen so that R is real. +!> +!> Below, wp=>sp stands for single precision from LA_CONSTANTS module. +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] F +!> \verbatim +!> F is COMPLEX(wp) +!> The first component of vector to be rotated. +!> \endverbatim +!> +!> \param[in] G +!> \verbatim +!> G is COMPLEX(wp) +!> The second component of vector to be rotated. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is REAL(wp) +!> The cosine of the rotation. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is COMPLEX(wp) +!> The sine of the rotation. +!> \endverbatim +!> +!> \param[out] R +!> \verbatim +!> R is COMPLEX(wp) +!> The nonzero component of the rotated vector. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \date August 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 +!> +!> \endverbatim +! +subroutine CLARTG( f, g, c, s, r ) + use LA_CONSTANTS, & + only: wp=>sp, zero=>szero, one=>sone, two=>stwo, czero, & + rtmin=>srtmin, rtmax=>srtmax, safmin=>ssafmin, safmax=>ssafmax +! +! -- LAPACK auxiliary routine (version 3.10.0) -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! February 2021 +! +! .. Scalar Arguments .. + real(wp) c + complex(wp) f, g, r, s +! .. +! .. Local Scalars .. + real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w + complex(wp) :: fs, gs, t +! .. +! .. Intrinsic Functions .. + intrinsic :: abs, aimag, conjg, max, min, real, sqrt +! .. +! .. Statement Functions .. + real(wp) :: ABSSQ +! .. +! .. Statement Function definitions .. + ABSSQ( t ) = real( t )**2 + aimag( t )**2 +! .. +! .. Executable Statements .. +! + if( g == czero ) then + c = one + s = czero + r = f + else if( f == czero ) then + c = zero + g1 = max( abs(real(g)), abs(aimag(g)) ) + if( g1 > rtmin .and. g1 < rtmax ) then +! +! Use unscaled algorithm +! + g2 = ABSSQ( g ) + d = sqrt( g2 ) + s = conjg( g ) / d + r = d + else +! +! Use scaled algorithm +! + u = min( safmax, max( safmin, g1 ) ) + uu = one / u + gs = g*uu + g2 = ABSSQ( gs ) + d = sqrt( g2 ) + s = conjg( gs ) / d + r = d*u + end if + else + f1 = max( abs(real(f)), abs(aimag(f)) ) + g1 = max( abs(real(g)), abs(aimag(g)) ) + if( f1 > rtmin .and. f1 < rtmax .and. & + g1 > rtmin .and. g1 < rtmax ) then +! +! Use unscaled algorithm +! + f2 = ABSSQ( f ) + g2 = ABSSQ( g ) + h2 = f2 + g2 + if( f2 > rtmin .and. h2 < rtmax ) then + d = sqrt( f2*h2 ) + else + d = sqrt( f2 )*sqrt( h2 ) + end if + p = 1 / d + c = f2*p + s = conjg( g )*( f*p ) + r = f*( h2*p ) + else +! +! Use scaled algorithm +! + u = min( safmax, max( safmin, f1, g1 ) ) + uu = one / u + gs = g*uu + g2 = ABSSQ( gs ) + if( f1*uu < rtmin ) then +! +! f is not well-scaled when scaled by g1. +! Use a different scaling for f. +! + v = min( safmax, max( safmin, f1 ) ) + vv = one / v + w = v * uu + fs = f*vv + f2 = ABSSQ( fs ) + h2 = f2*w**2 + g2 + else +! +! Otherwise use the same scaling for f and g. +! + w = one + fs = f*uu + f2 = ABSSQ( fs ) + h2 = f2 + g2 + end if + if( f2 > rtmin .and. h2 < rtmax ) then + d = sqrt( f2*h2 ) + else + d = sqrt( f2 )*sqrt( h2 ) + end if + p = 1 / d + c = ( f2*p )*w + s = conjg( gs )*( fs*p ) + r = ( fs*( h2*p ) )*u + end if + end if + return +end subroutine diff --git a/SRC/dlartg.f b/SRC/dlartg.f deleted file mode 100644 index 453bbe78ad..0000000000 --- a/SRC/dlartg.f +++ /dev/null @@ -1,201 +0,0 @@ -*> \brief \b DLARTG generates a plane rotation with real cosine and real sine. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download DLARTG + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE DLARTG( F, G, CS, SN, R ) -* -* .. Scalar Arguments .. -* DOUBLE PRECISION CS, F, G, R, SN -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> DLARTG generate a plane rotation so that -*> -*> [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. -*> [ -SN CS ] [ G ] [ 0 ] -*> -*> This is a slower, more accurate version of the BLAS1 routine DROTG, -*> with the following other differences: -*> F and G are unchanged on return. -*> If G=0, then CS=1 and SN=0. -*> If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any -*> floating point operations (saves work in DBDSQR when -*> there are zeros on the diagonal). -*> -*> If F exceeds G in magnitude, CS will be positive. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] F -*> \verbatim -*> F is DOUBLE PRECISION -*> The first component of vector to be rotated. -*> \endverbatim -*> -*> \param[in] G -*> \verbatim -*> G is DOUBLE PRECISION -*> The second component of vector to be rotated. -*> \endverbatim -*> -*> \param[out] CS -*> \verbatim -*> CS is DOUBLE PRECISION -*> The cosine of the rotation. -*> \endverbatim -*> -*> \param[out] SN -*> \verbatim -*> SN is DOUBLE PRECISION -*> The sine of the rotation. -*> \endverbatim -*> -*> \param[out] R -*> \verbatim -*> R is DOUBLE PRECISION -*> The nonzero component of the rotated vector. -*> -*> This version has a few statements commented out for thread safety -*> (machine parameters are computed on each entry). 10 feb 03, SJH. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup OTHERauxiliary -* -* ===================================================================== - SUBROUTINE DLARTG( F, G, CS, SN, R ) -* -* -- 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 .. - DOUBLE PRECISION CS, F, G, R, SN -* .. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) - DOUBLE PRECISION ONE - PARAMETER ( ONE = 1.0D0 ) - DOUBLE PRECISION TWO - PARAMETER ( TWO = 2.0D0 ) -* .. -* .. Local Scalars .. -* LOGICAL FIRST - INTEGER COUNT, I - DOUBLE PRECISION EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE -* .. -* .. External Functions .. - DOUBLE PRECISION DLAMCH - EXTERNAL DLAMCH -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, INT, LOG, MAX, SQRT -* .. -* .. Save statement .. -* SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 -* .. -* .. Data statements .. -* DATA FIRST / .TRUE. / -* .. -* .. Executable Statements .. -* -* IF( FIRST ) THEN - SAFMIN = DLAMCH( 'S' ) - EPS = DLAMCH( 'E' ) - SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / - $ LOG( DLAMCH( 'B' ) ) / TWO ) - SAFMX2 = ONE / SAFMN2 -* FIRST = .FALSE. -* END IF - IF( G.EQ.ZERO ) THEN - CS = ONE - SN = ZERO - R = F - ELSE IF( F.EQ.ZERO ) THEN - CS = ZERO - SN = ONE - R = G - ELSE - F1 = F - G1 = G - SCALE = MAX( ABS( F1 ), ABS( G1 ) ) - IF( SCALE.GE.SAFMX2 ) THEN - COUNT = 0 - 10 CONTINUE - COUNT = COUNT + 1 - F1 = F1*SAFMN2 - G1 = G1*SAFMN2 - SCALE = MAX( ABS( F1 ), ABS( G1 ) ) - IF( SCALE.GE.SAFMX2 .AND. COUNT .LT. 20) - $ GO TO 10 - R = SQRT( F1**2+G1**2 ) - CS = F1 / R - SN = G1 / R - DO 20 I = 1, COUNT - R = R*SAFMX2 - 20 CONTINUE - ELSE IF( SCALE.LE.SAFMN2 ) THEN - COUNT = 0 - 30 CONTINUE - COUNT = COUNT + 1 - F1 = F1*SAFMX2 - G1 = G1*SAFMX2 - SCALE = MAX( ABS( F1 ), ABS( G1 ) ) - IF( SCALE.LE.SAFMN2 ) - $ GO TO 30 - R = SQRT( F1**2+G1**2 ) - CS = F1 / R - SN = G1 / R - DO 40 I = 1, COUNT - R = R*SAFMN2 - 40 CONTINUE - ELSE - R = SQRT( F1**2+G1**2 ) - CS = F1 / R - SN = G1 / R - END IF - IF( ABS( F ).GT.ABS( G ) .AND. CS.LT.ZERO ) THEN - CS = -CS - SN = -SN - R = -R - END IF - END IF - RETURN -* -* End of DLARTG -* - END diff --git a/SRC/dlartg.f90 b/SRC/dlartg.f90 new file mode 100644 index 0000000000..03a708f863 --- /dev/null +++ b/SRC/dlartg.f90 @@ -0,0 +1,162 @@ +!> \brief \b DLARTG generates a plane rotation with real cosine and real sine. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! SUBROUTINE DLARTG( F, G, C, S, R ) +! +! .. Scalar Arguments .. +! REAL(wp) C, F, G, R, S +! .. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> DLARTG generates a plane rotation so that +!> +!> [ C S ] . [ F ] = [ R ] +!> [ -S C ] [ G ] [ 0 ] +!> +!> where C**2 + S**2 = 1. +!> +!> The mathematical formulas used for C and S are +!> R = sign(F) * sqrt(F**2 + G**2) +!> C = F / R +!> S = G / R +!> Hence C >= 0. The algorithm used to compute these quantities +!> incorporates scaling to avoid overflow or underflow in computing the +!> square root of the sum of squares. +!> +!> This version is discontinuous in R at F = 0 but it returns the same +!> C and S as ZLARTG for complex inputs (F,0) and (G,0). +!> +!> This is a more accurate version of the BLAS1 routine DROTG, +!> with the following other differences: +!> F and G are unchanged on return. +!> If G=0, then C=1 and S=0. +!> If F=0 and (G .ne. 0), then C=0 and S=sign(1,G) without doing any +!> floating point operations (saves work in DBDSQR when +!> there are zeros on the diagonal). +!> +!> If F exceeds G in magnitude, C will be positive. +!> +!> Below, wp=>dp stands for double precision from LA_CONSTANTS module. +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] F +!> \verbatim +!> F is REAL(wp) +!> The first component of vector to be rotated. +!> \endverbatim +!> +!> \param[in] G +!> \verbatim +!> G is REAL(wp) +!> The second component of vector to be rotated. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is REAL(wp) +!> The cosine of the rotation. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is REAL(wp) +!> The sine of the rotation. +!> \endverbatim +!> +!> \param[out] R +!> \verbatim +!> R is REAL(wp) +!> The nonzero component of the rotated vector. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \date July 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 +!> +!> \endverbatim +! +subroutine DLARTG( f, g, c, s, r ) + use LA_CONSTANTS, & + only: wp=>dp, zero=>dzero, half=>dhalf, one=>done, & + rtmin=>drtmin, rtmax=>drtmax, safmin=>dsafmin, safmax=>dsafmax +! +! -- LAPACK auxiliary routine (version 3.10.0) -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! February 2021 +! +! .. Scalar Arguments .. + real(wp) :: c, f, g, r, s +! .. +! .. Local Scalars .. + real(wp) :: d, f1, fs, g1, gs, p, u, uu +! .. +! .. Intrinsic Functions .. + intrinsic :: abs, sign, sqrt +! .. +! .. Executable Statements .. +! + f1 = abs( f ) + g1 = abs( g ) + if( g == zero ) then + c = one + s = zero + r = f + else if( f == zero ) then + c = zero + s = sign( one, g ) + r = g1 + else if( f1 > rtmin .and. f1 < rtmax .and. & + g1 > rtmin .and. g1 < rtmax ) then + d = sqrt( f*f + g*g ) + p = one / d + c = f1*p + s = g*sign( p, f ) + r = sign( d, f ) + else + u = min( safmax, max( safmin, f1, g1 ) ) + uu = one / u + fs = f*uu + gs = g*uu + d = sqrt( fs*fs + gs*gs ) + p = one / d + c = abs( fs )*p + s = gs*sign( p, f ) + r = sign( d, f )*u + end if + return +end subroutine diff --git a/SRC/slartg.f b/SRC/slartg.f deleted file mode 100644 index 6c23d57cc2..0000000000 --- a/SRC/slartg.f +++ /dev/null @@ -1,201 +0,0 @@ -*> \brief \b SLARTG generates a plane rotation with real cosine and real sine. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download SLARTG + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE SLARTG( F, G, CS, SN, R ) -* -* .. Scalar Arguments .. -* REAL CS, F, G, R, SN -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> SLARTG generate a plane rotation so that -*> -*> [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. -*> [ -SN CS ] [ G ] [ 0 ] -*> -*> This is a slower, more accurate version of the BLAS1 routine SROTG, -*> with the following other differences: -*> F and G are unchanged on return. -*> If G=0, then CS=1 and SN=0. -*> If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any -*> floating point operations (saves work in SBDSQR when -*> there are zeros on the diagonal). -*> -*> If F exceeds G in magnitude, CS will be positive. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] F -*> \verbatim -*> F is REAL -*> The first component of vector to be rotated. -*> \endverbatim -*> -*> \param[in] G -*> \verbatim -*> G is REAL -*> The second component of vector to be rotated. -*> \endverbatim -*> -*> \param[out] CS -*> \verbatim -*> CS is REAL -*> The cosine of the rotation. -*> \endverbatim -*> -*> \param[out] SN -*> \verbatim -*> SN is REAL -*> The sine of the rotation. -*> \endverbatim -*> -*> \param[out] R -*> \verbatim -*> R is REAL -*> The nonzero component of the rotated vector. -*> -*> This version has a few statements commented out for thread safety -*> (machine parameters are computed on each entry). 10 feb 03, SJH. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup OTHERauxiliary -* -* ===================================================================== - SUBROUTINE SLARTG( F, G, CS, SN, R ) -* -* -- 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 .. - REAL CS, F, G, R, SN -* .. -* -* ===================================================================== -* -* .. Parameters .. - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) - REAL ONE - PARAMETER ( ONE = 1.0E0 ) - REAL TWO - PARAMETER ( TWO = 2.0E0 ) -* .. -* .. Local Scalars .. -* LOGICAL FIRST - INTEGER COUNT, I - REAL EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE -* .. -* .. External Functions .. - REAL SLAMCH - EXTERNAL SLAMCH -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, INT, LOG, MAX, SQRT -* .. -* .. Save statement .. -* SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 -* .. -* .. Data statements .. -* DATA FIRST / .TRUE. / -* .. -* .. Executable Statements .. -* -* IF( FIRST ) THEN - SAFMIN = SLAMCH( 'S' ) - EPS = SLAMCH( 'E' ) - SAFMN2 = SLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / - $ LOG( SLAMCH( 'B' ) ) / TWO ) - SAFMX2 = ONE / SAFMN2 -* FIRST = .FALSE. -* END IF - IF( G.EQ.ZERO ) THEN - CS = ONE - SN = ZERO - R = F - ELSE IF( F.EQ.ZERO ) THEN - CS = ZERO - SN = ONE - R = G - ELSE - F1 = F - G1 = G - SCALE = MAX( ABS( F1 ), ABS( G1 ) ) - IF( SCALE.GE.SAFMX2 ) THEN - COUNT = 0 - 10 CONTINUE - COUNT = COUNT + 1 - F1 = F1*SAFMN2 - G1 = G1*SAFMN2 - SCALE = MAX( ABS( F1 ), ABS( G1 ) ) - IF( SCALE.GE.SAFMX2 .AND. COUNT .LT. 20) - $ GO TO 10 - R = SQRT( F1**2+G1**2 ) - CS = F1 / R - SN = G1 / R - DO 20 I = 1, COUNT - R = R*SAFMX2 - 20 CONTINUE - ELSE IF( SCALE.LE.SAFMN2 ) THEN - COUNT = 0 - 30 CONTINUE - COUNT = COUNT + 1 - F1 = F1*SAFMX2 - G1 = G1*SAFMX2 - SCALE = MAX( ABS( F1 ), ABS( G1 ) ) - IF( SCALE.LE.SAFMN2 ) - $ GO TO 30 - R = SQRT( F1**2+G1**2 ) - CS = F1 / R - SN = G1 / R - DO 40 I = 1, COUNT - R = R*SAFMN2 - 40 CONTINUE - ELSE - R = SQRT( F1**2+G1**2 ) - CS = F1 / R - SN = G1 / R - END IF - IF( ABS( F ).GT.ABS( G ) .AND. CS.LT.ZERO ) THEN - CS = -CS - SN = -SN - R = -R - END IF - END IF - RETURN -* -* End of SLARTG -* - END diff --git a/SRC/slartg.f90 b/SRC/slartg.f90 new file mode 100644 index 0000000000..2a936a919f --- /dev/null +++ b/SRC/slartg.f90 @@ -0,0 +1,162 @@ +!> \brief \b SLARTG generates a plane rotation with real cosine and real sine. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! SUBROUTINE SLARTG( F, G, C, S, R ) +! +! .. Scalar Arguments .. +! REAL(wp) C, F, G, R, S +! .. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> SLARTG generates a plane rotation so that +!> +!> [ C S ] . [ F ] = [ R ] +!> [ -S C ] [ G ] [ 0 ] +!> +!> where C**2 + S**2 = 1. +!> +!> The mathematical formulas used for C and S are +!> R = sign(F) * sqrt(F**2 + G**2) +!> C = F / R +!> S = G / R +!> Hence C >= 0. The algorithm used to compute these quantities +!> incorporates scaling to avoid overflow or underflow in computing the +!> square root of the sum of squares. +!> +!> This version is discontinuous in R at F = 0 but it returns the same +!> C and S as SLARTG for complex inputs (F,0) and (G,0). +!> +!> This is a more accurate version of the BLAS1 routine SROTG, +!> with the following other differences: +!> F and G are unchanged on return. +!> If G=0, then C=1 and S=0. +!> If F=0 and (G .ne. 0), then C=0 and S=sign(1,G) without doing any +!> floating point operations (saves work in SBDSQR when +!> there are zeros on the diagonal). +!> +!> If F exceeds G in magnitude, C will be positive. +!> +!> Below, wp=>sp stands for single precision from LA_CONSTANTS module. +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] F +!> \verbatim +!> F is REAL(wp) +!> The first component of vector to be rotated. +!> \endverbatim +!> +!> \param[in] G +!> \verbatim +!> G is REAL(wp) +!> The second component of vector to be rotated. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is REAL(wp) +!> The cosine of the rotation. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is REAL(wp) +!> The sine of the rotation. +!> \endverbatim +!> +!> \param[out] R +!> \verbatim +!> R is REAL(wp) +!> The nonzero component of the rotated vector. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \date July 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 +!> +!> \endverbatim +! +subroutine SLARTG( f, g, c, s, r ) + use LA_CONSTANTS, & + only: wp=>sp, zero=>szero, half=>shalf, one=>sone, & + rtmin=>srtmin, rtmax=>srtmax, safmin=>ssafmin, safmax=>ssafmax +! +! -- LAPACK auxiliary routine (version 3.10.0) -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! February 2021 +! +! .. Scalar Arguments .. + real(wp) :: c, f, g, r, s +! .. +! .. Local Scalars .. + real(wp) :: d, f1, fs, g1, gs, p, u, uu +! .. +! .. Intrinsic Functions .. + intrinsic :: abs, sign, sqrt +! .. +! .. Executable Statements .. +! + f1 = abs( f ) + g1 = abs( g ) + if( g == zero ) then + c = one + s = zero + r = f + else if( f == zero ) then + c = zero + s = sign( one, g ) + r = g1 + else if( f1 > rtmin .and. f1 < rtmax .and. & + g1 > rtmin .and. g1 < rtmax ) then + d = sqrt( f*f + g*g ) + p = one / d + c = f1*p + s = g*sign( p, f ) + r = sign( d, f ) + else + u = min( safmax, max( safmin, f1, g1 ) ) + uu = one / u + fs = f*uu + gs = g*uu + d = sqrt( fs*fs + gs*gs ) + p = one / d + c = abs( fs )*p + s = gs*sign( p, f ) + r = sign( d, f )*u + end if + return +end subroutine diff --git a/SRC/zlartg.f b/SRC/zlartg.f deleted file mode 100644 index b89761aaf5..0000000000 --- a/SRC/zlartg.f +++ /dev/null @@ -1,247 +0,0 @@ -*> \brief \b ZLARTG generates a plane rotation with real cosine and complex sine. -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download ZLARTG + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE ZLARTG( F, G, CS, SN, R ) -* -* .. Scalar Arguments .. -* DOUBLE PRECISION CS -* COMPLEX*16 F, G, R, SN -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> ZLARTG generates a plane rotation so that -*> -*> [ CS SN ] [ F ] [ R ] -*> [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1. -*> [ -SN CS ] [ G ] [ 0 ] -*> -*> This is a faster version of the BLAS1 routine ZROTG, except for -*> the following differences: -*> F and G are unchanged on return. -*> If G=0, then CS=1 and SN=0. -*> If F=0, then CS=0 and SN is chosen so that R is real. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in] F -*> \verbatim -*> F is COMPLEX*16 -*> The first component of vector to be rotated. -*> \endverbatim -*> -*> \param[in] G -*> \verbatim -*> G is COMPLEX*16 -*> The second component of vector to be rotated. -*> \endverbatim -*> -*> \param[out] CS -*> \verbatim -*> CS is DOUBLE PRECISION -*> The cosine of the rotation. -*> \endverbatim -*> -*> \param[out] SN -*> \verbatim -*> SN is COMPLEX*16 -*> The sine of the rotation. -*> \endverbatim -*> -*> \param[out] R -*> \verbatim -*> R is COMPLEX*16 -*> The nonzero component of the rotated vector. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup complex16OTHERauxiliary -* -*> \par Further Details: -* ===================== -*> -*> \verbatim -*> -*> 3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel -*> -*> This version has a few statements commented out for thread safety -*> (machine parameters are computed on each entry). 10 feb 03, SJH. -*> \endverbatim -*> -* ===================================================================== - SUBROUTINE ZLARTG( F, G, CS, SN, R ) -* -* -- 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 .. - DOUBLE PRECISION CS - COMPLEX*16 F, G, R, SN -* .. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION TWO, ONE, ZERO - PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) - COMPLEX*16 CZERO - PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) -* .. -* .. Local Scalars .. -* LOGICAL FIRST - INTEGER COUNT, I - DOUBLE PRECISION D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN, - $ SAFMN2, SAFMX2, SCALE - COMPLEX*16 FF, FS, GS -* .. -* .. External Functions .. - DOUBLE PRECISION DLAMCH, DLAPY2 - LOGICAL DISNAN - EXTERNAL DLAMCH, DLAPY2, DISNAN -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG, - $ MAX, SQRT -* .. -* .. Statement Functions .. - DOUBLE PRECISION ABS1, ABSSQ -* .. -* .. Statement Function definitions .. - ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) ) - ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2 -* .. -* .. Executable Statements .. -* - SAFMIN = DLAMCH( 'S' ) - EPS = DLAMCH( 'E' ) - SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / - $ LOG( DLAMCH( 'B' ) ) / TWO ) - SAFMX2 = ONE / SAFMN2 - SCALE = MAX( ABS1( F ), ABS1( G ) ) - FS = F - GS = G - COUNT = 0 - IF( SCALE.GE.SAFMX2 ) THEN - 10 CONTINUE - COUNT = COUNT + 1 - FS = FS*SAFMN2 - GS = GS*SAFMN2 - SCALE = SCALE*SAFMN2 - IF( SCALE.GE.SAFMX2 .AND. COUNT .LT. 20 ) - $ GO TO 10 - ELSE IF( SCALE.LE.SAFMN2 ) THEN - IF( G.EQ.CZERO.OR.DISNAN( ABS( G ) ) ) THEN - CS = ONE - SN = CZERO - R = F - RETURN - END IF - 20 CONTINUE - COUNT = COUNT - 1 - FS = FS*SAFMX2 - GS = GS*SAFMX2 - SCALE = SCALE*SAFMX2 - IF( SCALE.LE.SAFMN2 ) - $ GO TO 20 - END IF - F2 = ABSSQ( FS ) - G2 = ABSSQ( GS ) - IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN -* -* This is a rare case: F is very small. -* - IF( F.EQ.CZERO ) THEN - CS = ZERO - R = DLAPY2( DBLE( G ), DIMAG( G ) ) -* Do complex/real division explicitly with two real divisions - D = DLAPY2( DBLE( GS ), DIMAG( GS ) ) - SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D ) - RETURN - END IF - F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) ) -* G2 and G2S are accurate -* G2 is at least SAFMIN, and G2S is at least SAFMN2 - G2S = SQRT( G2 ) -* Error in CS from underflow in F2S is at most -* UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS -* If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN, -* and so CS .lt. sqrt(SAFMIN) -* If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN -* and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS) -* Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S - CS = F2S / G2S -* Make sure abs(FF) = 1 -* Do complex/real division explicitly with 2 real divisions - IF( ABS1( F ).GT.ONE ) THEN - D = DLAPY2( DBLE( F ), DIMAG( F ) ) - FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D ) - ELSE - DR = SAFMX2*DBLE( F ) - DI = SAFMX2*DIMAG( F ) - D = DLAPY2( DR, DI ) - FF = DCMPLX( DR / D, DI / D ) - END IF - SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S ) - R = CS*F + SN*G - ELSE -* -* This is the most common case. -* Neither F2 nor F2/G2 are less than SAFMIN -* F2S cannot overflow, and it is accurate -* - F2S = SQRT( ONE+G2 / F2 ) -* Do the F2S(real)*FS(complex) multiply with two real multiplies - R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) ) - CS = ONE / F2S - D = F2 + G2 -* Do complex/real division explicitly with two real divisions - SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D ) - SN = SN*DCONJG( GS ) - IF( COUNT.NE.0 ) THEN - IF( COUNT.GT.0 ) THEN - DO 30 I = 1, COUNT - R = R*SAFMX2 - 30 CONTINUE - ELSE - DO 40 I = 1, -COUNT - R = R*SAFMN2 - 40 CONTINUE - END IF - END IF - END IF - RETURN -* -* End of ZLARTG -* - END diff --git a/SRC/zlartg.f90 b/SRC/zlartg.f90 new file mode 100644 index 0000000000..e509898a1c --- /dev/null +++ b/SRC/zlartg.f90 @@ -0,0 +1,233 @@ +!> \brief \b ZLARTG generates a plane rotation with real cosine and complex sine. +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! SUBROUTINE ZLARTG( F, G, C, S, R ) +! +! .. Scalar Arguments .. +! REAL(wp) C +! COMPLEX(wp) F, G, R, S +! .. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> ZLARTG generates a plane rotation so that +!> +!> [ C S ] . [ F ] = [ R ] +!> [ -conjg(S) C ] [ G ] [ 0 ] +!> +!> where C is real and C**2 + |S|**2 = 1. +!> +!> The mathematical formulas used for C and S are +!> +!> sgn(x) = { x / |x|, x != 0 +!> { 1, x = 0 +!> +!> R = sgn(F) * sqrt(|F|**2 + |G|**2) +!> +!> C = |F| / sqrt(|F|**2 + |G|**2) +!> +!> S = sgn(F) * conjg(G) / sqrt(|F|**2 + |G|**2) +!> +!> When F and G are real, the formulas simplify to C = F/R and +!> S = G/R, and the returned values of C, S, and R should be +!> identical to those returned by DLARTG. +!> +!> The algorithm used to compute these quantities incorporates scaling +!> to avoid overflow or underflow in computing the square root of the +!> sum of squares. +!> +!> This is a faster version of the BLAS1 routine ZROTG, except for +!> the following differences: +!> F and G are unchanged on return. +!> If G=0, then C=1 and S=0. +!> If F=0, then C=0 and S is chosen so that R is real. +!> +!> Below, wp=>dp stands for double precision from LA_CONSTANTS module. +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in] F +!> \verbatim +!> F is COMPLEX(wp) +!> The first component of vector to be rotated. +!> \endverbatim +!> +!> \param[in] G +!> \verbatim +!> G is COMPLEX(wp) +!> The second component of vector to be rotated. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is REAL(wp) +!> The cosine of the rotation. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is COMPLEX(wp) +!> The sine of the rotation. +!> \endverbatim +!> +!> \param[out] R +!> \verbatim +!> R is COMPLEX(wp) +!> The nonzero component of the rotated vector. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \date August 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 +!> +!> \endverbatim +! +subroutine ZLARTG( f, g, c, s, r ) + use LA_CONSTANTS, & + only: wp=>dp, zero=>dzero, one=>done, two=>dtwo, czero=>zzero, & + rtmin=>drtmin, rtmax=>drtmax, safmin=>dsafmin, safmax=>dsafmax +! +! -- LAPACK auxiliary routine (version 3.10.0) -- +! -- LAPACK is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! February 2021 +! +! .. Scalar Arguments .. + real(wp) c + complex(wp) f, g, r, s +! .. +! .. Local Scalars .. + real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w + complex(wp) :: fs, gs, t +! .. +! .. Intrinsic Functions .. + intrinsic :: abs, aimag, conjg, max, min, real, sqrt +! .. +! .. Statement Functions .. + real(wp) :: ABSSQ +! .. +! .. Statement Function definitions .. + ABSSQ( t ) = real( t )**2 + aimag( t )**2 +! .. +! .. Executable Statements .. +! + if( g == czero ) then + c = one + s = czero + r = f + else if( f == czero ) then + c = zero + g1 = max( abs(real(g)), abs(aimag(g)) ) + if( g1 > rtmin .and. g1 < rtmax ) then +! +! Use unscaled algorithm +! + g2 = ABSSQ( g ) + d = sqrt( g2 ) + s = conjg( g ) / d + r = d + else +! +! Use scaled algorithm +! + u = min( safmax, max( safmin, g1 ) ) + uu = one / u + gs = g*uu + g2 = ABSSQ( gs ) + d = sqrt( g2 ) + s = conjg( gs ) / d + r = d*u + end if + else + f1 = max( abs(real(f)), abs(aimag(f)) ) + g1 = max( abs(real(g)), abs(aimag(g)) ) + if( f1 > rtmin .and. f1 < rtmax .and. & + g1 > rtmin .and. g1 < rtmax ) then +! +! Use unscaled algorithm +! + f2 = ABSSQ( f ) + g2 = ABSSQ( g ) + h2 = f2 + g2 + if( f2 > rtmin .and. h2 < rtmax ) then + d = sqrt( f2*h2 ) + else + d = sqrt( f2 )*sqrt( h2 ) + end if + p = 1 / d + c = f2*p + s = conjg( g )*( f*p ) + r = f*( h2*p ) + else +! +! Use scaled algorithm +! + u = min( safmax, max( safmin, f1, g1 ) ) + uu = one / u + gs = g*uu + g2 = ABSSQ( gs ) + if( f1*uu < rtmin ) then +! +! f is not well-scaled when scaled by g1. +! Use a different scaling for f. +! + v = min( safmax, max( safmin, f1 ) ) + vv = one / v + w = v * uu + fs = f*vv + f2 = ABSSQ( fs ) + h2 = f2*w**2 + g2 + else +! +! Otherwise use the same scaling for f and g. +! + w = one + fs = f*uu + f2 = ABSSQ( fs ) + h2 = f2 + g2 + end if + if( f2 > rtmin .and. h2 < rtmax ) then + d = sqrt( f2*h2 ) + else + d = sqrt( f2 )*sqrt( h2 ) + end if + p = 1 / d + c = ( f2*p )*w + s = conjg( gs )*( fs*p ) + r = ( fs*( h2*p ) )*u + end if + end if + return +end subroutine