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