|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 00002 $ WORK, INFO ) 00003 * 00004 * -- LAPACK PROTOTYPE routine (version 3.3.0) -- 00005 * 00006 * -- Written by Julie Langou of the Univ. of TN -- 00007 * November 2010 00008 * 00009 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00010 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00011 * 00012 * .. Scalar Arguments .. 00013 CHARACTER UPLO 00014 INTEGER INFO, LDA, LDB, N, NRHS 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IPIV( * ) 00018 DOUBLE COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * ZHETRS2 solves a system of linear equations A*X = B with a real 00025 * Hermitian matrix A using the factorization A = U*D*U**T or 00026 * A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * UPLO (input) CHARACTER*1 00032 * Specifies whether the details of the factorization are stored 00033 * as an upper or lower triangular matrix. 00034 * = 'U': Upper triangular, form is A = U*D*U**H; 00035 * = 'L': Lower triangular, form is A = L*D*L**H. 00036 * 00037 * N (input) INTEGER 00038 * The order of the matrix A. N >= 0. 00039 * 00040 * NRHS (input) INTEGER 00041 * The number of right hand sides, i.e., the number of columns 00042 * of the matrix B. NRHS >= 0. 00043 * 00044 * A (input) DOUBLE COMPLEX array, dimension (LDA,N) 00045 * The block diagonal matrix D and the multipliers used to 00046 * obtain the factor U or L as computed by ZHETRF. 00047 * 00048 * LDA (input) INTEGER 00049 * The leading dimension of the array A. LDA >= max(1,N). 00050 * 00051 * IPIV (input) INTEGER array, dimension (N) 00052 * Details of the interchanges and the block structure of D 00053 * as determined by ZHETRF. 00054 * 00055 * B (input/output) DOUBLE COMPLEX array, dimension (LDB,NRHS) 00056 * On entry, the right hand side matrix B. 00057 * On exit, the solution matrix X. 00058 * 00059 * LDB (input) INTEGER 00060 * The leading dimension of the array B. LDB >= max(1,N). 00061 * 00062 * WORK (workspace) REAL array, dimension (N) 00063 * 00064 * INFO (output) INTEGER 00065 * = 0: successful exit 00066 * < 0: if INFO = -i, the i-th argument had an illegal value 00067 * 00068 * ===================================================================== 00069 * 00070 * .. Parameters .. 00071 DOUBLE COMPLEX ONE 00072 PARAMETER ( ONE = (1.0D+0,0.0D+0) ) 00073 * .. 00074 * .. Local Scalars .. 00075 LOGICAL UPPER 00076 INTEGER I, IINFO, J, K, KP 00077 DOUBLE PRECISION S 00078 DOUBLE COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM 00079 * .. 00080 * .. External Functions .. 00081 LOGICAL LSAME 00082 EXTERNAL LSAME 00083 * .. 00084 * .. External Subroutines .. 00085 EXTERNAL ZLACGV, ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA 00086 * .. 00087 * .. Intrinsic Functions .. 00088 INTRINSIC DBLE, DCONJG, MAX 00089 * .. 00090 * .. Executable Statements .. 00091 * 00092 INFO = 0 00093 UPPER = LSAME( UPLO, 'U' ) 00094 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00095 INFO = -1 00096 ELSE IF( N.LT.0 ) THEN 00097 INFO = -2 00098 ELSE IF( NRHS.LT.0 ) THEN 00099 INFO = -3 00100 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00101 INFO = -5 00102 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00103 INFO = -8 00104 END IF 00105 IF( INFO.NE.0 ) THEN 00106 CALL XERBLA( 'ZHETRS2', -INFO ) 00107 RETURN 00108 END IF 00109 * 00110 * Quick return if possible 00111 * 00112 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00113 $ RETURN 00114 * 00115 * Convert A 00116 * 00117 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 00118 * 00119 IF( UPPER ) THEN 00120 * 00121 * Solve A*X = B, where A = U*D*U'. 00122 * 00123 * P' * B 00124 K=N 00125 DO WHILE ( K .GE. 1 ) 00126 IF( IPIV( K ).GT.0 ) THEN 00127 * 1 x 1 diagonal block 00128 * Interchange rows K and IPIV(K). 00129 KP = IPIV( K ) 00130 IF( KP.NE.K ) 00131 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00132 K=K-1 00133 ELSE 00134 * 2 x 2 diagonal block 00135 * Interchange rows K-1 and -IPIV(K). 00136 KP = -IPIV( K ) 00137 IF( KP.EQ.-IPIV( K-1 ) ) 00138 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 00139 K=K-2 00140 END IF 00141 END DO 00142 * 00143 * Compute (U \P' * B) -> B [ (U \P' * B) ] 00144 * 00145 CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N) 00146 * 00147 * Compute D \ B -> B [ D \ (U \P' * B) ] 00148 * 00149 I=N 00150 DO WHILE ( I .GE. 1 ) 00151 IF( IPIV(I) .GT. 0 ) THEN 00152 S = DBLE( ONE ) / DBLE( A( I, I ) ) 00153 CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 00154 ELSEIF ( I .GT. 1) THEN 00155 IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN 00156 AKM1K = WORK(I) 00157 AKM1 = A( I-1, I-1 ) / AKM1K 00158 AK = A( I, I ) / DCONJG( AKM1K ) 00159 DENOM = AKM1*AK - ONE 00160 DO 15 J = 1, NRHS 00161 BKM1 = B( I-1, J ) / AKM1K 00162 BK = B( I, J ) / DCONJG( AKM1K ) 00163 B( I-1, J ) = ( AK*BKM1-BK ) / DENOM 00164 B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM 00165 15 CONTINUE 00166 I = I - 1 00167 ENDIF 00168 ENDIF 00169 I = I - 1 00170 END DO 00171 * 00172 * Compute (U' \ B) -> B [ U' \ (D \ (U \P' * B) ) ] 00173 * 00174 CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,N,B,N) 00175 * 00176 * P * B [ P * (U' \ (D \ (U \P' * B) )) ] 00177 * 00178 K=1 00179 DO WHILE ( K .LE. N ) 00180 IF( IPIV( K ).GT.0 ) THEN 00181 * 1 x 1 diagonal block 00182 * Interchange rows K and IPIV(K). 00183 KP = IPIV( K ) 00184 IF( KP.NE.K ) 00185 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00186 K=K+1 00187 ELSE 00188 * 2 x 2 diagonal block 00189 * Interchange rows K-1 and -IPIV(K). 00190 KP = -IPIV( K ) 00191 IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) ) 00192 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00193 K=K+2 00194 ENDIF 00195 END DO 00196 * 00197 ELSE 00198 * 00199 * Solve A*X = B, where A = L*D*L'. 00200 * 00201 * P' * B 00202 K=1 00203 DO WHILE ( K .LE. N ) 00204 IF( IPIV( K ).GT.0 ) THEN 00205 * 1 x 1 diagonal block 00206 * Interchange rows K and IPIV(K). 00207 KP = IPIV( K ) 00208 IF( KP.NE.K ) 00209 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00210 K=K+1 00211 ELSE 00212 * 2 x 2 diagonal block 00213 * Interchange rows K and -IPIV(K+1). 00214 KP = -IPIV( K+1 ) 00215 IF( KP.EQ.-IPIV( K ) ) 00216 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 00217 K=K+2 00218 ENDIF 00219 END DO 00220 * 00221 * Compute (L \P' * B) -> B [ (L \P' * B) ] 00222 * 00223 CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N) 00224 * 00225 * Compute D \ B -> B [ D \ (L \P' * B) ] 00226 * 00227 I=1 00228 DO WHILE ( I .LE. N ) 00229 IF( IPIV(I) .GT. 0 ) THEN 00230 S = DBLE( ONE ) / DBLE( A( I, I ) ) 00231 CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 00232 ELSE 00233 AKM1K = WORK(I) 00234 AKM1 = A( I, I ) / DCONJG( AKM1K ) 00235 AK = A( I+1, I+1 ) / AKM1K 00236 DENOM = AKM1*AK - ONE 00237 DO 25 J = 1, NRHS 00238 BKM1 = B( I, J ) / DCONJG( AKM1K ) 00239 BK = B( I+1, J ) / AKM1K 00240 B( I, J ) = ( AK*BKM1-BK ) / DENOM 00241 B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 00242 25 CONTINUE 00243 I = I + 1 00244 ENDIF 00245 I = I + 1 00246 END DO 00247 * 00248 * Compute (L' \ B) -> B [ L' \ (D \ (L \P' * B) ) ] 00249 * 00250 CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,N,B,N) 00251 * 00252 * P * B [ P * (L' \ (D \ (L \P' * B) )) ] 00253 * 00254 K=N 00255 DO WHILE ( K .GE. 1 ) 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.GT.1 .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 END IF 00274 * 00275 * Revert A 00276 * 00277 CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) 00278 * 00279 RETURN 00280 * 00281 * End of ZHETRS2 00282 * 00283 END
1.7.3 
