001:       SUBROUTINE DLA_SYAMV( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y,
002:      $                      INCY )
003: *
004: *     -- LAPACK routine (version 3.2)                                 --
005: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
006: *     -- Jason Riedy of Univ. of California Berkeley.                 --
007: *     -- November 2008                                                --
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:       DOUBLE PRECISION   ALPHA, BETA
016:       INTEGER            INCX, INCY, LDA, N, UPLO
017: *     ..
018: *     .. Array Arguments ..
019:       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  DLA_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: *  Parameters
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  - DOUBLE PRECISION   .
063: *           On entry, ALPHA specifies the scalar alpha.
064: *           Unchanged on exit.
065: *
066: *  A      - DOUBLE PRECISION   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      - DOUBLE PRECISION   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   - DOUBLE PRECISION   .
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      - DOUBLE PRECISION   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: *
105: *  Level 2 Blas routine.
106: *
107: *  -- Written on 22-October-1986.
108: *     Jack Dongarra, Argonne National Lab.
109: *     Jeremy Du Croz, Nag Central Office.
110: *     Sven Hammarling, Nag Central Office.
111: *     Richard Hanson, Sandia National Labs.
112: *  -- Modified for the absolute-value product, April 2006
113: *     Jason Riedy, UC Berkeley
114: *
115: *     ..
116: *     .. Parameters ..
117:       DOUBLE PRECISION   ONE, ZERO
118:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
119: *     ..
120: *     .. Local Scalars ..
121:       LOGICAL            SYMB_ZERO
122:       DOUBLE PRECISION   TEMP, SAFE1
123:       INTEGER            I, INFO, IY, J, JX, KX, KY
124: *     ..
125: *     .. External Subroutines ..
126:       EXTERNAL           XERBLA, DLAMCH
127:       DOUBLE PRECISION   DLAMCH
128: *     ..
129: *     .. External Functions ..
130:       EXTERNAL           ILAUPLO
131:       INTEGER            ILAUPLO
132: *     ..
133: *     .. Intrinsic Functions ..
134:       INTRINSIC          MAX, ABS, SIGN
135: *     ..
136: *     .. Executable Statements ..
137: *
138: *     Test the input parameters.
139: *
140:       INFO = 0
141:       IF     ( UPLO.NE.ILAUPLO( 'U' ) .AND.
142:      $         UPLO.NE.ILAUPLO( 'L' ) ) THEN
143:          INFO = 1
144:       ELSE IF( N.LT.0 )THEN
145:          INFO = 2
146:       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
147:          INFO = 5
148:       ELSE IF( INCX.EQ.0 )THEN
149:          INFO = 7
150:       ELSE IF( INCY.EQ.0 )THEN
151:          INFO = 10
152:       END IF
153:       IF( INFO.NE.0 )THEN
154:          CALL XERBLA( 'DSYMV ', INFO )
155:          RETURN
156:       END IF
157: *
158: *     Quick return if possible.
159: *
160:       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
161:      $   RETURN
162: *
163: *     Set up the start points in  X  and  Y.
164: *
165:       IF( INCX.GT.0 )THEN
166:          KX = 1
167:       ELSE
168:          KX = 1 - ( N - 1 )*INCX
169:       END IF
170:       IF( INCY.GT.0 )THEN
171:          KY = 1
172:       ELSE
173:          KY = 1 - ( N - 1 )*INCY
174:       END IF
175: *
176: *     Set SAFE1 essentially to be the underflow threshold times the
177: *     number of additions in each row.
178: *
179:       SAFE1 = DLAMCH( 'Safe minimum' )
180:       SAFE1 = (N+1)*SAFE1
181: *
182: *     Form  y := alpha*abs(A)*abs(x) + beta*abs(y).
183: *
184: *     The O(N^2) SYMB_ZERO tests could be replaced by O(N) queries to
185: *     the inexact flag.  Still doesn't help change the iteration order
186: *     to per-column.
187: *
188:       IY = KY
189:       IF ( INCX.EQ.1 ) THEN
190:          DO I = 1, N
191:             IF ( BETA .EQ. ZERO ) THEN
192:                SYMB_ZERO = .TRUE.
193:                Y( IY ) = 0.0D+0
194:             ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
195:                SYMB_ZERO = .TRUE.
196:             ELSE
197:                SYMB_ZERO = .FALSE.
198:                Y( IY ) = BETA * ABS( Y( IY ) )
199:             END IF
200:             IF ( ALPHA .NE. ZERO ) THEN
201:                DO J = 1, N
202:                   IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
203:                      IF ( I .LE. J ) THEN
204:                         TEMP = ABS( A( I, J ) )
205:                      ELSE
206:                         TEMP = ABS( A( J, I ) )
207:                      END IF
208:                   ELSE
209:                      IF ( I .GE. J ) THEN
210:                         TEMP = ABS( A( I, J ) )
211:                      ELSE
212:                         TEMP = ABS( A( J, I ) )
213:                      END IF
214:                   END IF
215: 
216:                   SYMB_ZERO = SYMB_ZERO .AND.
217:      $                 ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
218: 
219:                   Y( IY ) = Y( IY ) + ALPHA*ABS( X( J ) )*TEMP
220:                END DO
221:             END IF
222: 
223:             IF ( .NOT.SYMB_ZERO )
224:      $           Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
225: 
226:             IY = IY + INCY
227:          END DO
228:       ELSE
229:          DO I = 1, N
230:             IF ( BETA .EQ. ZERO ) THEN
231:                SYMB_ZERO = .TRUE.
232:                Y( IY ) = 0.0D+0
233:             ELSE IF ( Y( IY ) .EQ. ZERO ) THEN
234:                SYMB_ZERO = .TRUE.
235:             ELSE
236:                SYMB_ZERO = .FALSE.
237:                Y( IY ) = BETA * ABS( Y( IY ) )
238:             END IF
239:             JX = KX
240:             IF ( ALPHA .NE. ZERO ) THEN
241:                DO J = 1, N
242:                   IF ( UPLO .EQ. ILAUPLO( 'U' ) ) THEN
243:                      IF ( I .LE. J ) THEN
244:                         TEMP = ABS( A( I, J ) )
245:                      ELSE
246:                         TEMP = ABS( A( J, I ) )
247:                      END IF
248:                   ELSE
249:                      IF ( I .GE. J ) THEN
250:                         TEMP = ABS( A( I, J ) )
251:                      ELSE
252:                         TEMP = ABS( A( J, I ) )
253:                      END IF
254:                   END IF
255: 
256:                   SYMB_ZERO = SYMB_ZERO .AND.
257:      $                 ( X( J ) .EQ. ZERO .OR. TEMP .EQ. ZERO )
258: 
259:                   Y( IY ) = Y( IY ) + ALPHA*ABS( X( JX ) )*TEMP
260:                   JX = JX + INCX
261:                END DO
262:             END IF
263: 
264:             IF ( .NOT.SYMB_ZERO )
265:      $           Y( IY ) = Y( IY ) + SIGN( SAFE1, Y( IY ) )
266: 
267:             IY = IY + INCY
268:          END DO
269:       END IF
270: *
271:       RETURN
272: *
273: *     End of DLA_SYAMV
274: *
275:       END
276: