![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHETRS2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHETRS2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrs2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrs2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrs2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 00022 * WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO 00026 * INTEGER INFO, LDA, LDB, N, NRHS 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IPIV( * ) 00030 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> ZHETRS2 solves a system of linear equations A*X = B with a complex 00040 *> Hermitian matrix A using the factorization A = U*D*U**H or 00041 *> A = L*D*L**H computed by ZHETRF and converted by ZSYCONV. 00042 *> \endverbatim 00043 * 00044 * Arguments: 00045 * ========== 00046 * 00047 *> \param[in] UPLO 00048 *> \verbatim 00049 *> UPLO is CHARACTER*1 00050 *> Specifies whether the details of the factorization are stored 00051 *> as an upper or lower triangular matrix. 00052 *> = 'U': Upper triangular, form is A = U*D*U**H; 00053 *> = 'L': Lower triangular, form is A = L*D*L**H. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] N 00057 *> \verbatim 00058 *> N is INTEGER 00059 *> The order of the matrix A. N >= 0. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] NRHS 00063 *> \verbatim 00064 *> NRHS is INTEGER 00065 *> The number of right hand sides, i.e., the number of columns 00066 *> of the matrix B. NRHS >= 0. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] A 00070 *> \verbatim 00071 *> A is COMPLEX*16 array, dimension (LDA,N) 00072 *> The block diagonal matrix D and the multipliers used to 00073 *> obtain the factor U or L as computed by ZHETRF. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] LDA 00077 *> \verbatim 00078 *> LDA is INTEGER 00079 *> The leading dimension of the array A. LDA >= max(1,N). 00080 *> \endverbatim 00081 *> 00082 *> \param[in] IPIV 00083 *> \verbatim 00084 *> IPIV is INTEGER array, dimension (N) 00085 *> Details of the interchanges and the block structure of D 00086 *> as determined by ZHETRF. 00087 *> \endverbatim 00088 *> 00089 *> \param[in,out] B 00090 *> \verbatim 00091 *> B is COMPLEX*16 array, dimension (LDB,NRHS) 00092 *> On entry, the right hand side matrix B. 00093 *> On exit, the solution matrix X. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] LDB 00097 *> \verbatim 00098 *> LDB is INTEGER 00099 *> The leading dimension of the array B. LDB >= max(1,N). 00100 *> \endverbatim 00101 *> 00102 *> \param[out] WORK 00103 *> \verbatim 00104 *> WORK is REAL array, dimension (N) 00105 *> \endverbatim 00106 *> 00107 *> \param[out] INFO 00108 *> \verbatim 00109 *> INFO is INTEGER 00110 *> = 0: successful exit 00111 *> < 0: if INFO = -i, the i-th argument had an illegal value 00112 *> \endverbatim 00113 * 00114 * Authors: 00115 * ======== 00116 * 00117 *> \author Univ. of Tennessee 00118 *> \author Univ. of California Berkeley 00119 *> \author Univ. of Colorado Denver 00120 *> \author NAG Ltd. 00121 * 00122 *> \date November 2011 00123 * 00124 *> \ingroup complex16HEcomputational 00125 * 00126 * ===================================================================== 00127 SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 00128 $ WORK, INFO ) 00129 * 00130 * -- LAPACK computational routine (version 3.4.0) -- 00131 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00133 * November 2011 00134 * 00135 * .. Scalar Arguments .. 00136 CHARACTER UPLO 00137 INTEGER INFO, LDA, LDB, N, NRHS 00138 * .. 00139 * .. Array Arguments .. 00140 INTEGER IPIV( * ) 00141 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00142 * .. 00143 * 00144 * ===================================================================== 00145 * 00146 * .. Parameters .. 00147 COMPLEX*16 ONE 00148 PARAMETER ( ONE = (1.0D+0,0.0D+0) ) 00149 * .. 00150 * .. Local Scalars .. 00151 LOGICAL UPPER 00152 INTEGER I, IINFO, J, K, KP 00153 DOUBLE PRECISION S 00154 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM 00155 * .. 00156 * .. External Functions .. 00157 LOGICAL LSAME 00158 EXTERNAL LSAME 00159 * .. 00160 * .. External Subroutines .. 00161 EXTERNAL ZLACGV, ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA 00162 * .. 00163 * .. Intrinsic Functions .. 00164 INTRINSIC DBLE, DCONJG, MAX 00165 * .. 00166 * .. Executable Statements .. 00167 * 00168 INFO = 0 00169 UPPER = LSAME( UPLO, 'U' ) 00170 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00171 INFO = -1 00172 ELSE IF( N.LT.0 ) THEN 00173 INFO = -2 00174 ELSE IF( NRHS.LT.0 ) THEN 00175 INFO = -3 00176 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00177 INFO = -5 00178 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00179 INFO = -8 00180 END IF 00181 IF( INFO.NE.0 ) THEN 00182 CALL XERBLA( 'ZHETRS2', -INFO ) 00183 RETURN 00184 END IF 00185 * 00186 * Quick return if possible 00187 * 00188 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00189 $ RETURN 00190 * 00191 * Convert A 00192 * 00193 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 00194 * 00195 IF( UPPER ) THEN 00196 * 00197 * Solve A*X = B, where A = U*D*U**H. 00198 * 00199 * P**T * B 00200 K=N 00201 DO WHILE ( K .GE. 1 ) 00202 IF( IPIV( K ).GT.0 ) THEN 00203 * 1 x 1 diagonal block 00204 * Interchange rows K and IPIV(K). 00205 KP = IPIV( K ) 00206 IF( KP.NE.K ) 00207 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00208 K=K-1 00209 ELSE 00210 * 2 x 2 diagonal block 00211 * Interchange rows K-1 and -IPIV(K). 00212 KP = -IPIV( K ) 00213 IF( KP.EQ.-IPIV( K-1 ) ) 00214 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 00215 K=K-2 00216 END IF 00217 END DO 00218 * 00219 * Compute (U \P**T * B) -> B [ (U \P**T * B) ] 00220 * 00221 CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB) 00222 * 00223 * Compute D \ B -> B [ D \ (U \P**T * B) ] 00224 * 00225 I=N 00226 DO WHILE ( I .GE. 1 ) 00227 IF( IPIV(I) .GT. 0 ) THEN 00228 S = DBLE( ONE ) / DBLE( A( I, I ) ) 00229 CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 00230 ELSEIF ( I .GT. 1) THEN 00231 IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN 00232 AKM1K = WORK(I) 00233 AKM1 = A( I-1, I-1 ) / AKM1K 00234 AK = A( I, I ) / DCONJG( AKM1K ) 00235 DENOM = AKM1*AK - ONE 00236 DO 15 J = 1, NRHS 00237 BKM1 = B( I-1, J ) / AKM1K 00238 BK = B( I, J ) / DCONJG( AKM1K ) 00239 B( I-1, J ) = ( AK*BKM1-BK ) / DENOM 00240 B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM 00241 15 CONTINUE 00242 I = I - 1 00243 ENDIF 00244 ENDIF 00245 I = I - 1 00246 END DO 00247 * 00248 * Compute (U**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ] 00249 * 00250 CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,LDA,B,LDB) 00251 * 00252 * P * B [ P * (U**H \ (D \ (U \P**T * B) )) ] 00253 * 00254 K=1 00255 DO WHILE ( K .LE. N ) 00256 IF( IPIV( K ).GT.0 ) THEN 00257 * 1 x 1 diagonal block 00258 * Interchange rows K and IPIV(K). 00259 KP = IPIV( K ) 00260 IF( KP.NE.K ) 00261 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00262 K=K+1 00263 ELSE 00264 * 2 x 2 diagonal block 00265 * Interchange rows K-1 and -IPIV(K). 00266 KP = -IPIV( K ) 00267 IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) ) 00268 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00269 K=K+2 00270 ENDIF 00271 END DO 00272 * 00273 ELSE 00274 * 00275 * Solve A*X = B, where A = L*D*L**H. 00276 * 00277 * P**T * B 00278 K=1 00279 DO WHILE ( K .LE. N ) 00280 IF( IPIV( K ).GT.0 ) THEN 00281 * 1 x 1 diagonal block 00282 * Interchange rows K and IPIV(K). 00283 KP = IPIV( K ) 00284 IF( KP.NE.K ) 00285 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00286 K=K+1 00287 ELSE 00288 * 2 x 2 diagonal block 00289 * Interchange rows K and -IPIV(K+1). 00290 KP = -IPIV( K+1 ) 00291 IF( KP.EQ.-IPIV( K ) ) 00292 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 00293 K=K+2 00294 ENDIF 00295 END DO 00296 * 00297 * Compute (L \P**T * B) -> B [ (L \P**T * B) ] 00298 * 00299 CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB) 00300 * 00301 * Compute D \ B -> B [ D \ (L \P**T * B) ] 00302 * 00303 I=1 00304 DO WHILE ( I .LE. N ) 00305 IF( IPIV(I) .GT. 0 ) THEN 00306 S = DBLE( ONE ) / DBLE( A( I, I ) ) 00307 CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 00308 ELSE 00309 AKM1K = WORK(I) 00310 AKM1 = A( I, I ) / DCONJG( AKM1K ) 00311 AK = A( I+1, I+1 ) / AKM1K 00312 DENOM = AKM1*AK - ONE 00313 DO 25 J = 1, NRHS 00314 BKM1 = B( I, J ) / DCONJG( AKM1K ) 00315 BK = B( I+1, J ) / AKM1K 00316 B( I, J ) = ( AK*BKM1-BK ) / DENOM 00317 B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 00318 25 CONTINUE 00319 I = I + 1 00320 ENDIF 00321 I = I + 1 00322 END DO 00323 * 00324 * Compute (L**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ] 00325 * 00326 CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,LDA,B,LDB) 00327 * 00328 * P * B [ P * (L**H \ (D \ (L \P**T * B) )) ] 00329 * 00330 K=N 00331 DO WHILE ( K .GE. 1 ) 00332 IF( IPIV( K ).GT.0 ) THEN 00333 * 1 x 1 diagonal block 00334 * Interchange rows K and IPIV(K). 00335 KP = IPIV( K ) 00336 IF( KP.NE.K ) 00337 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00338 K=K-1 00339 ELSE 00340 * 2 x 2 diagonal block 00341 * Interchange rows K-1 and -IPIV(K). 00342 KP = -IPIV( K ) 00343 IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) ) 00344 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00345 K=K-2 00346 ENDIF 00347 END DO 00348 * 00349 END IF 00350 * 00351 * Revert A 00352 * 00353 CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) 00354 * 00355 RETURN 00356 * 00357 * End of ZHETRS2 00358 * 00359 END