SRC/ctpttf.f

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