001:       SUBROUTINE ZSPR( UPLO, N, ALPHA, X, INCX, AP )
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, N
010:       COMPLEX*16         ALPHA
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX*16         AP( * ), X( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZSPR    performs the symmetric rank 1 operation
020: *
021: *     A := alpha*x*conjg( x' ) + A,
022: *
023: *  where alpha is a complex scalar, x is an n element vector and A is an
024: *  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: *  X        (input) COMPLEX*16 array, dimension at least
052: *           ( 1 + ( N - 1 )*abs( INCX ) ).
053: *           Before entry, the incremented array X must contain the N-
054: *           element vector x.
055: *           Unchanged on exit.
056: *
057: *  INCX     (input) INTEGER
058: *           On entry, INCX specifies the increment for the elements of
059: *           X. INCX must not be zero.
060: *           Unchanged on exit.
061: *
062: *  AP       (input/output) COMPLEX*16 array, dimension at least
063: *           ( ( N*( N + 1 ) )/2 ).
064: *           Before entry, with  UPLO = 'U' or 'u', the array AP must
065: *           contain the upper triangular part of the symmetric matrix
066: *           packed sequentially, column by column, so that AP( 1 )
067: *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
068: *           and a( 2, 2 ) respectively, and so on. On exit, the array
069: *           AP is overwritten by the upper triangular part of the
070: *           updated matrix.
071: *           Before entry, with UPLO = 'L' or 'l', the array AP must
072: *           contain the lower triangular part of the symmetric matrix
073: *           packed sequentially, column by column, so that AP( 1 )
074: *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
075: *           and a( 3, 1 ) respectively, and so on. On exit, the array
076: *           AP is overwritten by the lower triangular part of the
077: *           updated matrix.
078: *           Note that the imaginary parts of the diagonal elements need
079: *           not be set, they are assumed to be zero, and on exit they
080: *           are set to zero.
081: *
082: * =====================================================================
083: *
084: *     .. Parameters ..
085:       COMPLEX*16         ZERO
086:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
087: *     ..
088: *     .. Local Scalars ..
089:       INTEGER            I, INFO, IX, J, JX, K, KK, KX
090:       COMPLEX*16         TEMP
091: *     ..
092: *     .. External Functions ..
093:       LOGICAL            LSAME
094:       EXTERNAL           LSAME
095: *     ..
096: *     .. External Subroutines ..
097:       EXTERNAL           XERBLA
098: *     ..
099: *     .. Executable Statements ..
100: *
101: *     Test the input parameters.
102: *
103:       INFO = 0
104:       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
105:          INFO = 1
106:       ELSE IF( N.LT.0 ) THEN
107:          INFO = 2
108:       ELSE IF( INCX.EQ.0 ) THEN
109:          INFO = 5
110:       END IF
111:       IF( INFO.NE.0 ) THEN
112:          CALL XERBLA( 'ZSPR  ', INFO )
113:          RETURN
114:       END IF
115: *
116: *     Quick return if possible.
117: *
118:       IF( ( N.EQ.0 ) .OR. ( ALPHA.EQ.ZERO ) )
119:      $   RETURN
120: *
121: *     Set the start point in X if the increment is not unity.
122: *
123:       IF( INCX.LE.0 ) THEN
124:          KX = 1 - ( N-1 )*INCX
125:       ELSE IF( INCX.NE.1 ) THEN
126:          KX = 1
127:       END IF
128: *
129: *     Start the operations. In this version the elements of the array AP
130: *     are accessed sequentially with one pass through AP.
131: *
132:       KK = 1
133:       IF( LSAME( UPLO, 'U' ) ) THEN
134: *
135: *        Form  A  when upper triangle is stored in AP.
136: *
137:          IF( INCX.EQ.1 ) THEN
138:             DO 20 J = 1, N
139:                IF( X( J ).NE.ZERO ) THEN
140:                   TEMP = ALPHA*X( J )
141:                   K = KK
142:                   DO 10 I = 1, J - 1
143:                      AP( K ) = AP( K ) + X( I )*TEMP
144:                      K = K + 1
145:    10             CONTINUE
146:                   AP( KK+J-1 ) = AP( KK+J-1 ) + X( J )*TEMP
147:                ELSE
148:                   AP( KK+J-1 ) = AP( KK+J-1 )
149:                END IF
150:                KK = KK + J
151:    20       CONTINUE
152:          ELSE
153:             JX = KX
154:             DO 40 J = 1, N
155:                IF( X( JX ).NE.ZERO ) THEN
156:                   TEMP = ALPHA*X( JX )
157:                   IX = KX
158:                   DO 30 K = KK, KK + J - 2
159:                      AP( K ) = AP( K ) + X( IX )*TEMP
160:                      IX = IX + INCX
161:    30             CONTINUE
162:                   AP( KK+J-1 ) = AP( KK+J-1 ) + X( JX )*TEMP
163:                ELSE
164:                   AP( KK+J-1 ) = AP( KK+J-1 )
165:                END IF
166:                JX = JX + INCX
167:                KK = KK + J
168:    40       CONTINUE
169:          END IF
170:       ELSE
171: *
172: *        Form  A  when lower triangle is stored in AP.
173: *
174:          IF( INCX.EQ.1 ) THEN
175:             DO 60 J = 1, N
176:                IF( X( J ).NE.ZERO ) THEN
177:                   TEMP = ALPHA*X( J )
178:                   AP( KK ) = AP( KK ) + TEMP*X( J )
179:                   K = KK + 1
180:                   DO 50 I = J + 1, N
181:                      AP( K ) = AP( K ) + X( I )*TEMP
182:                      K = K + 1
183:    50             CONTINUE
184:                ELSE
185:                   AP( KK ) = AP( KK )
186:                END IF
187:                KK = KK + N - J + 1
188:    60       CONTINUE
189:          ELSE
190:             JX = KX
191:             DO 80 J = 1, N
192:                IF( X( JX ).NE.ZERO ) THEN
193:                   TEMP = ALPHA*X( JX )
194:                   AP( KK ) = AP( KK ) + TEMP*X( JX )
195:                   IX = JX
196:                   DO 70 K = KK + 1, KK + N - J
197:                      IX = IX + INCX
198:                      AP( K ) = AP( K ) + X( IX )*TEMP
199:    70             CONTINUE
200:                ELSE
201:                   AP( KK ) = AP( KK )
202:                END IF
203:                JX = JX + INCX
204:                KK = KK + N - J + 1
205:    80       CONTINUE
206:          END IF
207:       END IF
208: *
209:       RETURN
210: *
211: *     End of ZSPR
212: *
213:       END
214: