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