File:  [local] / rpl / lapack / lapack / ztpttf.f
Revision 1.15: download - view: text, annotated - select for diffs - revision graph
Sat Jun 17 11:07:05 2017 UTC (6 years, 10 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_27, rpl-4_1_26, HEAD
Cohérence.

    1: *> \brief \b ZTPTTF 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 ZTPTTF + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztpttf.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztpttf.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztpttf.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZTPTTF( TRANSR, UPLO, N, AP, ARF, 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: *> \par Purpose:
   32: *  =============
   33: *>
   34: *> \verbatim
   35: *>
   36: *> ZTPTTF 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*16 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*16 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 December 2016
   96: *
   97: *> \ingroup complex16OTHERcomputational
   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 ZTPTTF( TRANSR, UPLO, N, AP, ARF, INFO )
  209: *
  210: *  -- LAPACK computational routine (version 3.7.0) --
  211: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  212: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  213: *     December 2016
  214: *
  215: *     .. Scalar Arguments ..
  216:       CHARACTER          TRANSR, UPLO
  217:       INTEGER            INFO, N
  218: *     ..
  219: *     .. Array Arguments ..
  220:       COMPLEX*16         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          DCONJG, 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( 'ZTPTTF', -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 ) = DCONJG( 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 ) = DCONJG( 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 ) = DCONJG( 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 ) = DCONJG( 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 ) = DCONJG( 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 ) = DCONJG( 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 ) = DCONJG( 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 ) = DCONJG( 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 ) = DCONJG( 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 ZTPTTF
  538: *
  539:       END

CVSweb interface <joel.bertrand@systella.fr>