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