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