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