LAPACK 3.3.0

ctfttr.f

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