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