Annotation of rpl/lapack/lapack/dtfttp.f, revision 1.1

1.1     ! bertrand    1:       SUBROUTINE DTFTTP( TRANSR, UPLO, N, ARF, AP, INFO )
        !             2: *
        !             3: *  -- LAPACK routine (version 3.2.2)                                    --
        !             4: *
        !             5: *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
        !             6: *  -- June 2010                                                       --
        !             7: *
        !             8: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
        !             9: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
        !            10: *
        !            11: *     ..
        !            12: *     .. Scalar Arguments ..
        !            13:       CHARACTER          TRANSR, UPLO
        !            14:       INTEGER            INFO, N
        !            15: *     ..
        !            16: *     .. Array Arguments ..
        !            17:       DOUBLE PRECISION   AP( 0: * ), ARF( 0: * )
        !            18: *     ..
        !            19: *
        !            20: *  Purpose
        !            21: *  =======
        !            22: *
        !            23: *  DTFTTP copies a triangular matrix A from rectangular full packed
        !            24: *  format (TF) to standard packed format (TP).
        !            25: *
        !            26: *  Arguments
        !            27: *  =========
        !            28: *
        !            29: *  TRANSR  (input) CHARACTER
        !            30: *          = 'N':  ARF is in Normal format;
        !            31: *          = 'T':  ARF is in Transpose format;
        !            32: *
        !            33: *  UPLO    (input) CHARACTER
        !            34: *          = 'U':  A is upper triangular;
        !            35: *          = 'L':  A is lower triangular.
        !            36: *
        !            37: *  N       (input) INTEGER
        !            38: *          The order of the matrix A. N >= 0.
        !            39: *
        !            40: *  ARF     (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ),
        !            41: *          On entry, the upper or lower triangular matrix A stored in
        !            42: *          RFP format. For a further discussion see Notes below.
        !            43: *
        !            44: *  AP      (output) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ),
        !            45: *          On exit, the upper or lower triangular matrix A, packed
        !            46: *          columnwise in a linear array. The j-th column of A is stored
        !            47: *          in the array AP as follows:
        !            48: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
        !            49: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
        !            50: *
        !            51: *  INFO    (output) INTEGER
        !            52: *          = 0:  successful exit
        !            53: *          < 0:  if INFO = -i, the i-th argument had an illegal value
        !            54: *
        !            55: *  Further Details
        !            56: *  ===============
        !            57: *
        !            58: *  We first consider Rectangular Full Packed (RFP) Format when N is
        !            59: *  even. We give an example where N = 6.
        !            60: *
        !            61: *      AP is Upper             AP is Lower
        !            62: *
        !            63: *   00 01 02 03 04 05       00
        !            64: *      11 12 13 14 15       10 11
        !            65: *         22 23 24 25       20 21 22
        !            66: *            33 34 35       30 31 32 33
        !            67: *               44 45       40 41 42 43 44
        !            68: *                  55       50 51 52 53 54 55
        !            69: *
        !            70: *
        !            71: *  Let TRANSR = 'N'. RFP holds AP as follows:
        !            72: *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
        !            73: *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
        !            74: *  the transpose of the first three columns of AP upper.
        !            75: *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
        !            76: *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
        !            77: *  the transpose of the last three columns of AP lower.
        !            78: *  This covers the case N even and TRANSR = 'N'.
        !            79: *
        !            80: *         RFP A                   RFP A
        !            81: *
        !            82: *        03 04 05                33 43 53
        !            83: *        13 14 15                00 44 54
        !            84: *        23 24 25                10 11 55
        !            85: *        33 34 35                20 21 22
        !            86: *        00 44 45                30 31 32
        !            87: *        01 11 55                40 41 42
        !            88: *        02 12 22                50 51 52
        !            89: *
        !            90: *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
        !            91: *  transpose of RFP A above. One therefore gets:
        !            92: *
        !            93: *
        !            94: *           RFP A                   RFP A
        !            95: *
        !            96: *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
        !            97: *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
        !            98: *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
        !            99: *
        !           100: *
        !           101: *  We then consider Rectangular Full Packed (RFP) Format when N is
        !           102: *  odd. We give an example where N = 5.
        !           103: *
        !           104: *     AP is Upper                 AP is Lower
        !           105: *
        !           106: *   00 01 02 03 04              00
        !           107: *      11 12 13 14              10 11
        !           108: *         22 23 24              20 21 22
        !           109: *            33 34              30 31 32 33
        !           110: *               44              40 41 42 43 44
        !           111: *
        !           112: *
        !           113: *  Let TRANSR = 'N'. RFP holds AP as follows:
        !           114: *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
        !           115: *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
        !           116: *  the transpose of the first two columns of AP upper.
        !           117: *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
        !           118: *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
        !           119: *  the transpose of the last two columns of AP lower.
        !           120: *  This covers the case N odd and TRANSR = 'N'.
        !           121: *
        !           122: *         RFP A                   RFP A
        !           123: *
        !           124: *        02 03 04                00 33 43
        !           125: *        12 13 14                10 11 44
        !           126: *        22 23 24                20 21 22
        !           127: *        00 33 34                30 31 32
        !           128: *        01 11 44                40 41 42
        !           129: *
        !           130: *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
        !           131: *  transpose of RFP A above. One therefore gets:
        !           132: *
        !           133: *           RFP A                   RFP A
        !           134: *
        !           135: *     02 12 22 00 01             00 10 20 30 40 50
        !           136: *     03 13 23 33 11             33 11 21 31 41 51
        !           137: *     04 14 24 34 44             43 44 22 32 42 52
        !           138: *
        !           139: *  =====================================================================
        !           140: *
        !           141: *     .. Parameters ..
        !           142: *     ..
        !           143: *     .. Local Scalars ..
        !           144:       LOGICAL            LOWER, NISODD, NORMALTRANSR
        !           145:       INTEGER            N1, N2, K, NT
        !           146:       INTEGER            I, J, IJ
        !           147:       INTEGER            IJP, JP, LDA, JS
        !           148: *     ..
        !           149: *     .. External Functions ..
        !           150:       LOGICAL            LSAME
        !           151:       EXTERNAL           LSAME
        !           152: *     ..
        !           153: *     .. External Subroutines ..
        !           154:       EXTERNAL           XERBLA
        !           155: *     ..
        !           156: *     .. Executable Statements ..
        !           157: *
        !           158: *     Test the input parameters.
        !           159: *
        !           160:       INFO = 0
        !           161:       NORMALTRANSR = LSAME( TRANSR, 'N' )
        !           162:       LOWER = LSAME( UPLO, 'L' )
        !           163:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
        !           164:          INFO = -1
        !           165:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
        !           166:          INFO = -2
        !           167:       ELSE IF( N.LT.0 ) THEN
        !           168:          INFO = -3
        !           169:       END IF
        !           170:       IF( INFO.NE.0 ) THEN
        !           171:          CALL XERBLA( 'DTFTTP', -INFO )
        !           172:          RETURN
        !           173:       END IF
        !           174: *
        !           175: *     Quick return if possible
        !           176: *
        !           177:       IF( N.EQ.0 )
        !           178:      +   RETURN
        !           179: *
        !           180:       IF( N.EQ.1 ) THEN
        !           181:          IF( NORMALTRANSR ) THEN
        !           182:             AP( 0 ) = ARF( 0 )
        !           183:          ELSE
        !           184:             AP( 0 ) = ARF( 0 )
        !           185:          END IF
        !           186:          RETURN
        !           187:       END IF
        !           188: *
        !           189: *     Size of array ARF(0:NT-1)
        !           190: *
        !           191:       NT = N*( N+1 ) / 2
        !           192: *
        !           193: *     Set N1 and N2 depending on LOWER
        !           194: *
        !           195:       IF( LOWER ) THEN
        !           196:          N2 = N / 2
        !           197:          N1 = N - N2
        !           198:       ELSE
        !           199:          N1 = N / 2
        !           200:          N2 = N - N1
        !           201:       END IF
        !           202: *
        !           203: *     If N is odd, set NISODD = .TRUE.
        !           204: *     If N is even, set K = N/2 and NISODD = .FALSE.
        !           205: *
        !           206: *     set lda of ARF^C; ARF^C is (0:(N+1)/2-1,0:N-noe)
        !           207: *     where noe = 0 if n is even, noe = 1 if n is odd
        !           208: *
        !           209:       IF( MOD( N, 2 ).EQ.0 ) THEN
        !           210:          K = N / 2
        !           211:          NISODD = .FALSE.
        !           212:          LDA = N + 1
        !           213:       ELSE
        !           214:          NISODD = .TRUE.
        !           215:          LDA = N
        !           216:       END IF
        !           217: *
        !           218: *     ARF^C has lda rows and n+1-noe cols
        !           219: *
        !           220:       IF( .NOT.NORMALTRANSR )
        !           221:      +   LDA = ( N+1 ) / 2
        !           222: *
        !           223: *     start execution: there are eight cases
        !           224: *
        !           225:       IF( NISODD ) THEN
        !           226: *
        !           227: *        N is odd
        !           228: *
        !           229:          IF( NORMALTRANSR ) THEN
        !           230: *
        !           231: *           N is odd and TRANSR = 'N'
        !           232: *
        !           233:             IF( LOWER ) THEN
        !           234: *
        !           235: *             SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
        !           236: *             T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
        !           237: *             T1 -> a(0), T2 -> a(n), S -> a(n1); lda = n
        !           238: *
        !           239:                IJP = 0
        !           240:                JP = 0
        !           241:                DO J = 0, N2
        !           242:                   DO I = J, N - 1
        !           243:                      IJ = I + JP
        !           244:                      AP( IJP ) = ARF( IJ )
        !           245:                      IJP = IJP + 1
        !           246:                   END DO
        !           247:                   JP = JP + LDA
        !           248:                END DO
        !           249:                DO I = 0, N2 - 1
        !           250:                   DO J = 1 + I, N2
        !           251:                      IJ = I + J*LDA
        !           252:                      AP( IJP ) = ARF( IJ )
        !           253:                      IJP = IJP + 1
        !           254:                   END DO
        !           255:                END DO
        !           256: *
        !           257:             ELSE
        !           258: *
        !           259: *             SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
        !           260: *             T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
        !           261: *             T1 -> a(n2), T2 -> a(n1), S -> a(0)
        !           262: *
        !           263:                IJP = 0
        !           264:                DO J = 0, N1 - 1
        !           265:                   IJ = N2 + J
        !           266:                   DO I = 0, J
        !           267:                      AP( IJP ) = ARF( IJ )
        !           268:                      IJP = IJP + 1
        !           269:                      IJ = IJ + LDA
        !           270:                   END DO
        !           271:                END DO
        !           272:                JS = 0
        !           273:                DO J = N1, N - 1
        !           274:                   IJ = JS
        !           275:                   DO IJ = JS, JS + J
        !           276:                      AP( IJP ) = ARF( IJ )
        !           277:                      IJP = IJP + 1
        !           278:                   END DO
        !           279:                   JS = JS + LDA
        !           280:                END DO
        !           281: *
        !           282:             END IF
        !           283: *
        !           284:          ELSE
        !           285: *
        !           286: *           N is odd and TRANSR = 'T'
        !           287: *
        !           288:             IF( LOWER ) THEN
        !           289: *
        !           290: *              SRPA for LOWER, TRANSPOSE and N is odd
        !           291: *              T1 -> A(0,0) , T2 -> A(1,0) , S -> A(0,n1)
        !           292: *              T1 -> a(0+0) , T2 -> a(1+0) , S -> a(0+n1*n1); lda=n1
        !           293: *
        !           294:                IJP = 0
        !           295:                DO I = 0, N2
        !           296:                   DO IJ = I*( LDA+1 ), N*LDA - 1, LDA
        !           297:                      AP( IJP ) = ARF( IJ )
        !           298:                      IJP = IJP + 1
        !           299:                   END DO
        !           300:                END DO
        !           301:                JS = 1
        !           302:                DO J = 0, N2 - 1
        !           303:                   DO IJ = JS, JS + N2 - J - 1
        !           304:                      AP( IJP ) = ARF( IJ )
        !           305:                      IJP = IJP + 1
        !           306:                   END DO
        !           307:                   JS = JS + LDA + 1
        !           308:                END DO
        !           309: *
        !           310:             ELSE
        !           311: *
        !           312: *              SRPA for UPPER, TRANSPOSE and N is odd
        !           313: *              T1 -> A(0,n1+1), T2 -> A(0,n1), S -> A(0,0)
        !           314: *              T1 -> a(n2*n2), T2 -> a(n1*n2), S -> a(0); lda = n2
        !           315: *
        !           316:                IJP = 0
        !           317:                JS = N2*LDA
        !           318:                DO J = 0, N1 - 1
        !           319:                   DO IJ = JS, JS + J
        !           320:                      AP( IJP ) = ARF( IJ )
        !           321:                      IJP = IJP + 1
        !           322:                   END DO
        !           323:                   JS = JS + LDA
        !           324:                END DO
        !           325:                DO I = 0, N1
        !           326:                   DO IJ = I, I + ( N1+I )*LDA, LDA
        !           327:                      AP( IJP ) = ARF( IJ )
        !           328:                      IJP = IJP + 1
        !           329:                   END DO
        !           330:                END DO
        !           331: *
        !           332:             END IF
        !           333: *
        !           334:          END IF
        !           335: *
        !           336:       ELSE
        !           337: *
        !           338: *        N is even
        !           339: *
        !           340:          IF( NORMALTRANSR ) THEN
        !           341: *
        !           342: *           N is even and TRANSR = 'N'
        !           343: *
        !           344:             IF( LOWER ) THEN
        !           345: *
        !           346: *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
        !           347: *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
        !           348: *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
        !           349: *
        !           350:                IJP = 0
        !           351:                JP = 0
        !           352:                DO J = 0, K - 1
        !           353:                   DO I = J, N - 1
        !           354:                      IJ = 1 + I + JP
        !           355:                      AP( IJP ) = ARF( IJ )
        !           356:                      IJP = IJP + 1
        !           357:                   END DO
        !           358:                   JP = JP + LDA
        !           359:                END DO
        !           360:                DO I = 0, K - 1
        !           361:                   DO J = I, K - 1
        !           362:                      IJ = I + J*LDA
        !           363:                      AP( IJP ) = ARF( IJ )
        !           364:                      IJP = IJP + 1
        !           365:                   END DO
        !           366:                END DO
        !           367: *
        !           368:             ELSE
        !           369: *
        !           370: *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
        !           371: *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
        !           372: *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
        !           373: *
        !           374:                IJP = 0
        !           375:                DO J = 0, K - 1
        !           376:                   IJ = K + 1 + J
        !           377:                   DO I = 0, J
        !           378:                      AP( IJP ) = ARF( IJ )
        !           379:                      IJP = IJP + 1
        !           380:                      IJ = IJ + LDA
        !           381:                   END DO
        !           382:                END DO
        !           383:                JS = 0
        !           384:                DO J = K, N - 1
        !           385:                   IJ = JS
        !           386:                   DO IJ = JS, JS + J
        !           387:                      AP( IJP ) = ARF( IJ )
        !           388:                      IJP = IJP + 1
        !           389:                   END DO
        !           390:                   JS = JS + LDA
        !           391:                END DO
        !           392: *
        !           393:             END IF
        !           394: *
        !           395:          ELSE
        !           396: *
        !           397: *           N is even and TRANSR = 'T'
        !           398: *
        !           399:             IF( LOWER ) THEN
        !           400: *
        !           401: *              SRPA for LOWER, TRANSPOSE and N is even (see paper)
        !           402: *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
        !           403: *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
        !           404: *
        !           405:                IJP = 0
        !           406:                DO I = 0, K - 1
        !           407:                   DO IJ = I + ( I+1 )*LDA, ( N+1 )*LDA - 1, LDA
        !           408:                      AP( IJP ) = ARF( IJ )
        !           409:                      IJP = IJP + 1
        !           410:                   END DO
        !           411:                END DO
        !           412:                JS = 0
        !           413:                DO J = 0, K - 1
        !           414:                   DO IJ = JS, JS + K - J - 1
        !           415:                      AP( IJP ) = ARF( IJ )
        !           416:                      IJP = IJP + 1
        !           417:                   END DO
        !           418:                   JS = JS + LDA + 1
        !           419:                END DO
        !           420: *
        !           421:             ELSE
        !           422: *
        !           423: *              SRPA for UPPER, TRANSPOSE and N is even (see paper)
        !           424: *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0)
        !           425: *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
        !           426: *
        !           427:                IJP = 0
        !           428:                JS = ( K+1 )*LDA
        !           429:                DO J = 0, K - 1
        !           430:                   DO IJ = JS, JS + J
        !           431:                      AP( IJP ) = ARF( IJ )
        !           432:                      IJP = IJP + 1
        !           433:                   END DO
        !           434:                   JS = JS + LDA
        !           435:                END DO
        !           436:                DO I = 0, K - 1
        !           437:                   DO IJ = I, I + ( K+I )*LDA, LDA
        !           438:                      AP( IJP ) = ARF( IJ )
        !           439:                      IJP = IJP + 1
        !           440:                   END DO
        !           441:                END DO
        !           442: *
        !           443:             END IF
        !           444: *
        !           445:          END IF
        !           446: *
        !           447:       END IF
        !           448: *
        !           449:       RETURN
        !           450: *
        !           451: *     End of DTFTTP
        !           452: *
        !           453:       END

CVSweb interface <joel.bertrand@systella.fr>