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