File:  [local] / rpl / lapack / lapack / dtrttf.f
Revision 1.15: download - view: text, annotated - select for diffs - revision graph
Sat Jun 17 11:06:38 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 DTRTTF copies a triangular matrix from the standard full format (TR) 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 DTRTTF + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrttf.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrttf.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrttf.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DTRTTF( TRANSR, UPLO, N, A, LDA, ARF, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       CHARACTER          TRANSR, UPLO
   25: *       INTEGER            INFO, N, LDA
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       DOUBLE PRECISION   A( 0: LDA-1, 0: * ), ARF( 0: * )
   29: *       ..
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> DTRTTF copies a triangular matrix A from standard full format (TR)
   38: *> to rectangular full packed format (TF) .
   39: *> \endverbatim
   40: *
   41: *  Arguments:
   42: *  ==========
   43: *
   44: *> \param[in] TRANSR
   45: *> \verbatim
   46: *>          TRANSR is CHARACTER*1
   47: *>          = 'N':  ARF in Normal form is wanted;
   48: *>          = 'T':  ARF in Transpose form is wanted.
   49: *> \endverbatim
   50: *>
   51: *> \param[in] UPLO
   52: *> \verbatim
   53: *>          UPLO is CHARACTER*1
   54: *>          = 'U':  Upper triangle of A is stored;
   55: *>          = 'L':  Lower triangle of A is stored.
   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] A
   65: *> \verbatim
   66: *>          A is DOUBLE PRECISION array, dimension (LDA,N).
   67: *>          On entry, the triangular matrix A.  If UPLO = 'U', the
   68: *>          leading N-by-N upper triangular part of the array A contains
   69: *>          the upper triangular matrix, and the strictly lower
   70: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   71: *>          leading N-by-N lower triangular part of the array A contains
   72: *>          the lower triangular matrix, and the strictly upper
   73: *>          triangular part of A is not referenced.
   74: *> \endverbatim
   75: *>
   76: *> \param[in] LDA
   77: *> \verbatim
   78: *>          LDA is INTEGER
   79: *>          The leading dimension of the matrix A. LDA >= max(1,N).
   80: *> \endverbatim
   81: *>
   82: *> \param[out] ARF
   83: *> \verbatim
   84: *>          ARF is DOUBLE PRECISION array, dimension (NT).
   85: *>          NT=N*(N+1)/2. On exit, the triangular matrix A in RFP format.
   86: *> \endverbatim
   87: *>
   88: *> \param[out] INFO
   89: *> \verbatim
   90: *>          INFO is INTEGER
   91: *>          = 0:  successful exit
   92: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
   93: *> \endverbatim
   94: *
   95: *  Authors:
   96: *  ========
   97: *
   98: *> \author Univ. of Tennessee
   99: *> \author Univ. of California Berkeley
  100: *> \author Univ. of Colorado Denver
  101: *> \author NAG Ltd.
  102: *
  103: *> \date December 2016
  104: *
  105: *> \ingroup doubleOTHERcomputational
  106: *
  107: *> \par Further Details:
  108: *  =====================
  109: *>
  110: *> \verbatim
  111: *>
  112: *>  We first consider Rectangular Full Packed (RFP) Format when N is
  113: *>  even. We give an example where N = 6.
  114: *>
  115: *>      AP is Upper             AP is Lower
  116: *>
  117: *>   00 01 02 03 04 05       00
  118: *>      11 12 13 14 15       10 11
  119: *>         22 23 24 25       20 21 22
  120: *>            33 34 35       30 31 32 33
  121: *>               44 45       40 41 42 43 44
  122: *>                  55       50 51 52 53 54 55
  123: *>
  124: *>
  125: *>  Let TRANSR = 'N'. RFP holds AP as follows:
  126: *>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  127: *>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  128: *>  the transpose of the first three columns of AP upper.
  129: *>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  130: *>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  131: *>  the transpose of the last three columns of AP lower.
  132: *>  This covers the case N even and TRANSR = 'N'.
  133: *>
  134: *>         RFP A                   RFP A
  135: *>
  136: *>        03 04 05                33 43 53
  137: *>        13 14 15                00 44 54
  138: *>        23 24 25                10 11 55
  139: *>        33 34 35                20 21 22
  140: *>        00 44 45                30 31 32
  141: *>        01 11 55                40 41 42
  142: *>        02 12 22                50 51 52
  143: *>
  144: *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  145: *>  transpose of RFP A above. One therefore gets:
  146: *>
  147: *>
  148: *>           RFP A                   RFP A
  149: *>
  150: *>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
  151: *>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
  152: *>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
  153: *>
  154: *>
  155: *>  We then consider Rectangular Full Packed (RFP) Format when N is
  156: *>  odd. We give an example where N = 5.
  157: *>
  158: *>     AP is Upper                 AP is Lower
  159: *>
  160: *>   00 01 02 03 04              00
  161: *>      11 12 13 14              10 11
  162: *>         22 23 24              20 21 22
  163: *>            33 34              30 31 32 33
  164: *>               44              40 41 42 43 44
  165: *>
  166: *>
  167: *>  Let TRANSR = 'N'. RFP holds AP as follows:
  168: *>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  169: *>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  170: *>  the transpose of the first two columns of AP upper.
  171: *>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  172: *>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  173: *>  the transpose of the last two columns of AP lower.
  174: *>  This covers the case N odd and TRANSR = 'N'.
  175: *>
  176: *>         RFP A                   RFP A
  177: *>
  178: *>        02 03 04                00 33 43
  179: *>        12 13 14                10 11 44
  180: *>        22 23 24                20 21 22
  181: *>        00 33 34                30 31 32
  182: *>        01 11 44                40 41 42
  183: *>
  184: *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  185: *>  transpose of RFP A above. One therefore gets:
  186: *>
  187: *>           RFP A                   RFP A
  188: *>
  189: *>     02 12 22 00 01             00 10 20 30 40 50
  190: *>     03 13 23 33 11             33 11 21 31 41 51
  191: *>     04 14 24 34 44             43 44 22 32 42 52
  192: *> \endverbatim
  193: *
  194: *  =====================================================================
  195:       SUBROUTINE DTRTTF( TRANSR, UPLO, N, A, LDA, ARF, INFO )
  196: *
  197: *  -- LAPACK computational routine (version 3.7.0) --
  198: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  199: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  200: *     December 2016
  201: *
  202: *     .. Scalar Arguments ..
  203:       CHARACTER          TRANSR, UPLO
  204:       INTEGER            INFO, N, LDA
  205: *     ..
  206: *     .. Array Arguments ..
  207:       DOUBLE PRECISION   A( 0: LDA-1, 0: * ), ARF( 0: * )
  208: *     ..
  209: *
  210: *  =====================================================================
  211: *
  212: *     ..
  213: *     .. Local Scalars ..
  214:       LOGICAL            LOWER, NISODD, NORMALTRANSR
  215:       INTEGER            I, IJ, J, K, L, N1, N2, NT, NX2, NP1X2
  216: *     ..
  217: *     .. External Functions ..
  218:       LOGICAL            LSAME
  219:       EXTERNAL           LSAME
  220: *     ..
  221: *     .. External Subroutines ..
  222:       EXTERNAL           XERBLA
  223: *     ..
  224: *     .. Intrinsic Functions ..
  225:       INTRINSIC          MAX, MOD
  226: *     ..
  227: *     .. Executable Statements ..
  228: *
  229: *     Test the input parameters.
  230: *
  231:       INFO = 0
  232:       NORMALTRANSR = LSAME( TRANSR, 'N' )
  233:       LOWER = LSAME( UPLO, 'L' )
  234:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
  235:          INFO = -1
  236:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
  237:          INFO = -2
  238:       ELSE IF( N.LT.0 ) THEN
  239:          INFO = -3
  240:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  241:          INFO = -5
  242:       END IF
  243:       IF( INFO.NE.0 ) THEN
  244:          CALL XERBLA( 'DTRTTF', -INFO )
  245:          RETURN
  246:       END IF
  247: *
  248: *     Quick return if possible
  249: *
  250:       IF( N.LE.1 ) THEN
  251:          IF( N.EQ.1 ) THEN
  252:             ARF( 0 ) = A( 0, 0 )
  253:          END IF
  254:          RETURN
  255:       END IF
  256: *
  257: *     Size of array ARF(0:nt-1)
  258: *
  259:       NT = N*( N+1 ) / 2
  260: *
  261: *     Set N1 and N2 depending on LOWER: for N even N1=N2=K
  262: *
  263:       IF( LOWER ) THEN
  264:          N2 = N / 2
  265:          N1 = N - N2
  266:       ELSE
  267:          N1 = N / 2
  268:          N2 = N - N1
  269:       END IF
  270: *
  271: *     If N is odd, set NISODD = .TRUE., LDA=N+1 and A is (N+1)--by--K2.
  272: *     If N is even, set K = N/2 and NISODD = .FALSE., LDA=N and A is
  273: *     N--by--(N+1)/2.
  274: *
  275:       IF( MOD( N, 2 ).EQ.0 ) THEN
  276:          K = N / 2
  277:          NISODD = .FALSE.
  278:          IF( .NOT.LOWER )
  279:      $      NP1X2 = N + N + 2
  280:       ELSE
  281:          NISODD = .TRUE.
  282:          IF( .NOT.LOWER )
  283:      $      NX2 = N + N
  284:       END IF
  285: *
  286:       IF( NISODD ) THEN
  287: *
  288: *        N is odd
  289: *
  290:          IF( NORMALTRANSR ) THEN
  291: *
  292: *           N is odd and TRANSR = 'N'
  293: *
  294:             IF( LOWER ) THEN
  295: *
  296: *              N is odd, TRANSR = 'N', and UPLO = 'L'
  297: *
  298:                IJ = 0
  299:                DO J = 0, N2
  300:                   DO I = N1, N2 + J
  301:                      ARF( IJ ) = A( N2+J, I )
  302:                      IJ = IJ + 1
  303:                   END DO
  304:                   DO I = J, N - 1
  305:                      ARF( IJ ) = A( I, J )
  306:                      IJ = IJ + 1
  307:                   END DO
  308:                END DO
  309: *
  310:             ELSE
  311: *
  312: *              N is odd, TRANSR = 'N', and UPLO = 'U'
  313: *
  314:                IJ = NT - N
  315:                DO J = N - 1, N1, -1
  316:                   DO I = 0, J
  317:                      ARF( IJ ) = A( I, J )
  318:                      IJ = IJ + 1
  319:                   END DO
  320:                   DO L = J - N1, N1 - 1
  321:                      ARF( IJ ) = A( J-N1, L )
  322:                      IJ = IJ + 1
  323:                   END DO
  324:                   IJ = IJ - NX2
  325:                END DO
  326: *
  327:             END IF
  328: *
  329:          ELSE
  330: *
  331: *           N is odd and TRANSR = 'T'
  332: *
  333:             IF( LOWER ) THEN
  334: *
  335: *              N is odd, TRANSR = 'T', and UPLO = 'L'
  336: *
  337:                IJ = 0
  338:                DO J = 0, N2 - 1
  339:                   DO I = 0, J
  340:                      ARF( IJ ) = A( J, I )
  341:                      IJ = IJ + 1
  342:                   END DO
  343:                   DO I = N1 + J, N - 1
  344:                      ARF( IJ ) = A( I, N1+J )
  345:                      IJ = IJ + 1
  346:                   END DO
  347:                END DO
  348:                DO J = N2, N - 1
  349:                   DO I = 0, N1 - 1
  350:                      ARF( IJ ) = A( J, I )
  351:                      IJ = IJ + 1
  352:                   END DO
  353:                END DO
  354: *
  355:             ELSE
  356: *
  357: *              N is odd, TRANSR = 'T', and UPLO = 'U'
  358: *
  359:                IJ = 0
  360:                DO J = 0, N1
  361:                   DO I = N1, N - 1
  362:                      ARF( IJ ) = A( J, I )
  363:                      IJ = IJ + 1
  364:                   END DO
  365:                END DO
  366:                DO J = 0, N1 - 1
  367:                   DO I = 0, J
  368:                      ARF( IJ ) = A( I, J )
  369:                      IJ = IJ + 1
  370:                   END DO
  371:                   DO L = N2 + J, N - 1
  372:                      ARF( IJ ) = A( N2+J, L )
  373:                      IJ = IJ + 1
  374:                   END DO
  375:                END DO
  376: *
  377:             END IF
  378: *
  379:          END IF
  380: *
  381:       ELSE
  382: *
  383: *        N is even
  384: *
  385:          IF( NORMALTRANSR ) THEN
  386: *
  387: *           N is even and TRANSR = 'N'
  388: *
  389:             IF( LOWER ) THEN
  390: *
  391: *              N is even, TRANSR = 'N', and UPLO = 'L'
  392: *
  393:                IJ = 0
  394:                DO J = 0, K - 1
  395:                   DO I = K, K + J
  396:                      ARF( IJ ) = A( K+J, I )
  397:                      IJ = IJ + 1
  398:                   END DO
  399:                   DO I = J, N - 1
  400:                      ARF( IJ ) = A( I, J )
  401:                      IJ = IJ + 1
  402:                   END DO
  403:                END DO
  404: *
  405:             ELSE
  406: *
  407: *              N is even, TRANSR = 'N', and UPLO = 'U'
  408: *
  409:                IJ = NT - N - 1
  410:                DO J = N - 1, K, -1
  411:                   DO I = 0, J
  412:                      ARF( IJ ) = A( I, J )
  413:                      IJ = IJ + 1
  414:                   END DO
  415:                   DO L = J - K, K - 1
  416:                      ARF( IJ ) = A( J-K, L )
  417:                      IJ = IJ + 1
  418:                   END DO
  419:                   IJ = IJ - NP1X2
  420:                END DO
  421: *
  422:             END IF
  423: *
  424:          ELSE
  425: *
  426: *           N is even and TRANSR = 'T'
  427: *
  428:             IF( LOWER ) THEN
  429: *
  430: *              N is even, TRANSR = 'T', and UPLO = 'L'
  431: *
  432:                IJ = 0
  433:                J = K
  434:                DO I = K, N - 1
  435:                   ARF( IJ ) = A( I, J )
  436:                   IJ = IJ + 1
  437:                END DO
  438:                DO J = 0, K - 2
  439:                   DO I = 0, J
  440:                      ARF( IJ ) = A( J, I )
  441:                      IJ = IJ + 1
  442:                   END DO
  443:                   DO I = K + 1 + J, N - 1
  444:                      ARF( IJ ) = A( I, K+1+J )
  445:                      IJ = IJ + 1
  446:                   END DO
  447:                END DO
  448:                DO J = K - 1, N - 1
  449:                   DO I = 0, K - 1
  450:                      ARF( IJ ) = A( J, I )
  451:                      IJ = IJ + 1
  452:                   END DO
  453:                END DO
  454: *
  455:             ELSE
  456: *
  457: *              N is even, TRANSR = 'T', and UPLO = 'U'
  458: *
  459:                IJ = 0
  460:                DO J = 0, K
  461:                   DO I = K, N - 1
  462:                      ARF( IJ ) = A( J, I )
  463:                      IJ = IJ + 1
  464:                   END DO
  465:                END DO
  466:                DO J = 0, K - 2
  467:                   DO I = 0, J
  468:                      ARF( IJ ) = A( I, J )
  469:                      IJ = IJ + 1
  470:                   END DO
  471:                   DO L = K + 1 + J, N - 1
  472:                      ARF( IJ ) = A( K+1+J, L )
  473:                      IJ = IJ + 1
  474:                   END DO
  475:                END DO
  476: *              Note that here, on exit of the loop, J = K-1
  477:                DO I = 0, J
  478:                   ARF( IJ ) = A( I, J )
  479:                   IJ = IJ + 1
  480:                END DO
  481: *
  482:             END IF
  483: *
  484:          END IF
  485: *
  486:       END IF
  487: *
  488:       RETURN
  489: *
  490: *     End of DTRTTF
  491: *
  492:       END

CVSweb interface <joel.bertrand@systella.fr>