SRC/ztfttp.f

   1 *> \brief \b ZTFTTP copies a triangular matrix from the rectangular full packed format (TF) to the standard packed format (TP).
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZTFTTP + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztfttp.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztfttp.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztfttp.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          TRANSR, UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       COMPLEX*16         AP( 0: * ), ARF( 0: * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZTFTTP copies a triangular matrix A from rectangular full packed
  38 *> format (TF) to standard packed format (TP).
  39 *> \endverbatim
  40 *
  41 *  Arguments:
  42 *  ==========
  43 *
  44 *> \param[in] TRANSR
  45 *> \verbatim
  46 *>          TRANSR is CHARACTER*1
  47 *>          = 'N':  ARF is in Normal format;
  48 *>          = 'C':  ARF is in Conjugate-transpose format;
  49 *> \endverbatim
  50 *>
  51 *> \param[in] UPLO
  52 *> \verbatim
  53 *>          UPLO is CHARACTER*1
  54 *>          = 'U':  A is upper triangular;
  55 *>          = 'L':  A is lower triangular.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] N
  59 *> \verbatim
  60 *>          N is INTEGER
  61 *>          The order of the matrix A. N >= 0.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] ARF
  65 *> \verbatim
  66 *>          ARF is COMPLEX*16 array, dimension ( N*(N+1)/2 ),
  67 *>          On entry, the upper or lower triangular matrix A stored in
  68 *>          RFP format. For a further discussion see Notes below.
  69 *> \endverbatim
  70 *>
  71 *> \param[out] AP
  72 *> \verbatim
  73 *>          AP is COMPLEX*16 array, dimension ( N*(N+1)/2 ),
  74 *>          On exit, the upper or lower triangular matrix A, packed
  75 *>          columnwise in a linear array. The j-th column of A is stored
  76 *>          in the array AP as follows:
  77 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  78 *>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  79 *> \endverbatim
  80 *>
  81 *> \param[out] INFO
  82 *> \verbatim
  83 *>          INFO is INTEGER
  84 *>          = 0:  successful exit
  85 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  86 *> \endverbatim
  87 *
  88 *  Authors:
  89 *  ========
  90 *
  91 *> \author Univ. of Tennessee
  92 *> \author Univ. of California Berkeley
  93 *> \author Univ. of Colorado Denver
  94 *> \author NAG Ltd.
  95 *
  96 *> \date September 2012
  97 *
  98 *> \ingroup complex16OTHERcomputational
  99 *
 100 *> \par Further Details:
 101 *  =====================
 102 *>
 103 *> \verbatim
 104 *>
 105 *>  We first consider Standard Packed Format when N is even.
 106 *>  We give an example where N = 6.
 107 *>
 108 *>      AP is Upper             AP is Lower
 109 *>
 110 *>   00 01 02 03 04 05       00
 111 *>      11 12 13 14 15       10 11
 112 *>         22 23 24 25       20 21 22
 113 *>            33 34 35       30 31 32 33
 114 *>               44 45       40 41 42 43 44
 115 *>                  55       50 51 52 53 54 55
 116 *>
 117 *>
 118 *>  Let TRANSR = 'N'. RFP holds AP as follows:
 119 *>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
 120 *>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
 121 *>  conjugate-transpose of the first three columns of AP upper.
 122 *>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
 123 *>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
 124 *>  conjugate-transpose of the last three columns of AP lower.
 125 *>  To denote conjugate we place -- above the element. This covers the
 126 *>  case N even and TRANSR = 'N'.
 127 *>
 128 *>         RFP A                   RFP A
 129 *>
 130 *>                                -- -- --
 131 *>        03 04 05                33 43 53
 132 *>                                   -- --
 133 *>        13 14 15                00 44 54
 134 *>                                      --
 135 *>        23 24 25                10 11 55
 136 *>
 137 *>        33 34 35                20 21 22
 138 *>        --
 139 *>        00 44 45                30 31 32
 140 *>        -- --
 141 *>        01 11 55                40 41 42
 142 *>        -- -- --
 143 *>        02 12 22                50 51 52
 144 *>
 145 *>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
 146 *>  transpose of RFP A above. One therefore gets:
 147 *>
 148 *>
 149 *>           RFP A                   RFP A
 150 *>
 151 *>     -- -- -- --                -- -- -- -- -- --
 152 *>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
 153 *>     -- -- -- -- --                -- -- -- -- --
 154 *>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
 155 *>     -- -- -- -- -- --                -- -- -- --
 156 *>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
 157 *>
 158 *>
 159 *>  We next consider Standard Packed Format when N is odd.
 160 *>  We give an example where N = 5.
 161 *>
 162 *>     AP is Upper                 AP is Lower
 163 *>
 164 *>   00 01 02 03 04              00
 165 *>      11 12 13 14              10 11
 166 *>         22 23 24              20 21 22
 167 *>            33 34              30 31 32 33
 168 *>               44              40 41 42 43 44
 169 *>
 170 *>
 171 *>  Let TRANSR = 'N'. RFP holds AP as follows:
 172 *>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
 173 *>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
 174 *>  conjugate-transpose of the first two   columns of AP upper.
 175 *>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
 176 *>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
 177 *>  conjugate-transpose of the last two   columns of AP lower.
 178 *>  To denote conjugate we place -- above the element. This covers the
 179 *>  case N odd  and TRANSR = 'N'.
 180 *>
 181 *>         RFP A                   RFP A
 182 *>
 183 *>                                   -- --
 184 *>        02 03 04                00 33 43
 185 *>                                      --
 186 *>        12 13 14                10 11 44
 187 *>
 188 *>        22 23 24                20 21 22
 189 *>        --
 190 *>        00 33 34                30 31 32
 191 *>        -- --
 192 *>        01 11 44                40 41 42
 193 *>
 194 *>  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-
 195 *>  transpose of RFP A above. One therefore gets:
 196 *>
 197 *>
 198 *>           RFP A                   RFP A
 199 *>
 200 *>     -- -- --                   -- -- -- -- -- --
 201 *>     02 12 22 00 01             00 10 20 30 40 50
 202 *>     -- -- -- --                   -- -- -- -- --
 203 *>     03 13 23 33 11             33 11 21 31 41 51
 204 *>     -- -- -- -- --                   -- -- -- --
 205 *>     04 14 24 34 44             43 44 22 32 42 52
 206 *> \endverbatim
 207 *>
 208 *  =====================================================================
 209       SUBROUTINE ZTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
 210 *
 211 *  -- LAPACK computational routine (version 3.4.2) --
 212 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 213 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 214 *     September 2012
 215 *
 216 *     .. Scalar Arguments ..
 217       CHARACTER          TRANSR, UPLO
 218       INTEGER            INFO, N
 219 *     ..
 220 *     .. Array Arguments ..
 221       COMPLEX*16         AP( 0: * ), ARF( 0: * )
 222 *     ..
 223 *
 224 *  =====================================================================
 225 *
 226 *     .. Parameters ..
 227 *     ..
 228 *     .. Local Scalars ..
 229       LOGICAL            LOWER, NISODD, NORMALTRANSR
 230       INTEGER            N1, N2, K, NT
 231       INTEGER            I, J, IJ
 232       INTEGER            IJP, JP, LDA, JS
 233 *     ..
 234 *     .. External Functions ..
 235       LOGICAL            LSAME
 236       EXTERNAL           LSAME
 237 *     ..
 238 *     .. External Subroutines ..
 239       EXTERNAL           XERBLA
 240 *     ..
 241 *     .. Intrinsic Functions ..
 242       INTRINSIC          DCONJG
 243 *     ..
 244 *     .. Intrinsic Functions ..
 245 *     ..
 246 *     .. Executable Statements ..
 247 *
 248 *     Test the input parameters.
 249 *
 250       INFO = 0
 251       NORMALTRANSR = LSAME( TRANSR, 'N' )
 252       LOWER = LSAME( UPLO, 'L' )
 253       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN
 254          INFO = -1
 255       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
 256          INFO = -2
 257       ELSE IF( N.LT.0 ) THEN
 258          INFO = -3
 259       END IF
 260       IF( INFO.NE.0 ) THEN
 261          CALL XERBLA( 'ZTFTTP', -INFO )
 262          RETURN
 263       END IF
 264 *
 265 *     Quick return if possible
 266 *
 267       IF( N.EQ.0 )
 268      $   RETURN
 269 *
 270       IF( N.EQ.1 ) THEN
 271          IF( NORMALTRANSR ) THEN
 272             AP( 0 ) = ARF( 0 )
 273          ELSE
 274             AP( 0 ) = DCONJG( ARF( 0 ) )
 275          END IF
 276          RETURN
 277       END IF
 278 *
 279 *     Size of array ARF(0:NT-1)
 280 *
 281       NT = N*( N+1 ) / 2
 282 *
 283 *     Set N1 and N2 depending on LOWER
 284 *
 285       IF( LOWER ) THEN
 286          N2 = N / 2
 287          N1 = N - N2
 288       ELSE
 289          N1 = N / 2
 290          N2 = N - N1
 291       END IF
 292 *
 293 *     If N is odd, set NISODD = .TRUE.
 294 *     If N is even, set K = N/2 and NISODD = .FALSE.
 295 *
 296 *     set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
 297 *     where noe = 0 if n is even, noe = 1 if n is odd
 298 *
 299       IF( MOD( N, 2 ).EQ.0 ) THEN
 300          K = N / 2
 301          NISODD = .FALSE.
 302          LDA = N + 1
 303       ELSE
 304          NISODD = .TRUE.
 305          LDA = N
 306       END IF
 307 *
 308 *     ARF^C has lda rows and n+1-noe cols
 309 *
 310       IF( .NOT.NORMALTRANSR )
 311      $   LDA = ( N+1 ) / 2
 312 *
 313 *     start execution: there are eight cases
 314 *
 315       IF( NISODD ) THEN
 316 *
 317 *        N is odd
 318 *
 319          IF( NORMALTRANSR ) THEN
 320 *
 321 *           N is odd and TRANSR = 'N'
 322 *
 323             IF( LOWER ) THEN
 324 *
 325 *             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
 326 *             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
 327 *             T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n
 328 *
 329                IJP = 0
 330                JP = 0
 331                DO J = 0, N2
 332                   DO I = J, N - 1
 333                      IJ = I + JP
 334                      AP( IJP ) = ARF( IJ )
 335                      IJP = IJP + 1
 336                   END DO
 337                   JP = JP + LDA
 338                END DO
 339                DO I = 0, N2 - 1
 340                   DO J = 1 + I, N2
 341                      IJ = I + J*LDA
 342                      AP( IJP ) = DCONJG( ARF( IJ ) )
 343                      IJP = IJP + 1
 344                   END DO
 345                END DO
 346 *
 347             ELSE
 348 *
 349 *             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
 350 *             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
 351 *             T1 -> a(n2), T2 -> a(n1), S -> a(0)
 352 *
 353                IJP = 0
 354                DO J = 0, N1 - 1
 355                   IJ = N2 + J
 356                   DO I = 0, J
 357                      AP( IJP ) = DCONJG( ARF( IJ ) )
 358                      IJP = IJP + 1
 359                      IJ = IJ + LDA
 360                   END DO
 361                END DO
 362                JS = 0
 363                DO J = N1, N - 1
 364                   IJ = JS
 365                   DO IJ = JS, JS + J
 366                      AP( IJP ) = ARF( IJ )
 367                      IJP = IJP + 1
 368                   END DO
 369                   JS = JS + LDA
 370                END DO
 371 *
 372             END IF
 373 *
 374          ELSE
 375 *
 376 *           N is odd and TRANSR = 'C'
 377 *
 378             IF( LOWER ) THEN
 379 *
 380 *              SRPA for LOWER, TRANSPOSE and N is odd
 381 *              T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
 382 *              T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
 383 *
 384                IJP = 0
 385                DO I = 0, N2
 386                   DO IJ = I*( LDA+1 ), N*LDA - 1, LDA
 387                      AP( IJP ) = DCONJG( ARF( IJ ) )
 388                      IJP = IJP + 1
 389                   END DO
 390                END DO
 391                JS = 1
 392                DO J = 0, N2 - 1
 393                   DO IJ = JS, JS + N2 - J - 1
 394                      AP( IJP ) = ARF( IJ )
 395                      IJP = IJP + 1
 396                   END DO
 397                   JS = JS + LDA + 1
 398                END DO
 399 *
 400             ELSE
 401 *
 402 *              SRPA for UPPER, TRANSPOSE and N is odd
 403 *              T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
 404 *              T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
 405 *
 406                IJP = 0
 407                JS = N2*LDA
 408                DO J = 0, N1 - 1
 409                   DO IJ = JS, JS + J
 410                      AP( IJP ) = ARF( IJ )
 411                      IJP = IJP + 1
 412                   END DO
 413                   JS = JS + LDA
 414                END DO
 415                DO I = 0, N1
 416                   DO IJ = I, I + ( N1+I )*LDA, LDA
 417                      AP( IJP ) = DCONJG( ARF( IJ ) )
 418                      IJP = IJP + 1
 419                   END DO
 420                END DO
 421 *
 422             END IF
 423 *
 424          END IF
 425 *
 426       ELSE
 427 *
 428 *        N is even
 429 *
 430          IF( NORMALTRANSR ) THEN
 431 *
 432 *           N is even and TRANSR = 'N'
 433 *
 434             IF( LOWER ) THEN
 435 *
 436 *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
 437 *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
 438 *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
 439 *
 440                IJP = 0
 441                JP = 0
 442                DO J = 0, K - 1
 443                   DO I = J, N - 1
 444                      IJ = 1 + I + JP
 445                      AP( IJP ) = ARF( IJ )
 446                      IJP = IJP + 1
 447                   END DO
 448                   JP = JP + LDA
 449                END DO
 450                DO I = 0, K - 1
 451                   DO J = I, K - 1
 452                      IJ = I + J*LDA
 453                      AP( IJP ) = DCONJG( ARF( IJ ) )
 454                      IJP = IJP + 1
 455                   END DO
 456                END DO
 457 *
 458             ELSE
 459 *
 460 *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
 461 *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
 462 *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
 463 *
 464                IJP = 0
 465                DO J = 0, K - 1
 466                   IJ = K + 1 + J
 467                   DO I = 0, J
 468                      AP( IJP ) = DCONJG( ARF( IJ ) )
 469                      IJP = IJP + 1
 470                      IJ = IJ + LDA
 471                   END DO
 472                END DO
 473                JS = 0
 474                DO J = K, N - 1
 475                   IJ = JS
 476                   DO IJ = JS, JS + J
 477                      AP( IJP ) = ARF( IJ )
 478                      IJP = IJP + 1
 479                   END DO
 480                   JS = JS + LDA
 481                END DO
 482 *
 483             END IF
 484 *
 485          ELSE
 486 *
 487 *           N is even and TRANSR = 'C'
 488 *
 489             IF( LOWER ) THEN
 490 *
 491 *              SRPA for LOWER, TRANSPOSE and N is even (see paper)
 492 *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
 493 *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
 494 *
 495                IJP = 0
 496                DO I = 0, K - 1
 497                   DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA
 498                      AP( IJP ) = DCONJG( ARF( IJ ) )
 499                      IJP = IJP + 1
 500                   END DO
 501                END DO
 502                JS = 0
 503                DO J = 0, K - 1
 504                   DO IJ = JS, JS + K - J - 1
 505                      AP( IJP ) = ARF( IJ )
 506                      IJP = IJP + 1
 507                   END DO
 508                   JS = JS + LDA + 1
 509                END DO
 510 *
 511             ELSE
 512 *
 513 *              SRPA for UPPER, TRANSPOSE and N is even (see paper)
 514 *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
 515 *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
 516 *
 517                IJP = 0
 518                JS = ( K+1 )*LDA
 519                DO J = 0, K - 1
 520                   DO IJ = JS, JS + J
 521                      AP( IJP ) = ARF( IJ )
 522                      IJP = IJP + 1
 523                   END DO
 524                   JS = JS + LDA
 525                END DO
 526                DO I = 0, K - 1
 527                   DO IJ = I, I + ( K+I )*LDA, LDA
 528                      AP( IJP ) = DCONJG( ARF( IJ ) )
 529                      IJP = IJP + 1
 530                   END DO
 531                END DO
 532 *
 533             END IF
 534 *
 535          END IF
 536 *
 537       END IF
 538 *
 539       RETURN
 540 *
 541 *     End of ZTFTTP
 542 *
 543       END