LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dla_syamv.f
Go to the documentation of this file.
1*> \brief \b DLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bounds.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLA_SYAMV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_syamv.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_syamv.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_syamv.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
22* INCY )
23*
24* .. Scalar Arguments ..
25* DOUBLE PRECISION ALPHA, BETA
26* INTEGER INCX, INCY, LDA, N, UPLO
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DLA_SYAMV performs the matrix-vector operation
39*>
40*> y := alpha*abs(A)*abs(x) + beta*abs(y),
41*>
42*> where alpha and beta are scalars, x and y are vectors and A is an
43*> n by n symmetric matrix.
44*>
45*> This function is primarily used in calculating error bounds.
46*> To protect against underflow during evaluation, components in
47*> the resulting vector are perturbed away from zero by (N+1)
48*> times the underflow threshold. To prevent unnecessarily large
49*> errors for block-structure embedded in general matrices,
50*> "symbolically" zero components are not perturbed. A zero
51*> entry is considered "symbolic" if all multiplications involved
52*> in computing that entry have at least one zero multiplicand.
53*> \endverbatim
54*
55* Arguments:
56* ==========
57*
58*> \param[in] UPLO
59*> \verbatim
60*> UPLO is INTEGER
61*> On entry, UPLO specifies whether the upper or lower
62*> triangular part of the array A is to be referenced as
63*> follows:
64*>
65*> UPLO = BLAS_UPPER Only the upper triangular part of A
66*> is to be referenced.
67*>
68*> UPLO = BLAS_LOWER Only the lower triangular part of A
69*> is to be referenced.
70*>
71*> Unchanged on exit.
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*> N is INTEGER
77*> On entry, N specifies the number of columns of the matrix A.
78*> N must be at least zero.
79*> Unchanged on exit.
80*> \endverbatim
81*>
82*> \param[in] ALPHA
83*> \verbatim
84*> ALPHA is DOUBLE PRECISION .
85*> On entry, ALPHA specifies the scalar alpha.
86*> Unchanged on exit.
87*> \endverbatim
88*>
89*> \param[in] A
90*> \verbatim
91*> A is DOUBLE PRECISION array, dimension ( LDA, n ).
92*> Before entry, the leading m by n part of the array A must
93*> contain the matrix of coefficients.
94*> Unchanged on exit.
95*> \endverbatim
96*>
97*> \param[in] LDA
98*> \verbatim
99*> LDA is INTEGER
100*> On entry, LDA specifies the first dimension of A as declared
101*> in the calling (sub) program. LDA must be at least
102*> max( 1, n ).
103*> Unchanged on exit.
104*> \endverbatim
105*>
106*> \param[in] X
107*> \verbatim
108*> X is DOUBLE PRECISION array, dimension
109*> ( 1 + ( n - 1 )*abs( INCX ) )
110*> Before entry, the incremented array X must contain the
111*> vector x.
112*> Unchanged on exit.
113*> \endverbatim
114*>
115*> \param[in] INCX
116*> \verbatim
117*> INCX is INTEGER
118*> On entry, INCX specifies the increment for the elements of
119*> X. INCX must not be zero.
120*> Unchanged on exit.
121*> \endverbatim
122*>
123*> \param[in] BETA
124*> \verbatim
125*> BETA is DOUBLE PRECISION .
126*> On entry, BETA specifies the scalar beta. When BETA is
127*> supplied as zero then Y need not be set on input.
128*> Unchanged on exit.
129*> \endverbatim
130*>
131*> \param[in,out] Y
132*> \verbatim
133*> Y is DOUBLE PRECISION array, dimension
134*> ( 1 + ( n - 1 )*abs( INCY ) )
135*> Before entry with BETA non-zero, the incremented array Y
136*> must contain the vector y. On exit, Y is overwritten by the
137*> updated vector y.
138*> \endverbatim
139*>
140*> \param[in] INCY
141*> \verbatim
142*> INCY is INTEGER
143*> On entry, INCY specifies the increment for the elements of
144*> Y. INCY must not be zero.
145*> Unchanged on exit.
146*> \endverbatim
147*
148* Authors:
149* ========
150*
151*> \author Univ. of Tennessee
152*> \author Univ. of California Berkeley
153*> \author Univ. of Colorado Denver
154*> \author NAG Ltd.
155*
156*> \ingroup doubleSYcomputational
157*
158*> \par Further Details:
159* =====================
160*>
161*> \verbatim
162*>
163*> Level 2 Blas routine.
164*>
165*> -- Written on 22-October-1986.
166*> Jack Dongarra, Argonne National Lab.
167*> Jeremy Du Croz, Nag Central Office.
168*> Sven Hammarling, Nag Central Office.
169*> Richard Hanson, Sandia National Labs.
170*> -- Modified for the absolute-value product, April 2006
171*> Jason Riedy, UC Berkeley
172*> \endverbatim
173*>
174* =====================================================================
175 SUBROUTINE dla_syamv( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
176 $ INCY )
177*
178* -- LAPACK computational routine --
179* -- LAPACK is a software package provided by Univ. of Tennessee, --
180* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181*
182* .. Scalar Arguments ..
183 DOUBLE PRECISION ALPHA, BETA
184 INTEGER INCX, INCY, LDA, N, UPLO
185* ..
186* .. Array Arguments ..
187 DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
188* ..
189*
190* =====================================================================
191*
192* .. Parameters ..
193 DOUBLE PRECISION ONE, ZERO
194 parameter( one = 1.0d+0, zero = 0.0d+0 )
195* ..
196* .. Local Scalars ..
197 LOGICAL SYMB_ZERO
198 DOUBLE PRECISION TEMP, SAFE1
199 INTEGER I, INFO, IY, J, JX, KX, KY
200* ..
201* .. External Subroutines ..
202 EXTERNAL xerbla, dlamch
203 DOUBLE PRECISION DLAMCH
204* ..
205* .. External Functions ..
206 EXTERNAL ilauplo
207 INTEGER ILAUPLO
208* ..
209* .. Intrinsic Functions ..
210 INTRINSIC max, abs, sign
211* ..
212* .. Executable Statements ..
213*
214* Test the input parameters.
215*
216 info = 0
217 IF ( uplo.NE.ilauplo( 'U' ) .AND.
218 $ uplo.NE.ilauplo( 'L' ) ) THEN
219 info = 1
220 ELSE IF( n.LT.0 )THEN
221 info = 2
222 ELSE IF( lda.LT.max( 1, n ) )THEN
223 info = 5
224 ELSE IF( incx.EQ.0 )THEN
225 info = 7
226 ELSE IF( incy.EQ.0 )THEN
227 info = 10
228 END IF
229 IF( info.NE.0 )THEN
230 CALL xerbla( 'DLA_SYAMV', info )
231 RETURN
232 END IF
233*
234* Quick return if possible.
235*
236 IF( ( n.EQ.0 ).OR.( ( alpha.EQ.zero ).AND.( beta.EQ.one ) ) )
237 $ RETURN
238*
239* Set up the start points in X and Y.
240*
241 IF( incx.GT.0 )THEN
242 kx = 1
243 ELSE
244 kx = 1 - ( n - 1 )*incx
245 END IF
246 IF( incy.GT.0 )THEN
247 ky = 1
248 ELSE
249 ky = 1 - ( n - 1 )*incy
250 END IF
251*
252* Set SAFE1 essentially to be the underflow threshold times the
253* number of additions in each row.
254*
255 safe1 = dlamch( 'Safe minimum' )
256 safe1 = (n+1)*safe1
257*
258* Form y := alpha*abs(A)*abs(x) + beta*abs(y).
259*
260* The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
261* the inexact flag. Still doesn't help change the iteration order
262* to per-column.
263*
264 iy = ky
265 IF ( incx.EQ.1 ) THEN
266 IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
267 DO i = 1, n
268 IF ( beta .EQ. zero ) THEN
269 symb_zero = .true.
270 y( iy ) = 0.0d+0
271 ELSE IF ( y( iy ) .EQ. zero ) THEN
272 symb_zero = .true.
273 ELSE
274 symb_zero = .false.
275 y( iy ) = beta * abs( y( iy ) )
276 END IF
277 IF ( alpha .NE. zero ) THEN
278 DO j = 1, i
279 temp = abs( a( j, i ) )
280 symb_zero = symb_zero .AND.
281 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
282
283 y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
284 END DO
285 DO j = i+1, n
286 temp = abs( a( i, j ) )
287 symb_zero = symb_zero .AND.
288 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
289
290 y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
291 END DO
292 END IF
293
294 IF ( .NOT.symb_zero )
295 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
296
297 iy = iy + incy
298 END DO
299 ELSE
300 DO i = 1, n
301 IF ( beta .EQ. zero ) THEN
302 symb_zero = .true.
303 y( iy ) = 0.0d+0
304 ELSE IF ( y( iy ) .EQ. zero ) THEN
305 symb_zero = .true.
306 ELSE
307 symb_zero = .false.
308 y( iy ) = beta * abs( y( iy ) )
309 END IF
310 IF ( alpha .NE. zero ) THEN
311 DO j = 1, i
312 temp = abs( a( i, j ) )
313 symb_zero = symb_zero .AND.
314 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
315
316 y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
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( j ) )*temp
324 END DO
325 END IF
326
327 IF ( .NOT.symb_zero )
328 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
329
330 iy = iy + incy
331 END DO
332 END IF
333 ELSE
334 IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
335 DO i = 1, n
336 IF ( beta .EQ. zero ) THEN
337 symb_zero = .true.
338 y( iy ) = 0.0d+0
339 ELSE IF ( y( iy ) .EQ. zero ) THEN
340 symb_zero = .true.
341 ELSE
342 symb_zero = .false.
343 y( iy ) = beta * abs( y( iy ) )
344 END IF
345 jx = kx
346 IF ( alpha .NE. zero ) THEN
347 DO j = 1, i
348 temp = abs( a( j, i ) )
349 symb_zero = symb_zero .AND.
350 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
351
352 y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
353 jx = jx + incx
354 END DO
355 DO j = i+1, n
356 temp = abs( a( i, j ) )
357 symb_zero = symb_zero .AND.
358 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
359
360 y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
361 jx = jx + incx
362 END DO
363 END IF
364
365 IF ( .NOT.symb_zero )
366 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
367
368 iy = iy + incy
369 END DO
370 ELSE
371 DO i = 1, n
372 IF ( beta .EQ. zero ) THEN
373 symb_zero = .true.
374 y( iy ) = 0.0d+0
375 ELSE IF ( y( iy ) .EQ. zero ) THEN
376 symb_zero = .true.
377 ELSE
378 symb_zero = .false.
379 y( iy ) = beta * abs( y( iy ) )
380 END IF
381 jx = kx
382 IF ( alpha .NE. zero ) THEN
383 DO j = 1, i
384 temp = abs( a( i, j ) )
385 symb_zero = symb_zero .AND.
386 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
387
388 y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
389 jx = jx + incx
390 END DO
391 DO j = i+1, n
392 temp = abs( a( j, i ) )
393 symb_zero = symb_zero .AND.
394 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
395
396 y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
397 jx = jx + incx
398 END DO
399 END IF
400
401 IF ( .NOT.symb_zero )
402 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
403
404 iy = iy + incy
405 END DO
406 END IF
407
408 END IF
409*
410 RETURN
411*
412* End of DLA_SYAMV
413*
414 END
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
integer function ilauplo(UPLO)
ILAUPLO
Definition: ilauplo.f:58
subroutine dla_syamv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DLA_SYAMV computes a matrix-vector product using a symmetric indefinite matrix to calculate error bou...
Definition: dla_syamv.f:177