 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZPFTRI( TRANSR, UPLO, N, A, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * 00005 * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 00006 * -- April 2011 -- 00007 * 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER TRANSR, UPLO 00013 INTEGER INFO, N 00014 * .. Array Arguments .. 00015 COMPLEX*16 A( 0: * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZPFTRI computes the inverse of a complex Hermitian positive definite 00022 * matrix A using the Cholesky factorization A = U**H*U or A = L*L**H 00023 * computed by ZPFTRF. 00024 * 00025 * Arguments 00026 * ========= 00027 * 00028 * TRANSR (input) CHARACTER*1 00029 * = 'N': The Normal TRANSR of RFP A is stored; 00030 * = 'C': The Conjugate-transpose TRANSR of RFP A is stored. 00031 * 00032 * UPLO (input) CHARACTER*1 00033 * = 'U': Upper triangle of A is stored; 00034 * = 'L': Lower triangle of A is stored. 00035 * 00036 * N (input) INTEGER 00037 * The order of the matrix A. N >= 0. 00038 * 00039 * A (input/output) COMPLEX*16 array, dimension ( N*(N+1)/2 ); 00040 * On entry, the Hermitian matrix A in RFP format. RFP format is 00041 * described by TRANSR, UPLO, and N as follows: If TRANSR = 'N' 00042 * then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is 00043 * (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is 00044 * the Conjugate-transpose of RFP A as defined when 00045 * TRANSR = 'N'. The contents of RFP A are defined by UPLO as 00046 * follows: If UPLO = 'U' the RFP A contains the nt elements of 00047 * upper packed A. If UPLO = 'L' the RFP A contains the elements 00048 * of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR = 00049 * 'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N 00050 * is odd. See the Note below for more details. 00051 * 00052 * On exit, the Hermitian inverse of the original matrix, in the 00053 * same storage format. 00054 * 00055 * INFO (output) INTEGER 00056 * = 0: successful exit 00057 * < 0: if INFO = -i, the i-th argument had an illegal value 00058 * > 0: if INFO = i, the (i,i) element of the factor U or L is 00059 * zero, and the inverse could not be computed. 00060 * 00061 * Further Details 00062 * =============== 00063 * 00064 * We first consider Standard Packed Format when N is even. 00065 * We give an example where N = 6. 00066 * 00067 * AP is Upper AP is Lower 00068 * 00069 * 00 01 02 03 04 05 00 00070 * 11 12 13 14 15 10 11 00071 * 22 23 24 25 20 21 22 00072 * 33 34 35 30 31 32 33 00073 * 44 45 40 41 42 43 44 00074 * 55 50 51 52 53 54 55 00075 * 00076 * 00077 * Let TRANSR = 'N'. RFP holds AP as follows: 00078 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00079 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00080 * conjugate-transpose of the first three columns of AP upper. 00081 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00082 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00083 * conjugate-transpose of the last three columns of AP lower. 00084 * To denote conjugate we place -- above the element. This covers the 00085 * case N even and TRANSR = 'N'. 00086 * 00087 * RFP A RFP A 00088 * 00089 * -- -- -- 00090 * 03 04 05 33 43 53 00091 * -- -- 00092 * 13 14 15 00 44 54 00093 * -- 00094 * 23 24 25 10 11 55 00095 * 00096 * 33 34 35 20 21 22 00097 * -- 00098 * 00 44 45 30 31 32 00099 * -- -- 00100 * 01 11 55 40 41 42 00101 * -- -- -- 00102 * 02 12 22 50 51 52 00103 * 00104 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00105 * transpose of RFP A above. One therefore gets: 00106 * 00107 * 00108 * RFP A RFP A 00109 * 00110 * -- -- -- -- -- -- -- -- -- -- 00111 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00112 * -- -- -- -- -- -- -- -- -- -- 00113 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00114 * -- -- -- -- -- -- -- -- -- -- 00115 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00116 * 00117 * 00118 * We next consider Standard Packed Format when N is odd. 00119 * We give an example where N = 5. 00120 * 00121 * AP is Upper AP is Lower 00122 * 00123 * 00 01 02 03 04 00 00124 * 11 12 13 14 10 11 00125 * 22 23 24 20 21 22 00126 * 33 34 30 31 32 33 00127 * 44 40 41 42 43 44 00128 * 00129 * 00130 * Let TRANSR = 'N'. RFP holds AP as follows: 00131 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00132 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00133 * conjugate-transpose of the first two columns of AP upper. 00134 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00135 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00136 * conjugate-transpose of the last two columns of AP lower. 00137 * To denote conjugate we place -- above the element. This covers the 00138 * case N odd and TRANSR = 'N'. 00139 * 00140 * RFP A RFP A 00141 * 00142 * -- -- 00143 * 02 03 04 00 33 43 00144 * -- 00145 * 12 13 14 10 11 44 00146 * 00147 * 22 23 24 20 21 22 00148 * -- 00149 * 00 33 34 30 31 32 00150 * -- -- 00151 * 01 11 44 40 41 42 00152 * 00153 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00154 * transpose of RFP A above. One therefore gets: 00155 * 00156 * 00157 * RFP A RFP A 00158 * 00159 * -- -- -- -- -- -- -- -- -- 00160 * 02 12 22 00 01 00 10 20 30 40 50 00161 * -- -- -- -- -- -- -- -- -- 00162 * 03 13 23 33 11 33 11 21 31 41 51 00163 * -- -- -- -- -- -- -- -- -- 00164 * 04 14 24 34 44 43 44 22 32 42 52 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 DOUBLE PRECISION ONE 00170 COMPLEX*16 CONE 00171 PARAMETER ( ONE = 1.D0, CONE = ( 1.D0, 0.D0 ) ) 00172 * .. 00173 * .. Local Scalars .. 00174 LOGICAL LOWER, NISODD, NORMALTRANSR 00175 INTEGER N1, N2, K 00176 * .. 00177 * .. External Functions .. 00178 LOGICAL LSAME 00179 EXTERNAL LSAME 00180 * .. 00181 * .. External Subroutines .. 00182 EXTERNAL XERBLA, ZTFTRI, ZLAUUM, ZTRMM, ZHERK 00183 * .. 00184 * .. Intrinsic Functions .. 00185 INTRINSIC MOD 00186 * .. 00187 * .. Executable Statements .. 00188 * 00189 * Test the input parameters. 00190 * 00191 INFO = 0 00192 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00193 LOWER = LSAME( UPLO, 'L' ) 00194 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 00195 INFO = -1 00196 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00197 INFO = -2 00198 ELSE IF( N.LT.0 ) THEN 00199 INFO = -3 00200 END IF 00201 IF( INFO.NE.0 ) THEN 00202 CALL XERBLA( 'ZPFTRI', -INFO ) 00203 RETURN 00204 END IF 00205 * 00206 * Quick return if possible 00207 * 00208 IF( N.EQ.0 ) 00209 $ RETURN 00210 * 00211 * Invert the triangular Cholesky factor U or L. 00212 * 00213 CALL ZTFTRI( TRANSR, UPLO, 'N', N, A, INFO ) 00214 IF( INFO.GT.0 ) 00215 $ RETURN 00216 * 00217 * If N is odd, set NISODD = .TRUE. 00218 * If N is even, set K = N/2 and NISODD = .FALSE. 00219 * 00220 IF( MOD( N, 2 ).EQ.0 ) THEN 00221 K = N / 2 00222 NISODD = .FALSE. 00223 ELSE 00224 NISODD = .TRUE. 00225 END IF 00226 * 00227 * Set N1 and N2 depending on LOWER 00228 * 00229 IF( LOWER ) THEN 00230 N2 = N / 2 00231 N1 = N - N2 00232 ELSE 00233 N1 = N / 2 00234 N2 = N - N1 00235 END IF 00236 * 00237 * Start execution of triangular matrix multiply: inv(U)*inv(U)^C or 00238 * inv(L)^C*inv(L). There are eight cases. 00239 * 00240 IF( NISODD ) THEN 00241 * 00242 * N is odd 00243 * 00244 IF( NORMALTRANSR ) THEN 00245 * 00246 * N is odd and TRANSR = 'N' 00247 * 00248 IF( LOWER ) THEN 00249 * 00250 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) ) 00251 * T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0) 00252 * T1 -> a(0), T2 -> a(n), S -> a(N1) 00253 * 00254 CALL ZLAUUM( 'L', N1, A( 0 ), N, INFO ) 00255 CALL ZHERK( 'L', 'C', N1, N2, ONE, A( N1 ), N, ONE, 00256 $ A( 0 ), N ) 00257 CALL ZTRMM( 'L', 'U', 'N', 'N', N2, N1, CONE, A( N ), N, 00258 $ A( N1 ), N ) 00259 CALL ZLAUUM( 'U', N2, A( N ), N, INFO ) 00260 * 00261 ELSE 00262 * 00263 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1) 00264 * T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0) 00265 * T1 -> a(N2), T2 -> a(N1), S -> a(0) 00266 * 00267 CALL ZLAUUM( 'L', N1, A( N2 ), N, INFO ) 00268 CALL ZHERK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE, 00269 $ A( N2 ), N ) 00270 CALL ZTRMM( 'R', 'U', 'C', 'N', N1, N2, CONE, A( N1 ), N, 00271 $ A( 0 ), N ) 00272 CALL ZLAUUM( 'U', N2, A( N1 ), N, INFO ) 00273 * 00274 END IF 00275 * 00276 ELSE 00277 * 00278 * N is odd and TRANSR = 'C' 00279 * 00280 IF( LOWER ) THEN 00281 * 00282 * SRPA for LOWER, TRANSPOSE, and N is odd 00283 * T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1) 00284 * 00285 CALL ZLAUUM( 'U', N1, A( 0 ), N1, INFO ) 00286 CALL ZHERK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE, 00287 $ A( 0 ), N1 ) 00288 CALL ZTRMM( 'R', 'L', 'N', 'N', N1, N2, CONE, A( 1 ), N1, 00289 $ A( N1*N1 ), N1 ) 00290 CALL ZLAUUM( 'L', N2, A( 1 ), N1, INFO ) 00291 * 00292 ELSE 00293 * 00294 * SRPA for UPPER, TRANSPOSE, and N is odd 00295 * T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0) 00296 * 00297 CALL ZLAUUM( 'U', N1, A( N2*N2 ), N2, INFO ) 00298 CALL ZHERK( 'U', 'C', N1, N2, ONE, A( 0 ), N2, ONE, 00299 $ A( N2*N2 ), N2 ) 00300 CALL ZTRMM( 'L', 'L', 'C', 'N', N2, N1, CONE, A( N1*N2 ), 00301 $ N2, A( 0 ), N2 ) 00302 CALL ZLAUUM( 'L', N2, A( N1*N2 ), N2, INFO ) 00303 * 00304 END IF 00305 * 00306 END IF 00307 * 00308 ELSE 00309 * 00310 * N is even 00311 * 00312 IF( NORMALTRANSR ) THEN 00313 * 00314 * N is even and TRANSR = 'N' 00315 * 00316 IF( LOWER ) THEN 00317 * 00318 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00319 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 00320 * T1 -> a(1), T2 -> a(0), S -> a(k+1) 00321 * 00322 CALL ZLAUUM( 'L', K, A( 1 ), N+1, INFO ) 00323 CALL ZHERK( 'L', 'C', K, K, ONE, A( K+1 ), N+1, ONE, 00324 $ A( 1 ), N+1 ) 00325 CALL ZTRMM( 'L', 'U', 'N', 'N', K, K, CONE, A( 0 ), N+1, 00326 $ A( K+1 ), N+1 ) 00327 CALL ZLAUUM( 'U', K, A( 0 ), N+1, INFO ) 00328 * 00329 ELSE 00330 * 00331 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00332 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 00333 * T1 -> a(k+1), T2 -> a(k), S -> a(0) 00334 * 00335 CALL ZLAUUM( 'L', K, A( K+1 ), N+1, INFO ) 00336 CALL ZHERK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE, 00337 $ A( K+1 ), N+1 ) 00338 CALL ZTRMM( 'R', 'U', 'C', 'N', K, K, CONE, A( K ), N+1, 00339 $ A( 0 ), N+1 ) 00340 CALL ZLAUUM( 'U', K, A( K ), N+1, INFO ) 00341 * 00342 END IF 00343 * 00344 ELSE 00345 * 00346 * N is even and TRANSR = 'C' 00347 * 00348 IF( LOWER ) THEN 00349 * 00350 * SRPA for LOWER, TRANSPOSE, and N is even (see paper) 00351 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1), 00352 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 00353 * 00354 CALL ZLAUUM( 'U', K, A( K ), K, INFO ) 00355 CALL ZHERK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE, 00356 $ A( K ), K ) 00357 CALL ZTRMM( 'R', 'L', 'N', 'N', K, K, CONE, A( 0 ), K, 00358 $ A( K*( K+1 ) ), K ) 00359 CALL ZLAUUM( 'L', K, A( 0 ), K, INFO ) 00360 * 00361 ELSE 00362 * 00363 * SRPA for UPPER, TRANSPOSE, and N is even (see paper) 00364 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0), 00365 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 00366 * 00367 CALL ZLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO ) 00368 CALL ZHERK( 'U', 'C', K, K, ONE, A( 0 ), K, ONE, 00369 $ A( K*( K+1 ) ), K ) 00370 CALL ZTRMM( 'L', 'L', 'C', 'N', K, K, CONE, A( K*K ), K, 00371 $ A( 0 ), K ) 00372 CALL ZLAUUM( 'L', K, A( K*K ), K, INFO ) 00373 * 00374 END IF 00375 * 00376 END IF 00377 * 00378 END IF 00379 * 00380 RETURN 00381 * 00382 * End of ZPFTRI 00383 * 00384 END
1.7.3 
