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

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

CVSweb interface <joel.bertrand@systella.fr>