Annotation of rpl/lapack/lapack/dla_gbamv.f, revision 1.5
1.1 bertrand 1: SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
2: $ INCX, BETA, Y, INCY )
3: *
1.5 ! bertrand 4: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 5: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
6: * -- Jason Riedy of Univ. of California Berkeley. --
7: * -- June 2010 --
8: *
9: * -- LAPACK is a software package provided by Univ. of Tennessee, --
10: * -- Univ. of California Berkeley and NAG Ltd. --
11: *
12: IMPLICIT NONE
13: * ..
14: * .. Scalar Arguments ..
15: DOUBLE PRECISION ALPHA, BETA
16: INTEGER INCX, INCY, LDAB, M, N, KL, KU, TRANS
17: * ..
18: * .. Array Arguments ..
19: DOUBLE PRECISION AB( LDAB, * ), X( * ), Y( * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * DLA_GBAMV performs one of the matrix-vector operations
26: *
27: * y := alpha*abs(A)*abs(x) + beta*abs(y),
1.5 ! bertrand 28: * or y := alpha*abs(A)**T*abs(x) + beta*abs(y),
1.1 bertrand 29: *
30: * where alpha and beta are scalars, x and y are vectors and A is an
31: * m by n matrix.
32: *
33: * This function is primarily used in calculating error bounds.
34: * To protect against underflow during evaluation, components in
35: * the resulting vector are perturbed away from zero by (N+1)
36: * times the underflow threshold. To prevent unnecessarily large
37: * errors for block-structure embedded in general matrices,
38: * "symbolically" zero components are not perturbed. A zero
39: * entry is considered "symbolic" if all multiplications involved
40: * in computing that entry have at least one zero multiplicand.
41: *
42: * Arguments
43: * ==========
44: *
45: * TRANS (input) INTEGER
46: * On entry, TRANS specifies the operation to be performed as
47: * follows:
48: *
49: * BLAS_NO_TRANS y := alpha*abs(A)*abs(x) + beta*abs(y)
1.5 ! bertrand 50: * BLAS_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
! 51: * BLAS_CONJ_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
1.1 bertrand 52: *
53: * Unchanged on exit.
54: *
1.5 ! bertrand 55: * M (input) INTEGER
1.1 bertrand 56: * On entry, M specifies the number of rows of the matrix A.
57: * M must be at least zero.
58: * Unchanged on exit.
59: *
1.5 ! bertrand 60: * N (input) INTEGER
1.1 bertrand 61: * On entry, N specifies the number of columns of the matrix A.
62: * N must be at least zero.
63: * Unchanged on exit.
64: *
1.5 ! bertrand 65: * KL (input) INTEGER
1.1 bertrand 66: * The number of subdiagonals within the band of A. KL >= 0.
67: *
1.5 ! bertrand 68: * KU (input) INTEGER
1.1 bertrand 69: * The number of superdiagonals within the band of A. KU >= 0.
70: *
1.5 ! bertrand 71: * ALPHA (input) DOUBLE PRECISION
1.1 bertrand 72: * On entry, ALPHA specifies the scalar alpha.
73: * Unchanged on exit.
74: *
1.5 ! bertrand 75: * AB (input) DOUBLE PRECISION array of DIMENSION ( LDAB, n )
! 76: * Before entry, the leading m by n part of the array AB must
1.1 bertrand 77: * contain the matrix of coefficients.
78: * Unchanged on exit.
79: *
1.5 ! bertrand 80: * LDAB (input) INTEGER
! 81: * On entry, LDA specifies the first dimension of AB as declared
! 82: * in the calling (sub) program. LDAB must be at least
1.1 bertrand 83: * max( 1, m ).
84: * Unchanged on exit.
85: *
1.5 ! bertrand 86: * X (input) DOUBLE PRECISION array, dimension
1.1 bertrand 87: * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
88: * and at least
89: * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
90: * Before entry, the incremented array X must contain the
91: * vector x.
92: * Unchanged on exit.
93: *
1.5 ! bertrand 94: * INCX (input) INTEGER
1.1 bertrand 95: * On entry, INCX specifies the increment for the elements of
96: * X. INCX must not be zero.
97: * Unchanged on exit.
98: *
1.5 ! bertrand 99: * BETA (input) DOUBLE PRECISION
1.1 bertrand 100: * On entry, BETA specifies the scalar beta. When BETA is
101: * supplied as zero then Y need not be set on input.
102: * Unchanged on exit.
103: *
1.5 ! bertrand 104: * Y (input/output) DOUBLE PRECISION array, dimension
1.1 bertrand 105: * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
106: * and at least
107: * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
108: * Before entry with BETA non-zero, the incremented array Y
109: * must contain the vector y. On exit, Y is overwritten by the
110: * updated vector y.
111: *
1.5 ! bertrand 112: * INCY (input) INTEGER
1.1 bertrand 113: * On entry, INCY specifies the increment for the elements of
114: * Y. INCY must not be zero.
115: * Unchanged on exit.
116: *
117: *
118: * Level 2 Blas routine.
119: *
120: * =====================================================================
1.5 ! bertrand 121: *
1.1 bertrand 122: * .. Parameters ..
123: DOUBLE PRECISION ONE, ZERO
124: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
125: * ..
126: * .. Local Scalars ..
127: LOGICAL SYMB_ZERO
128: DOUBLE PRECISION TEMP, SAFE1
129: INTEGER I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE
130: * ..
131: * .. External Subroutines ..
132: EXTERNAL XERBLA, DLAMCH
133: DOUBLE PRECISION DLAMCH
134: * ..
135: * .. External Functions ..
136: EXTERNAL ILATRANS
137: INTEGER ILATRANS
138: * ..
139: * .. Intrinsic Functions ..
140: INTRINSIC MAX, ABS, SIGN
141: * ..
142: * .. Executable Statements ..
143: *
144: * Test the input parameters.
145: *
146: INFO = 0
147: IF ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
148: $ .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
149: $ .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN
150: INFO = 1
151: ELSE IF( M.LT.0 )THEN
152: INFO = 2
153: ELSE IF( N.LT.0 )THEN
154: INFO = 3
155: ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN
156: INFO = 4
157: ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
158: INFO = 5
159: ELSE IF( LDAB.LT.KL+KU+1 )THEN
160: INFO = 6
161: ELSE IF( INCX.EQ.0 )THEN
162: INFO = 8
163: ELSE IF( INCY.EQ.0 )THEN
164: INFO = 11
165: END IF
166: IF( INFO.NE.0 )THEN
167: CALL XERBLA( 'DLA_GBAMV ', INFO )
168: RETURN
169: END IF
170: *
171: * Quick return if possible.
172: *
173: IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
174: $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
175: $ RETURN
176: *
177: * Set LENX and LENY, the lengths of the vectors x and y, and set
178: * up the start points in X and Y.
179: *
180: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
181: LENX = N
182: LENY = M
183: ELSE
184: LENX = M
185: LENY = N
186: END IF
187: IF( INCX.GT.0 )THEN
188: KX = 1
189: ELSE
190: KX = 1 - ( LENX - 1 )*INCX
191: END IF
192: IF( INCY.GT.0 )THEN
193: KY = 1
194: ELSE
195: KY = 1 - ( LENY - 1 )*INCY
196: END IF
197: *
198: * Set SAFE1 essentially to be the underflow threshold times the
199: * number of additions in each row.
200: *
201: SAFE1 = DLAMCH( 'Safe minimum' )
202: SAFE1 = (N+1)*SAFE1
203: *
204: * Form y := alpha*abs(A)*abs(x) + beta*abs(y).
205: *
206: * The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
207: * the inexact flag. Still doesn't help change the iteration order
208: * to per-column.
209: *
210: KD = KU + 1
211: KE = KL + 1
212: IY = KY
213: IF ( INCX.EQ.1 ) THEN
214: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
215: DO I = 1, LENY
216: IF ( BETA .EQ. ZERO ) THEN
217: SYMB_ZERO = .TRUE.
218: Y( IY ) = 0.0D+0
219: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
220: SYMB_ZERO = .TRUE.
221: ELSE
222: SYMB_ZERO = .FALSE.
223: Y( IY ) = BETA * ABS( Y( IY ) )
224: END IF
225: IF ( ALPHA .NE. ZERO ) THEN
226: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
227: TEMP = ABS( AB( KD+I-J, J ) )
228: SYMB_ZERO = SYMB_ZERO .AND.
229: $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
230:
231: Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
232: END DO
233: END IF
234:
235: IF ( .NOT.SYMB_ZERO )
236: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
237: IY = IY + INCY
238: END DO
239: ELSE
240: DO I = 1, LENY
241: IF ( BETA .EQ. ZERO ) THEN
242: SYMB_ZERO = .TRUE.
243: Y( IY ) = 0.0D+0
244: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
245: SYMB_ZERO = .TRUE.
246: ELSE
247: SYMB_ZERO = .FALSE.
248: Y( IY ) = BETA * ABS( Y( IY ) )
249: END IF
250: IF ( ALPHA .NE. ZERO ) THEN
251: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
252: TEMP = ABS( AB( KE-I+J, I ) )
253: SYMB_ZERO = SYMB_ZERO .AND.
254: $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
255:
256: Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
257: END DO
258: END IF
259:
260: IF ( .NOT.SYMB_ZERO )
261: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
262: IY = IY + INCY
263: END DO
264: END IF
265: ELSE
266: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
267: DO I = 1, LENY
268: IF ( BETA .EQ. ZERO ) THEN
269: SYMB_ZERO = .TRUE.
270: Y( IY ) = 0.0D+0
271: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
272: SYMB_ZERO = .TRUE.
273: ELSE
274: SYMB_ZERO = .FALSE.
275: Y( IY ) = BETA * ABS( Y( IY ) )
276: END IF
277: IF ( ALPHA .NE. ZERO ) THEN
278: JX = KX
279: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
280: TEMP = ABS( AB( KD+I-J, J ) )
281: SYMB_ZERO = SYMB_ZERO .AND.
282: $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
283:
284: Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
285: JX = JX + INCX
286: END DO
287: END IF
288:
289: IF ( .NOT.SYMB_ZERO )
290: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
291:
292: IY = IY + INCY
293: END DO
294: ELSE
295: DO I = 1, LENY
296: IF ( BETA .EQ. ZERO ) THEN
297: SYMB_ZERO = .TRUE.
298: Y( IY ) = 0.0D+0
299: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
300: SYMB_ZERO = .TRUE.
301: ELSE
302: SYMB_ZERO = .FALSE.
303: Y( IY ) = BETA * ABS( Y( IY ) )
304: END IF
305: IF ( ALPHA .NE. ZERO ) THEN
306: JX = KX
307: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
308: TEMP = ABS( AB( KE-I+J, I ) )
309: SYMB_ZERO = SYMB_ZERO .AND.
310: $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
311:
312: Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
313: JX = JX + INCX
314: END DO
315: END IF
316:
317: IF ( .NOT.SYMB_ZERO )
318: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
319:
320: IY = IY + INCY
321: END DO
322: END IF
323:
324: END IF
325: *
326: RETURN
327: *
328: * End of DLA_GBAMV
329: *
330: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_gbamv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack