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