1       SUBROUTINE ZLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
  2      $                      INCY )
  3 *
  4 *     -- LAPACK routine (version 3.2.2)                                 --
  5 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
  6 *     -- Jason Riedy of Univ. of California Berkeley.                 --
  7 *     -- June 2010                                                    --
  8 *
  9 *     -- LAPACK is a software package provided by Univ. of Tennessee, --
 10 *     -- Univ. of California Berkeley and NAG Ltd.                    --
 11 *
 12       IMPLICIT NONE
 13 *     ..
 14 *     .. Scalar Arguments ..
 15       DOUBLE PRECISION   ALPHA, BETA
 16       INTEGER            INCX, INCY, LDA, N
 17       INTEGER            UPLO
 18 *     ..
 19 *     .. Array Arguments ..
 20       COMPLEX*16         A( LDA, * ), X( * )
 21       DOUBLE PRECISION   Y( * )
 22 *     ..
 23 *
 24 *  Purpose
 25 *  =======
 26 *
 27 *  ZLA_SYAMV  performs the matrix-vector operation
 28 *
 29 *          y := alpha*abs(A)*abs(x) + beta*abs(y),
 30 *
 31 *  where alpha and beta are scalars, x and y are vectors and A is an
 32 *  n by n symmetric matrix.
 33 *
 34 *  This function is primarily used in calculating error bounds.
 35 *  To protect against underflow during evaluation, components in
 36 *  the resulting vector are perturbed away from zero by (N+1)
 37 *  times the underflow threshold.  To prevent unnecessarily large
 38 *  errors for block-structure embedded in general matrices,
 39 *  "symbolically" zero components are not perturbed.  A zero
 40 *  entry is considered "symbolic" if all multiplications involved
 41 *  in computing that entry have at least one zero multiplicand.
 42 *
 43 *  Arguments
 44 *  ==========
 45 *
 46 *  UPLO    (input) INTEGER
 47 *           On entry, UPLO specifies whether the upper or lower
 48 *           triangular part of the array A is to be referenced as
 49 *           follows:
 50 *
 51 *              UPLO = BLAS_UPPER   Only the upper triangular part of A
 52 *                                  is to be referenced.
 53 *
 54 *              UPLO = BLAS_LOWER   Only the lower triangular part of A
 55 *                                  is to be referenced.
 56 *
 57 *           Unchanged on exit.
 58 *
 59 *  N       (input) INTEGER
 60 *           On entry, N specifies the number of columns of the matrix A.
 61 *           N must be at least zero.
 62 *           Unchanged on exit.
 63 *
 64 *  ALPHA  - DOUBLE PRECISION   .
 65 *           On entry, ALPHA specifies the scalar alpha.
 66 *           Unchanged on exit.
 67 *
 68 *  A      - COMPLEX*16         array of DIMENSION ( LDA, n ).
 69 *           Before entry, the leading m by n part of the array A must
 70 *           contain the matrix of coefficients.
 71 *           Unchanged on exit.
 72 *
 73 *  LDA     (input) INTEGER
 74 *           On entry, LDA specifies the first dimension of A as declared
 75 *           in the calling (sub) program. LDA must be at least
 76 *           max( 1, n ).
 77 *           Unchanged on exit.
 78 *
 79 *  X      - COMPLEX*16         array of DIMENSION at least
 80 *           ( 1 + ( n - 1 )*abs( INCX ) )
 81 *           Before entry, the incremented array X must contain the
 82 *           vector x.
 83 *           Unchanged on exit.
 84 *
 85 *  INCX    (input) INTEGER
 86 *           On entry, INCX specifies the increment for the elements of
 87 *           X. INCX must not be zero.
 88 *           Unchanged on exit.
 89 *
 90 *  BETA   - DOUBLE PRECISION   .
 91 *           On entry, BETA specifies the scalar beta. When BETA is
 92 *           supplied as zero then Y need not be set on input.
 93 *           Unchanged on exit.
 94 *
 95 *  Y       (input/output) DOUBLE PRECISION  array, dimension
 96 *           ( 1 + ( n - 1 )*abs( INCY ) )
 97 *           Before entry with BETA non-zero, the incremented array Y
 98 *           must contain the vector y. On exit, Y is overwritten by the
 99 *           updated vector y.
100 *
101 *  INCY    (input) INTEGER
102 *           On entry, INCY specifies the increment for the elements of
103 *           Y. INCY must not be zero.
104 *           Unchanged on exit.
105 *
106 *  Further Details
107 *  ===============
108 *
109 *  Level 2 Blas routine.
110 *
111 *  -- Written on 22-October-1986.
112 *     Jack Dongarra, Argonne National Lab.
113 *     Jeremy Du Croz, Nag Central Office.
114 *     Sven Hammarling, Nag Central Office.
115 *     Richard Hanson, Sandia National Labs.
116 *  -- Modified for the absolute-value product, April 2006
117 *     Jason Riedy, UC Berkeley
118 *
119 *  =====================================================================
120 *
121 *     .. Parameters ..
122       DOUBLE PRECISION   ONE, ZERO
123       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
124 *     ..
125 *     .. Local Scalars ..
126       LOGICAL            SYMB_ZERO
127       DOUBLE PRECISION   TEMP, SAFE1
128       INTEGER            I, INFO, IY, J, JX, KX, KY
129       COMPLEX*16         ZDUM
130 *     ..
131 *     .. External Subroutines ..
132       EXTERNAL           XERBLA, DLAMCH
133       DOUBLE PRECISION   DLAMCH
134 *     ..
135 *     .. External Functions ..
136       EXTERNAL           ILAUPLO
137       INTEGER            ILAUPLO
138 *     ..
139 *     .. Intrinsic Functions ..
140       INTRINSIC          MAXABSSIGN, REAL, DIMAG
141 *     ..
142 *     .. Statement Functions ..
143       DOUBLE PRECISION   CABS1
144 *     ..
145 *     .. Statement Function Definitions ..
146       CABS1( ZDUM ) = ABSDBLE ( ZDUM ) ) + ABSDIMAG ( ZDUM ) )
147 *     ..
148 *     .. Executable Statements ..
149 *
150 *     Test the input parameters.
151 *
152       INFO = 0
153       IF     ( UPLO.NE.ILAUPLO( 'U' ) .AND.
154      $         UPLO.NE.ILAUPLO( 'L' ) )THEN
155          INFO = 1
156       ELSE IF( N.LT.0 )THEN
157          INFO = 2
158       ELSE IF( LDA.LT.MAX1, N ) )THEN
159          INFO = 5
160       ELSE IF( INCX.EQ.0 )THEN
161          INFO = 7
162       ELSE IF( INCY.EQ.0 )THEN
163          INFO = 10
164       END IF
165       IF( INFO.NE.0 )THEN
166          CALL XERBLA( 'DSYMV ', INFO )
167          RETURN
168       END IF
169 *
170 *     Quick return if possible.
171 *
172       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
173      $   RETURN
174 *
175 *     Set up the start points in  X  and  Y.
176 *
177       IF( INCX.GT.0 )THEN
178          KX = 1
179       ELSE
180          KX = 1 - ( N - 1 )*INCX
181       END IF
182       IF( INCY.GT.0 )THEN
183          KY = 1
184       ELSE
185          KY = 1 - ( N - 1 )*INCY
186       END IF
187 *
188 *     Set SAFE1 essentially to be the underflow threshold times the
189 *     number of additions in each row.
190 *
191       SAFE1 = DLAMCH( 'Safe minimum' )
192       SAFE1 = (N+1)*SAFE1
193 *
194 *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
195 *
196 *     The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
197 *     the inexact flag.  Still doesn't help change the iteration order
198 *     to per-column.
199 *
200       IY = KY
201       IF ( INCX.EQ.1 ) THEN
202          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
203             DO I = 1, N
204                IF ( BETA .EQ. ZERO ) THEN
205                   SYMB_ZERO = .TRUE.
206                   Y( IY ) = 0.0D+0
207                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
208                   SYMB_ZERO = .TRUE.
209                ELSE
210                   SYMB_ZERO = .FALSE.
211                   Y( IY ) = BETA * ABS( Y( IY ) )
212                END IF
213                IF ( ALPHA .NE. ZERO ) THEN
214                   DO J = 1, I
215                      TEMP = CABS1( A( J, I ) )
216                      SYMB_ZERO = SYMB_ZERO .AND.
217      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
218 
219                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
220                   END DO
221                   DO J = I+1, N
222                      TEMP = CABS1( A( I, J ) )
223                      SYMB_ZERO = SYMB_ZERO .AND.
224      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
225 
226                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
227                   END DO
228                END IF
229 
230                IF ( .NOT.SYMB_ZERO )
231      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
232 
233                IY = IY + INCY
234             END DO
235          ELSE
236             DO I = 1, N
237                IF ( BETA .EQ. ZERO ) THEN
238                   SYMB_ZERO = .TRUE.
239                   Y( IY ) = 0.0D+0
240                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
241                   SYMB_ZERO = .TRUE.
242                ELSE
243                   SYMB_ZERO = .FALSE.
244                   Y( IY ) = BETA * ABS( Y( IY ) )
245                END IF
246                IF ( ALPHA .NE. ZERO ) THEN
247                   DO J = 1, I
248                      TEMP = CABS1( A( I, J ) )
249                      SYMB_ZERO = SYMB_ZERO .AND.
250      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
251 
252                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
253                   END DO
254                   DO J = I+1, N
255                      TEMP = CABS1( A( J, I ) )
256                      SYMB_ZERO = SYMB_ZERO .AND.
257      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
258 
259                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( J ) )*TEMP
260                   END DO
261                END IF
262 
263                IF ( .NOT.SYMB_ZERO )
264      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
265 
266                IY = IY + INCY
267             END DO
268          END IF
269       ELSE
270          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
271             DO I = 1, N
272                IF ( BETA .EQ. ZERO ) THEN
273                   SYMB_ZERO = .TRUE.
274                   Y( IY ) = 0.0D+0
275                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
276                   SYMB_ZERO = .TRUE.
277                ELSE
278                   SYMB_ZERO = .FALSE.
279                   Y( IY ) = BETA * ABS( Y( IY ) )
280                END IF
281                JX = KX
282                IF ( ALPHA .NE. ZERO ) THEN
283                   DO J = 1, I
284                      TEMP = CABS1( A( J, I ) )
285                      SYMB_ZERO = SYMB_ZERO .AND.
286      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
287 
288                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
289                      JX = JX + INCX
290                   END DO
291                   DO J = I+1, N
292                      TEMP = CABS1( A( I, J ) )
293                      SYMB_ZERO = SYMB_ZERO .AND.
294      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
295 
296                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
297                      JX = JX + INCX
298                   END DO
299                END IF
300 
301                IF ( .NOT.SYMB_ZERO )
302      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
303 
304                IY = IY + INCY
305             END DO
306          ELSE
307             DO I = 1, N
308                IF ( BETA .EQ. ZERO ) THEN
309                   SYMB_ZERO = .TRUE.
310                   Y( IY ) = 0.0D+0
311                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
312                   SYMB_ZERO = .TRUE.
313                ELSE
314                   SYMB_ZERO = .FALSE.
315                   Y( IY ) = BETA * ABS( Y( IY ) )
316                END IF
317                JX = KX
318                IF ( ALPHA .NE. ZERO ) THEN
319                   DO J = 1, I
320                      TEMP = CABS1( A( I, J ) )
321                      SYMB_ZERO = SYMB_ZERO .AND.
322      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
323 
324                      Y( IY ) = Y( IY ) + ALPHA*CABS1( X( JX ) )*TEMP
325                      JX = JX + INCX
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( JX ) )*TEMP
333                      JX = JX + INCX
334                   END DO
335                END IF
336 
337                IF ( .NOT.SYMB_ZERO )
338      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
339 
340                IY = IY + INCY
341             END DO
342          END IF
343 
344       END IF
345 *
346       RETURN
347 *
348 *     End of ZLA_SYAMV
349 *
350       END