LAPACK 3.3.0
|
00001 SUBROUTINE ZLA_GEAMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, 00002 $ Y, INCY ) 00003 * 00004 * -- LAPACK routine (version 3.2.2) -- 00005 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 00006 * -- Jason Riedy of Univ. of California Berkeley. -- 00007 * -- June 2010 -- 00008 * 00009 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00010 * -- Univ. of California Berkeley and NAG Ltd. -- 00011 * 00012 IMPLICIT NONE 00013 * .. 00014 * .. Scalar Arguments .. 00015 DOUBLE PRECISION ALPHA, BETA 00016 INTEGER INCX, INCY, LDA, M, N 00017 INTEGER TRANS 00018 * .. 00019 * .. Array Arguments .. 00020 COMPLEX*16 A( LDA, * ), X( * ) 00021 DOUBLE PRECISION Y( * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * ZLA_GEAMV performs one of the matrix-vector operations 00028 * 00029 * y := alpha*abs(A)*abs(x) + beta*abs(y), 00030 * or y := alpha*abs(A)'*abs(x) + beta*abs(y), 00031 * 00032 * where alpha and beta are scalars, x and y are vectors and A is an 00033 * m by n matrix. 00034 * 00035 * This function is primarily used in calculating error bounds. 00036 * To protect against underflow during evaluation, components in 00037 * the resulting vector are perturbed away from zero by (N+1) 00038 * times the underflow threshold. To prevent unnecessarily large 00039 * errors for block-structure embedded in general matrices, 00040 * "symbolically" zero components are not perturbed. A zero 00041 * entry is considered "symbolic" if all multiplications involved 00042 * in computing that entry have at least one zero multiplicand. 00043 * 00044 * Arguments 00045 * ========== 00046 * 00047 * TRANS (input) INTEGER 00048 * On entry, TRANS specifies the operation to be performed as 00049 * follows: 00050 * 00051 * BLAS_NO_TRANS y := alpha*abs(A)*abs(x) + beta*abs(y) 00052 * BLAS_TRANS y := alpha*abs(A')*abs(x) + beta*abs(y) 00053 * BLAS_CONJ_TRANS y := alpha*abs(A')*abs(x) + beta*abs(y) 00054 * 00055 * Unchanged on exit. 00056 * 00057 * M (input) INTEGER 00058 * On entry, M specifies the number of rows of the matrix A. 00059 * M must be at least zero. 00060 * Unchanged on exit. 00061 * 00062 * N (input) INTEGER 00063 * On entry, N specifies the number of columns of the matrix A. 00064 * N must be at least zero. 00065 * Unchanged on exit. 00066 * 00067 * ALPHA - DOUBLE PRECISION 00068 * On entry, ALPHA specifies the scalar alpha. 00069 * Unchanged on exit. 00070 * 00071 * A - COMPLEX*16 array of DIMENSION ( LDA, n ) 00072 * Before entry, the leading m by n part of the array A must 00073 * contain the matrix of coefficients. 00074 * Unchanged on exit. 00075 * 00076 * LDA (input) INTEGER 00077 * On entry, LDA specifies the first dimension of A as declared 00078 * in the calling (sub) program. LDA must be at least 00079 * max( 1, m ). 00080 * Unchanged on exit. 00081 * 00082 * X - COMPLEX*16 array of DIMENSION at least 00083 * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' 00084 * and at least 00085 * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. 00086 * Before entry, the incremented array X must contain the 00087 * vector x. 00088 * Unchanged on exit. 00089 * 00090 * INCX (input) INTEGER 00091 * On entry, INCX specifies the increment for the elements of 00092 * X. INCX must not be zero. 00093 * Unchanged on exit. 00094 * 00095 * BETA - DOUBLE PRECISION 00096 * On entry, BETA specifies the scalar beta. When BETA is 00097 * supplied as zero then Y need not be set on input. 00098 * Unchanged on exit. 00099 * 00100 * Y (input/output) DOUBLE PRECISION array, dimension 00101 * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' 00102 * and at least 00103 * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. 00104 * Before entry with BETA non-zero, the incremented array Y 00105 * must contain the vector y. On exit, Y is overwritten by the 00106 * updated vector y. 00107 * 00108 * INCY (input) INTEGER 00109 * On entry, INCY specifies the increment for the elements of 00110 * Y. INCY must not be zero. 00111 * Unchanged on exit. 00112 * 00113 * 00114 * Level 2 Blas routine. 00115 * 00116 * ===================================================================== 00117 * 00118 * .. Parameters .. 00119 COMPLEX*16 ONE, ZERO 00120 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00121 * .. 00122 * .. Local Scalars .. 00123 LOGICAL SYMB_ZERO 00124 DOUBLE PRECISION TEMP, SAFE1 00125 INTEGER I, INFO, IY, J, JX, KX, KY, LENX, LENY 00126 COMPLEX*16 CDUM 00127 * .. 00128 * .. External Subroutines .. 00129 EXTERNAL XERBLA, DLAMCH 00130 DOUBLE PRECISION DLAMCH 00131 * .. 00132 * .. External Functions .. 00133 EXTERNAL ILATRANS 00134 INTEGER ILATRANS 00135 * .. 00136 * .. Intrinsic Functions .. 00137 INTRINSIC MAX, ABS, REAL, DIMAG, SIGN 00138 * .. 00139 * .. Statement Functions .. 00140 DOUBLE PRECISION CABS1 00141 * .. 00142 * .. Statement Function Definitions .. 00143 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 00144 * .. 00145 * .. Executable Statements .. 00146 * 00147 * Test the input parameters. 00148 * 00149 INFO = 0 00150 IF ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) ) 00151 $ .OR. ( TRANS.EQ.ILATRANS( 'T' ) ) 00152 $ .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN 00153 INFO = 1 00154 ELSE IF( M.LT.0 )THEN 00155 INFO = 2 00156 ELSE IF( N.LT.0 )THEN 00157 INFO = 3 00158 ELSE IF( LDA.LT.MAX( 1, M ) )THEN 00159 INFO = 6 00160 ELSE IF( INCX.EQ.0 )THEN 00161 INFO = 8 00162 ELSE IF( INCY.EQ.0 )THEN 00163 INFO = 11 00164 END IF 00165 IF( INFO.NE.0 )THEN 00166 CALL XERBLA( 'ZLA_GEAMV ', INFO ) 00167 RETURN 00168 END IF 00169 * 00170 * Quick return if possible. 00171 * 00172 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. 00173 $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) 00174 $ RETURN 00175 * 00176 * Set LENX and LENY, the lengths of the vectors x and y, and set 00177 * up the start points in X and Y. 00178 * 00179 IF( TRANS.EQ.ILATRANS( 'N' ) )THEN 00180 LENX = N 00181 LENY = M 00182 ELSE 00183 LENX = M 00184 LENY = N 00185 END IF 00186 IF( INCX.GT.0 )THEN 00187 KX = 1 00188 ELSE 00189 KX = 1 - ( LENX - 1 )*INCX 00190 END IF 00191 IF( INCY.GT.0 )THEN 00192 KY = 1 00193 ELSE 00194 KY = 1 - ( LENY - 1 )*INCY 00195 END IF 00196 * 00197 * Set SAFE1 essentially to be the underflow threshold times the 00198 * number of additions in each row. 00199 * 00200 SAFE1 = DLAMCH( 'Safe minimum' ) 00201 SAFE1 = (N+1)*SAFE1 00202 * 00203 * Form y := alpha*abs(A)*abs(x) + beta*abs(y). 00204 * 00205 * The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to 00206 * the inexact flag. Still doesn't help change the iteration order 00207 * to per-column. 00208 * 00209 IY = KY 00210 IF ( INCX.EQ.1 ) THEN 00211 IF( TRANS.EQ.ILATRANS( 'N' ) )THEN 00212 DO I = 1, LENY 00213 IF ( BETA .EQ. 0.0D+0 ) THEN 00214 SYMB_ZERO = .TRUE. 00215 Y( IY ) = 0.0D+0 00216 ELSE IF ( Y( IY ) .EQ. 0.0D+0 ) THEN 00217 SYMB_ZERO = .TRUE. 00218 ELSE 00219 SYMB_ZERO = .FALSE. 00220 Y( IY ) = BETA * ABS( Y( IY ) ) 00221 END IF 00222 IF ( ALPHA .NE. 0.0D+0 ) THEN 00223 DO J = 1, LENX 00224 TEMP = CABS1( A( I, J ) ) 00225 SYMB_ZERO = SYMB_ZERO .AND. 00226 $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00227 00228 Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP 00229 END DO 00230 END IF 00231 00232 IF ( .NOT.SYMB_ZERO ) Y( IY ) = 00233 $ Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00234 00235 IY = IY + INCY 00236 END DO 00237 ELSE 00238 DO I = 1, LENY 00239 IF ( BETA .EQ. 0.0D+0 ) THEN 00240 SYMB_ZERO = .TRUE. 00241 Y( IY ) = 0.0D+0 00242 ELSE IF ( Y( IY ) .EQ. 0.0D+0 ) THEN 00243 SYMB_ZERO = .TRUE. 00244 ELSE 00245 SYMB_ZERO = .FALSE. 00246 Y( IY ) = BETA * ABS( Y( IY ) ) 00247 END IF 00248 IF ( ALPHA .NE. 0.0D+0 ) THEN 00249 DO J = 1, LENX 00250 TEMP = CABS1( A( J, I ) ) 00251 SYMB_ZERO = SYMB_ZERO .AND. 00252 $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00253 00254 Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP 00255 END DO 00256 END IF 00257 00258 IF ( .NOT.SYMB_ZERO ) Y( IY ) = 00259 $ Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00260 00261 IY = IY + INCY 00262 END DO 00263 END IF 00264 ELSE 00265 IF( TRANS.EQ.ILATRANS( 'N' ) )THEN 00266 DO I = 1, LENY 00267 IF ( BETA .EQ. 0.0D+0 ) THEN 00268 SYMB_ZERO = .TRUE. 00269 Y( IY ) = 0.0D+0 00270 ELSE IF ( Y( IY ) .EQ. 0.0D+0 ) THEN 00271 SYMB_ZERO = .TRUE. 00272 ELSE 00273 SYMB_ZERO = .FALSE. 00274 Y( IY ) = BETA * ABS( Y( IY ) ) 00275 END IF 00276 IF ( ALPHA .NE. 0.0D+0 ) THEN 00277 JX = KX 00278 DO J = 1, LENX 00279 TEMP = CABS1( A( I, J ) ) 00280 SYMB_ZERO = SYMB_ZERO .AND. 00281 $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00282 00283 Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP 00284 JX = JX + INCX 00285 END DO 00286 END IF 00287 00288 IF ( .NOT.SYMB_ZERO ) Y( IY ) = 00289 $ Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00290 00291 IY = IY + INCY 00292 END DO 00293 ELSE 00294 DO I = 1, LENY 00295 IF ( BETA .EQ. 0.0D+0 ) THEN 00296 SYMB_ZERO = .TRUE. 00297 Y( IY ) = 0.0D+0 00298 ELSE IF ( Y( IY ) .EQ. 0.0D+0 ) THEN 00299 SYMB_ZERO = .TRUE. 00300 ELSE 00301 SYMB_ZERO = .FALSE. 00302 Y( IY ) = BETA * ABS( Y( IY ) ) 00303 END IF 00304 IF ( ALPHA .NE. 0.0D+0 ) THEN 00305 JX = KX 00306 DO J = 1, LENX 00307 TEMP = CABS1( A( J, I ) ) 00308 SYMB_ZERO = SYMB_ZERO .AND. 00309 $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00310 00311 Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP 00312 JX = JX + INCX 00313 END DO 00314 END IF 00315 00316 IF ( .NOT.SYMB_ZERO ) Y( IY ) = 00317 $ Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00318 00319 IY = IY + INCY 00320 END DO 00321 END IF 00322 00323 END IF 00324 * 00325 RETURN 00326 * 00327 * End of ZLA_GEAMV 00328 * 00329 END