1: *> \brief \b DLA_GBAMV performs a matrix-vector operation to calculate error bounds.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLA_GBAMV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gbamv.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gbamv.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gbamv.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
22: * INCX, BETA, Y, INCY )
23: *
24: * .. Scalar Arguments ..
25: * DOUBLE PRECISION ALPHA, BETA
26: * INTEGER INCX, INCY, LDAB, M, N, KL, KU, TRANS
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION AB( LDAB, * ), X( * ), Y( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLA_GBAMV performs one of the matrix-vector operations
39: *>
40: *> y := alpha*abs(A)*abs(x) + beta*abs(y),
41: *> or y := alpha*abs(A)**T*abs(x) + beta*abs(y),
42: *>
43: *> where alpha and beta are scalars, x and y are vectors and A is an
44: *> m by n matrix.
45: *>
46: *> This function is primarily used in calculating error bounds.
47: *> To protect against underflow during evaluation, components in
48: *> the resulting vector are perturbed away from zero by (N+1)
49: *> times the underflow threshold. To prevent unnecessarily large
50: *> errors for block-structure embedded in general matrices,
51: *> "symbolically" zero components are not perturbed. A zero
52: *> entry is considered "symbolic" if all multiplications involved
53: *> in computing that entry have at least one zero multiplicand.
54: *> \endverbatim
55: *
56: * Arguments:
57: * ==========
58: *
59: *> \param[in] TRANS
60: *> \verbatim
61: *> TRANS is INTEGER
62: *> On entry, TRANS specifies the operation to be performed as
63: *> follows:
64: *>
65: *> BLAS_NO_TRANS y := alpha*abs(A)*abs(x) + beta*abs(y)
66: *> BLAS_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
67: *> BLAS_CONJ_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
68: *>
69: *> Unchanged on exit.
70: *> \endverbatim
71: *>
72: *> \param[in] M
73: *> \verbatim
74: *> M is INTEGER
75: *> On entry, M specifies the number of rows of the matrix A.
76: *> M must be at least zero.
77: *> Unchanged on exit.
78: *> \endverbatim
79: *>
80: *> \param[in] N
81: *> \verbatim
82: *> N is INTEGER
83: *> On entry, N specifies the number of columns of the matrix A.
84: *> N must be at least zero.
85: *> Unchanged on exit.
86: *> \endverbatim
87: *>
88: *> \param[in] KL
89: *> \verbatim
90: *> KL is INTEGER
91: *> The number of subdiagonals within the band of A. KL >= 0.
92: *> \endverbatim
93: *>
94: *> \param[in] KU
95: *> \verbatim
96: *> KU is INTEGER
97: *> The number of superdiagonals within the band of A. KU >= 0.
98: *> \endverbatim
99: *>
100: *> \param[in] ALPHA
101: *> \verbatim
102: *> ALPHA is DOUBLE PRECISION
103: *> On entry, ALPHA specifies the scalar alpha.
104: *> Unchanged on exit.
105: *> \endverbatim
106: *>
107: *> \param[in] AB
108: *> \verbatim
109: *> AB is DOUBLE PRECISION array, dimension ( LDAB, n )
110: *> Before entry, the leading m by n part of the array AB must
111: *> contain the matrix of coefficients.
112: *> Unchanged on exit.
113: *> \endverbatim
114: *>
115: *> \param[in] LDAB
116: *> \verbatim
117: *> LDAB is INTEGER
118: *> On entry, LDA specifies the first dimension of AB as declared
119: *> in the calling (sub) program. LDAB must be at least
120: *> max( 1, m ).
121: *> Unchanged on exit.
122: *> \endverbatim
123: *>
124: *> \param[in] X
125: *> \verbatim
126: *> X is DOUBLE PRECISION array, dimension
127: *> ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
128: *> and at least
129: *> ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
130: *> Before entry, the incremented array X must contain the
131: *> vector x.
132: *> Unchanged on exit.
133: *> \endverbatim
134: *>
135: *> \param[in] INCX
136: *> \verbatim
137: *> INCX is INTEGER
138: *> On entry, INCX specifies the increment for the elements of
139: *> X. INCX must not be zero.
140: *> Unchanged on exit.
141: *> \endverbatim
142: *>
143: *> \param[in] BETA
144: *> \verbatim
145: *> BETA is DOUBLE PRECISION
146: *> On entry, BETA specifies the scalar beta. When BETA is
147: *> supplied as zero then Y need not be set on input.
148: *> Unchanged on exit.
149: *> \endverbatim
150: *>
151: *> \param[in,out] Y
152: *> \verbatim
153: *> Y is DOUBLE PRECISION array, dimension
154: *> ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
155: *> and at least
156: *> ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
157: *> Before entry with BETA non-zero, the incremented array Y
158: *> must contain the vector y. On exit, Y is overwritten by the
159: *> updated vector y.
160: *> \endverbatim
161: *>
162: *> \param[in] INCY
163: *> \verbatim
164: *> INCY is INTEGER
165: *> On entry, INCY specifies the increment for the elements of
166: *> Y. INCY must not be zero.
167: *> Unchanged on exit.
168: *>
169: *> Level 2 Blas routine.
170: *> \endverbatim
171: *
172: * Authors:
173: * ========
174: *
175: *> \author Univ. of Tennessee
176: *> \author Univ. of California Berkeley
177: *> \author Univ. of Colorado Denver
178: *> \author NAG Ltd.
179: *
180: *> \ingroup doubleGBcomputational
181: *
182: * =====================================================================
183: SUBROUTINE DLA_GBAMV( TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X,
184: $ INCX, BETA, Y, INCY )
185: *
186: * -- LAPACK computational routine --
187: * -- LAPACK is a software package provided by Univ. of Tennessee, --
188: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189: *
190: * .. Scalar Arguments ..
191: DOUBLE PRECISION ALPHA, BETA
192: INTEGER INCX, INCY, LDAB, M, N, KL, KU, TRANS
193: * ..
194: * .. Array Arguments ..
195: DOUBLE PRECISION AB( LDAB, * ), X( * ), Y( * )
196: * ..
197: *
198: * =====================================================================
199: *
200: * .. Parameters ..
201: DOUBLE PRECISION ONE, ZERO
202: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
203: * ..
204: * .. Local Scalars ..
205: LOGICAL SYMB_ZERO
206: DOUBLE PRECISION TEMP, SAFE1
207: INTEGER I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE
208: * ..
209: * .. External Subroutines ..
210: EXTERNAL XERBLA, DLAMCH
211: DOUBLE PRECISION DLAMCH
212: * ..
213: * .. External Functions ..
214: EXTERNAL ILATRANS
215: INTEGER ILATRANS
216: * ..
217: * .. Intrinsic Functions ..
218: INTRINSIC MAX, ABS, SIGN
219: * ..
220: * .. Executable Statements ..
221: *
222: * Test the input parameters.
223: *
224: INFO = 0
225: IF ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
226: $ .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
227: $ .OR. ( TRANS.EQ.ILATRANS( 'C' ) ) ) ) THEN
228: INFO = 1
229: ELSE IF( M.LT.0 )THEN
230: INFO = 2
231: ELSE IF( N.LT.0 )THEN
232: INFO = 3
233: ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN
234: INFO = 4
235: ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
236: INFO = 5
237: ELSE IF( LDAB.LT.KL+KU+1 )THEN
238: INFO = 6
239: ELSE IF( INCX.EQ.0 )THEN
240: INFO = 8
241: ELSE IF( INCY.EQ.0 )THEN
242: INFO = 11
243: END IF
244: IF( INFO.NE.0 )THEN
245: CALL XERBLA( 'DLA_GBAMV ', INFO )
246: RETURN
247: END IF
248: *
249: * Quick return if possible.
250: *
251: IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
252: $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
253: $ RETURN
254: *
255: * Set LENX and LENY, the lengths of the vectors x and y, and set
256: * up the start points in X and Y.
257: *
258: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
259: LENX = N
260: LENY = M
261: ELSE
262: LENX = M
263: LENY = N
264: END IF
265: IF( INCX.GT.0 )THEN
266: KX = 1
267: ELSE
268: KX = 1 - ( LENX - 1 )*INCX
269: END IF
270: IF( INCY.GT.0 )THEN
271: KY = 1
272: ELSE
273: KY = 1 - ( LENY - 1 )*INCY
274: END IF
275: *
276: * Set SAFE1 essentially to be the underflow threshold times the
277: * number of additions in each row.
278: *
279: SAFE1 = DLAMCH( 'Safe minimum' )
280: SAFE1 = (N+1)*SAFE1
281: *
282: * Form y := alpha*abs(A)*abs(x) + beta*abs(y).
283: *
284: * The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
285: * the inexact flag. Still doesn't help change the iteration order
286: * to per-column.
287: *
288: KD = KU + 1
289: KE = KL + 1
290: IY = KY
291: IF ( INCX.EQ.1 ) THEN
292: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
293: DO I = 1, LENY
294: IF ( BETA .EQ. ZERO ) THEN
295: SYMB_ZERO = .TRUE.
296: Y( IY ) = 0.0D+0
297: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
298: SYMB_ZERO = .TRUE.
299: ELSE
300: SYMB_ZERO = .FALSE.
301: Y( IY ) = BETA * ABS( Y( IY ) )
302: END IF
303: IF ( ALPHA .NE. ZERO ) THEN
304: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
305: TEMP = ABS( AB( KD+I-J, J ) )
306: SYMB_ZERO = SYMB_ZERO .AND.
307: $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
308:
309: Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
310: END DO
311: END IF
312:
313: IF ( .NOT.SYMB_ZERO )
314: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
315: IY = IY + INCY
316: END DO
317: ELSE
318: DO I = 1, LENY
319: IF ( BETA .EQ. ZERO ) THEN
320: SYMB_ZERO = .TRUE.
321: Y( IY ) = 0.0D+0
322: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
323: SYMB_ZERO = .TRUE.
324: ELSE
325: SYMB_ZERO = .FALSE.
326: Y( IY ) = BETA * ABS( Y( IY ) )
327: END IF
328: IF ( ALPHA .NE. ZERO ) THEN
329: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
330: TEMP = ABS( AB( KE-I+J, I ) )
331: SYMB_ZERO = SYMB_ZERO .AND.
332: $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
333:
334: Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
335: END DO
336: END IF
337:
338: IF ( .NOT.SYMB_ZERO )
339: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
340: IY = IY + INCY
341: END DO
342: END IF
343: ELSE
344: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
345: DO I = 1, LENY
346: IF ( BETA .EQ. ZERO ) THEN
347: SYMB_ZERO = .TRUE.
348: Y( IY ) = 0.0D+0
349: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
350: SYMB_ZERO = .TRUE.
351: ELSE
352: SYMB_ZERO = .FALSE.
353: Y( IY ) = BETA * ABS( Y( IY ) )
354: END IF
355: IF ( ALPHA .NE. ZERO ) THEN
356: JX = KX
357: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
358: TEMP = ABS( AB( KD+I-J, J ) )
359: SYMB_ZERO = SYMB_ZERO .AND.
360: $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
361:
362: Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
363: JX = JX + INCX
364: END DO
365: END IF
366:
367: IF ( .NOT.SYMB_ZERO )
368: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
369:
370: IY = IY + INCY
371: END DO
372: ELSE
373: DO I = 1, LENY
374: IF ( BETA .EQ. ZERO ) THEN
375: SYMB_ZERO = .TRUE.
376: Y( IY ) = 0.0D+0
377: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
378: SYMB_ZERO = .TRUE.
379: ELSE
380: SYMB_ZERO = .FALSE.
381: Y( IY ) = BETA * ABS( Y( IY ) )
382: END IF
383: IF ( ALPHA .NE. ZERO ) THEN
384: JX = KX
385: DO J = MAX( I-KL, 1 ), MIN( I+KU, LENX )
386: TEMP = ABS( AB( KE-I+J, I ) )
387: SYMB_ZERO = SYMB_ZERO .AND.
388: $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
389:
390: Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
391: JX = JX + INCX
392: END DO
393: END IF
394:
395: IF ( .NOT.SYMB_ZERO )
396: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
397:
398: IY = IY + INCY
399: END DO
400: END IF
401:
402: END IF
403: *
404: RETURN
405: *
406: * End of DLA_GBAMV
407: *
408: 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