floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #4126 from martin-frbg/lapack839 

Add C/ZRSCL for reciprocal scaling of a complex vector (Reference-LAPACK PR 839)
This commit is contained in:
Martin Kroeker 2023-07-10 01:49:33 +02:00 committed by GitHub
parent 351645b8af e3277e134c
commit 22ad23abb1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 420 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
cmake/lapack.cmake
View File

@ -187,7 +187,7 @@ set(CLASRC
cposv.f cposvx.f cpotrf2.f cpotri.f cpstrf.f cpstf2.f
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 csprfs.f cspsv.f
crot.f crscl.f cspcon.f csprfs.f cspsv.f
cspsvx.f csptrf.f csptri.f csptrs.f csrscl.f cstedc.f
cstegr.f cstein.f csteqr.f csycon.f
csyrfs.f csysv.f csysvx.f csytf2.f csytrf.f csytri.f
@ -381,7 +381,7 @@ set(ZLASRC
zposv.f zposvx.f zpotrf2.f zpotri.f zpotrs.f zpstrf.f zpstf2.f
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 zsprfs.f zspsv.f
zrot.f zrscl.f zspcon.f zsprfs.f zspsv.f
zspsvx.f zsptrf.f zsptri.f zsptrs.f zdrscl.f zstedc.f
zstegr.f zstein.f zsteqr.f zsycon.f
zsyrfs.f zsysv.f zsysvx.f zsytf2.f zsytrf.f zsytri.f

4
lapack-netlib/SRC/Makefile
View File

@ -280,7 +280,7 @@ CLASRC_O = \
cposv.o cposvx.o cpotf2.o cpotri.o cpstrf.o cpstf2.o \
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 \
crot.o crscl.o cspcon.o cspmv.o cspr.o csprfs.o cspsv.o \
cspsvx.o csptrf.o csptri.o csptrs.o csrscl.o cstedc.o \
cstegr.o cstein.o csteqr.o \
csycon.o csymv.o \
@ -488,7 +488,7 @@ ZLASRC_O = \
zposv.o zposvx.o zpotf2.o zpotrf.o zpotri.o zpotrs.o zpstrf.o zpstf2.o \
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 \
zrot.o zrscl.o zspcon.o zspmv.o zspr.o zsprfs.o zspsv.o \
zspsvx.o zsptrf.o zsptri.o zsptrs.o zdrscl.o zstedc.o \
zstegr.o zstein.o zsteqr.o \
zsycon.o zsymv.o \

25
lapack-netlib/SRC/cgetf2.f
View File

@ -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
*

202
lapack-netlib/SRC/crscl.f Normal file
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/crscl.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/crscl.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/crscl.f">
*> [TXT]</a>
*> \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

17
lapack-netlib/SRC/zgetf2.f
View File

@ -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
*

203
lapack-netlib/SRC/zrscl.f Normal file
View File

@ -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
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zdrscl.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zdrscl.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zdrscl.f">
*> [TXT]</a>
*> \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