001:       SUBROUTINE ZSPMV( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          UPLO
009:       INTEGER            INCX, INCY, N
010:       COMPLEX*16         ALPHA, BETA
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX*16         AP( * ), X( * ), Y( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZSPMV  performs the matrix-vector operation
020: *
021: *     y := alpha*A*x + beta*y,
022: *
023: *  where alpha and beta are scalars, x and y are n element vectors and
024: *  A is an n by n symmetric matrix, supplied in packed form.
025: *
026: *  Arguments
027: *  ==========
028: *
029: *  UPLO     (input) CHARACTER*1
030: *           On entry, UPLO specifies whether the upper or lower
031: *           triangular part of the matrix A is supplied in the packed
032: *           array AP as follows:
033: *
034: *              UPLO = 'U' or 'u'   The upper triangular part of A is
035: *                                  supplied in AP.
036: *
037: *              UPLO = 'L' or 'l'   The lower triangular part of A is
038: *                                  supplied in AP.
039: *
040: *           Unchanged on exit.
041: *
042: *  N        (input) INTEGER
043: *           On entry, N specifies the order of the matrix A.
044: *           N must be at least zero.
045: *           Unchanged on exit.
046: *
047: *  ALPHA    (input) COMPLEX*16
048: *           On entry, ALPHA specifies the scalar alpha.
049: *           Unchanged on exit.
050: *
051: *  AP       (input) COMPLEX*16 array, dimension at least
052: *           ( ( N*( N + 1 ) )/2 ).
053: *           Before entry, with UPLO = 'U' or 'u', the array AP must
054: *           contain the upper triangular part of the symmetric matrix
055: *           packed sequentially, column by column, so that AP( 1 )
056: *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
057: *           and a( 2, 2 ) respectively, and so on.
058: *           Before entry, with UPLO = 'L' or 'l', the array AP must
059: *           contain the lower triangular part of the symmetric matrix
060: *           packed sequentially, column by column, so that AP( 1 )
061: *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
062: *           and a( 3, 1 ) respectively, and so on.
063: *           Unchanged on exit.
064: *
065: *  X        (input) COMPLEX*16 array, dimension at least
066: *           ( 1 + ( N - 1 )*abs( INCX ) ).
067: *           Before entry, the incremented array X must contain the N-
068: *           element vector x.
069: *           Unchanged on exit.
070: *
071: *  INCX     (input) INTEGER
072: *           On entry, INCX specifies the increment for the elements of
073: *           X. INCX must not be zero.
074: *           Unchanged on exit.
075: *
076: *  BETA     (input) COMPLEX*16
077: *           On entry, BETA specifies the scalar beta. When BETA is
078: *           supplied as zero then Y need not be set on input.
079: *           Unchanged on exit.
080: *
081: *  Y        (input/output) COMPLEX*16 array, dimension at least
082: *           ( 1 + ( N - 1 )*abs( INCY ) ).
083: *           Before entry, the incremented array Y must contain the n
084: *           element vector y. On exit, Y is overwritten by the updated
085: *           vector y.
086: *
087: *  INCY     (input) INTEGER
088: *           On entry, INCY specifies the increment for the elements of
089: *           Y. INCY must not be zero.
090: *           Unchanged on exit.
091: *
092: * =====================================================================
093: *
094: *     .. Parameters ..
095:       COMPLEX*16         ONE
096:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
097:       COMPLEX*16         ZERO
098:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
099: *     ..
100: *     .. Local Scalars ..
101:       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
102:       COMPLEX*16         TEMP1, TEMP2
103: *     ..
104: *     .. External Functions ..
105:       LOGICAL            LSAME
106:       EXTERNAL           LSAME
107: *     ..
108: *     .. External Subroutines ..
109:       EXTERNAL           XERBLA
110: *     ..
111: *     .. Executable Statements ..
112: *
113: *     Test the input parameters.
114: *
115:       INFO = 0
116:       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
117:          INFO = 1
118:       ELSE IF( N.LT.0 ) THEN
119:          INFO = 2
120:       ELSE IF( INCX.EQ.0 ) THEN
121:          INFO = 6
122:       ELSE IF( INCY.EQ.0 ) THEN
123:          INFO = 9
124:       END IF
125:       IF( INFO.NE.0 ) THEN
126:          CALL XERBLA( 'ZSPMV ', INFO )
127:          RETURN
128:       END IF
129: *
130: *     Quick return if possible.
131: *
132:       IF( ( N.EQ.0 ) .OR. ( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ONE ) ) )
133:      $   RETURN
134: *
135: *     Set up the start points in  X  and  Y.
136: *
137:       IF( INCX.GT.0 ) THEN
138:          KX = 1
139:       ELSE
140:          KX = 1 - ( N-1 )*INCX
141:       END IF
142:       IF( INCY.GT.0 ) THEN
143:          KY = 1
144:       ELSE
145:          KY = 1 - ( N-1 )*INCY
146:       END IF
147: *
148: *     Start the operations. In this version the elements of the array AP
149: *     are accessed sequentially with one pass through AP.
150: *
151: *     First form  y := beta*y.
152: *
153:       IF( BETA.NE.ONE ) THEN
154:          IF( INCY.EQ.1 ) THEN
155:             IF( BETA.EQ.ZERO ) THEN
156:                DO 10 I = 1, N
157:                   Y( I ) = ZERO
158:    10          CONTINUE
159:             ELSE
160:                DO 20 I = 1, N
161:                   Y( I ) = BETA*Y( I )
162:    20          CONTINUE
163:             END IF
164:          ELSE
165:             IY = KY
166:             IF( BETA.EQ.ZERO ) THEN
167:                DO 30 I = 1, N
168:                   Y( IY ) = ZERO
169:                   IY = IY + INCY
170:    30          CONTINUE
171:             ELSE
172:                DO 40 I = 1, N
173:                   Y( IY ) = BETA*Y( IY )
174:                   IY = IY + INCY
175:    40          CONTINUE
176:             END IF
177:          END IF
178:       END IF
179:       IF( ALPHA.EQ.ZERO )
180:      $   RETURN
181:       KK = 1
182:       IF( LSAME( UPLO, 'U' ) ) THEN
183: *
184: *        Form  y  when AP contains the upper triangle.
185: *
186:          IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN
187:             DO 60 J = 1, N
188:                TEMP1 = ALPHA*X( J )
189:                TEMP2 = ZERO
190:                K = KK
191:                DO 50 I = 1, J - 1
192:                   Y( I ) = Y( I ) + TEMP1*AP( K )
193:                   TEMP2 = TEMP2 + AP( K )*X( I )
194:                   K = K + 1
195:    50          CONTINUE
196:                Y( J ) = Y( J ) + TEMP1*AP( KK+J-1 ) + ALPHA*TEMP2
197:                KK = KK + J
198:    60       CONTINUE
199:          ELSE
200:             JX = KX
201:             JY = KY
202:             DO 80 J = 1, N
203:                TEMP1 = ALPHA*X( JX )
204:                TEMP2 = ZERO
205:                IX = KX
206:                IY = KY
207:                DO 70 K = KK, KK + J - 2
208:                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
209:                   TEMP2 = TEMP2 + AP( K )*X( IX )
210:                   IX = IX + INCX
211:                   IY = IY + INCY
212:    70          CONTINUE
213:                Y( JY ) = Y( JY ) + TEMP1*AP( KK+J-1 ) + ALPHA*TEMP2
214:                JX = JX + INCX
215:                JY = JY + INCY
216:                KK = KK + J
217:    80       CONTINUE
218:          END IF
219:       ELSE
220: *
221: *        Form  y  when AP contains the lower triangle.
222: *
223:          IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN
224:             DO 100 J = 1, N
225:                TEMP1 = ALPHA*X( J )
226:                TEMP2 = ZERO
227:                Y( J ) = Y( J ) + TEMP1*AP( KK )
228:                K = KK + 1
229:                DO 90 I = J + 1, N
230:                   Y( I ) = Y( I ) + TEMP1*AP( K )
231:                   TEMP2 = TEMP2 + AP( K )*X( I )
232:                   K = K + 1
233:    90          CONTINUE
234:                Y( J ) = Y( J ) + ALPHA*TEMP2
235:                KK = KK + ( N-J+1 )
236:   100       CONTINUE
237:          ELSE
238:             JX = KX
239:             JY = KY
240:             DO 120 J = 1, N
241:                TEMP1 = ALPHA*X( JX )
242:                TEMP2 = ZERO
243:                Y( JY ) = Y( JY ) + TEMP1*AP( KK )
244:                IX = JX
245:                IY = JY
246:                DO 110 K = KK + 1, KK + N - J
247:                   IX = IX + INCX
248:                   IY = IY + INCY
249:                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
250:                   TEMP2 = TEMP2 + AP( K )*X( IX )
251:   110          CONTINUE
252:                Y( JY ) = Y( JY ) + ALPHA*TEMP2
253:                JX = JX + INCX
254:                JY = JY + INCY
255:                KK = KK + ( N-J+1 )
256:   120       CONTINUE
257:          END IF
258:       END IF
259: *
260:       RETURN
261: *
262: *     End of ZSPMV
263: *
264:       END
265: