1*> \brief \b CLA_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 CLA_SYAMV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cla_syamv.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_syamv.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_syamv.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
22*                             INCY )
23*
24*       .. Scalar Arguments ..
25*       REAL               ALPHA, BETA
26*       INTEGER            INCX, INCY, LDA, N
27*       INTEGER            UPLO
28*       ..
29*       .. Array Arguments ..
30*       COMPLEX            A( LDA, * ), X( * )
31*       REAL               Y( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> CLA_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 REAL .
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 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 array, dimension
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 REAL .
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 REAL 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 complexSYcomputational
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 CLA_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      REAL               ALPHA, BETA
186      INTEGER            INCX, INCY, LDA, N
187      INTEGER            UPLO
188*     ..
189*     .. Array Arguments ..
190      COMPLEX            A( LDA, * ), X( * )
191      REAL               Y( * )
192*     ..
193*
194*  =====================================================================
195*
196*     .. Parameters ..
197      REAL               ONE, ZERO
198      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
199*     ..
200*     .. Local Scalars ..
201      LOGICAL            SYMB_ZERO
202      REAL               TEMP, SAFE1
203      INTEGER            I, INFO, IY, J, JX, KX, KY
204      COMPLEX            ZDUM
205*     ..
206*     .. External Subroutines ..
207      EXTERNAL           XERBLA, SLAMCH
208      REAL               SLAMCH
209*     ..
210*     .. External Functions ..
211      EXTERNAL           ILAUPLO
212      INTEGER            ILAUPLO
213*     ..
214*     .. Intrinsic Functions ..
215      INTRINSIC          MAX, ABS, SIGN, REAL, AIMAG
216*     ..
217*     .. Statement Functions ..
218      REAL               CABS1
219*     ..
220*     .. Statement Function Definitions ..
221      CABS1( ZDUM ) = ABS( REAL ( ZDUM ) ) + ABS( AIMAG ( 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( 'CLA_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 = SLAMCH( '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.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.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.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.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 CLA_SYAMV
424*
425      END
426