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