1       SUBROUTINE ZSPMV( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
  2 *
  3 *  -- LAPACK auxiliary routine (version 3.2) --
  4 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  5 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          UPLO
 10       INTEGER            INCX, INCY, N
 11       COMPLEX*16         ALPHA, BETA
 12 *     ..
 13 *     .. Array Arguments ..
 14       COMPLEX*16         AP( * ), X( * ), Y( * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  ZSPMV  performs the matrix-vector operation
 21 *
 22 *     y := alpha*A*x + beta*y,
 23 *
 24 *  where alpha and beta are scalars, x and y are n element vectors and
 25 *  A is an n by n symmetric matrix, supplied in packed form.
 26 *
 27 *  Arguments
 28 *  ==========
 29 *
 30 *  UPLO     (input) CHARACTER*1
 31 *           On entry, UPLO specifies whether the upper or lower
 32 *           triangular part of the matrix A is supplied in the packed
 33 *           array AP as follows:
 34 *
 35 *              UPLO = 'U' or 'u'   The upper triangular part of A is
 36 *                                  supplied in AP.
 37 *
 38 *              UPLO = 'L' or 'l'   The lower triangular part of A is
 39 *                                  supplied in AP.
 40 *
 41 *           Unchanged on exit.
 42 *
 43 *  N        (input) INTEGER
 44 *           On entry, N specifies the order of the matrix A.
 45 *           N must be at least zero.
 46 *           Unchanged on exit.
 47 *
 48 *  ALPHA    (input) COMPLEX*16
 49 *           On entry, ALPHA specifies the scalar alpha.
 50 *           Unchanged on exit.
 51 *
 52 *  AP       (input) COMPLEX*16 array, dimension at least
 53 *           ( ( N*( N + 1 ) )/2 ).
 54 *           Before entry, with UPLO = 'U' or 'u', the array AP must
 55 *           contain the upper triangular part of the symmetric matrix
 56 *           packed sequentially, column by column, so that AP( 1 )
 57 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
 58 *           and a( 2, 2 ) respectively, and so on.
 59 *           Before entry, with UPLO = 'L' or 'l', the array AP must
 60 *           contain the lower triangular part of the symmetric matrix
 61 *           packed sequentially, column by column, so that AP( 1 )
 62 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
 63 *           and a( 3, 1 ) respectively, and so on.
 64 *           Unchanged on exit.
 65 *
 66 *  X        (input) COMPLEX*16 array, dimension at least
 67 *           ( 1 + ( N - 1 )*abs( INCX ) ).
 68 *           Before entry, the incremented array X must contain the N-
 69 *           element vector x.
 70 *           Unchanged on exit.
 71 *
 72 *  INCX     (input) INTEGER
 73 *           On entry, INCX specifies the increment for the elements of
 74 *           X. INCX must not be zero.
 75 *           Unchanged on exit.
 76 *
 77 *  BETA     (input) COMPLEX*16
 78 *           On entry, BETA specifies the scalar beta. When BETA is
 79 *           supplied as zero then Y need not be set on input.
 80 *           Unchanged on exit.
 81 *
 82 *  Y        (input/output) COMPLEX*16 array, dimension at least
 83 *           ( 1 + ( N - 1 )*abs( INCY ) ).
 84 *           Before entry, the incremented array Y must contain the n
 85 *           element vector y. On exit, Y is overwritten by the updated
 86 *           vector y.
 87 *
 88 *  INCY     (input) INTEGER
 89 *           On entry, INCY specifies the increment for the elements of
 90 *           Y. INCY must not be zero.
 91 *           Unchanged on exit.
 92 *
 93 * =====================================================================
 94 *
 95 *     .. Parameters ..
 96       COMPLEX*16         ONE
 97       PARAMETER          ( ONE = ( 1.0D+00.0D+0 ) )
 98       COMPLEX*16         ZERO
 99       PARAMETER          ( ZERO = ( 0.0D+00.0D+0 ) )
100 *     ..
101 *     .. Local Scalars ..
102       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
103       COMPLEX*16         TEMP1, TEMP2
104 *     ..
105 *     .. External Functions ..
106       LOGICAL            LSAME
107       EXTERNAL           LSAME
108 *     ..
109 *     .. External Subroutines ..
110       EXTERNAL           XERBLA
111 *     ..
112 *     .. Executable Statements ..
113 *
114 *     Test the input parameters.
115 *
116       INFO = 0
117       IF.NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
118          INFO = 1
119       ELSE IF( N.LT.0 ) THEN
120          INFO = 2
121       ELSE IF( INCX.EQ.0 ) THEN
122          INFO = 6
123       ELSE IF( INCY.EQ.0 ) THEN
124          INFO = 9
125       END IF
126       IF( INFO.NE.0 ) THEN
127          CALL XERBLA( 'ZSPMV ', INFO )
128          RETURN
129       END IF
130 *
131 *     Quick return if possible.
132 *
133       IF( ( N.EQ.0 ) .OR. ( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ONE ) ) )
134      $   RETURN
135 *
136 *     Set up the start points in  X  and  Y.
137 *
138       IF( INCX.GT.0 ) THEN
139          KX = 1
140       ELSE
141          KX = 1 - ( N-1 )*INCX
142       END IF
143       IF( INCY.GT.0 ) THEN
144          KY = 1
145       ELSE
146          KY = 1 - ( N-1 )*INCY
147       END IF
148 *
149 *     Start the operations. In this version the elements of the array AP
150 *     are accessed sequentially with one pass through AP.
151 *
152 *     First form  y := beta*y.
153 *
154       IF( BETA.NE.ONE ) THEN
155          IF( INCY.EQ.1 ) THEN
156             IF( BETA.EQ.ZERO ) THEN
157                DO 10 I = 1, N
158                   Y( I ) = ZERO
159    10          CONTINUE
160             ELSE
161                DO 20 I = 1, N
162                   Y( I ) = BETA*Y( I )
163    20          CONTINUE
164             END IF
165          ELSE
166             IY = KY
167             IF( BETA.EQ.ZERO ) THEN
168                DO 30 I = 1, N
169                   Y( IY ) = ZERO
170                   IY = IY + INCY
171    30          CONTINUE
172             ELSE
173                DO 40 I = 1, N
174                   Y( IY ) = BETA*Y( IY )
175                   IY = IY + INCY
176    40          CONTINUE
177             END IF
178          END IF
179       END IF
180       IF( ALPHA.EQ.ZERO )
181      $   RETURN
182       KK = 1
183       IF( LSAME( UPLO, 'U' ) ) THEN
184 *
185 *        Form  y  when AP contains the upper triangle.
186 *
187          IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN
188             DO 60 J = 1, N
189                TEMP1 = ALPHA*X( J )
190                TEMP2 = ZERO
191                K = KK
192                DO 50 I = 1, J - 1
193                   Y( I ) = Y( I ) + TEMP1*AP( K )
194                   TEMP2 = TEMP2 + AP( K )*X( I )
195                   K = K + 1
196    50          CONTINUE
197                Y( J ) = Y( J ) + TEMP1*AP( KK+J-1 ) + ALPHA*TEMP2
198                KK = KK + J
199    60       CONTINUE
200          ELSE
201             JX = KX
202             JY = KY
203             DO 80 J = 1, N
204                TEMP1 = ALPHA*X( JX )
205                TEMP2 = ZERO
206                IX = KX
207                IY = KY
208                DO 70 K = KK, KK + J - 2
209                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
210                   TEMP2 = TEMP2 + AP( K )*X( IX )
211                   IX = IX + INCX
212                   IY = IY + INCY
213    70          CONTINUE
214                Y( JY ) = Y( JY ) + TEMP1*AP( KK+J-1 ) + ALPHA*TEMP2
215                JX = JX + INCX
216                JY = JY + INCY
217                KK = KK + J
218    80       CONTINUE
219          END IF
220       ELSE
221 *
222 *        Form  y  when AP contains the lower triangle.
223 *
224          IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN
225             DO 100 J = 1, N
226                TEMP1 = ALPHA*X( J )
227                TEMP2 = ZERO
228                Y( J ) = Y( J ) + TEMP1*AP( KK )
229                K = KK + 1
230                DO 90 I = J + 1, N
231                   Y( I ) = Y( I ) + TEMP1*AP( K )
232                   TEMP2 = TEMP2 + AP( K )*X( I )
233                   K = K + 1
234    90          CONTINUE
235                Y( J ) = Y( J ) + ALPHA*TEMP2
236                KK = KK + ( N-J+1 )
237   100       CONTINUE
238          ELSE
239             JX = KX
240             JY = KY
241             DO 120 J = 1, N
242                TEMP1 = ALPHA*X( JX )
243                TEMP2 = ZERO
244                Y( JY ) = Y( JY ) + TEMP1*AP( KK )
245                IX = JX
246                IY = JY
247                DO 110 K = KK + 1, KK + N - J
248                   IX = IX + INCX
249                   IY = IY + INCY
250                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
251                   TEMP2 = TEMP2 + AP( K )*X( IX )
252   110          CONTINUE
253                Y( JY ) = Y( JY ) + ALPHA*TEMP2
254                JX = JX + INCX
255                JY = JY + INCY
256                KK = KK + ( N-J+1 )
257   120       CONTINUE
258          END IF
259       END IF
260 *
261       RETURN
262 *
263 *     End of ZSPMV
264 *
265       END