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