Return to dtfttp.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    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