Return to dpftri.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence Lapack.

    1: *> \brief \b DPFTRI
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    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 )
  193: *
  194: *  -- LAPACK computational routine (version 3.4.0) --
  195: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  196: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  197: *     November 2011
  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 )
  248:      $   RETURN
  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 )
  254:      $   RETURN
  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,
  295:      $                     A( 0 ), N )
  296:                CALL DTRMM( 'L', 'U', 'N', 'N', N2, N1, ONE, A( N ), N,
  297:      $                     A( N1 ), N )
  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,
  308:      $                     A( N2 ), N )
  309:                CALL DTRMM( 'R', 'U', 'T', 'N', N1, N2, ONE, A( N1 ), N,
  310:      $                     A( 0 ), N )
  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,
  326:      $                     A( 0 ), N1 )
  327:                CALL DTRMM( 'R', 'L', 'N', 'N', N1, N2, ONE, A( 1 ), N1,
  328:      $                     A( N1*N1 ), N1 )
  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,
  338:      $                     A( N2*N2 ), N2 )
  339:                CALL DTRMM( 'L', 'L', 'T', 'N', N2, N1, ONE, A( N1*N2 ),
  340:      $                     N2, A( 0 ), N2 )
  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,
  363:      $                     A( 1 ), N+1 )
  364:                CALL DTRMM( 'L', 'U', 'N', 'N', K, K, ONE, A( 0 ), N+1,
  365:      $                     A( K+1 ), N+1 )
  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,
  376:      $                     A( K+1 ), N+1 )
  377:                CALL DTRMM( 'R', 'U', 'T', 'N', K, K, ONE, A( K ), N+1,
  378:      $                     A( 0 ), N+1 )
  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,
  395:      $                     A( K ), K )
  396:                CALL DTRMM( 'R', 'L', 'N', 'N', K, K, ONE, A( 0 ), K,
  397:      $                     A( K*( K+1 ) ), K )
  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,
  408:      $                     A( K*( K+1 ) ), K )
  409:                CALL DTRMM( 'L', 'L', 'T', 'N', K, K, ONE, A( K*K ), K,
  410:      $                     A( 0 ), K )
  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