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