SRC/dtfttp.f

   1 *> \brief \b DTFTTP 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 DTFTTP + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtfttp.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtfttp.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtfttp.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          TRANSR, UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   AP( 0: * ), ARF( 0: * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DTFTTP 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 *>          = 'T':  ARF is in 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERcomputational
  99 *
 100 *> \par Further Details:
 101 *  =====================
 102 *>
 103 *> \verbatim
 104 *>
 105 *>  We first consider Rectangular Full Packed (RFP) Format when N is
 106 *>  even. 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 *>  the 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 *>  the transpose of the last three columns of AP lower.
 125 *>  This covers the case N even and TRANSR = 'N'.
 126 *>
 127 *>         RFP A                   RFP A
 128 *>
 129 *>        03 04 05                33 43 53
 130 *>        13 14 15                00 44 54
 131 *>        23 24 25                10 11 55
 132 *>        33 34 35                20 21 22
 133 *>        00 44 45                30 31 32
 134 *>        01 11 55                40 41 42
 135 *>        02 12 22                50 51 52
 136 *>
 137 *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
 138 *>  transpose of RFP A above. One therefore gets:
 139 *>
 140 *>
 141 *>           RFP A                   RFP A
 142 *>
 143 *>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
 144 *>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
 145 *>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
 146 *>
 147 *>
 148 *>  We then consider Rectangular Full Packed (RFP) Format when N is
 149 *>  odd. We give an example where N = 5.
 150 *>
 151 *>     AP is Upper                 AP is Lower
 152 *>
 153 *>   00 01 02 03 04              00
 154 *>      11 12 13 14              10 11
 155 *>         22 23 24              20 21 22
 156 *>            33 34              30 31 32 33
 157 *>               44              40 41 42 43 44
 158 *>
 159 *>
 160 *>  Let TRANSR = 'N'. RFP holds AP as follows:
 161 *>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
 162 *>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
 163 *>  the transpose of the first two columns of AP upper.
 164 *>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
 165 *>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
 166 *>  the transpose of the last two columns of AP lower.
 167 *>  This covers the case N odd and TRANSR = 'N'.
 168 *>
 169 *>         RFP A                   RFP A
 170 *>
 171 *>        02 03 04                00 33 43
 172 *>        12 13 14                10 11 44
 173 *>        22 23 24                20 21 22
 174 *>        00 33 34                30 31 32
 175 *>        01 11 44                40 41 42
 176 *>
 177 *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
 178 *>  transpose of RFP A above. One therefore gets:
 179 *>
 180 *>           RFP A                   RFP A
 181 *>
 182 *>     02 12 22 00 01             00 10 20 30 40 50
 183 *>     03 13 23 33 11             33 11 21 31 41 51
 184 *>     04 14 24 34 44             43 44 22 32 42 52
 185 *> \endverbatim
 186 *>
 187 *  =====================================================================
 188       SUBROUTINE DTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
 189 *
 190 *  -- LAPACK computational routine (version 3.4.2) --
 191 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 192 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 193 *     September 2012
 194 *
 195 *     .. Scalar Arguments ..
 196       CHARACTER          TRANSR, UPLO
 197       INTEGER            INFO, N
 198 *     ..
 199 *     .. Array Arguments ..
 200       DOUBLE PRECISION   AP( 0: * ), ARF( 0: * )
 201 *     ..
 202 *
 203 *  =====================================================================
 204 *
 205 *     .. Parameters ..
 206 *     ..
 207 *     .. Local Scalars ..
 208       LOGICAL            LOWER, NISODD, NORMALTRANSR
 209       INTEGER            N1, N2, K, NT
 210       INTEGER            I, J, IJ
 211       INTEGER            IJP, JP, LDA, JS
 212 *     ..
 213 *     .. External Functions ..
 214       LOGICAL            LSAME
 215       EXTERNAL           LSAME
 216 *     ..
 217 *     .. External Subroutines ..
 218       EXTERNAL           XERBLA
 219 *     ..
 220 *     .. Executable Statements ..
 221 *
 222 *     Test the input parameters.
 223 *
 224       INFO = 0
 225       NORMALTRANSR = LSAME( TRANSR, 'N' )
 226       LOWER = LSAME( UPLO, 'L' )
 227       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
 228          INFO = -1
 229       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
 230          INFO = -2
 231       ELSE IF( N.LT.0 ) THEN
 232          INFO = -3
 233       END IF
 234       IF( INFO.NE.0 ) THEN
 235          CALL XERBLA( 'DTFTTP', -INFO )
 236          RETURN
 237       END IF
 238 *
 239 *     Quick return if possible
 240 *
 241       IF( N.EQ.0 )
 242      $   RETURN
 243 *
 244       IF( N.EQ.1 ) THEN
 245          IF( NORMALTRANSR ) THEN
 246             AP( 0 ) = ARF( 0 )
 247          ELSE
 248             AP( 0 ) = ARF( 0 )
 249          END IF
 250          RETURN
 251       END IF
 252 *
 253 *     Size of array ARF(0:NT-1)
 254 *
 255       NT = N*( N+1 ) / 2
 256 *
 257 *     Set N1 and N2 depending on LOWER
 258 *
 259       IF( LOWER ) THEN
 260          N2 = N / 2
 261          N1 = N - N2
 262       ELSE
 263          N1 = N / 2
 264          N2 = N - N1
 265       END IF
 266 *
 267 *     If N is odd, set NISODD = .TRUE.
 268 *     If N is even, set K = N/2 and NISODD = .FALSE.
 269 *
 270 *     set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
 271 *     where noe = 0 if n is even, noe = 1 if n is odd
 272 *
 273       IF( MOD( N, 2 ).EQ.0 ) THEN
 274          K = N / 2
 275          NISODD = .FALSE.
 276          LDA = N + 1
 277       ELSE
 278          NISODD = .TRUE.
 279          LDA = N
 280       END IF
 281 *
 282 *     ARF^C has lda rows and n+1-noe cols
 283 *
 284       IF( .NOT.NORMALTRANSR )
 285      $   LDA = ( N+1 ) / 2
 286 *
 287 *     start execution: there are eight cases
 288 *
 289       IF( NISODD ) THEN
 290 *
 291 *        N is odd
 292 *
 293          IF( NORMALTRANSR ) THEN
 294 *
 295 *           N is odd and TRANSR = 'N'
 296 *
 297             IF( LOWER ) THEN
 298 *
 299 *             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
 300 *             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
 301 *             T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n
 302 *
 303                IJP = 0
 304                JP = 0
 305                DO J = 0, N2
 306                   DO I = J, N - 1
 307                      IJ = I + JP
 308                      AP( IJP ) = ARF( IJ )
 309                      IJP = IJP + 1
 310                   END DO
 311                   JP = JP + LDA
 312                END DO
 313                DO I = 0, N2 - 1
 314                   DO J = 1 + I, N2
 315                      IJ = I + J*LDA
 316                      AP( IJP ) = ARF( IJ )
 317                      IJP = IJP + 1
 318                   END DO
 319                END DO
 320 *
 321             ELSE
 322 *
 323 *             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
 324 *             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
 325 *             T1 -> a(n2), T2 -> a(n1), S -> a(0)
 326 *
 327                IJP = 0
 328                DO J = 0, N1 - 1
 329                   IJ = N2 + J
 330                   DO I = 0, J
 331                      AP( IJP ) = ARF( IJ )
 332                      IJP = IJP + 1
 333                      IJ = IJ + LDA
 334                   END DO
 335                END DO
 336                JS = 0
 337                DO J = N1, N - 1
 338                   IJ = JS
 339                   DO IJ = JS, JS + J
 340                      AP( IJP ) = ARF( IJ )
 341                      IJP = IJP + 1
 342                   END DO
 343                   JS = JS + LDA
 344                END DO
 345 *
 346             END IF
 347 *
 348          ELSE
 349 *
 350 *           N is odd and TRANSR = 'T'
 351 *
 352             IF( LOWER ) THEN
 353 *
 354 *              SRPA for LOWER, TRANSPOSE and N is odd
 355 *              T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
 356 *              T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
 357 *
 358                IJP = 0
 359                DO I = 0, N2
 360                   DO IJ = I*( LDA+1 ), N*LDA - 1, LDA
 361                      AP( IJP ) = ARF( IJ )
 362                      IJP = IJP + 1
 363                   END DO
 364                END DO
 365                JS = 1
 366                DO J = 0, N2 - 1
 367                   DO IJ = JS, JS + N2 - J - 1
 368                      AP( IJP ) = ARF( IJ )
 369                      IJP = IJP + 1
 370                   END DO
 371                   JS = JS + LDA + 1
 372                END DO
 373 *
 374             ELSE
 375 *
 376 *              SRPA for UPPER, TRANSPOSE and N is odd
 377 *              T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
 378 *              T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
 379 *
 380                IJP = 0
 381                JS = N2*LDA
 382                DO J = 0, N1 - 1
 383                   DO IJ = JS, JS + J
 384                      AP( IJP ) = ARF( IJ )
 385                      IJP = IJP + 1
 386                   END DO
 387                   JS = JS + LDA
 388                END DO
 389                DO I = 0, N1
 390                   DO IJ = I, I + ( N1+I )*LDA, LDA
 391                      AP( IJP ) = ARF( IJ )
 392                      IJP = IJP + 1
 393                   END DO
 394                END DO
 395 *
 396             END IF
 397 *
 398          END IF
 399 *
 400       ELSE
 401 *
 402 *        N is even
 403 *
 404          IF( NORMALTRANSR ) THEN
 405 *
 406 *           N is even and TRANSR = 'N'
 407 *
 408             IF( LOWER ) THEN
 409 *
 410 *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
 411 *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
 412 *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
 413 *
 414                IJP = 0
 415                JP = 0
 416                DO J = 0, K - 1
 417                   DO I = J, N - 1
 418                      IJ = 1 + I + JP
 419                      AP( IJP ) = ARF( IJ )
 420                      IJP = IJP + 1
 421                   END DO
 422                   JP = JP + LDA
 423                END DO
 424                DO I = 0, K - 1
 425                   DO J = I, K - 1
 426                      IJ = I + J*LDA
 427                      AP( IJP ) = ARF( IJ )
 428                      IJP = IJP + 1
 429                   END DO
 430                END DO
 431 *
 432             ELSE
 433 *
 434 *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
 435 *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
 436 *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
 437 *
 438                IJP = 0
 439                DO J = 0, K - 1
 440                   IJ = K + 1 + J
 441                   DO I = 0, J
 442                      AP( IJP ) = ARF( IJ )
 443                      IJP = IJP + 1
 444                      IJ = IJ + LDA
 445                   END DO
 446                END DO
 447                JS = 0
 448                DO J = K, N - 1
 449                   IJ = JS
 450                   DO IJ = JS, JS + J
 451                      AP( IJP ) = ARF( IJ )
 452                      IJP = IJP + 1
 453                   END DO
 454                   JS = JS + LDA
 455                END DO
 456 *
 457             END IF
 458 *
 459          ELSE
 460 *
 461 *           N is even and TRANSR = 'T'
 462 *
 463             IF( LOWER ) THEN
 464 *
 465 *              SRPA for LOWER, TRANSPOSE and N is even (see paper)
 466 *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
 467 *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
 468 *
 469                IJP = 0
 470                DO I = 0, K - 1
 471                   DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA
 472                      AP( IJP ) = ARF( IJ )
 473                      IJP = IJP + 1
 474                   END DO
 475                END DO
 476                JS = 0
 477                DO J = 0, K - 1
 478                   DO IJ = JS, JS + K - J - 1
 479                      AP( IJP ) = ARF( IJ )
 480                      IJP = IJP + 1
 481                   END DO
 482                   JS = JS + LDA + 1
 483                END DO
 484 *
 485             ELSE
 486 *
 487 *              SRPA for UPPER, TRANSPOSE and N is even (see paper)
 488 *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
 489 *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
 490 *
 491                IJP = 0
 492                JS = ( K+1 )*LDA
 493                DO J = 0, K - 1
 494                   DO IJ = JS, JS + J
 495                      AP( IJP ) = ARF( IJ )
 496                      IJP = IJP + 1
 497                   END DO
 498                   JS = JS + LDA
 499                END DO
 500                DO I = 0, K - 1
 501                   DO IJ = I, I + ( K+I )*LDA, LDA
 502                      AP( IJP ) = ARF( IJ )
 503                      IJP = IJP + 1
 504                   END DO
 505                END DO
 506 *
 507             END IF
 508 *
 509          END IF
 510 *
 511       END IF
 512 *
 513       RETURN
 514 *
 515 *     End of DTFTTP
 516 *
 517       END