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