|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CHPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO ) 00002 * 00003 * -- LAPACK 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 INFO, LDB, N, NRHS 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IPIV( * ) 00014 COMPLEX AP( * ), B( LDB, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CHPTRS solves a system of linear equations A*X = B with a complex 00021 * Hermitian matrix A stored in packed format using the factorization 00022 * A = U*D*U**H or A = L*D*L**H computed by CHPTRF. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * UPLO (input) CHARACTER*1 00028 * Specifies whether the details of the factorization are stored 00029 * as an upper or lower triangular matrix. 00030 * = 'U': Upper triangular, form is A = U*D*U**H; 00031 * = 'L': Lower triangular, form is A = L*D*L**H. 00032 * 00033 * N (input) INTEGER 00034 * The order of the matrix A. N >= 0. 00035 * 00036 * NRHS (input) INTEGER 00037 * The number of right hand sides, i.e., the number of columns 00038 * of the matrix B. NRHS >= 0. 00039 * 00040 * AP (input) COMPLEX array, dimension (N*(N+1)/2) 00041 * The block diagonal matrix D and the multipliers used to 00042 * obtain the factor U or L as computed by CHPTRF, stored as a 00043 * packed triangular matrix. 00044 * 00045 * IPIV (input) INTEGER array, dimension (N) 00046 * Details of the interchanges and the block structure of D 00047 * as determined by CHPTRF. 00048 * 00049 * B (input/output) COMPLEX array, dimension (LDB,NRHS) 00050 * On entry, the right hand side matrix B. 00051 * On exit, the solution matrix X. 00052 * 00053 * LDB (input) INTEGER 00054 * The leading dimension of the array B. LDB >= max(1,N). 00055 * 00056 * INFO (output) INTEGER 00057 * = 0: successful exit 00058 * < 0: if INFO = -i, the i-th argument had an illegal value 00059 * 00060 * ===================================================================== 00061 * 00062 * .. Parameters .. 00063 COMPLEX ONE 00064 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) 00065 * .. 00066 * .. Local Scalars .. 00067 LOGICAL UPPER 00068 INTEGER J, K, KC, KP 00069 REAL S 00070 COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM 00071 * .. 00072 * .. External Functions .. 00073 LOGICAL LSAME 00074 EXTERNAL LSAME 00075 * .. 00076 * .. External Subroutines .. 00077 EXTERNAL CGEMV, CGERU, CLACGV, CSSCAL, CSWAP, XERBLA 00078 * .. 00079 * .. Intrinsic Functions .. 00080 INTRINSIC CONJG, MAX, REAL 00081 * .. 00082 * .. Executable Statements .. 00083 * 00084 INFO = 0 00085 UPPER = LSAME( UPLO, 'U' ) 00086 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00087 INFO = -1 00088 ELSE IF( N.LT.0 ) THEN 00089 INFO = -2 00090 ELSE IF( NRHS.LT.0 ) THEN 00091 INFO = -3 00092 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00093 INFO = -7 00094 END IF 00095 IF( INFO.NE.0 ) THEN 00096 CALL XERBLA( 'CHPTRS', -INFO ) 00097 RETURN 00098 END IF 00099 * 00100 * Quick return if possible 00101 * 00102 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00103 $ RETURN 00104 * 00105 IF( UPPER ) THEN 00106 * 00107 * Solve A*X = B, where A = U*D*U'. 00108 * 00109 * First solve U*D*X = B, overwriting B with X. 00110 * 00111 * K is the main loop index, decreasing from N to 1 in steps of 00112 * 1 or 2, depending on the size of the diagonal blocks. 00113 * 00114 K = N 00115 KC = N*( N+1 ) / 2 + 1 00116 10 CONTINUE 00117 * 00118 * If K < 1, exit from loop. 00119 * 00120 IF( K.LT.1 ) 00121 $ GO TO 30 00122 * 00123 KC = KC - K 00124 IF( IPIV( K ).GT.0 ) THEN 00125 * 00126 * 1 x 1 diagonal block 00127 * 00128 * Interchange rows K and IPIV(K). 00129 * 00130 KP = IPIV( K ) 00131 IF( KP.NE.K ) 00132 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00133 * 00134 * Multiply by inv(U(K)), where U(K) is the transformation 00135 * stored in column K of A. 00136 * 00137 CALL CGERU( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB, 00138 $ B( 1, 1 ), LDB ) 00139 * 00140 * Multiply by the inverse of the diagonal block. 00141 * 00142 S = REAL( ONE ) / REAL( AP( KC+K-1 ) ) 00143 CALL CSSCAL( NRHS, S, B( K, 1 ), LDB ) 00144 K = K - 1 00145 ELSE 00146 * 00147 * 2 x 2 diagonal block 00148 * 00149 * Interchange rows K-1 and -IPIV(K). 00150 * 00151 KP = -IPIV( K ) 00152 IF( KP.NE.K-1 ) 00153 $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 00154 * 00155 * Multiply by inv(U(K)), where U(K) is the transformation 00156 * stored in columns K-1 and K of A. 00157 * 00158 CALL CGERU( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB, 00159 $ B( 1, 1 ), LDB ) 00160 CALL CGERU( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1, 00161 $ B( K-1, 1 ), LDB, B( 1, 1 ), LDB ) 00162 * 00163 * Multiply by the inverse of the diagonal block. 00164 * 00165 AKM1K = AP( KC+K-2 ) 00166 AKM1 = AP( KC-1 ) / AKM1K 00167 AK = AP( KC+K-1 ) / CONJG( AKM1K ) 00168 DENOM = AKM1*AK - ONE 00169 DO 20 J = 1, NRHS 00170 BKM1 = B( K-1, J ) / AKM1K 00171 BK = B( K, J ) / CONJG( AKM1K ) 00172 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM 00173 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM 00174 20 CONTINUE 00175 KC = KC - K + 1 00176 K = K - 2 00177 END IF 00178 * 00179 GO TO 10 00180 30 CONTINUE 00181 * 00182 * Next solve U'*X = B, overwriting B with X. 00183 * 00184 * K is the main loop index, increasing from 1 to N in steps of 00185 * 1 or 2, depending on the size of the diagonal blocks. 00186 * 00187 K = 1 00188 KC = 1 00189 40 CONTINUE 00190 * 00191 * If K > N, exit from loop. 00192 * 00193 IF( K.GT.N ) 00194 $ GO TO 50 00195 * 00196 IF( IPIV( K ).GT.0 ) THEN 00197 * 00198 * 1 x 1 diagonal block 00199 * 00200 * Multiply by inv(U'(K)), where U(K) is the transformation 00201 * stored in column K of A. 00202 * 00203 IF( K.GT.1 ) THEN 00204 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00205 CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B, 00206 $ LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB ) 00207 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00208 END IF 00209 * 00210 * Interchange rows K and IPIV(K). 00211 * 00212 KP = IPIV( K ) 00213 IF( KP.NE.K ) 00214 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00215 KC = KC + K 00216 K = K + 1 00217 ELSE 00218 * 00219 * 2 x 2 diagonal block 00220 * 00221 * Multiply by inv(U'(K+1)), where U(K+1) is the transformation 00222 * stored in columns K and K+1 of A. 00223 * 00224 IF( K.GT.1 ) THEN 00225 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00226 CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B, 00227 $ LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB ) 00228 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00229 * 00230 CALL CLACGV( NRHS, B( K+1, 1 ), LDB ) 00231 CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B, 00232 $ LDB, AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB ) 00233 CALL CLACGV( NRHS, B( K+1, 1 ), LDB ) 00234 END IF 00235 * 00236 * Interchange rows K and -IPIV(K). 00237 * 00238 KP = -IPIV( K ) 00239 IF( KP.NE.K ) 00240 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00241 KC = KC + 2*K + 1 00242 K = K + 2 00243 END IF 00244 * 00245 GO TO 40 00246 50 CONTINUE 00247 * 00248 ELSE 00249 * 00250 * Solve A*X = B, where A = L*D*L'. 00251 * 00252 * First solve L*D*X = B, overwriting B with X. 00253 * 00254 * K is the main loop index, increasing from 1 to N in steps of 00255 * 1 or 2, depending on the size of the diagonal blocks. 00256 * 00257 K = 1 00258 KC = 1 00259 60 CONTINUE 00260 * 00261 * If K > N, exit from loop. 00262 * 00263 IF( K.GT.N ) 00264 $ GO TO 80 00265 * 00266 IF( IPIV( K ).GT.0 ) THEN 00267 * 00268 * 1 x 1 diagonal block 00269 * 00270 * Interchange rows K and IPIV(K). 00271 * 00272 KP = IPIV( K ) 00273 IF( KP.NE.K ) 00274 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00275 * 00276 * Multiply by inv(L(K)), where L(K) is the transformation 00277 * stored in column K of A. 00278 * 00279 IF( K.LT.N ) 00280 $ CALL CGERU( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ), 00281 $ LDB, B( K+1, 1 ), LDB ) 00282 * 00283 * Multiply by the inverse of the diagonal block. 00284 * 00285 S = REAL( ONE ) / REAL( AP( KC ) ) 00286 CALL CSSCAL( NRHS, S, B( K, 1 ), LDB ) 00287 KC = KC + N - K + 1 00288 K = K + 1 00289 ELSE 00290 * 00291 * 2 x 2 diagonal block 00292 * 00293 * Interchange rows K+1 and -IPIV(K). 00294 * 00295 KP = -IPIV( K ) 00296 IF( KP.NE.K+1 ) 00297 $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 00298 * 00299 * Multiply by inv(L(K)), where L(K) is the transformation 00300 * stored in columns K and K+1 of A. 00301 * 00302 IF( K.LT.N-1 ) THEN 00303 CALL CGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ), 00304 $ LDB, B( K+2, 1 ), LDB ) 00305 CALL CGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1, 00306 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB ) 00307 END IF 00308 * 00309 * Multiply by the inverse of the diagonal block. 00310 * 00311 AKM1K = AP( KC+1 ) 00312 AKM1 = AP( KC ) / CONJG( AKM1K ) 00313 AK = AP( KC+N-K+1 ) / AKM1K 00314 DENOM = AKM1*AK - ONE 00315 DO 70 J = 1, NRHS 00316 BKM1 = B( K, J ) / CONJG( AKM1K ) 00317 BK = B( K+1, J ) / AKM1K 00318 B( K, J ) = ( AK*BKM1-BK ) / DENOM 00319 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 00320 70 CONTINUE 00321 KC = KC + 2*( N-K ) + 1 00322 K = K + 2 00323 END IF 00324 * 00325 GO TO 60 00326 80 CONTINUE 00327 * 00328 * Next solve L'*X = B, overwriting B with X. 00329 * 00330 * K is the main loop index, decreasing from N to 1 in steps of 00331 * 1 or 2, depending on the size of the diagonal blocks. 00332 * 00333 K = N 00334 KC = N*( N+1 ) / 2 + 1 00335 90 CONTINUE 00336 * 00337 * If K < 1, exit from loop. 00338 * 00339 IF( K.LT.1 ) 00340 $ GO TO 100 00341 * 00342 KC = KC - ( N-K+1 ) 00343 IF( IPIV( K ).GT.0 ) THEN 00344 * 00345 * 1 x 1 diagonal block 00346 * 00347 * Multiply by inv(L'(K)), where L(K) is the transformation 00348 * stored in column K of A. 00349 * 00350 IF( K.LT.N ) THEN 00351 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00352 CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE, 00353 $ B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE, 00354 $ B( K, 1 ), LDB ) 00355 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00356 END IF 00357 * 00358 * Interchange rows K and IPIV(K). 00359 * 00360 KP = IPIV( K ) 00361 IF( KP.NE.K ) 00362 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00363 K = K - 1 00364 ELSE 00365 * 00366 * 2 x 2 diagonal block 00367 * 00368 * Multiply by inv(L'(K-1)), where L(K-1) is the transformation 00369 * stored in columns K-1 and K of A. 00370 * 00371 IF( K.LT.N ) THEN 00372 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00373 CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE, 00374 $ B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE, 00375 $ B( K, 1 ), LDB ) 00376 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 00377 * 00378 CALL CLACGV( NRHS, B( K-1, 1 ), LDB ) 00379 CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE, 00380 $ B( K+1, 1 ), LDB, AP( KC-( N-K ) ), 1, ONE, 00381 $ B( K-1, 1 ), LDB ) 00382 CALL CLACGV( NRHS, B( K-1, 1 ), LDB ) 00383 END IF 00384 * 00385 * Interchange rows K and -IPIV(K). 00386 * 00387 KP = -IPIV( K ) 00388 IF( KP.NE.K ) 00389 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00390 KC = KC - ( N-K+2 ) 00391 K = K - 2 00392 END IF 00393 * 00394 GO TO 90 00395 100 CONTINUE 00396 END IF 00397 * 00398 RETURN 00399 * 00400 * End of CHPTRS 00401 * 00402 END
1.7.3 
