diff --git a/lapack-netlib/SRC/cgetf2.f b/lapack-netlib/SRC/cgetf2.f index aac989970..995ee40ec 100644 --- a/lapack-netlib/SRC/cgetf2.f +++ b/lapack-netlib/SRC/cgetf2.f @@ -101,7 +101,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexGEcomputational +*> \ingroup getf2 * * ===================================================================== SUBROUTINE CGETF2( M, N, A, LDA, IPIV, INFO ) @@ -126,16 +126,14 @@ $ 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 @@ * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN -* -* Compute machine safe minimum -* - SFMIN = SLAMCH('S') * DO 10 J = 1, MIN( M, N ) * @@ -181,15 +175,8 @@ * * 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 - 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/lapack-netlib/SRC/crscl.f b/lapack-netlib/SRC/crscl.f new file mode 100644 index 000000000..22919cd62 --- /dev/null +++ b/lapack-netlib/SRC/crscl.f @@ -0,0 +1,202 @@ +*> \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 SAFMAX, SAFMIN, OV, AR, AI, ABSR, ABSI, UR + % , UI +* .. +* .. External Functions .. + REAL SLAMCH + COMPLEX CLADIV + EXTERNAL SLAMCH, CLADIV +* .. +* .. External Subroutines .. + EXTERNAL CSCAL, CSSCAL, CSRSCL +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.LE.0 ) + $ RETURN +* +* Get machine parameters +* + SAFMIN = SLAMCH( 'S' ) + SAFMAX = ONE / SAFMIN + OV = SLAMCH( 'O' ) +* +* Initialize constants related to A. +* + AR = REAL( A ) + AI = AIMAG( A ) + ABSR = ABS( AR ) + ABSI = ABS( AI ) +* + IF( AI.EQ.ZERO ) THEN +* If alpha is real, then we can use csrscl + CALL CSRSCL( 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 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 + CALL CSCAL( N, CMPLX( 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 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 + 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 + ELSE + CALL CSCAL( N, CMPLX( ONE / UR, -ONE / UI ), X, INCX ) + END IF + END IF +* + RETURN +* +* End of CRSCL +* + END diff --git a/lapack-netlib/SRC/zgetf2.f b/lapack-netlib/SRC/zgetf2.f index c247f8645..7c63dbbee 100644 --- a/lapack-netlib/SRC/zgetf2.f +++ b/lapack-netlib/SRC/zgetf2.f @@ -101,7 +101,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16GEcomputational +*> \ingroup getf2 * * ===================================================================== SUBROUTINE ZGETF2( M, N, A, LDA, IPIV, INFO ) @@ -127,7 +127,7 @@ * .. * .. Local Scalars .. DOUBLE PRECISION SFMIN - INTEGER I, J, JP + INTEGER J, JP * .. * .. External Functions .. DOUBLE PRECISION DLAMCH @@ -135,7 +135,7 @@ EXTERNAL DLAMCH, IZAMAX * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZGERU, ZSCAL, ZSWAP + EXTERNAL XERBLA, ZGERU, ZRSCL, ZSWAP * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -181,15 +181,8 @@ * * 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/lapack-netlib/SRC/zrscl.f b/lapack-netlib/SRC/zrscl.f new file mode 100644 index 000000000..970f6de75 --- /dev/null +++ b/lapack-netlib/SRC/zrscl.f @@ -0,0 +1,203 @@ +*> \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*16 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