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