|
| 1 | +*> \brief \b CRSCL multiplies a vector by the reciprocal of a real scalar. |
| 2 | +* |
| 3 | +* =========== DOCUMENTATION =========== |
| 4 | +* |
| 5 | +* Online html documentation available at |
| 6 | +* http://www.netlib.org/lapack/explore-html/ |
| 7 | +* |
| 8 | +*> \htmlonly |
| 9 | +*> Download CRSCL + dependencies |
| 10 | +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/crscl.f"> |
| 11 | +*> [TGZ]</a> |
| 12 | +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/crscl.f"> |
| 13 | +*> [ZIP]</a> |
| 14 | +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/crscl.f"> |
| 15 | +*> [TXT]</a> |
| 16 | +*> \endhtmlonly |
| 17 | +* |
| 18 | +* Definition: |
| 19 | +* =========== |
| 20 | +* |
| 21 | +* SUBROUTINE CRSCL( N, A, X, INCX ) |
| 22 | +* |
| 23 | +* .. Scalar Arguments .. |
| 24 | +* INTEGER INCX, N |
| 25 | +* COMPLEX A |
| 26 | +* .. |
| 27 | +* .. Array Arguments .. |
| 28 | +* COMPLEX X( * ) |
| 29 | +* .. |
| 30 | +* |
| 31 | +* |
| 32 | +*> \par Purpose: |
| 33 | +* ============= |
| 34 | +*> |
| 35 | +*> \verbatim |
| 36 | +*> |
| 37 | +*> CRSCL multiplies an n-element complex vector x by the complex scalar |
| 38 | +*> 1/a. This is done without overflow or underflow as long as |
| 39 | +*> the final result x/a does not overflow or underflow. |
| 40 | +*> \endverbatim |
| 41 | +* |
| 42 | +* Arguments: |
| 43 | +* ========== |
| 44 | +* |
| 45 | +*> \param[in] N |
| 46 | +*> \verbatim |
| 47 | +*> N is INTEGER |
| 48 | +*> The number of components of the vector x. |
| 49 | +*> \endverbatim |
| 50 | +*> |
| 51 | +*> \param[in] A |
| 52 | +*> \verbatim |
| 53 | +*> A is COMPLEX |
| 54 | +*> The scalar a which is used to divide each component of x. |
| 55 | +*> A must not be 0, or the subroutine will divide by zero. |
| 56 | +*> \endverbatim |
| 57 | +*> |
| 58 | +*> \param[in,out] X |
| 59 | +*> \verbatim |
| 60 | +*> X is COMPLEX array, dimension |
| 61 | +*> (1+(N-1)*abs(INCX)) |
| 62 | +*> The n-element vector x. |
| 63 | +*> \endverbatim |
| 64 | +*> |
| 65 | +*> \param[in] INCX |
| 66 | +*> \verbatim |
| 67 | +*> INCX is INTEGER |
| 68 | +*> The increment between successive values of the vector X. |
| 69 | +*> > 0: X(1) = X(1) and X(1+(i-1)*INCX) = x(i), 1< i<= n |
| 70 | +*> \endverbatim |
| 71 | +* |
| 72 | +* Authors: |
| 73 | +* ======== |
| 74 | +* |
| 75 | +*> \author Univ. of Tennessee |
| 76 | +*> \author Univ. of California Berkeley |
| 77 | +*> \author Univ. of Colorado Denver |
| 78 | +*> \author NAG Ltd. |
| 79 | +* |
| 80 | +*> \ingroup complexOTHERauxiliary |
| 81 | +* |
| 82 | +* ===================================================================== |
| 83 | + SUBROUTINE CRSCL( N, A, X, INCX ) |
| 84 | +* |
| 85 | +* -- LAPACK auxiliary routine -- |
| 86 | +* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
| 87 | +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
| 88 | +* |
| 89 | +* .. Scalar Arguments .. |
| 90 | + INTEGER INCX, N |
| 91 | + COMPLEX A |
| 92 | +* .. |
| 93 | +* .. Array Arguments .. |
| 94 | + COMPLEX X( * ) |
| 95 | +* .. |
| 96 | +* |
| 97 | +* ===================================================================== |
| 98 | +* |
| 99 | +* .. Parameters .. |
| 100 | + REAL ZERO, ONE |
| 101 | + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) |
| 102 | +* .. |
| 103 | +* .. Local Scalars .. |
| 104 | + REAL SAFMAX, SAFMIN, OV, AR, AI, ABSR, ABSI, UR |
| 105 | + % , UI |
| 106 | +* .. |
| 107 | +* .. External Functions .. |
| 108 | + REAL SLAMCH |
| 109 | + COMPLEX CLADIV |
| 110 | + EXTERNAL SLAMCH, CLADIV |
| 111 | +* .. |
| 112 | +* .. External Subroutines .. |
| 113 | + EXTERNAL CSCAL, CSSCAL, CSRSCL |
| 114 | +* .. |
| 115 | +* .. Intrinsic Functions .. |
| 116 | + INTRINSIC ABS |
| 117 | +* .. |
| 118 | +* .. Executable Statements .. |
| 119 | +* |
| 120 | +* Quick return if possible |
| 121 | +* |
| 122 | + IF( N.LE.0 ) |
| 123 | + $ RETURN |
| 124 | +* |
| 125 | +* Get machine parameters |
| 126 | +* |
| 127 | + SAFMIN = SLAMCH( 'S' ) |
| 128 | + SAFMAX = ONE / SAFMIN |
| 129 | + OV = SLAMCH( 'O' ) |
| 130 | +* |
| 131 | +* Initialize constants related to A. |
| 132 | +* |
| 133 | + AR = REAL( A ) |
| 134 | + AI = AIMAG( A ) |
| 135 | + ABSR = ABS( AR ) |
| 136 | + ABSI = ABS( AI ) |
| 137 | +* |
| 138 | + IF( AI.EQ.ZERO ) THEN |
| 139 | +* If alpha is real, then we can use csrscl |
| 140 | + CALL CSRSCL( N, AR, X, INCX ) |
| 141 | +* |
| 142 | + ELSE IF( AR.EQ.ZERO ) THEN |
| 143 | +* If alpha has a zero real part, then we follow the same rules as if |
| 144 | +* alpha were real. |
| 145 | + IF( ABSI.GT.SAFMAX ) THEN |
| 146 | + CALL CSSCAL( N, SAFMIN, X, INCX ) |
| 147 | + CALL CSCAL( N, CMPLX( ZERO, -SAFMAX / AI ), X, INCX ) |
| 148 | + ELSE IF( ABSI.LT.SAFMIN ) THEN |
| 149 | + CALL CSCAL( N, CMPLX( ZERO, -SAFMIN / AI ), X, INCX ) |
| 150 | + CALL CSSCAL( N, SAFMAX, X, INCX ) |
| 151 | + ELSE |
| 152 | + CALL CSCAL( N, CMPLX( ZERO, -ONE / AI ), X, INCX ) |
| 153 | + END IF |
| 154 | +* |
| 155 | + ELSE |
| 156 | +* The following numbers can be computed. |
| 157 | +* They are the inverse of the real and imaginary parts of 1/alpha. |
| 158 | +* Note that a and b are always different from zero. |
| 159 | +* NaNs are only possible if either: |
| 160 | +* 1. alphaR or alphaI is NaN. |
| 161 | +* 2. alphaR and alphaI are both infinite, in which case it makes sense |
| 162 | +* to propagate a NaN. |
| 163 | + UR = AR + AI * ( AI / AR ) |
| 164 | + UI = AI + AR * ( AR / AI ) |
| 165 | +* |
| 166 | + IF( (ABS( UR ).LT.SAFMIN).OR.(ABS( UI ).LT.SAFMIN) ) THEN |
| 167 | +* This means that both alphaR and alphaI are very small. |
| 168 | + CALL CSCAL( N, CMPLX( SAFMIN / UR, -SAFMIN / UI ), X, INCX ) |
| 169 | + CALL CSSCAL( N, SAFMAX, X, INCX ) |
| 170 | + ELSE IF( (ABS( UR ).GT.SAFMAX).OR.(ABS( UI ).GT.SAFMAX) ) THEN |
| 171 | + IF( (ABSR.GT.OV).OR.(ABSI.GT.OV) ) THEN |
| 172 | +* This means that a and b are both Inf. No need for scaling. |
| 173 | + CALL CSCAL( N, CMPLX( ONE / UR, -ONE / UI ), X, INCX ) |
| 174 | + ELSE |
| 175 | + CALL CSSCAL( N, SAFMIN, X, INCX ) |
| 176 | + IF( (ABS( UR ).GT.OV).OR.(ABS( UI ).GT.OV) ) THEN |
| 177 | +* Infs were generated. We do proper scaling to avoid them. |
| 178 | + IF( ABSR.GE.ABSI ) THEN |
| 179 | +* ABS( UR ) <= ABS( UI ) |
| 180 | + UR = (SAFMIN * AR) + SAFMIN * (AI * ( AI / AR )) |
| 181 | + UI = (SAFMIN * AI) + AR * ( (SAFMIN * AR) / AI ) |
| 182 | + ELSE |
| 183 | +* ABS( UR ) > ABS( UI ) |
| 184 | + UR = (SAFMIN * AR) + AI * ( (SAFMIN * AI) / AR ) |
| 185 | + UI = (SAFMIN * AI) + SAFMIN * (AR * ( AR / AI )) |
| 186 | + END IF |
| 187 | + CALL CSCAL( N, CMPLX( ONE / UR, -ONE / UI ), X, INCX ) |
| 188 | + ELSE |
| 189 | + CALL CSCAL( N, CMPLX( SAFMAX / UR, -SAFMAX / UI ), |
| 190 | + $ X, INCX ) |
| 191 | + END IF |
| 192 | + END IF |
| 193 | + ELSE |
| 194 | + CALL CSCAL( N, CMPLX( ONE / UR, -ONE / UI ), X, INCX ) |
| 195 | + END IF |
| 196 | + END IF |
| 197 | +* |
| 198 | + RETURN |
| 199 | +* |
| 200 | +* End of CRSCL |
| 201 | +* |
| 202 | + END |
0 commit comments