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