|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CSPR( UPLO, N, ALPHA, X, INCX, AP ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INCX, N 00011 COMPLEX ALPHA 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX AP( * ), X( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CSPR performs the symmetric rank 1 operation 00021 * 00022 * A := alpha*x*conjg( x' ) + A, 00023 * 00024 * where alpha is a complex scalar, x is an n element vector and A is an 00025 * n by n symmetric matrix, supplied in packed form. 00026 * 00027 * Arguments 00028 * ========== 00029 * 00030 * UPLO (input) CHARACTER*1 00031 * On entry, UPLO specifies whether the upper or lower 00032 * triangular part of the matrix A is supplied in the packed 00033 * array AP as follows: 00034 * 00035 * UPLO = 'U' or 'u' The upper triangular part of A is 00036 * supplied in AP. 00037 * 00038 * UPLO = 'L' or 'l' The lower triangular part of A is 00039 * supplied in AP. 00040 * 00041 * Unchanged on exit. 00042 * 00043 * N (input) INTEGER 00044 * On entry, N specifies the order of the matrix A. 00045 * N must be at least zero. 00046 * Unchanged on exit. 00047 * 00048 * ALPHA (input) COMPLEX 00049 * On entry, ALPHA specifies the scalar alpha. 00050 * Unchanged on exit. 00051 * 00052 * X (input) COMPLEX array, dimension at least 00053 * ( 1 + ( N - 1 )*abs( INCX ) ). 00054 * Before entry, the incremented array X must contain the N- 00055 * element vector x. 00056 * Unchanged on exit. 00057 * 00058 * INCX (input) INTEGER 00059 * On entry, INCX specifies the increment for the elements of 00060 * X. INCX must not be zero. 00061 * Unchanged on exit. 00062 * 00063 * AP (input/output) COMPLEX array, dimension at least 00064 * ( ( N*( N + 1 ) )/2 ). 00065 * Before entry, with UPLO = 'U' or 'u', the array AP must 00066 * contain the upper triangular part of the symmetric matrix 00067 * packed sequentially, column by column, so that AP( 1 ) 00068 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) 00069 * and a( 2, 2 ) respectively, and so on. On exit, the array 00070 * AP is overwritten by the upper triangular part of the 00071 * updated matrix. 00072 * Before entry, with UPLO = 'L' or 'l', the array AP must 00073 * contain the lower triangular part of the symmetric matrix 00074 * packed sequentially, column by column, so that AP( 1 ) 00075 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) 00076 * and a( 3, 1 ) respectively, and so on. On exit, the array 00077 * AP is overwritten by the lower triangular part of the 00078 * updated matrix. 00079 * Note that the imaginary parts of the diagonal elements need 00080 * not be set, they are assumed to be zero, and on exit they 00081 * are set to zero. 00082 * 00083 * ===================================================================== 00084 * 00085 * .. Parameters .. 00086 COMPLEX ZERO 00087 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) 00088 * .. 00089 * .. Local Scalars .. 00090 INTEGER I, INFO, IX, J, JX, K, KK, KX 00091 COMPLEX TEMP 00092 * .. 00093 * .. External Functions .. 00094 LOGICAL LSAME 00095 EXTERNAL LSAME 00096 * .. 00097 * .. External Subroutines .. 00098 EXTERNAL XERBLA 00099 * .. 00100 * .. Executable Statements .. 00101 * 00102 * Test the input parameters. 00103 * 00104 INFO = 0 00105 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00106 INFO = 1 00107 ELSE IF( N.LT.0 ) THEN 00108 INFO = 2 00109 ELSE IF( INCX.EQ.0 ) THEN 00110 INFO = 5 00111 END IF 00112 IF( INFO.NE.0 ) THEN 00113 CALL XERBLA( 'CSPR ', INFO ) 00114 RETURN 00115 END IF 00116 * 00117 * Quick return if possible. 00118 * 00119 IF( ( N.EQ.0 ) .OR. ( ALPHA.EQ.ZERO ) ) 00120 $ RETURN 00121 * 00122 * Set the start point in X if the increment is not unity. 00123 * 00124 IF( INCX.LE.0 ) THEN 00125 KX = 1 - ( N-1 )*INCX 00126 ELSE IF( INCX.NE.1 ) THEN 00127 KX = 1 00128 END IF 00129 * 00130 * Start the operations. In this version the elements of the array AP 00131 * are accessed sequentially with one pass through AP. 00132 * 00133 KK = 1 00134 IF( LSAME( UPLO, 'U' ) ) THEN 00135 * 00136 * Form A when upper triangle is stored in AP. 00137 * 00138 IF( INCX.EQ.1 ) THEN 00139 DO 20 J = 1, N 00140 IF( X( J ).NE.ZERO ) THEN 00141 TEMP = ALPHA*X( J ) 00142 K = KK 00143 DO 10 I = 1, J - 1 00144 AP( K ) = AP( K ) + X( I )*TEMP 00145 K = K + 1 00146 10 CONTINUE 00147 AP( KK+J-1 ) = AP( KK+J-1 ) + X( J )*TEMP 00148 ELSE 00149 AP( KK+J-1 ) = AP( KK+J-1 ) 00150 END IF 00151 KK = KK + J 00152 20 CONTINUE 00153 ELSE 00154 JX = KX 00155 DO 40 J = 1, N 00156 IF( X( JX ).NE.ZERO ) THEN 00157 TEMP = ALPHA*X( JX ) 00158 IX = KX 00159 DO 30 K = KK, KK + J - 2 00160 AP( K ) = AP( K ) + X( IX )*TEMP 00161 IX = IX + INCX 00162 30 CONTINUE 00163 AP( KK+J-1 ) = AP( KK+J-1 ) + X( JX )*TEMP 00164 ELSE 00165 AP( KK+J-1 ) = AP( KK+J-1 ) 00166 END IF 00167 JX = JX + INCX 00168 KK = KK + J 00169 40 CONTINUE 00170 END IF 00171 ELSE 00172 * 00173 * Form A when lower triangle is stored in AP. 00174 * 00175 IF( INCX.EQ.1 ) THEN 00176 DO 60 J = 1, N 00177 IF( X( J ).NE.ZERO ) THEN 00178 TEMP = ALPHA*X( J ) 00179 AP( KK ) = AP( KK ) + TEMP*X( J ) 00180 K = KK + 1 00181 DO 50 I = J + 1, N 00182 AP( K ) = AP( K ) + X( I )*TEMP 00183 K = K + 1 00184 50 CONTINUE 00185 ELSE 00186 AP( KK ) = AP( KK ) 00187 END IF 00188 KK = KK + N - J + 1 00189 60 CONTINUE 00190 ELSE 00191 JX = KX 00192 DO 80 J = 1, N 00193 IF( X( JX ).NE.ZERO ) THEN 00194 TEMP = ALPHA*X( JX ) 00195 AP( KK ) = AP( KK ) + TEMP*X( JX ) 00196 IX = JX 00197 DO 70 K = KK + 1, KK + N - J 00198 IX = IX + INCX 00199 AP( K ) = AP( K ) + X( IX )*TEMP 00200 70 CONTINUE 00201 ELSE 00202 AP( KK ) = AP( KK ) 00203 END IF 00204 JX = JX + INCX 00205 KK = KK + N - J + 1 00206 80 CONTINUE 00207 END IF 00208 END IF 00209 * 00210 RETURN 00211 * 00212 * End of CSPR 00213 * 00214 END
1.7.3 
