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