LAPACK 3.3.0

ctrttf.f

Go to the documentation of this file.
00001       SUBROUTINE CTRTTF( TRANSR, UPLO, N, A, LDA, ARF, 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 *  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
 All Files Functions