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