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