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