 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CTPTTF( TRANSR, UPLO, N, AP, ARF, 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 * Purpose 00020 * ======= 00021 * 00022 * CTPTTF copies a triangular matrix A from standard packed format (TP) 00023 * to rectangular full packed format (TF). 00024 * 00025 * Arguments 00026 * ========= 00027 * 00028 * TRANSR (input) CHARACTER*1 00029 * = 'N': ARF in Normal format is wanted; 00030 * = 'C': ARF in Conjugate-transpose format is wanted. 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 * AP (input) COMPLEX array, dimension ( N*(N+1)/2 ), 00040 * On entry, the upper or lower triangular matrix A, packed 00041 * columnwise in a linear array. The j-th column of A is stored 00042 * in the array AP as follows: 00043 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00044 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 00045 * 00046 * ARF (output) COMPLEX array, dimension ( N*(N+1)/2 ), 00047 * On exit, the upper or lower triangular matrix A stored in 00048 * RFP format. For a further discussion see Notes below. 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 CONJG, MOD 00178 * .. 00179 * .. Executable Statements .. 00180 * 00181 * Test the input parameters. 00182 * 00183 INFO = 0 00184 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00185 LOWER = LSAME( UPLO, 'L' ) 00186 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 00187 INFO = -1 00188 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00189 INFO = -2 00190 ELSE IF( N.LT.0 ) THEN 00191 INFO = -3 00192 END IF 00193 IF( INFO.NE.0 ) THEN 00194 CALL XERBLA( 'CTPTTF', -INFO ) 00195 RETURN 00196 END IF 00197 * 00198 * Quick return if possible 00199 * 00200 IF( N.EQ.0 ) 00201 $ RETURN 00202 * 00203 IF( N.EQ.1 ) THEN 00204 IF( NORMALTRANSR ) THEN 00205 ARF( 0 ) = AP( 0 ) 00206 ELSE 00207 ARF( 0 ) = CONJG( AP( 0 ) ) 00208 END IF 00209 RETURN 00210 END IF 00211 * 00212 * Size of array ARF(0:NT-1) 00213 * 00214 NT = N*( N+1 ) / 2 00215 * 00216 * Set N1 and N2 depending on LOWER 00217 * 00218 IF( LOWER ) THEN 00219 N2 = N / 2 00220 N1 = N - N2 00221 ELSE 00222 N1 = N / 2 00223 N2 = N - N1 00224 END IF 00225 * 00226 * If N is odd, set NISODD = .TRUE. 00227 * If N is even, set K = N/2 and NISODD = .FALSE. 00228 * 00229 * set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe) 00230 * where noe = 0 if n is even, noe = 1 if n is odd 00231 * 00232 IF( MOD( N, 2 ).EQ.0 ) THEN 00233 K = N / 2 00234 NISODD = .FALSE. 00235 LDA = N + 1 00236 ELSE 00237 NISODD = .TRUE. 00238 LDA = N 00239 END IF 00240 * 00241 * ARF^C has lda rows and n+1-noe cols 00242 * 00243 IF( .NOT.NORMALTRANSR ) 00244 $ LDA = ( N+1 ) / 2 00245 * 00246 * start execution: there are eight cases 00247 * 00248 IF( NISODD ) THEN 00249 * 00250 * N is odd 00251 * 00252 IF( NORMALTRANSR ) THEN 00253 * 00254 * N is odd and TRANSR = 'N' 00255 * 00256 IF( LOWER ) THEN 00257 * 00258 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) ) 00259 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0) 00260 * T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n 00261 * 00262 IJP = 0 00263 JP = 0 00264 DO J = 0, N2 00265 DO I = J, N - 1 00266 IJ = I + JP 00267 ARF( IJ ) = AP( IJP ) 00268 IJP = IJP + 1 00269 END DO 00270 JP = JP + LDA 00271 END DO 00272 DO I = 0, N2 - 1 00273 DO J = 1 + I, N2 00274 IJ = I + J*LDA 00275 ARF( IJ ) = CONJG( AP( IJP ) ) 00276 IJP = IJP + 1 00277 END DO 00278 END DO 00279 * 00280 ELSE 00281 * 00282 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1) 00283 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0) 00284 * T1 -> a(n2), T2 -> a(n1), S -> a(0) 00285 * 00286 IJP = 0 00287 DO J = 0, N1 - 1 00288 IJ = N2 + J 00289 DO I = 0, J 00290 ARF( IJ ) = CONJG( AP( IJP ) ) 00291 IJP = IJP + 1 00292 IJ = IJ + LDA 00293 END DO 00294 END DO 00295 JS = 0 00296 DO J = N1, N - 1 00297 IJ = JS 00298 DO IJ = JS, JS + J 00299 ARF( IJ ) = AP( IJP ) 00300 IJP = IJP + 1 00301 END DO 00302 JS = JS + LDA 00303 END DO 00304 * 00305 END IF 00306 * 00307 ELSE 00308 * 00309 * N is odd and TRANSR = 'C' 00310 * 00311 IF( LOWER ) THEN 00312 * 00313 * SRPA for LOWER, TRANSPOSE and N is odd 00314 * T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1) 00315 * T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1 00316 * 00317 IJP = 0 00318 DO I = 0, N2 00319 DO IJ = I*( LDA+1 ), N*LDA - 1, LDA 00320 ARF( IJ ) = CONJG( AP( IJP ) ) 00321 IJP = IJP + 1 00322 END DO 00323 END DO 00324 JS = 1 00325 DO J = 0, N2 - 1 00326 DO IJ = JS, JS + N2 - J - 1 00327 ARF( IJ ) = AP( IJP ) 00328 IJP = IJP + 1 00329 END DO 00330 JS = JS + LDA + 1 00331 END DO 00332 * 00333 ELSE 00334 * 00335 * SRPA for UPPER, TRANSPOSE and N is odd 00336 * T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0) 00337 * T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2 00338 * 00339 IJP = 0 00340 JS = N2*LDA 00341 DO J = 0, N1 - 1 00342 DO IJ = JS, JS + J 00343 ARF( IJ ) = AP( IJP ) 00344 IJP = IJP + 1 00345 END DO 00346 JS = JS + LDA 00347 END DO 00348 DO I = 0, N1 00349 DO IJ = I, I + ( N1+I )*LDA, LDA 00350 ARF( IJ ) = CONJG( AP( IJP ) ) 00351 IJP = IJP + 1 00352 END DO 00353 END DO 00354 * 00355 END IF 00356 * 00357 END IF 00358 * 00359 ELSE 00360 * 00361 * N is even 00362 * 00363 IF( NORMALTRANSR ) THEN 00364 * 00365 * N is even and TRANSR = 'N' 00366 * 00367 IF( LOWER ) THEN 00368 * 00369 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00370 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0) 00371 * T1 -> a(1), T2 -> a(0), S -> a(k+1) 00372 * 00373 IJP = 0 00374 JP = 0 00375 DO J = 0, K - 1 00376 DO I = J, N - 1 00377 IJ = 1 + I + JP 00378 ARF( IJ ) = AP( IJP ) 00379 IJP = IJP + 1 00380 END DO 00381 JP = JP + LDA 00382 END DO 00383 DO I = 0, K - 1 00384 DO J = I, K - 1 00385 IJ = I + J*LDA 00386 ARF( IJ ) = CONJG( AP( IJP ) ) 00387 IJP = IJP + 1 00388 END DO 00389 END DO 00390 * 00391 ELSE 00392 * 00393 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) ) 00394 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0) 00395 * T1 -> a(k+1), T2 -> a(k), S -> a(0) 00396 * 00397 IJP = 0 00398 DO J = 0, K - 1 00399 IJ = K + 1 + J 00400 DO I = 0, J 00401 ARF( IJ ) = CONJG( AP( IJP ) ) 00402 IJP = IJP + 1 00403 IJ = IJ + LDA 00404 END DO 00405 END DO 00406 JS = 0 00407 DO J = K, N - 1 00408 IJ = JS 00409 DO IJ = JS, JS + J 00410 ARF( IJ ) = AP( IJP ) 00411 IJP = IJP + 1 00412 END DO 00413 JS = JS + LDA 00414 END DO 00415 * 00416 END IF 00417 * 00418 ELSE 00419 * 00420 * N is even and TRANSR = 'C' 00421 * 00422 IF( LOWER ) THEN 00423 * 00424 * SRPA for LOWER, TRANSPOSE and N is even (see paper) 00425 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1) 00426 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k 00427 * 00428 IJP = 0 00429 DO I = 0, K - 1 00430 DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA 00431 ARF( IJ ) = CONJG( AP( IJP ) ) 00432 IJP = IJP + 1 00433 END DO 00434 END DO 00435 JS = 0 00436 DO J = 0, K - 1 00437 DO IJ = JS, JS + K - J - 1 00438 ARF( IJ ) = AP( IJP ) 00439 IJP = IJP + 1 00440 END DO 00441 JS = JS + LDA + 1 00442 END DO 00443 * 00444 ELSE 00445 * 00446 * SRPA for UPPER, TRANSPOSE and N is even (see paper) 00447 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0) 00448 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k 00449 * 00450 IJP = 0 00451 JS = ( K+1 )*LDA 00452 DO J = 0, K - 1 00453 DO IJ = JS, JS + J 00454 ARF( IJ ) = AP( IJP ) 00455 IJP = IJP + 1 00456 END DO 00457 JS = JS + LDA 00458 END DO 00459 DO I = 0, K - 1 00460 DO IJ = I, I + ( K+I )*LDA, LDA 00461 ARF( IJ ) = CONJG( AP( IJP ) ) 00462 IJP = IJP + 1 00463 END DO 00464 END DO 00465 * 00466 END IF 00467 * 00468 END IF 00469 * 00470 END IF 00471 * 00472 RETURN 00473 * 00474 * End of CTPTTF 00475 * 00476 END
1.7.3 
