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