LAPACK 3.3.0

ctfttp.f

Go to the documentation of this file.
00001       SUBROUTINE CTFTTP( 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       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
 All Files Functions