File:  [local] / rpl / lapack / lapack / dtrttf.f
Revision 1.17: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:14 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

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

CVSweb interface <joel.bertrand@systella.fr>