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