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