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