diff --git a/BLAS/SRC/CMakeLists.txt b/BLAS/SRC/CMakeLists.txt index 0078dca400..da780b9885 100644 --- a/BLAS/SRC/CMakeLists.txt +++ b/BLAS/SRC/CMakeLists.txt @@ -29,16 +29,16 @@ # Level 1 BLAS #--------------------------------------------------------- set(SBLAS1 isamax.f sasum.f saxpy.f scopy.f sdot.f snrm2.f - srot.f srotg.f sscal.f sswap.f sdsdot.f srotmg.f srotm.f) + srot.f srotg.f90 sscal.f sswap.f sdsdot.f srotmg.f srotm.f) set(CBLAS1 scabs1.f scasum.f scnrm2.f icamax.f caxpy.f ccopy.f - cdotc.f cdotu.f csscal.f crotg.f cscal.f cswap.f csrot.f) + cdotc.f cdotu.f csscal.f crotg.f90 cscal.f cswap.f csrot.f) set(DBLAS1 idamax.f dasum.f daxpy.f dcopy.f ddot.f dnrm2.f - drot.f drotg.f dscal.f dsdot.f dswap.f drotmg.f drotm.f) + drot.f drotg.f90 dscal.f dsdot.f dswap.f drotmg.f drotm.f) set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f izamax.f zaxpy.f zcopy.f - zdotc.f zdotu.f zdscal.f zrotg.f zscal.f zswap.f zdrot.f) + zdotc.f zdotu.f zdscal.f zrotg.f90 zscal.f zswap.f zdrot.f) set(CB1AUX isamax.f sasum.f saxpy.f scopy.f snrm2.f sscal.f) diff --git a/BLAS/SRC/Makefile b/BLAS/SRC/Makefile index 66bb964218..9b0c8907d8 100644 --- a/BLAS/SRC/Makefile +++ b/BLAS/SRC/Makefile @@ -56,6 +56,10 @@ TOPSRCDIR = ../.. include $(TOPSRCDIR)/make.inc +.SUFFIXES: .f90 .o +.f90.o: + $(FC) $(FFLAGS) -c -o $@ $< + .PHONY: all all: $(BLASLIB) diff --git a/BLAS/SRC/crotg.f b/BLAS/SRC/crotg.f deleted file mode 100644 index 01be542770..0000000000 --- a/BLAS/SRC/crotg.f +++ /dev/null @@ -1,94 +0,0 @@ -*> \brief \b CROTG -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -* Definition: -* =========== -* -* SUBROUTINE CROTG(CA,CB,C,S) -* -* .. Scalar Arguments .. -* COMPLEX CA,CB,S -* REAL C -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> CROTG determines a complex Givens rotation. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in,out] CA -*> \verbatim -*> CA is COMPLEX -*> \endverbatim -*> -*> \param[in] CB -*> \verbatim -*> CB is COMPLEX -*> \endverbatim -*> -*> \param[out] C -*> \verbatim -*> C is REAL -*> \endverbatim -*> -*> \param[out] S -*> \verbatim -*> S is COMPLEX -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup complex_blas_level1 -* -* ===================================================================== - SUBROUTINE CROTG(CA,CB,C,S) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX CA,CB,S - REAL C -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - COMPLEX ALPHA - REAL NORM,SCALE -* .. -* .. Intrinsic Functions .. - INTRINSIC CABS,CONJG,SQRT -* .. - IF (CABS(CA).EQ.0.) THEN - C = 0. - S = (1.,0.) - CA = CB - ELSE - SCALE = CABS(CA) + CABS(CB) - NORM = SCALE*SQRT((CABS(CA/SCALE))**2+ (CABS(CB/SCALE))**2) - ALPHA = CA/CABS(CA) - C = CABS(CA)/NORM - S = ALPHA*CONJG(CB)/NORM - CA = ALPHA*NORM - END IF - RETURN - END diff --git a/BLAS/SRC/crotg.f90 b/BLAS/SRC/crotg.f90 new file mode 100644 index 0000000000..8cb1eb495e --- /dev/null +++ b/BLAS/SRC/crotg.f90 @@ -0,0 +1,229 @@ +!> \brief \b CROTG +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! CROTG constructs a plane rotation +! [ c s ] [ a ] = [ r ] +! [ -conjg(s) c ] [ b ] [ 0 ] +! where c is real, s ic complex, and c**2 + conjg(s)*s = 1. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> The computation uses the formulas +!> |x| = sqrt( Re(x)**2 + Im(x)**2 ) +!> sgn(x) = x / |x| if x /= 0 +!> = 1 if x = 0 +!> c = |a| / sqrt(|a|**2 + |b|**2) +!> s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2) +!> When a and b are real and r /= 0, the formulas simplify to +!> r = sgn(a)*sqrt(|a|**2 + |b|**2) +!> c = a / r +!> s = b / r +!> the same as in CROTG when |a| > |b|. When |b| >= |a|, the +!> sign of c and s will be different from those computed by CROTG +!> if the signs of a and b are not the same. +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in,out] A +!> \verbatim +!> A is COMPLEX +!> On entry, the scalar a. +!> On exit, the scalar r. +!> \endverbatim +!> +!> \param[in] B +!> \verbatim +!> B is COMPLEX +!> The scalar b. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is REAL +!> The scalar c. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is REAL +!> The scalar s. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +! +!> \ingroup single_blas_level1 +! +!> \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 CROTG( a, b, c, s ) + integer, parameter :: wp = kind(1.e0) +! +! -- Reference BLAS level1 routine -- +! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! +! .. Constants .. + real(wp), parameter :: zero = 0.0_wp + real(wp), parameter :: one = 1.0_wp + complex(wp), parameter :: czero = 0.0_wp +! .. +! .. Scaling constants .. + real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( & + minexponent(0._wp)-1, & + 1-maxexponent(0._wp) & + ) + real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( & + 1-minexponent(0._wp), & + maxexponent(0._wp)-1 & + ) + real(wp), parameter :: rtmin = sqrt( real(radix(0._wp),wp)**max( & + minexponent(0._wp)-1, & + 1-maxexponent(0._wp) & + ) / epsilon(0._wp) ) + real(wp), parameter :: rtmax = sqrt( real(radix(0._wp),wp)**max( & + 1-minexponent(0._wp), & + maxexponent(0._wp)-1 & + ) * epsilon(0._wp) ) +! .. +! .. Scalar Arguments .. + real(wp) :: c + complex(wp) :: a, b, s +! .. +! .. Local Scalars .. + real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w + complex(wp) :: f, fs, g, gs, r, 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 .. +! + f = a + g = b + 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 + a = r + return +end subroutine diff --git a/BLAS/SRC/drotg.f b/BLAS/SRC/drotg.f deleted file mode 100644 index bcaec27f5e..0000000000 --- a/BLAS/SRC/drotg.f +++ /dev/null @@ -1,106 +0,0 @@ -*> \brief \b DROTG -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -* Definition: -* =========== -* -* SUBROUTINE DROTG(DA,DB,C,S) -* -* .. Scalar Arguments .. -* DOUBLE PRECISION C,DA,DB,S -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> DROTG construct givens plane rotation. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in,out] DA -*> \verbatim -*> DA is DOUBLE PRECISION -*> \endverbatim -*> -*> \param[in,out] DB -*> \verbatim -*> DB is DOUBLE PRECISION -*> \endverbatim -*> -*> \param[out] C -*> \verbatim -*> C is DOUBLE PRECISION -*> \endverbatim -*> -*> \param[out] S -*> \verbatim -*> S is DOUBLE PRECISION -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup double_blas_level1 -* -*> \par Further Details: -* ===================== -*> -*> \verbatim -*> -*> jack dongarra, linpack, 3/11/78. -*> \endverbatim -*> -* ===================================================================== - SUBROUTINE DROTG(DA,DB,C,S) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - DOUBLE PRECISION C,DA,DB,S -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - DOUBLE PRECISION R,ROE,SCALE,Z -* .. -* .. Intrinsic Functions .. - INTRINSIC DABS,DSIGN,DSQRT -* .. - SCALE = DABS(DA) + DABS(DB) - IF (SCALE.EQ.0.0d0) THEN - C = 1.0d0 - S = 0.0d0 - R = 0.0d0 - Z = 0.0d0 - ELSE - ROE = DB - IF (DABS(DA).GT.DABS(DB)) ROE = DA - R = SCALE*DSQRT((DA/SCALE)**2+ (DB/SCALE)**2) - R = DSIGN(1.0d0,ROE)*R - C = DA/R - S = DB/R - Z = 1.0d0 - IF (DABS(DA).GT.DABS(DB)) Z = S - IF (DABS(DB).GE.DABS(DA) .AND. C.NE.0.0d0) Z = 1.0d0/C - END IF - DA = R - DB = Z - RETURN - END diff --git a/BLAS/SRC/drotg.f90 b/BLAS/SRC/drotg.f90 new file mode 100644 index 0000000000..fd309c0af9 --- /dev/null +++ b/BLAS/SRC/drotg.f90 @@ -0,0 +1,151 @@ +!> \brief \b DROTG +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! DROTG constructs a plane rotation +! [ c s ] [ a ] = [ r ] +! [ -s c ] [ b ] [ 0 ] +! satisfying c**2 + s**2 = 1. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> The computation uses the formulas +!> sigma = sgn(a) if |a| > |b| +!> = sgn(b) if |b| >= |a| +!> r = sigma*sqrt( a**2 + b**2 ) +!> c = 1; s = 0 if r = 0 +!> c = a/r; s = b/r if r != 0 +!> The subroutine also computes +!> z = s if |a| > |b|, +!> = 1/c if |b| >= |a| and c != 0 +!> = 1 if c = 0 +!> This allows c and s to be reconstructed from z as follows: +!> If z = 1, set c = 0, s = 1. +!> If |z| < 1, set c = sqrt(1 - z**2) and s = z. +!> If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2). +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in,out] A +!> \verbatim +!> A is DOUBLE PRECISION +!> On entry, the scalar a. +!> On exit, the scalar r. +!> \endverbatim +!> +!> \param[in,out] B +!> \verbatim +!> B is DOUBLE PRECISION +!> On entry, the scalar b. +!> On exit, the scalar z. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is DOUBLE PRECISION +!> The scalar c. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is DOUBLE PRECISION +!> The scalar s. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +! +!> \ingroup single_blas_level1 +! +!> \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 DROTG( a, b, c, s ) + integer, parameter :: wp = kind(1.d0) +! +! -- Reference BLAS level1 routine -- +! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! +! .. Constants .. + real(wp), parameter :: zero = 0.0_wp + real(wp), parameter :: one = 1.0_wp +! .. +! .. Scaling constants .. + real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( & + minexponent(0._wp)-1, & + 1-maxexponent(0._wp) & + ) + real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( & + 1-minexponent(0._wp), & + maxexponent(0._wp)-1 & + ) +! .. +! .. Scalar Arguments .. + real(wp) :: a, b, c, s +! .. +! .. Local Scalars .. + real(wp) :: anorm, bnorm, scl, sigma, r, z +! .. + anorm = abs(a) + bnorm = abs(b) + if( anorm == zero ) then + c = zero + s = one + a = b + b = one + else if( bnorm == zero ) then + c = one + s = zero + b = zero + else + scl = min( safmax, max( safmin, anorm, bnorm ) ) + if( anorm > bnorm ) then + sigma = sign(one,a) + else + sigma = sign(one,b) + end if + r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) ) + c = a/r + s = b/r + if( anorm > bnorm ) then + z = s + else if( c /= zero ) then + z = one/c + else + z = one + end if + a = r + b = z + end if + return +end subroutine diff --git a/BLAS/SRC/srotg.f b/BLAS/SRC/srotg.f deleted file mode 100644 index bd101edf25..0000000000 --- a/BLAS/SRC/srotg.f +++ /dev/null @@ -1,106 +0,0 @@ -*> \brief \b SROTG -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -* Definition: -* =========== -* -* SUBROUTINE SROTG(SA,SB,C,S) -* -* .. Scalar Arguments .. -* REAL C,S,SA,SB -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> SROTG construct givens plane rotation. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in,out] SA -*> \verbatim -*> SA is REAL -*> \endverbatim -*> -*> \param[in,out] SB -*> \verbatim -*> SB is REAL -*> \endverbatim -*> -*> \param[out] C -*> \verbatim -*> C is REAL -*> \endverbatim -*> -*> \param[out] S -*> \verbatim -*> S is REAL -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup single_blas_level1 -* -*> \par Further Details: -* ===================== -*> -*> \verbatim -*> -*> jack dongarra, linpack, 3/11/78. -*> \endverbatim -*> -* ===================================================================== - SUBROUTINE SROTG(SA,SB,C,S) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - REAL C,S,SA,SB -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - REAL R,ROE,SCALE,Z -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS,SIGN,SQRT -* .. - SCALE = ABS(SA) + ABS(SB) - IF (SCALE.EQ.0.0) THEN - C = 1.0 - S = 0.0 - R = 0.0 - Z = 0.0 - ELSE - ROE = SB - IF (ABS(SA).GT.ABS(SB)) ROE = SA - R = SCALE*SQRT((SA/SCALE)**2+ (SB/SCALE)**2) - R = SIGN(1.0,ROE)*R - C = SA/R - S = SB/R - Z = 1.0 - IF (ABS(SA).GT.ABS(SB)) Z = S - IF (ABS(SB).GE.ABS(SA) .AND. C.NE.0.0) Z = 1.0/C - END IF - SA = R - SB = Z - RETURN - END diff --git a/BLAS/SRC/srotg.f90 b/BLAS/SRC/srotg.f90 new file mode 100644 index 0000000000..7ca8f41779 --- /dev/null +++ b/BLAS/SRC/srotg.f90 @@ -0,0 +1,151 @@ +!> \brief \b SROTG +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! SROTG constructs a plane rotation +! [ c s ] [ a ] = [ r ] +! [ -s c ] [ b ] [ 0 ] +! satisfying c**2 + s**2 = 1. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> The computation uses the formulas +!> sigma = sgn(a) if |a| > |b| +!> = sgn(b) if |b| >= |a| +!> r = sigma*sqrt( a**2 + b**2 ) +!> c = 1; s = 0 if r = 0 +!> c = a/r; s = b/r if r != 0 +!> The subroutine also computes +!> z = s if |a| > |b|, +!> = 1/c if |b| >= |a| and c != 0 +!> = 1 if c = 0 +!> This allows c and s to be reconstructed from z as follows: +!> If z = 1, set c = 0, s = 1. +!> If |z| < 1, set c = sqrt(1 - z**2) and s = z. +!> If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2). +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in,out] A +!> \verbatim +!> A is REAL +!> On entry, the scalar a. +!> On exit, the scalar r. +!> \endverbatim +!> +!> \param[in,out] B +!> \verbatim +!> B is REAL +!> On entry, the scalar b. +!> On exit, the scalar z. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is REAL +!> The scalar c. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is REAL +!> The scalar s. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +! +!> \ingroup single_blas_level1 +! +!> \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 SROTG( a, b, c, s ) + integer, parameter :: wp = kind(1.e0) +! +! -- Reference BLAS level1 routine -- +! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! +! .. Constants .. + real(wp), parameter :: zero = 0.0_wp + real(wp), parameter :: one = 1.0_wp +! .. +! .. Scaling constants .. + real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( & + minexponent(0._wp)-1, & + 1-maxexponent(0._wp) & + ) + real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( & + 1-minexponent(0._wp), & + maxexponent(0._wp)-1 & + ) +! .. +! .. Scalar Arguments .. + real(wp) :: a, b, c, s +! .. +! .. Local Scalars .. + real(wp) :: anorm, bnorm, scl, sigma, r, z +! .. + anorm = abs(a) + bnorm = abs(b) + if( anorm == zero ) then + c = zero + s = one + a = b + b = one + else if( bnorm == zero ) then + c = one + s = zero + b = zero + else + scl = min( safmax, max( safmin, anorm, bnorm ) ) + if( anorm > bnorm ) then + sigma = sign(one,a) + else + sigma = sign(one,b) + end if + r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) ) + c = a/r + s = b/r + if( anorm > bnorm ) then + z = s + else if( c /= zero ) then + z = one/c + else + z = one + end if + a = r + b = z + end if + return +end subroutine diff --git a/BLAS/SRC/zrotg.f b/BLAS/SRC/zrotg.f deleted file mode 100644 index de30974409..0000000000 --- a/BLAS/SRC/zrotg.f +++ /dev/null @@ -1,95 +0,0 @@ -*> \brief \b ZROTG -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -* Definition: -* =========== -* -* SUBROUTINE ZROTG(CA,CB,C,S) -* -* .. Scalar Arguments .. -* COMPLEX*16 CA,CB,S -* DOUBLE PRECISION C -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> ZROTG determines a double complex Givens rotation. -*> \endverbatim -* -* Arguments: -* ========== -* -*> \param[in,out] CA -*> \verbatim -*> CA is COMPLEX*16 -*> \endverbatim -*> -*> \param[in] CB -*> \verbatim -*> CB is COMPLEX*16 -*> \endverbatim -*> -*> \param[out] C -*> \verbatim -*> C is DOUBLE PRECISION -*> \endverbatim -*> -*> \param[out] S -*> \verbatim -*> S is COMPLEX*16 -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup complex16_blas_level1 -* -* ===================================================================== - SUBROUTINE ZROTG(CA,CB,C,S) -* -* -- Reference BLAS level1 routine -- -* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - COMPLEX*16 CA,CB,S - DOUBLE PRECISION C -* .. -* -* ===================================================================== -* -* .. Local Scalars .. - COMPLEX*16 ALPHA - DOUBLE PRECISION NORM,SCALE -* .. -* .. Intrinsic Functions .. - INTRINSIC CDABS,DCMPLX,DCONJG,DSQRT -* .. - IF (CDABS(CA).EQ.0.0d0) THEN - C = 0.0d0 - S = (1.0d0,0.0d0) - CA = CB - ELSE - SCALE = CDABS(CA) + CDABS(CB) - NORM = SCALE*DSQRT((CDABS(CA/DCMPLX(SCALE,0.0d0)))**2+ - $ (CDABS(CB/DCMPLX(SCALE,0.0d0)))**2) - ALPHA = CA/CDABS(CA) - C = CDABS(CA)/NORM - S = ALPHA*DCONJG(CB)/NORM - CA = ALPHA*NORM - END IF - RETURN - END diff --git a/BLAS/SRC/zrotg.f90 b/BLAS/SRC/zrotg.f90 new file mode 100644 index 0000000000..b23f8d1353 --- /dev/null +++ b/BLAS/SRC/zrotg.f90 @@ -0,0 +1,229 @@ +!> \brief \b ZROTG +! +! =========== DOCUMENTATION =========== +! +! Online html documentation available at +! http://www.netlib.org/lapack/explore-html/ +! +! Definition: +! =========== +! +! ZROTG constructs a plane rotation +! [ c s ] [ a ] = [ r ] +! [ -conjg(s) c ] [ b ] [ 0 ] +! where c is real, s ic complex, and c**2 + conjg(s)*s = 1. +! +!> \par Purpose: +! ============= +!> +!> \verbatim +!> +!> The computation uses the formulas +!> |x| = sqrt( Re(x)**2 + Im(x)**2 ) +!> sgn(x) = x / |x| if x /= 0 +!> = 1 if x = 0 +!> c = |a| / sqrt(|a|**2 + |b|**2) +!> s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2) +!> When a and b are real and r /= 0, the formulas simplify to +!> r = sgn(a)*sqrt(|a|**2 + |b|**2) +!> c = a / r +!> s = b / r +!> the same as in ZROTG when |a| > |b|. When |b| >= |a|, the +!> sign of c and s will be different from those computed by ZROTG +!> if the signs of a and b are not the same. +!> +!> \endverbatim +! +! Arguments: +! ========== +! +!> \param[in,out] A +!> \verbatim +!> A is DOUBLE COMPLEX +!> On entry, the scalar a. +!> On exit, the scalar r. +!> \endverbatim +!> +!> \param[in] B +!> \verbatim +!> B is DOUBLE COMPLEX +!> The scalar b. +!> \endverbatim +!> +!> \param[out] C +!> \verbatim +!> C is DOUBLE PRECISION +!> The scalar c. +!> \endverbatim +!> +!> \param[out] S +!> \verbatim +!> S is DOUBLE PRECISION +!> The scalar s. +!> \endverbatim +! +! Authors: +! ======== +! +!> \author Edward Anderson, Lockheed Martin +! +!> \par Contributors: +! ================== +!> +!> Weslley Pereira, University of Colorado Denver, USA +! +!> \ingroup single_blas_level1 +! +!> \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 ZROTG( a, b, c, s ) + integer, parameter :: wp = kind(1.d0) +! +! -- Reference BLAS level1 routine -- +! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +! +! .. Constants .. + real(wp), parameter :: zero = 0.0_wp + real(wp), parameter :: one = 1.0_wp + complex(wp), parameter :: czero = 0.0_wp +! .. +! .. Scaling constants .. + real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( & + minexponent(0._wp)-1, & + 1-maxexponent(0._wp) & + ) + real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( & + 1-minexponent(0._wp), & + maxexponent(0._wp)-1 & + ) + real(wp), parameter :: rtmin = sqrt( real(radix(0._wp),wp)**max( & + minexponent(0._wp)-1, & + 1-maxexponent(0._wp) & + ) / epsilon(0._wp) ) + real(wp), parameter :: rtmax = sqrt( real(radix(0._wp),wp)**max( & + 1-minexponent(0._wp), & + maxexponent(0._wp)-1 & + ) * epsilon(0._wp) ) +! .. +! .. Scalar Arguments .. + real(wp) :: c + complex(wp) :: a, b, s +! .. +! .. Local Scalars .. + real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w + complex(wp) :: f, fs, g, gs, r, 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 .. +! + f = a + g = b + 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 + a = r + return +end subroutine diff --git a/CMakeLists.txt b/CMakeLists.txt index 73ab55f3d2..b36fec1540 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,9 +16,6 @@ set(CMAKE_MODULE_PATH "${LAPACK_SOURCE_DIR}/CMAKE" ${CMAKE_MODULE_PATH}) # Export all symbols on Windows when building shared libraries SET(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS TRUE) -# Tell CMake that our Fortran sources are written in fixed format. -set(CMAKE_Fortran_FORMAT FIXED) - # Set a default build type if none was specified if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) message(STATUS "Setting build type to 'Release' as none was specified.")