001:       SUBROUTINE SLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
002:      $                      INCY )
003: *
004: *     -- LAPACK routine (version 3.2.1)                                 --
005: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
006: *     -- Jason Riedy of Univ. of California Berkeley.                 --
007: *     -- April 2009                                                   --
008: *
009: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
010: *     -- Univ. of California Berkeley and NAG Ltd.                    --
011: *
012:       IMPLICIT NONE
013: *     ..
014: *     .. Scalar Arguments ..
015:       REAL               ALPHA, BETA
016:       INTEGER            INCX, INCY, LDA, N, UPLO
017: *     ..
018: *     .. Array Arguments ..
019:       REAL               A( LDA, * ), X( * ), Y( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  SLA_SYAMV  performs the matrix-vector operation
026: *
027: *          y := alpha*abs(A)*abs(x) + beta*abs(y),
028: *
029: *  where alpha and beta are scalars, x and y are vectors and A is an
030: *  n by n symmetric matrix.
031: *
032: *  This function is primarily used in calculating error bounds.
033: *  To protect against underflow during evaluation, components in
034: *  the resulting vector are perturbed away from zero by (N+1)
035: *  times the underflow threshold.  To prevent unnecessarily large
036: *  errors for block-structure embedded in general matrices,
037: *  "symbolically" zero components are not perturbed.  A zero
038: *  entry is considered "symbolic" if all multiplications involved
039: *  in computing that entry have at least one zero multiplicand.
040: *
041: *  Arguments
042: *  ==========
043: *
044: *  UPLO   - INTEGER
045: *           On entry, UPLO specifies whether the upper or lower
046: *           triangular part of the array A is to be referenced as
047: *           follows:
048: *
049: *              UPLO = BLAS_UPPER   Only the upper triangular part of A
050: *                                  is to be referenced.
051: *
052: *              UPLO = BLAS_LOWER   Only the lower triangular part of A
053: *                                  is to be referenced.
054: *
055: *           Unchanged on exit.
056: *
057: *  N      - INTEGER.
058: *           On entry, N specifies the number of columns of the matrix A.
059: *           N must be at least zero.
060: *           Unchanged on exit.
061: *
062: *  ALPHA  - REAL            .
063: *           On entry, ALPHA specifies the scalar alpha.
064: *           Unchanged on exit.
065: *
066: *  A      - REAL             array of DIMENSION ( LDA, n ).
067: *           Before entry, the leading m by n part of the array A must
068: *           contain the matrix of coefficients.
069: *           Unchanged on exit.
070: *
071: *  LDA    - INTEGER.
072: *           On entry, LDA specifies the first dimension of A as declared
073: *           in the calling (sub) program. LDA must be at least
074: *           max( 1, n ).
075: *           Unchanged on exit.
076: *
077: *  X      - REAL             array of DIMENSION at least
078: *           ( 1 + ( n - 1 )*abs( INCX ) )
079: *           Before entry, the incremented array X must contain the
080: *           vector x.
081: *           Unchanged on exit.
082: *
083: *  INCX   - INTEGER.
084: *           On entry, INCX specifies the increment for the elements of
085: *           X. INCX must not be zero.
086: *           Unchanged on exit.
087: *
088: *  BETA   - REAL            .
089: *           On entry, BETA specifies the scalar beta. When BETA is
090: *           supplied as zero then Y need not be set on input.
091: *           Unchanged on exit.
092: *
093: *  Y      - REAL             array of DIMENSION at least
094: *           ( 1 + ( n - 1 )*abs( INCY ) )
095: *           Before entry with BETA non-zero, the incremented array Y
096: *           must contain the vector y. On exit, Y is overwritten by the
097: *           updated vector y.
098: *
099: *  INCY   - INTEGER.
100: *           On entry, INCY specifies the increment for the elements of
101: *           Y. INCY must not be zero.
102: *           Unchanged on exit.
103: *
104: *  Further Details
105: *  ===============
106: *
107: *  Level 2 Blas routine.
108: *
109: *  -- Written on 22-October-1986.
110: *     Jack Dongarra, Argonne National Lab.
111: *     Jeremy Du Croz, Nag Central Office.
112: *     Sven Hammarling, Nag Central Office.
113: *     Richard Hanson, Sandia National Labs.
114: *  -- Modified for the absolute-value product, April 2006
115: *     Jason Riedy, UC Berkeley
116: *
117: *  =====================================================================
118: *
119: *     .. Parameters ..
120:       REAL               ONE, ZERO
121:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
122: *     ..
123: *     .. Local Scalars ..
124:       LOGICAL            SYMB_ZERO
125:       REAL               TEMP, SAFE1
126:       INTEGER            I, INFO, IY, J, JX, KX, KY
127: *     ..
128: *     .. External Subroutines ..
129:       EXTERNAL           XERBLA, SLAMCH
130:       REAL               SLAMCH
131: *     ..
132: *     .. External Functions ..
133:       EXTERNAL           ILAUPLO
134:       INTEGER            ILAUPLO
135: *     ..
136: *     .. Intrinsic Functions ..
137:       INTRINSIC          MAX, ABS, SIGN
138: *     ..
139: *     .. Executable Statements ..
140: *
141: *     Test the input parameters.
142: *
143:       INFO = 0
144:       IF     ( UPLO.NE.ILAUPLO( 'U' ) .AND.
145:      $         UPLO.NE.ILAUPLO( 'L' ) ) THEN
146:          INFO = 1
147:       ELSE IF( N.LT.0 )THEN
148:          INFO = 2
149:       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
150:          INFO = 5
151:       ELSE IF( INCX.EQ.0 )THEN
152:          INFO = 7
153:       ELSE IF( INCY.EQ.0 )THEN
154:          INFO = 10
155:       END IF
156:       IF( INFO.NE.0 )THEN
157:          CALL XERBLA( 'SSYMV ', INFO )
158:          RETURN
159:       END IF
160: *
161: *     Quick return if possible.
162: *
163:       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
164:      $   RETURN
165: *
166: *     Set up the start points in  X  and  Y.
167: *
168:       IF( INCX.GT.0 )THEN
169:          KX = 1
170:       ELSE
171:          KX = 1 - ( N - 1 )*INCX
172:       END IF
173:       IF( INCY.GT.0 )THEN
174:          KY = 1
175:       ELSE
176:          KY = 1 - ( N - 1 )*INCY
177:       END IF
178: *
179: *     Set SAFE1 essentially to be the underflow threshold times the
180: *     number of additions in each row.
181: *
182:       SAFE1 = SLAMCH( 'Safe minimum' )
183:       SAFE1 = (N+1)*SAFE1
184: *
185: *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
186: *
187: *     The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
188: *     the inexact flag.  Still doesn't help change the iteration order
189: *     to per-column.
190: *
191:       IY = KY
192:       IF ( INCX.EQ.1 ) THEN
193:          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
194:             DO I = 1, N
195:                IF ( BETA .EQ. ZERO ) THEN
196:                   SYMB_ZERO = .TRUE.
197:                   Y( IY ) = 0.0
198:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
199:                   SYMB_ZERO = .TRUE.
200:                ELSE
201:                   SYMB_ZERO = .FALSE.
202:                   Y( IY ) = BETA * ABS( Y( IY ) )
203:                END IF
204:                IF ( ALPHA .NE. ZERO ) THEN
205:                   DO J = 1, I
206:                      TEMP = ABS( A( J, I ) )
207:                      SYMB_ZERO = SYMB_ZERO .AND.
208:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
209: 
210:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
211:                   END DO
212:                   DO J = I+1, N
213:                      TEMP = ABS( A( I, J ) )
214:                      SYMB_ZERO = SYMB_ZERO .AND.
215:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
216: 
217:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
218:                   END DO
219:                END IF
220: 
221:                IF ( .NOT.SYMB_ZERO )
222:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
223: 
224:                IY = IY + INCY
225:             END DO
226:          ELSE
227:             DO I = 1, N
228:                IF ( BETA .EQ. ZERO ) THEN
229:                   SYMB_ZERO = .TRUE.
230:                   Y( IY ) = 0.0
231:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
232:                   SYMB_ZERO = .TRUE.
233:                ELSE
234:                   SYMB_ZERO = .FALSE.
235:                   Y( IY ) = BETA * ABS( Y( IY ) )
236:                END IF
237:                IF ( ALPHA .NE. ZERO ) THEN
238:                   DO J = 1, I
239:                      TEMP = ABS( A( I, J ) )
240:                      SYMB_ZERO = SYMB_ZERO .AND.
241:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
242: 
243:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
244:                   END DO
245:                   DO J = I+1, N
246:                      TEMP = ABS( A( J, I ) )
247:                      SYMB_ZERO = SYMB_ZERO .AND.
248:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
249: 
250:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
251:                   END DO
252:                END IF
253: 
254:                IF ( .NOT.SYMB_ZERO )
255:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
256: 
257:                IY = IY + INCY
258:             END DO
259:          END IF
260:       ELSE
261:          IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
262:             DO I = 1, N
263:                IF ( BETA .EQ. ZERO ) THEN
264:                   SYMB_ZERO = .TRUE.
265:                   Y( IY ) = 0.0
266:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
267:                   SYMB_ZERO = .TRUE.
268:                ELSE
269:                   SYMB_ZERO = .FALSE.
270:                   Y( IY ) = BETA * ABS( Y( IY ) )
271:                END IF
272:                JX = KX
273:                IF ( ALPHA .NE. ZERO ) THEN
274:                   DO J = 1, I
275:                      TEMP = ABS( A( J, I ) )
276:                      SYMB_ZERO = SYMB_ZERO .AND.
277:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
278: 
279:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
280:                      JX = JX + INCX
281:                   END DO
282:                   DO J = I+1, N
283:                      TEMP = ABS( A( I, J ) )
284:                      SYMB_ZERO = SYMB_ZERO .AND.
285:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
286: 
287:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
288:                      JX = JX + INCX
289:                   END DO
290:                END IF
291: 
292:                IF ( .NOT.SYMB_ZERO )
293:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
294: 
295:                IY = IY + INCY
296:             END DO
297:          ELSE
298:             DO I = 1, N
299:                IF ( BETA .EQ. ZERO ) THEN
300:                   SYMB_ZERO = .TRUE.
301:                   Y( IY ) = 0.0
302:                ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
303:                   SYMB_ZERO = .TRUE.
304:                ELSE
305:                   SYMB_ZERO = .FALSE.
306:                   Y( IY ) = BETA * ABS( Y( IY ) )
307:                END IF
308:                JX = KX
309:                IF ( ALPHA .NE. ZERO ) THEN
310:                   DO J = 1, I
311:                      TEMP = ABS( A( I, J ) )
312:                      SYMB_ZERO = SYMB_ZERO .AND.
313:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
314: 
315:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
316:                      JX = JX + INCX
317:                   END DO
318:                   DO J = I+1, N
319:                      TEMP = ABS( A( J, I ) )
320:                      SYMB_ZERO = SYMB_ZERO .AND.
321:      $                    ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
322: 
323:                      Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
324:                      JX = JX + INCX
325:                   END DO
326:                END IF
327: 
328:                IF ( .NOT.SYMB_ZERO )
329:      $              Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
330: 
331:                IY = IY + INCY
332:             END DO
333:          END IF
334: 
335:       END IF
336: *
337:       RETURN
338: *
339: *     End of SLA_SYAMV
340: *
341:       END
342: