LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \ingroup complex16HEcomputational
158 *
159 *> \par Further Details:
160 * =====================
161 *>
162 *> \verbatim
163 *>
164 *> Level 2 Blas routine.
165 *>
166 *> -- Written on 22-October-1986.
167 *> Jack Dongarra, Argonne National Lab.
168 *> Jeremy Du Croz, Nag Central Office.
169 *> Sven Hammarling, Nag Central Office.
170 *> Richard Hanson, Sandia National Labs.
171 *> -- Modified for the absolute-value product, April 2006
172 *> Jason Riedy, UC Berkeley
173 *> \endverbatim
174 *>
175 * =====================================================================
176  SUBROUTINE zla_heamv( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
177  $ INCY )
178 *
179 * -- LAPACK computational routine --
180 * -- LAPACK is a software package provided by Univ. of Tennessee, --
181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182 *
183 * .. Scalar Arguments ..
184  DOUBLE PRECISION ALPHA, BETA
185  INTEGER INCX, INCY, LDA, N, UPLO
186 * ..
187 * .. Array Arguments ..
188  COMPLEX*16 A( LDA, * ), X( * )
189  DOUBLE PRECISION Y( * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  DOUBLE PRECISION ONE, ZERO
196  parameter( one = 1.0d+0, zero = 0.0d+0 )
197 * ..
198 * .. Local Scalars ..
199  LOGICAL SYMB_ZERO
200  DOUBLE PRECISION TEMP, SAFE1
201  INTEGER I, INFO, IY, J, JX, KX, KY
202  COMPLEX*16 ZDUM
203 * ..
204 * .. External Subroutines ..
205  EXTERNAL xerbla, dlamch
206  DOUBLE PRECISION DLAMCH
207 * ..
208 * .. External Functions ..
209  EXTERNAL ilauplo
210  INTEGER ILAUPLO
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC max, abs, sign, real, dimag
214 * ..
215 * .. Statement Functions ..
216  DOUBLE PRECISION CABS1
217 * ..
218 * .. Statement Function Definitions ..
219  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
220 * ..
221 * .. Executable Statements ..
222 *
223 * Test the input parameters.
224 *
225  info = 0
226  IF ( uplo.NE.ilauplo( 'U' ) .AND.
227  $ uplo.NE.ilauplo( 'L' ) )THEN
228  info = 1
229  ELSE IF( n.LT.0 )THEN
230  info = 2
231  ELSE IF( lda.LT.max( 1, n ) )THEN
232  info = 5
233  ELSE IF( incx.EQ.0 )THEN
234  info = 7
235  ELSE IF( incy.EQ.0 )THEN
236  info = 10
237  END IF
238  IF( info.NE.0 )THEN
239  CALL xerbla( 'ZHEMV ', info )
240  RETURN
241  END IF
242 *
243 * Quick return if possible.
244 *
245  IF( ( n.EQ.0 ).OR.( ( alpha.EQ.zero ).AND.( beta.EQ.one ) ) )
246  $ RETURN
247 *
248 * Set up the start points in X and Y.
249 *
250  IF( incx.GT.0 )THEN
251  kx = 1
252  ELSE
253  kx = 1 - ( n - 1 )*incx
254  END IF
255  IF( incy.GT.0 )THEN
256  ky = 1
257  ELSE
258  ky = 1 - ( n - 1 )*incy
259  END IF
260 *
261 * Set SAFE1 essentially to be the underflow threshold times the
262 * number of additions in each row.
263 *
264  safe1 = dlamch( 'Safe minimum' )
265  safe1 = (n+1)*safe1
266 *
267 * Form y := alpha*abs(A)*abs(x) + beta*abs(y).
268 *
269 * The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
270 * the inexact flag. Still doesn't help change the iteration order
271 * to per-column.
272 *
273  iy = ky
274  IF ( incx.EQ.1 ) THEN
275  IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
276  DO i = 1, n
277  IF ( beta .EQ. zero ) THEN
278  symb_zero = .true.
279  y( iy ) = 0.0d+0
280  ELSE IF ( y( iy ) .EQ. zero ) THEN
281  symb_zero = .true.
282  ELSE
283  symb_zero = .false.
284  y( iy ) = beta * abs( y( iy ) )
285  END IF
286  IF ( alpha .NE. zero ) THEN
287  DO j = 1, i
288  temp = cabs1( a( j, i ) )
289  symb_zero = symb_zero .AND.
290  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
291 
292  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
293  END DO
294  DO j = i+1, n
295  temp = cabs1( a( i, j ) )
296  symb_zero = symb_zero .AND.
297  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
298 
299  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
300  END DO
301  END IF
302 
303  IF (.NOT.symb_zero)
304  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
305 
306  iy = iy + incy
307  END DO
308  ELSE
309  DO i = 1, n
310  IF ( beta .EQ. zero ) THEN
311  symb_zero = .true.
312  y( iy ) = 0.0d+0
313  ELSE IF ( y( iy ) .EQ. zero ) THEN
314  symb_zero = .true.
315  ELSE
316  symb_zero = .false.
317  y( iy ) = beta * abs( y( iy ) )
318  END IF
319  IF ( alpha .NE. zero ) THEN
320  DO j = 1, i
321  temp = cabs1( a( i, j ) )
322  symb_zero = symb_zero .AND.
323  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
324 
325  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
326  END DO
327  DO j = i+1, n
328  temp = cabs1( a( j, i ) )
329  symb_zero = symb_zero .AND.
330  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
331 
332  y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
333  END DO
334  END IF
335 
336  IF (.NOT.symb_zero)
337  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
338 
339  iy = iy + incy
340  END DO
341  END IF
342  ELSE
343  IF ( uplo .EQ. ilauplo( 'U' ) ) THEN
344  DO i = 1, n
345  IF ( beta .EQ. zero ) THEN
346  symb_zero = .true.
347  y( iy ) = 0.0d+0
348  ELSE IF ( y( iy ) .EQ. zero ) THEN
349  symb_zero = .true.
350  ELSE
351  symb_zero = .false.
352  y( iy ) = beta * abs( y( iy ) )
353  END IF
354  jx = kx
355  IF ( alpha .NE. zero ) THEN
356  DO j = 1, i
357  temp = cabs1( a( j, i ) )
358  symb_zero = symb_zero .AND.
359  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
360 
361  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
362  jx = jx + incx
363  END DO
364  DO j = i+1, n
365  temp = cabs1( a( i, j ) )
366  symb_zero = symb_zero .AND.
367  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
368 
369  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
370  jx = jx + incx
371  END DO
372  END IF
373 
374  IF ( .NOT.symb_zero )
375  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
376 
377  iy = iy + incy
378  END DO
379  ELSE
380  DO i = 1, n
381  IF ( beta .EQ. zero ) THEN
382  symb_zero = .true.
383  y( iy ) = 0.0d+0
384  ELSE IF ( y( iy ) .EQ. zero ) THEN
385  symb_zero = .true.
386  ELSE
387  symb_zero = .false.
388  y( iy ) = beta * abs( y( iy ) )
389  END IF
390  jx = kx
391  IF ( alpha .NE. zero ) THEN
392  DO j = 1, i
393  temp = cabs1( a( i, j ) )
394  symb_zero = symb_zero .AND.
395  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
396 
397  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
398  jx = jx + incx
399  END DO
400  DO j = i+1, n
401  temp = cabs1( a( j, i ) )
402  symb_zero = symb_zero .AND.
403  $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
404 
405  y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
406  jx = jx + incx
407  END DO
408  END IF
409 
410  IF ( .NOT.symb_zero )
411  $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
412 
413  iy = iy + incy
414  END DO
415  END IF
416 
417  END IF
418 *
419  RETURN
420 *
421 * End of ZLA_HEAMV
422 *
423  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 zla_heamv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZLA_HEAMV computes a matrix-vector product using a Hermitian indefinite matrix to calculate error bou...
Definition: zla_heamv.f:178