Annotation of rpl/lapack/lapack/dpftri.f, revision 1.10

1.7       bertrand    1: *> \brief \b DPFTRI
                      2: *
                      3: *  =========== DOCUMENTATION ===========
1.1       bertrand    4: *
1.7       bertrand    5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
1.1       bertrand    7: *
1.7       bertrand    8: *> \htmlonly
                      9: *> Download DPFTRI + dependencies 
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpftri.f"> 
                     11: *> [TGZ]</a> 
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpftri.f"> 
                     13: *> [ZIP]</a> 
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpftri.f"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DPFTRI( TRANSR, UPLO, N, A, INFO )
                     22: * 
                     23: *       .. Scalar Arguments ..
                     24: *       CHARACTER          TRANSR, UPLO
                     25: *       INTEGER            INFO, N
                     26: *       .. Array Arguments ..
                     27: *       DOUBLE PRECISION         A( 0: * )
                     28: *       ..
                     29: *  
                     30: *
                     31: *> \par Purpose:
                     32: *  =============
                     33: *>
                     34: *> \verbatim
                     35: *>
                     36: *> DPFTRI computes the inverse of a (real) symmetric positive definite
                     37: *> matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
                     38: *> computed by DPFTRF.
                     39: *> \endverbatim
                     40: *
                     41: *  Arguments:
                     42: *  ==========
                     43: *
                     44: *> \param[in] TRANSR
                     45: *> \verbatim
                     46: *>          TRANSR is CHARACTER*1
                     47: *>          = 'N':  The Normal TRANSR of RFP A is stored;
                     48: *>          = 'T':  The Transpose TRANSR of RFP A is stored.
                     49: *> \endverbatim
                     50: *>
                     51: *> \param[in] UPLO
                     52: *> \verbatim
                     53: *>          UPLO is CHARACTER*1
                     54: *>          = 'U':  Upper triangle of A is stored;
                     55: *>          = 'L':  Lower triangle of A is stored.
                     56: *> \endverbatim
                     57: *>
                     58: *> \param[in] N
                     59: *> \verbatim
                     60: *>          N is INTEGER
                     61: *>          The order of the matrix A.  N >= 0.
                     62: *> \endverbatim
                     63: *>
                     64: *> \param[in,out] A
                     65: *> \verbatim
                     66: *>          A is DOUBLE PRECISION array, dimension ( N*(N+1)/2 )
                     67: *>          On entry, the symmetric matrix A in RFP format. RFP format is
                     68: *>          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
                     69: *>          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
                     70: *>          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
                     71: *>          the transpose of RFP A as defined when
                     72: *>          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
                     73: *>          follows: If UPLO = 'U' the RFP A contains the nt elements of
                     74: *>          upper packed A. If UPLO = 'L' the RFP A contains the elements
                     75: *>          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
                     76: *>          'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
                     77: *>          is odd. See the Note below for more details.
                     78: *>
                     79: *>          On exit, the symmetric inverse of the original matrix, in the
                     80: *>          same storage format.
                     81: *> \endverbatim
                     82: *>
                     83: *> \param[out] INFO
                     84: *> \verbatim
                     85: *>          INFO is INTEGER
                     86: *>          = 0:  successful exit
                     87: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
                     88: *>          > 0:  if INFO = i, the (i,i) element of the factor U or L is
                     89: *>                zero, and the inverse could not be computed.
                     90: *> \endverbatim
                     91: *
                     92: *  Authors:
                     93: *  ========
                     94: *
                     95: *> \author Univ. of Tennessee 
                     96: *> \author Univ. of California Berkeley 
                     97: *> \author Univ. of Colorado Denver 
                     98: *> \author NAG Ltd. 
                     99: *
                    100: *> \date November 2011
                    101: *
                    102: *> \ingroup doubleOTHERcomputational
                    103: *
                    104: *> \par Further Details:
                    105: *  =====================
                    106: *>
                    107: *> \verbatim
                    108: *>
                    109: *>  We first consider Rectangular Full Packed (RFP) Format when N is
                    110: *>  even. We give an example where N = 6.
                    111: *>
                    112: *>      AP is Upper             AP is Lower
                    113: *>
                    114: *>   00 01 02 03 04 05       00
                    115: *>      11 12 13 14 15       10 11
                    116: *>         22 23 24 25       20 21 22
                    117: *>            33 34 35       30 31 32 33
                    118: *>               44 45       40 41 42 43 44
                    119: *>                  55       50 51 52 53 54 55
                    120: *>
                    121: *>
                    122: *>  Let TRANSR = 'N'. RFP holds AP as follows:
                    123: *>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
                    124: *>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
                    125: *>  the transpose of the first three columns of AP upper.
                    126: *>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
                    127: *>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
                    128: *>  the transpose of the last three columns of AP lower.
                    129: *>  This covers the case N even and TRANSR = 'N'.
                    130: *>
                    131: *>         RFP A                   RFP A
                    132: *>
                    133: *>        03 04 05                33 43 53
                    134: *>        13 14 15                00 44 54
                    135: *>        23 24 25                10 11 55
                    136: *>        33 34 35                20 21 22
                    137: *>        00 44 45                30 31 32
                    138: *>        01 11 55                40 41 42
                    139: *>        02 12 22                50 51 52
                    140: *>
                    141: *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
                    142: *>  transpose of RFP A above. One therefore gets:
                    143: *>
                    144: *>
                    145: *>           RFP A                   RFP A
                    146: *>
                    147: *>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
                    148: *>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
                    149: *>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
                    150: *>
                    151: *>
                    152: *>  We then consider Rectangular Full Packed (RFP) Format when N is
                    153: *>  odd. We give an example where N = 5.
                    154: *>
                    155: *>     AP is Upper                 AP is Lower
                    156: *>
                    157: *>   00 01 02 03 04              00
                    158: *>      11 12 13 14              10 11
                    159: *>         22 23 24              20 21 22
                    160: *>            33 34              30 31 32 33
                    161: *>               44              40 41 42 43 44
                    162: *>
                    163: *>
                    164: *>  Let TRANSR = 'N'. RFP holds AP as follows:
                    165: *>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
                    166: *>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
                    167: *>  the transpose of the first two columns of AP upper.
                    168: *>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
                    169: *>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
                    170: *>  the transpose of the last two columns of AP lower.
                    171: *>  This covers the case N odd and TRANSR = 'N'.
                    172: *>
                    173: *>         RFP A                   RFP A
                    174: *>
                    175: *>        02 03 04                00 33 43
                    176: *>        12 13 14                10 11 44
                    177: *>        22 23 24                20 21 22
                    178: *>        00 33 34                30 31 32
                    179: *>        01 11 44                40 41 42
                    180: *>
                    181: *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
                    182: *>  transpose of RFP A above. One therefore gets:
                    183: *>
                    184: *>           RFP A                   RFP A
                    185: *>
                    186: *>     02 12 22 00 01             00 10 20 30 40 50
                    187: *>     03 13 23 33 11             33 11 21 31 41 51
                    188: *>     04 14 24 34 44             43 44 22 32 42 52
                    189: *> \endverbatim
                    190: *>
                    191: *  =====================================================================
                    192:       SUBROUTINE DPFTRI( TRANSR, UPLO, N, A, INFO )
1.1       bertrand  193: *
1.7       bertrand  194: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  195: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    196: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.7       bertrand  197: *     November 2011
1.1       bertrand  198: *
                    199: *     .. Scalar Arguments ..
                    200:       CHARACTER          TRANSR, UPLO
                    201:       INTEGER            INFO, N
                    202: *     .. Array Arguments ..
                    203:       DOUBLE PRECISION         A( 0: * )
                    204: *     ..
                    205: *
                    206: *  =====================================================================
                    207: *
                    208: *     .. Parameters ..
                    209:       DOUBLE PRECISION   ONE
                    210:       PARAMETER          ( ONE = 1.0D+0 )
                    211: *     ..
                    212: *     .. Local Scalars ..
                    213:       LOGICAL            LOWER, NISODD, NORMALTRANSR
                    214:       INTEGER            N1, N2, K
                    215: *     ..
                    216: *     .. External Functions ..
                    217:       LOGICAL            LSAME
                    218:       EXTERNAL           LSAME
                    219: *     ..
                    220: *     .. External Subroutines ..
                    221:       EXTERNAL           XERBLA, DTFTRI, DLAUUM, DTRMM, DSYRK
                    222: *     ..
                    223: *     .. Intrinsic Functions ..
                    224:       INTRINSIC          MOD
                    225: *     ..
                    226: *     .. Executable Statements ..
                    227: *
                    228: *     Test the input parameters.
                    229: *
                    230:       INFO = 0
                    231:       NORMALTRANSR = LSAME( TRANSR, 'N' )
                    232:       LOWER = LSAME( UPLO, 'L' )
                    233:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
                    234:          INFO = -1
                    235:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
                    236:          INFO = -2
                    237:       ELSE IF( N.LT.0 ) THEN
                    238:          INFO = -3
                    239:       END IF
                    240:       IF( INFO.NE.0 ) THEN
                    241:          CALL XERBLA( 'DPFTRI', -INFO )
                    242:          RETURN
                    243:       END IF
                    244: *
                    245: *     Quick return if possible
                    246: *
                    247:       IF( N.EQ.0 )
1.6       bertrand  248:      $   RETURN
1.1       bertrand  249: *
                    250: *     Invert the triangular Cholesky factor U or L.
                    251: *
                    252:       CALL DTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
                    253:       IF( INFO.GT.0 )
1.6       bertrand  254:      $   RETURN
1.1       bertrand  255: *
                    256: *     If N is odd, set NISODD = .TRUE.
                    257: *     If N is even, set K = N/2 and NISODD = .FALSE.
                    258: *
                    259:       IF( MOD( N, 2 ).EQ.0 ) THEN
                    260:          K = N / 2
                    261:          NISODD = .FALSE.
                    262:       ELSE
                    263:          NISODD = .TRUE.
                    264:       END IF
                    265: *
                    266: *     Set N1 and N2 depending on LOWER
                    267: *
                    268:       IF( LOWER ) THEN
                    269:          N2 = N / 2
                    270:          N1 = N - N2
                    271:       ELSE
                    272:          N1 = N / 2
                    273:          N2 = N - N1
                    274:       END IF
                    275: *
                    276: *     Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
                    277: *     inv(L)^C*inv(L). There are eight cases.
                    278: *
                    279:       IF( NISODD ) THEN
                    280: *
                    281: *        N is odd
                    282: *
                    283:          IF( NORMALTRANSR ) THEN
                    284: *
                    285: *           N is odd and TRANSR = 'N'
                    286: *
                    287:             IF( LOWER ) THEN
                    288: *
                    289: *              SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
                    290: *              T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
                    291: *              T1 -> a(0), T2 -> a(n), S -> a(N1)
                    292: *
                    293:                CALL DLAUUM( 'L', N1, A( 0 ), N, INFO )
                    294:                CALL DSYRK( 'L', 'T', N1, N2, ONE, A( N1 ), N, ONE,
1.6       bertrand  295:      $                     A( 0 ), N )
1.1       bertrand  296:                CALL DTRMM( 'L', 'U', 'N', 'N', N2, N1, ONE, A( N ), N,
1.6       bertrand  297:      $                     A( N1 ), N )
1.1       bertrand  298:                CALL DLAUUM( 'U', N2, A( N ), N, INFO )
                    299: *
                    300:             ELSE
                    301: *
                    302: *              SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
                    303: *              T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
                    304: *              T1 -> a(N2), T2 -> a(N1), S -> a(0)
                    305: *
                    306:                CALL DLAUUM( 'L', N1, A( N2 ), N, INFO )
                    307:                CALL DSYRK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
1.6       bertrand  308:      $                     A( N2 ), N )
1.1       bertrand  309:                CALL DTRMM( 'R', 'U', 'T', 'N', N1, N2, ONE, A( N1 ), N,
1.6       bertrand  310:      $                     A( 0 ), N )
1.1       bertrand  311:                CALL DLAUUM( 'U', N2, A( N1 ), N, INFO )
                    312: *
                    313:             END IF
                    314: *
                    315:          ELSE
                    316: *
                    317: *           N is odd and TRANSR = 'T'
                    318: *
                    319:             IF( LOWER ) THEN
                    320: *
                    321: *              SRPA for LOWER, TRANSPOSE, and N is odd
                    322: *              T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
                    323: *
                    324:                CALL DLAUUM( 'U', N1, A( 0 ), N1, INFO )
                    325:                CALL DSYRK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
1.6       bertrand  326:      $                     A( 0 ), N1 )
1.1       bertrand  327:                CALL DTRMM( 'R', 'L', 'N', 'N', N1, N2, ONE, A( 1 ), N1,
1.6       bertrand  328:      $                     A( N1*N1 ), N1 )
1.1       bertrand  329:                CALL DLAUUM( 'L', N2, A( 1 ), N1, INFO )
                    330: *
                    331:             ELSE
                    332: *
                    333: *              SRPA for UPPER, TRANSPOSE, and N is odd
                    334: *              T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
                    335: *
                    336:                CALL DLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
                    337:                CALL DSYRK( 'U', 'T', N1, N2, ONE, A( 0 ), N2, ONE,
1.6       bertrand  338:      $                     A( N2*N2 ), N2 )
1.1       bertrand  339:                CALL DTRMM( 'L', 'L', 'T', 'N', N2, N1, ONE, A( N1*N2 ),
1.6       bertrand  340:      $                     N2, A( 0 ), N2 )
1.1       bertrand  341:                CALL DLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
                    342: *
                    343:             END IF
                    344: *
                    345:          END IF
                    346: *
                    347:       ELSE
                    348: *
                    349: *        N is even
                    350: *
                    351:          IF( NORMALTRANSR ) THEN
                    352: *
                    353: *           N is even and TRANSR = 'N'
                    354: *
                    355:             IF( LOWER ) THEN
                    356: *
                    357: *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
                    358: *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
                    359: *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
                    360: *
                    361:                CALL DLAUUM( 'L', K, A( 1 ), N+1, INFO )
                    362:                CALL DSYRK( 'L', 'T', K, K, ONE, A( K+1 ), N+1, ONE,
1.6       bertrand  363:      $                     A( 1 ), N+1 )
1.1       bertrand  364:                CALL DTRMM( 'L', 'U', 'N', 'N', K, K, ONE, A( 0 ), N+1,
1.6       bertrand  365:      $                     A( K+1 ), N+1 )
1.1       bertrand  366:                CALL DLAUUM( 'U', K, A( 0 ), N+1, INFO )
                    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:                CALL DLAUUM( 'L', K, A( K+1 ), N+1, INFO )
                    375:                CALL DSYRK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
1.6       bertrand  376:      $                     A( K+1 ), N+1 )
1.1       bertrand  377:                CALL DTRMM( 'R', 'U', 'T', 'N', K, K, ONE, A( K ), N+1,
1.6       bertrand  378:      $                     A( 0 ), N+1 )
1.1       bertrand  379:                CALL DLAUUM( 'U', K, A( K ), N+1, INFO )
                    380: *
                    381:             END IF
                    382: *
                    383:          ELSE
                    384: *
                    385: *           N is even and TRANSR = 'T'
                    386: *
                    387:             IF( LOWER ) THEN
                    388: *
                    389: *              SRPA for LOWER, TRANSPOSE, and N is even (see paper)
                    390: *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
                    391: *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
                    392: *
                    393:                CALL DLAUUM( 'U', K, A( K ), K, INFO )
                    394:                CALL DSYRK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
1.6       bertrand  395:      $                     A( K ), K )
1.1       bertrand  396:                CALL DTRMM( 'R', 'L', 'N', 'N', K, K, ONE, A( 0 ), K,
1.6       bertrand  397:      $                     A( K*( K+1 ) ), K )
1.1       bertrand  398:                CALL DLAUUM( 'L', K, A( 0 ), K, INFO )
                    399: *
                    400:             ELSE
                    401: *
                    402: *              SRPA for UPPER, TRANSPOSE, and N is even (see paper)
                    403: *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0),
                    404: *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
                    405: *
                    406:                CALL DLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
                    407:                CALL DSYRK( 'U', 'T', K, K, ONE, A( 0 ), K, ONE,
1.6       bertrand  408:      $                     A( K*( K+1 ) ), K )
1.1       bertrand  409:                CALL DTRMM( 'L', 'L', 'T', 'N', K, K, ONE, A( K*K ), K,
1.6       bertrand  410:      $                     A( 0 ), K )
1.1       bertrand  411:                CALL DLAUUM( 'L', K, A( K*K ), K, INFO )
                    412: *
                    413:             END IF
                    414: *
                    415:          END IF
                    416: *
                    417:       END IF
                    418: *
                    419:       RETURN
                    420: *
                    421: *     End of DPFTRI
                    422: *
                    423:       END

CVSweb interface <joel.bertrand@systella.fr>