Annotation of rpl/lapack/lapack/dla_geamv.f, revision 1.5
1.1 bertrand 1: SUBROUTINE DLA_GEAMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA,
2: $ Y, INCY )
3: *
1.5 ! bertrand 4: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 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, TRANS
17: * ..
18: * .. Array Arguments ..
19: DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * DLA_GEAMV performs one of the matrix-vector operations
26: *
27: * y := alpha*abs(A)*abs(x) + beta*abs(y),
1.5 ! bertrand 28: * or y := alpha*abs(A)**T*abs(x) + beta*abs(y),
1.1 bertrand 29: *
30: * where alpha and beta are scalars, x and y are vectors and A is an
31: * m by n 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: * TRANS (input) INTEGER
46: * On entry, TRANS specifies the operation to be performed as
47: * follows:
48: *
49: * BLAS_NO_TRANS y := alpha*abs(A)*abs(x) + beta*abs(y)
1.5 ! bertrand 50: * BLAS_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
! 51: * BLAS_CONJ_TRANS y := alpha*abs(A**T)*abs(x) + beta*abs(y)
1.1 bertrand 52: *
53: * Unchanged on exit.
54: *
1.5 ! bertrand 55: * M (input) INTEGER
1.1 bertrand 56: * On entry, M specifies the number of rows of the matrix A.
57: * M must be at least zero.
58: * Unchanged on exit.
59: *
1.5 ! bertrand 60: * N (input) INTEGER
1.1 bertrand 61: * On entry, N specifies the number of columns of the matrix A.
62: * N must be at least zero.
63: * Unchanged on exit.
64: *
1.5 ! bertrand 65: * ALPHA (input) DOUBLE PRECISION
1.1 bertrand 66: * On entry, ALPHA specifies the scalar alpha.
67: * Unchanged on exit.
68: *
1.5 ! bertrand 69: * A (input) DOUBLE PRECISION array of DIMENSION ( LDA, n )
1.1 bertrand 70: * Before entry, the leading m by n part of the array A must
71: * contain the matrix of coefficients.
72: * Unchanged on exit.
73: *
1.5 ! bertrand 74: * LDA (input) INTEGER
1.1 bertrand 75: * On entry, LDA specifies the first dimension of A as declared
76: * in the calling (sub) program. LDA must be at least
77: * max( 1, m ).
78: * Unchanged on exit.
79: *
1.5 ! bertrand 80: * X (input) DOUBLE PRECISION array, dimension
1.1 bertrand 81: * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
82: * and at least
83: * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
84: * Before entry, the incremented array X must contain the
85: * vector x.
86: * Unchanged on exit.
87: *
1.5 ! bertrand 88: * INCX (input) INTEGER
1.1 bertrand 89: * On entry, INCX specifies the increment for the elements of
90: * X. INCX must not be zero.
91: * Unchanged on exit.
92: *
1.5 ! bertrand 93: * BETA (input) DOUBLE PRECISION
1.1 bertrand 94: * On entry, BETA specifies the scalar beta. When BETA is
95: * supplied as zero then Y need not be set on input.
96: * Unchanged on exit.
97: *
1.5 ! bertrand 98: * Y (input/output) DOUBLE PRECISION
1.1 bertrand 99: * Array of DIMENSION at least
100: * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
101: * and at least
102: * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
103: * Before entry with BETA non-zero, the incremented array Y
104: * must contain the vector y. On exit, Y is overwritten by the
105: * updated vector y.
106: *
1.5 ! bertrand 107: * INCY (input) INTEGER
1.1 bertrand 108: * On entry, INCY specifies the increment for the elements of
109: * Y. INCY must not be zero.
110: * Unchanged on exit.
111: *
112: * Level 2 Blas routine.
113: *
114: * =====================================================================
115: *
116: * .. Parameters ..
117: DOUBLE PRECISION ONE, ZERO
118: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
119: * ..
120: * .. Local Scalars ..
121: LOGICAL SYMB_ZERO
122: DOUBLE PRECISION TEMP, SAFE1
123: INTEGER I, INFO, IY, J, JX, KX, KY, LENX, LENY
124: * ..
125: * .. External Subroutines ..
126: EXTERNAL XERBLA, DLAMCH
127: DOUBLE PRECISION DLAMCH
128: * ..
129: * .. External Functions ..
130: EXTERNAL ILATRANS
131: INTEGER ILATRANS
132: * ..
133: * .. Intrinsic Functions ..
134: INTRINSIC MAX, ABS, SIGN
135: * ..
136: * .. Executable Statements ..
137: *
138: * Test the input parameters.
139: *
140: INFO = 0
141: IF ( .NOT.( ( TRANS.EQ.ILATRANS( 'N' ) )
142: $ .OR. ( TRANS.EQ.ILATRANS( 'T' ) )
143: $ .OR. ( TRANS.EQ.ILATRANS( 'C' )) ) ) THEN
144: INFO = 1
145: ELSE IF( M.LT.0 )THEN
146: INFO = 2
147: ELSE IF( N.LT.0 )THEN
148: INFO = 3
149: ELSE IF( LDA.LT.MAX( 1, M ) )THEN
150: INFO = 6
151: ELSE IF( INCX.EQ.0 )THEN
152: INFO = 8
153: ELSE IF( INCY.EQ.0 )THEN
154: INFO = 11
155: END IF
156: IF( INFO.NE.0 )THEN
157: CALL XERBLA( 'DLA_GEAMV ', INFO )
158: RETURN
159: END IF
160: *
161: * Quick return if possible.
162: *
163: IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
164: $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
165: $ RETURN
166: *
167: * Set LENX and LENY, the lengths of the vectors x and y, and set
168: * up the start points in X and Y.
169: *
170: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
171: LENX = N
172: LENY = M
173: ELSE
174: LENX = M
175: LENY = N
176: END IF
177: IF( INCX.GT.0 )THEN
178: KX = 1
179: ELSE
180: KX = 1 - ( LENX - 1 )*INCX
181: END IF
182: IF( INCY.GT.0 )THEN
183: KY = 1
184: ELSE
185: KY = 1 - ( LENY - 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(M*N) 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( TRANS.EQ.ILATRANS( 'N' ) )THEN
203: DO I = 1, LENY
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, LENX
215: TEMP = ABS( A( I, J ) )
216: SYMB_ZERO = SYMB_ZERO .AND.
217: $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
218:
219: Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
220: END DO
221: END IF
222:
223: IF ( .NOT.SYMB_ZERO )
224: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
225:
226: IY = IY + INCY
227: END DO
228: ELSE
229: DO I = 1, LENY
230: IF ( BETA .EQ. ZERO ) THEN
231: SYMB_ZERO = .TRUE.
232: Y( IY ) = 0.0D+0
233: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
234: SYMB_ZERO = .TRUE.
235: ELSE
236: SYMB_ZERO = .FALSE.
237: Y( IY ) = BETA * ABS( Y( IY ) )
238: END IF
239: IF ( ALPHA .NE. ZERO ) THEN
240: DO J = 1, LENX
241: TEMP = ABS( A( J, I ) )
242: SYMB_ZERO = SYMB_ZERO .AND.
243: $ ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
244:
245: Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
246: END DO
247: END IF
248:
249: IF ( .NOT.SYMB_ZERO )
250: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
251:
252: IY = IY + INCY
253: END DO
254: END IF
255: ELSE
256: IF( TRANS.EQ.ILATRANS( 'N' ) )THEN
257: DO I = 1, LENY
258: IF ( BETA .EQ. ZERO ) THEN
259: SYMB_ZERO = .TRUE.
260: Y( IY ) = 0.0D+0
261: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
262: SYMB_ZERO = .TRUE.
263: ELSE
264: SYMB_ZERO = .FALSE.
265: Y( IY ) = BETA * ABS( Y( IY ) )
266: END IF
267: IF ( ALPHA .NE. ZERO ) THEN
268: JX = KX
269: DO J = 1, LENX
270: TEMP = ABS( A( I, J ) )
271: SYMB_ZERO = SYMB_ZERO .AND.
272: $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
273:
274: Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
275: JX = JX + INCX
276: END DO
277: END IF
278:
279: IF (.NOT.SYMB_ZERO)
280: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
281:
282: IY = IY + INCY
283: END DO
284: ELSE
285: DO I = 1, LENY
286: IF ( BETA .EQ. ZERO ) THEN
287: SYMB_ZERO = .TRUE.
288: Y( IY ) = 0.0D+0
289: ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
290: SYMB_ZERO = .TRUE.
291: ELSE
292: SYMB_ZERO = .FALSE.
293: Y( IY ) = BETA * ABS( Y( IY ) )
294: END IF
295: IF ( ALPHA .NE. ZERO ) THEN
296: JX = KX
297: DO J = 1, LENX
298: TEMP = ABS( A( J, I ) )
299: SYMB_ZERO = SYMB_ZERO .AND.
300: $ ( X( JX ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
301:
302: Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
303: JX = JX + INCX
304: END DO
305: END IF
306:
307: IF (.NOT.SYMB_ZERO)
308: $ Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
309:
310: IY = IY + INCY
311: END DO
312: END IF
313:
314: END IF
315: *
316: RETURN
317: *
318: * End of DLA_GEAMV
319: *
320: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_geamv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack