Annotation of rpl/lapack/lapack/dtpttf.f, revision 1.14

1.10      bertrand    1: *> \brief \b DTPTTF copies a triangular matrix from the standard packed format (TP) to the rectangular full packed format (TF).
1.1       bertrand    2: *
1.7       bertrand    3: *  =========== DOCUMENTATION ===========
1.1       bertrand    4: *
1.14    ! bertrand    5: * Online html documentation available at
        !             6: *            http://www.netlib.org/lapack/explore-html/
1.1       bertrand    7: *
1.7       bertrand    8: *> \htmlonly
1.14    ! bertrand    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">
1.7       bertrand   15: *> [TXT]</a>
1.14    ! bertrand   16: *> \endhtmlonly
1.7       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DTPTTF( TRANSR, UPLO, N, AP, ARF, INFO )
1.14    ! bertrand   22: *
1.7       bertrand   23: *       .. Scalar Arguments ..
                     24: *       CHARACTER          TRANSR, UPLO
                     25: *       INTEGER            INFO, N
                     26: *       ..
                     27: *       .. Array Arguments ..
                     28: *       DOUBLE PRECISION   AP( 0: * ), ARF( 0: * )
1.14    ! bertrand   29: *
1.7       bertrand   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: *
1.14    ! bertrand   90: *> \author Univ. of Tennessee
        !            91: *> \author Univ. of California Berkeley
        !            92: *> \author Univ. of Colorado Denver
        !            93: *> \author NAG Ltd.
1.7       bertrand   94: *
1.14    ! bertrand   95: *> \date December 2016
1.7       bertrand   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: *
1.14    ! bertrand  189: *  -- LAPACK computational routine (version 3.7.0) --
1.1       bertrand  190: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    191: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.14    ! bertrand  192: *     December 2016
1.1       bertrand  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 )
1.6       bertrand  243:      $   RETURN
1.1       bertrand  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 )
1.6       bertrand  286:      $   LDA = ( N+1 ) / 2
1.1       bertrand  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

CVSweb interface <joel.bertrand@systella.fr>