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