LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X, 00002 $ INCX, BETA, Y, INCY ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 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 REAL ALPHA, BETA 00016 INTEGER INCX, INCY, LDAB, M, N, KL, KU, TRANS 00017 * .. 00018 * .. Array Arguments .. 00019 REAL AB( LDAB, * ), X( * ), Y( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * SLA_GBAMV performs one of the matrix-vector operations 00026 * 00027 * y := alpha*abs(A)*abs(x) + beta*abs(y), 00028 * or y := alpha*abs(A)**T*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**T)*abs(x) + beta*abs(y) 00051 * BLAS_CONJ_TRANS y := alpha*abs(A**T)*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 * KL (input) INTEGER 00066 * The number of subdiagonals within the band of A. KL >= 0. 00067 * 00068 * KU (input) INTEGER 00069 * The number of superdiagonals within the band of A. KU >= 0. 00070 * 00071 * ALPHA (input) REAL 00072 * On entry, ALPHA specifies the scalar alpha. 00073 * Unchanged on exit. 00074 * 00075 * AB (input) REAL array of DIMENSION ( LDAB, n ) 00076 * Before entry, the leading m by n part of the array AB must 00077 * contain the matrix of coefficients. 00078 * Unchanged on exit. 00079 * 00080 * LDAB (input) INTEGER 00081 * On entry, LDA specifies the first dimension of AB as declared 00082 * in the calling (sub) program. LDAB must be at least 00083 * max( 1, m ). 00084 * Unchanged on exit. 00085 * 00086 * X (input) REAL array, dimension 00087 * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' 00088 * and at least 00089 * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. 00090 * Before entry, the incremented array X must contain the 00091 * vector x. 00092 * Unchanged on exit. 00093 * 00094 * INCX (input) INTEGER 00095 * On entry, INCX specifies the increment for the elements of 00096 * X. INCX must not be zero. 00097 * Unchanged on exit. 00098 * 00099 * BETA (input) REAL 00100 * On entry, BETA specifies the scalar beta. When BETA is 00101 * supplied as zero then Y need not be set on input. 00102 * Unchanged on exit. 00103 * 00104 * Y (input/output) REAL array, dimension 00105 * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' 00106 * and at least 00107 * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. 00108 * Before entry with BETA non-zero, the incremented array Y 00109 * must contain the vector y. On exit, Y is overwritten by the 00110 * updated vector y. 00111 * 00112 * INCY (input) INTEGER 00113 * On entry, INCY specifies the increment for the elements of 00114 * Y. INCY must not be zero. 00115 * Unchanged on exit. 00116 * 00117 * 00118 * Level 2 Blas routine. 00119 * 00120 * ===================================================================== 00121 00122 * .. Parameters .. 00123 REAL ONE, ZERO 00124 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00125 * .. 00126 * .. Local Scalars .. 00127 LOGICAL SYMB_ZERO 00128 REAL TEMP, SAFE1 00129 INTEGER I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE 00130 * .. 00131 * .. External Subroutines .. 00132 EXTERNAL XERBLA, SLAMCH 00133 REAL SLAMCH 00134 * .. 00135 * .. External Functions .. 00136 EXTERNAL ILATRANS 00137 INTEGER ILATRANS 00138 * .. 00139 * .. Intrinsic Functions .. 00140 INTRINSIC MAX, ABS, SIGN 00141 * .. 00142 * .. Executable Statements .. 00143 * 00144 * Test the input parameters. 00145 * 00146 INFO = 0 00147 IF ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) ) 00148 $ .OR. ( TRANS.EQ.ILATRANS( 'T' ) ) 00149 $ .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN 00150 INFO = 1 00151 ELSE IF( M.LT.0 )THEN 00152 INFO = 2 00153 ELSE IF( N.LT.0 )THEN 00154 INFO = 3 00155 ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN 00156 INFO = 4 00157 ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN 00158 INFO = 5 00159 ELSE IF( LDAB.LT.KL+KU+1 )THEN 00160 INFO = 6 00161 ELSE IF( INCX.EQ.0 )THEN 00162 INFO = 8 00163 ELSE IF( INCY.EQ.0 )THEN 00164 INFO = 11 00165 END IF 00166 IF( INFO.NE.0 )THEN 00167 CALL XERBLA( 'SLA_GBAMV ', INFO ) 00168 RETURN 00169 END IF 00170 * 00171 * Quick return if possible. 00172 * 00173 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. 00174 $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) 00175 $ RETURN 00176 * 00177 * Set LENX and LENY, the lengths of the vectors x and y, and set 00178 * up the start points in X and Y. 00179 * 00180 IF( TRANS.EQ.ILATRANS( 'N' ) )THEN 00181 LENX = N 00182 LENY = M 00183 ELSE 00184 LENX = M 00185 LENY = N 00186 END IF 00187 IF( INCX.GT.0 )THEN 00188 KX = 1 00189 ELSE 00190 KX = 1 - ( LENX - 1 )*INCX 00191 END IF 00192 IF( INCY.GT.0 )THEN 00193 KY = 1 00194 ELSE 00195 KY = 1 - ( LENY - 1 )*INCY 00196 END IF 00197 * 00198 * Set SAFE1 essentially to be the underflow threshold times the 00199 * number of additions in each row. 00200 * 00201 SAFE1 = SLAMCH( 'Safe minimum' ) 00202 SAFE1 = (N+1)*SAFE1 00203 * 00204 * Form y := alpha*abs(A)*abs(x) + beta*abs(y). 00205 * 00206 * The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to 00207 * the inexact flag. Still doesn't help change the iteration order 00208 * to per-column. 00209 * 00210 KD = KU + 1 00211 KE = KL + 1 00212 IY = KY 00213 IF ( INCX.EQ.1 ) THEN 00214 IF( TRANS.EQ.ILATRANS( 'N' ) )THEN 00215 DO I = 1, LENY 00216 IF ( BETA .EQ. ZERO ) THEN 00217 SYMB_ZERO = .TRUE. 00218 Y( IY ) = 0.0 00219 ELSE IF ( Y( IY ) .EQ. ZERO ) THEN 00220 SYMB_ZERO = .TRUE. 00221 ELSE 00222 SYMB_ZERO = .FALSE. 00223 Y( IY ) = BETA * ABS( Y( IY ) ) 00224 END IF 00225 IF ( ALPHA .NE. ZERO ) THEN 00226 DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX ) 00227 TEMP = ABS( AB( KD+I-J, J ) ) 00228 SYMB_ZERO = SYMB_ZERO .AND. 00229 $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00230 00231 Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP 00232 END DO 00233 END IF 00234 00235 IF ( .NOT.SYMB_ZERO ) 00236 $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00237 IY = IY + INCY 00238 END DO 00239 ELSE 00240 DO I = 1, LENY 00241 IF ( BETA .EQ. ZERO ) THEN 00242 SYMB_ZERO = .TRUE. 00243 Y( IY ) = 0.0 00244 ELSE IF ( Y( IY ) .EQ. ZERO ) THEN 00245 SYMB_ZERO = .TRUE. 00246 ELSE 00247 SYMB_ZERO = .FALSE. 00248 Y( IY ) = BETA * ABS( Y( IY ) ) 00249 END IF 00250 IF ( ALPHA .NE. ZERO ) THEN 00251 DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX ) 00252 TEMP = ABS( AB( KE-I+J, I ) ) 00253 SYMB_ZERO = SYMB_ZERO .AND. 00254 $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00255 00256 Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP 00257 END DO 00258 END IF 00259 00260 IF ( .NOT.SYMB_ZERO ) 00261 $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00262 IY = IY + INCY 00263 END DO 00264 END IF 00265 ELSE 00266 IF( TRANS.EQ.ILATRANS( 'N' ) )THEN 00267 DO I = 1, LENY 00268 IF ( BETA .EQ. ZERO ) THEN 00269 SYMB_ZERO = .TRUE. 00270 Y( IY ) = 0.0 00271 ELSE IF ( Y( IY ) .EQ. ZERO ) THEN 00272 SYMB_ZERO = .TRUE. 00273 ELSE 00274 SYMB_ZERO = .FALSE. 00275 Y( IY ) = BETA * ABS( Y( IY ) ) 00276 END IF 00277 IF ( ALPHA .NE. ZERO ) THEN 00278 JX = KX 00279 DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX ) 00280 TEMP = ABS( AB( KD+I-J, J ) ) 00281 SYMB_ZERO = SYMB_ZERO .AND. 00282 $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00283 00284 Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP 00285 JX = JX + INCX 00286 END DO 00287 END IF 00288 00289 IF ( .NOT.SYMB_ZERO ) 00290 $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00291 00292 IY = IY + INCY 00293 END DO 00294 ELSE 00295 DO I = 1, LENY 00296 IF ( BETA .EQ. ZERO ) THEN 00297 SYMB_ZERO = .TRUE. 00298 Y( IY ) = 0.0 00299 ELSE IF ( Y( IY ) .EQ. ZERO ) THEN 00300 SYMB_ZERO = .TRUE. 00301 ELSE 00302 SYMB_ZERO = .FALSE. 00303 Y( IY ) = BETA * ABS( Y( IY ) ) 00304 END IF 00305 IF ( ALPHA .NE. ZERO ) THEN 00306 JX = KX 00307 DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX ) 00308 TEMP = ABS( AB( KE-I+J, I ) ) 00309 SYMB_ZERO = SYMB_ZERO .AND. 00310 $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO ) 00311 00312 Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP 00313 JX = JX + INCX 00314 END DO 00315 END IF 00316 00317 IF ( .NOT.SYMB_ZERO ) 00318 $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) ) 00319 00320 IY = IY + INCY 00321 END DO 00322 END IF 00323 00324 END IF 00325 * 00326 RETURN 00327 * 00328 * End of SLA_GBAMV 00329 * 00330 END