![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTFTTP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTFTTP + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztfttp.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztfttp.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztfttp.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTFTTP( TRANSR, UPLO, N, ARF, AP, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER TRANSR, UPLO 00025 * INTEGER INFO, N 00026 * .. 00027 * .. Array Arguments .. 00028 * COMPLEX*16 AP( 0: * ), ARF( 0: * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> ZTFTTP copies a triangular matrix A from rectangular full packed 00038 *> format (TF) to standard packed format (TP). 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] TRANSR 00045 *> \verbatim 00046 *> TRANSR is CHARACTER*1 00047 *> = 'N': ARF is in Normal format; 00048 *> = 'C': ARF is in Conjugate-transpose format; 00049 *> \endverbatim 00050 *> 00051 *> \param[in] UPLO 00052 *> \verbatim 00053 *> UPLO is CHARACTER*1 00054 *> = 'U': A is upper triangular; 00055 *> = 'L': A is lower triangular. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] N 00059 *> \verbatim 00060 *> N is INTEGER 00061 *> The order of the matrix A. N >= 0. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] ARF 00065 *> \verbatim 00066 *> ARF is COMPLEX*16 array, dimension ( N*(N+1)/2 ), 00067 *> On entry, the upper or lower triangular matrix A stored in 00068 *> RFP format. For a further discussion see Notes below. 00069 *> \endverbatim 00070 *> 00071 *> \param[out] AP 00072 *> \verbatim 00073 *> AP is COMPLEX*16 array, dimension ( N*(N+1)/2 ), 00074 *> On exit, the upper or lower triangular matrix A, packed 00075 *> columnwise in a linear array. The j-th column of A is stored 00076 *> in the array AP as follows: 00077 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00078 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 00079 *> \endverbatim 00080 *> 00081 *> \param[out] INFO 00082 *> \verbatim 00083 *> INFO is INTEGER 00084 *> = 0: successful exit 00085 *> < 0: if INFO = -i, the i-th argument had an illegal value 00086 *> \endverbatim 00087 * 00088 * Authors: 00089 * ======== 00090 * 00091 *> \author Univ. of Tennessee 00092 *> \author Univ. of California Berkeley 00093 *> \author Univ. of Colorado Denver 00094 *> \author NAG Ltd. 00095 * 00096 *> \date November 2011 00097 * 00098 *> \ingroup complex16OTHERcomputational 00099 * 00100 *> \par Further Details: 00101 * ===================== 00102 *> 00103 *> \verbatim 00104 *> 00105 *> We first consider Standard Packed Format when N is even. 00106 *> We give an example where N = 6. 00107 *> 00108 *> AP is Upper AP is Lower 00109 *> 00110 *> 00 01 02 03 04 05 00 00111 *> 11 12 13 14 15 10 11 00112 *> 22 23 24 25 20 21 22 00113 *> 33 34 35 30 31 32 33 00114 *> 44 45 40 41 42 43 44 00115 *> 55 50 51 52 53 54 55 00116 *> 00117 *> 00118 *> Let TRANSR = 'N'. RFP holds AP as follows: 00119 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00120 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00121 *> conjugate-transpose of the first three columns of AP upper. 00122 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00123 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00124 *> conjugate-transpose of the last three columns of AP lower. 00125 *> To denote conjugate we place -- above the element. This covers the 00126 *> case N even and TRANSR = 'N'. 00127 *> 00128 *> RFP A RFP A 00129 *> 00130 *> -- -- -- 00131 *> 03 04 05 33 43 53 00132 *> -- -- 00133 *> 13 14 15 00 44 54 00134 *> -- 00135 *> 23 24 25 10 11 55 00136 *> 00137 *> 33 34 35 20 21 22 00138 *> -- 00139 *> 00 44 45 30 31 32 00140 *> -- -- 00141 *> 01 11 55 40 41 42 00142 *> -- -- -- 00143 *> 02 12 22 50 51 52 00144 *> 00145 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00146 *> transpose of RFP A above. One therefore gets: 00147 *> 00148 *> 00149 *> RFP A RFP A 00150 *> 00151 *> -- -- -- -- -- -- -- -- -- -- 00152 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00153 *> -- -- -- -- -- -- -- -- -- -- 00154 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00155 *> -- -- -- -- -- -- -- -- -- -- 00156 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00157 *> 00158 *> 00159 *> We next consider Standard Packed Format when N is odd. 00160 *> We give an example where N = 5. 00161 *> 00162 *> AP is Upper AP is Lower 00163 *> 00164 *> 00 01 02 03 04 00 00165 *> 11 12 13 14 10 11 00166 *> 22 23 24 20 21 22 00167 *> 33 34 30 31 32 33 00168 *> 44 40 41 42 43 44 00169 *> 00170 *> 00171 *> Let TRANSR = 'N'. RFP holds AP as follows: 00172 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00173 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00174 *> conjugate-transpose of the first two columns of AP upper. 00175 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00176 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00177 *> conjugate-transpose of the last two columns of AP lower. 00178 *> To denote conjugate we place -- above the element. This covers the 00179 *> case N odd and TRANSR = 'N'. 00180 *> 00181 *> RFP A RFP A 00182 *> 00183 *> -- -- 00184 *> 02 03 04 00 33 43 00185 *> -- 00186 *> 12 13 14 10 11 44 00187 *> 00188 *> 22 23 24 20 21 22 00189 *> -- 00190 *> 00 33 34 30 31 32 00191 *> -- -- 00192 *> 01 11 44 40 41 42 00193 *> 00194 *> Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00195 *> transpose of RFP A above. One therefore gets: 00196 *> 00197 *> 00198 *> RFP A RFP A 00199 *> 00200 *> -- -- -- -- -- -- -- -- -- 00201 *> 02 12 22 00 01 00 10 20 30 40 50 00202 *> -- -- -- -- -- -- -- -- -- 00203 *> 03 13 23 33 11 33 11 21 31 41 51 00204 *> -- -- -- -- -- -- -- -- -- 00205 *> 04 14 24 34 44 43 44 22 32 42 52 00206 *> \endverbatim 00207 *> 00208 * ===================================================================== 00209 SUBROUTINE ZTFTTP( TRANSR, UPLO, N, ARF, AP, INFO ) 00210 * 00211 * -- LAPACK computational routine (version 3.4.0) -- 00212 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00213 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00214 * November 2011 00215 * 00216 * .. Scalar Arguments .. 00217 CHARACTER TRANSR, UPLO 00218 INTEGER INFO, N 00219 * .. 00220 * .. Array Arguments .. 00221 COMPLEX*16 AP( 0: * ), ARF( 0: * ) 00222 * .. 00223 * 00224 * ===================================================================== 00225 * 00226 * .. Parameters .. 00227 * .. 00228 * .. Local Scalars .. 00229 LOGICAL LOWER, NISODD, NORMALTRANSR 00230 INTEGER N1, N2, K, NT 00231 INTEGER I, J, IJ 00232 INTEGER IJP, JP, LDA, JS 00233 * .. 00234 * .. External Functions .. 00235 LOGICAL LSAME 00236 EXTERNAL LSAME 00237 * .. 00238 * .. External Subroutines .. 00239 EXTERNAL XERBLA 00240 * .. 00241 * .. Intrinsic Functions .. 00242 INTRINSIC DCONJG 00243 * .. 00244 * .. Intrinsic Functions .. 00245 * .. 00246 * .. Executable Statements .. 00247 * 00248 * Test the input parameters. 00249 * 00250 INFO = 0 00251 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00252 LOWER = LSAME( UPLO, 'L' ) 00253 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 00254 INFO = -1 00255 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00256 INFO = -2 00257 ELSE IF( N.LT.0 ) THEN 00258 INFO = -3 00259 END IF 00260 IF( INFO.NE.0 ) THEN 00261 CALL XERBLA( 'ZTFTTP', -INFO ) 00262 RETURN 00263 END IF 00264 * 00265 * Quick return if possible 00266 * 00267 IF( N.EQ.0 ) 00268 $ RETURN 00269 * 00270 IF( N.EQ.1 ) THEN 00271 IF( NORMALTRANSR ) THEN 00272 AP( 0 ) = ARF( 0 ) 00273 ELSE 00274 AP( 0 ) = DCONJG( ARF( 0 ) ) 00275 END IF 00276 RETURN 00277 END IF 00278 * 00279 * Size of array ARF(0:NT-1) 00280 * 00281 NT = N*( N+1 ) / 2 00282 * 00283 * Set N1 and N2 depending on LOWER 00284 * 00285 IF( LOWER ) THEN 00286 N2 = N / 2 00287 N1 = N - N2 00288 ELSE 00289 N1 = N / 2 00290 N2 = N - N1 00291 END IF 00292 * 00293 * If N is odd, set NISODD = .TRUE. 00294 * If N is even, set K = N/2 and NISODD = .FALSE. 00295 * 00296 * set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe) 00297 * where noe = 0 if n is even, noe = 1 if n is odd 00298 * 00299 IF( MOD( N, 2 ).EQ.0 ) THEN 00300 K = N / 2 00301 NISODD = .FALSE. 00302 LDA = N + 1 00303 ELSE 00304 NISODD = .TRUE. 00305 LDA = N 00306 END IF 00307 * 00308 * ARF^C has lda rows and n+1-noe cols 00309 * 00310 IF( .NOT.NORMALTRANSR ) 00311 $ LDA = ( N+1 ) / 2 00312 * 00313 * start execution: there are eight cases 00314 * 00315 IF( NISODD ) THEN 00316 * 00317 * N is odd 00318 * 00319 IF( NORMALTRANSR ) THEN 00320 * 00321 * N is odd and TRANSR = 'N' 00322 * 00323 IF( LOWER ) THEN 00324 * 00325 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) ) 00326 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0) 00327 * T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n 00328 * 00329 IJP = 0 00330 JP = 0 00331 DO J = 0, N2 00332 DO I = J, N - 1 00333 IJ = I + JP 00334 AP( IJP ) = ARF( IJ ) 00335 IJP = IJP + 1 00336 END DO 00337 JP = JP + LDA 00338 END DO 00339 DO I = 0, N2 - 1 00340 DO J = 1 + I, N2 00341 IJ = I + J*LDA 00342 AP( IJP ) = DCONJG( ARF( IJ ) ) 00343 IJP = IJP + 1 00344 END DO 00345 END DO 00346 * 00347 ELSE 00348 * 00349 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1) 00350 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0) 00351 * T1 -> a(n2), T2 -> a(n1), S -> a(0) 00352 * 00353 IJP = 0 00354 DO J = 0, N1 - 1 00355 IJ = N2 + J 00356 DO I = 0, J 00357 AP( IJP ) = DCONJG( ARF( IJ ) ) 00358 IJP = IJP + 1 00359 IJ = IJ + LDA 00360 END DO 00361 END DO 00362 JS = 0 00363 DO J = N1, N - 1 00364 IJ = JS 00365 DO IJ = JS, JS + J 00366 AP( IJP ) = ARF( IJ ) 00367 IJP = IJP + 1 00368 END DO 00369 JS = JS + LDA 00370 END DO 00371 * 00372 END IF 00373 * 00374 ELSE 00375 * 00376 * N is odd and TRANSR = 'C' 00377 * 00378 IF( LOWER ) THEN 00379 * 00380 * SRPA for LOWER, TRANSPOSE and N is odd 00381 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1) 00382 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1 00383 * 00384 IJP = 0 00385 DO I = 0, N2 00386 DO IJ = I*( LDA+1 ), N*LDA - 1, LDA 00387 AP( IJP ) = DCONJG( ARF( IJ ) ) 00388 IJP = IJP + 1 00389 END DO 00390 END DO 00391 JS = 1 00392 DO J = 0, N2 - 1 00393 DO IJ = JS, JS + N2 - J - 1 00394 AP( IJP ) = ARF( IJ ) 00395 IJP = IJP + 1 00396 END DO 00397 JS = JS + LDA + 1 00398 END DO 00399 * 00400 ELSE 00401 * 00402 * SRPA for UPPER, TRANSPOSE and N is odd 00403 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0) 00404 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2 00405 * 00406 IJP = 0 00407 JS = N2*LDA 00408 DO J = 0, N1 - 1 00409 DO IJ = JS, JS + J 00410 AP( IJP ) = ARF( IJ ) 00411 IJP = IJP + 1 00412 END DO 00413 JS = JS + LDA 00414 END DO 00415 DO I = 0, N1 00416 DO IJ = I, I + ( N1+I )*LDA, LDA 00417 AP( IJP ) = DCONJG( ARF( IJ ) ) 00418 IJP = IJP + 1 00419 END DO 00420 END DO 00421 * 00422 END IF 00423 * 00424 END IF 00425 * 00426 ELSE 00427 * 00428 * N is even 00429 * 00430 IF( NORMALTRANSR ) THEN 00431 * 00432 * N is even and TRANSR = 'N' 00433 * 00434 IF( LOWER ) THEN 00435 * 00436 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00437 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 00438 * T1 -> a(1), T2 -> a(0), S -> a(k+1) 00439 * 00440 IJP = 0 00441 JP = 0 00442 DO J = 0, K - 1 00443 DO I = J, N - 1 00444 IJ = 1 + I + JP 00445 AP( IJP ) = ARF( IJ ) 00446 IJP = IJP + 1 00447 END DO 00448 JP = JP + LDA 00449 END DO 00450 DO I = 0, K - 1 00451 DO J = I, K - 1 00452 IJ = I + J*LDA 00453 AP( IJP ) = DCONJG( ARF( IJ ) ) 00454 IJP = IJP + 1 00455 END DO 00456 END DO 00457 * 00458 ELSE 00459 * 00460 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00461 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 00462 * T1 -> a(k+1), T2 -> a(k), S -> a(0) 00463 * 00464 IJP = 0 00465 DO J = 0, K - 1 00466 IJ = K + 1 + J 00467 DO I = 0, J 00468 AP( IJP ) = DCONJG( ARF( IJ ) ) 00469 IJP = IJP + 1 00470 IJ = IJ + LDA 00471 END DO 00472 END DO 00473 JS = 0 00474 DO J = K, N - 1 00475 IJ = JS 00476 DO IJ = JS, JS + J 00477 AP( IJP ) = ARF( IJ ) 00478 IJP = IJP + 1 00479 END DO 00480 JS = JS + LDA 00481 END DO 00482 * 00483 END IF 00484 * 00485 ELSE 00486 * 00487 * N is even and TRANSR = 'C' 00488 * 00489 IF( LOWER ) THEN 00490 * 00491 * SRPA for LOWER, TRANSPOSE and N is even (see paper) 00492 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1) 00493 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 00494 * 00495 IJP = 0 00496 DO I = 0, K - 1 00497 DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA 00498 AP( IJP ) = DCONJG( ARF( IJ ) ) 00499 IJP = IJP + 1 00500 END DO 00501 END DO 00502 JS = 0 00503 DO J = 0, K - 1 00504 DO IJ = JS, JS + K - J - 1 00505 AP( IJP ) = ARF( IJ ) 00506 IJP = IJP + 1 00507 END DO 00508 JS = JS + LDA + 1 00509 END DO 00510 * 00511 ELSE 00512 * 00513 * SRPA for UPPER, TRANSPOSE and N is even (see paper) 00514 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0) 00515 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 00516 * 00517 IJP = 0 00518 JS = ( K+1 )*LDA 00519 DO J = 0, K - 1 00520 DO IJ = JS, JS + J 00521 AP( IJP ) = ARF( IJ ) 00522 IJP = IJP + 1 00523 END DO 00524 JS = JS + LDA 00525 END DO 00526 DO I = 0, K - 1 00527 DO IJ = I, I + ( K+I )*LDA, LDA 00528 AP( IJP ) = DCONJG( ARF( IJ ) ) 00529 IJP = IJP + 1 00530 END DO 00531 END DO 00532 * 00533 END IF 00534 * 00535 END IF 00536 * 00537 END IF 00538 * 00539 RETURN 00540 * 00541 * End of ZTFTTP 00542 * 00543 END