SRC/dtpttf.f

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