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