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