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