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