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