File:  [local] / rpl / lapack / lapack / dpftri.f
Revision 1.16: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:03 2023 UTC (14 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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: *> \ingroup doubleOTHERcomputational
  101: *
  102: *> \par Further Details:
  103: *  =====================
  104: *>
  105: *> \verbatim
  106: *>
  107: *>  We first consider Rectangular Full Packed (RFP) Format when N is
  108: *>  even. We give an example where N = 6.
  109: *>
  110: *>      AP is Upper             AP is Lower
  111: *>
  112: *>   00 01 02 03 04 05       00
  113: *>      11 12 13 14 15       10 11
  114: *>         22 23 24 25       20 21 22
  115: *>            33 34 35       30 31 32 33
  116: *>               44 45       40 41 42 43 44
  117: *>                  55       50 51 52 53 54 55
  118: *>
  119: *>
  120: *>  Let TRANSR = 'N'. RFP holds AP as follows:
  121: *>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  122: *>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  123: *>  the transpose of the first three columns of AP upper.
  124: *>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  125: *>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  126: *>  the transpose of the last three columns of AP lower.
  127: *>  This covers the case N even and TRANSR = 'N'.
  128: *>
  129: *>         RFP A                   RFP A
  130: *>
  131: *>        03 04 05                33 43 53
  132: *>        13 14 15                00 44 54
  133: *>        23 24 25                10 11 55
  134: *>        33 34 35                20 21 22
  135: *>        00 44 45                30 31 32
  136: *>        01 11 55                40 41 42
  137: *>        02 12 22                50 51 52
  138: *>
  139: *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  140: *>  transpose of RFP A above. One therefore gets:
  141: *>
  142: *>
  143: *>           RFP A                   RFP A
  144: *>
  145: *>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
  146: *>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
  147: *>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
  148: *>
  149: *>
  150: *>  We then consider Rectangular Full Packed (RFP) Format when N is
  151: *>  odd. We give an example where N = 5.
  152: *>
  153: *>     AP is Upper                 AP is Lower
  154: *>
  155: *>   00 01 02 03 04              00
  156: *>      11 12 13 14              10 11
  157: *>         22 23 24              20 21 22
  158: *>            33 34              30 31 32 33
  159: *>               44              40 41 42 43 44
  160: *>
  161: *>
  162: *>  Let TRANSR = 'N'. RFP holds AP as follows:
  163: *>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  164: *>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  165: *>  the transpose of the first two columns of AP upper.
  166: *>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  167: *>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  168: *>  the transpose of the last two columns of AP lower.
  169: *>  This covers the case N odd and TRANSR = 'N'.
  170: *>
  171: *>         RFP A                   RFP A
  172: *>
  173: *>        02 03 04                00 33 43
  174: *>        12 13 14                10 11 44
  175: *>        22 23 24                20 21 22
  176: *>        00 33 34                30 31 32
  177: *>        01 11 44                40 41 42
  178: *>
  179: *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  180: *>  transpose of RFP A above. One therefore gets:
  181: *>
  182: *>           RFP A                   RFP A
  183: *>
  184: *>     02 12 22 00 01             00 10 20 30 40 50
  185: *>     03 13 23 33 11             33 11 21 31 41 51
  186: *>     04 14 24 34 44             43 44 22 32 42 52
  187: *> \endverbatim
  188: *>
  189: *  =====================================================================
  190:       SUBROUTINE DPFTRI( TRANSR, UPLO, N, A, INFO )
  191: *
  192: *  -- LAPACK computational routine --
  193: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  194: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  195: *
  196: *     .. Scalar Arguments ..
  197:       CHARACTER          TRANSR, UPLO
  198:       INTEGER            INFO, N
  199: *     .. Array Arguments ..
  200:       DOUBLE PRECISION         A( 0: * )
  201: *     ..
  202: *
  203: *  =====================================================================
  204: *
  205: *     .. Parameters ..
  206:       DOUBLE PRECISION   ONE
  207:       PARAMETER          ( ONE = 1.0D+0 )
  208: *     ..
  209: *     .. Local Scalars ..
  210:       LOGICAL            LOWER, NISODD, NORMALTRANSR
  211:       INTEGER            N1, N2, K
  212: *     ..
  213: *     .. External Functions ..
  214:       LOGICAL            LSAME
  215:       EXTERNAL           LSAME
  216: *     ..
  217: *     .. External Subroutines ..
  218:       EXTERNAL           XERBLA, DTFTRI, DLAUUM, DTRMM, DSYRK
  219: *     ..
  220: *     .. Intrinsic Functions ..
  221:       INTRINSIC          MOD
  222: *     ..
  223: *     .. Executable Statements ..
  224: *
  225: *     Test the input parameters.
  226: *
  227:       INFO = 0
  228:       NORMALTRANSR = LSAME( TRANSR, 'N' )
  229:       LOWER = LSAME( UPLO, 'L' )
  230:       IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
  231:          INFO = -1
  232:       ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
  233:          INFO = -2
  234:       ELSE IF( N.LT.0 ) THEN
  235:          INFO = -3
  236:       END IF
  237:       IF( INFO.NE.0 ) THEN
  238:          CALL XERBLA( 'DPFTRI', -INFO )
  239:          RETURN
  240:       END IF
  241: *
  242: *     Quick return if possible
  243: *
  244:       IF( N.EQ.0 )
  245:      $   RETURN
  246: *
  247: *     Invert the triangular Cholesky factor U or L.
  248: *
  249:       CALL DTFTRI( TRANSR, UPLO, 'N', N, A, INFO )
  250:       IF( INFO.GT.0 )
  251:      $   RETURN
  252: *
  253: *     If N is odd, set NISODD = .TRUE.
  254: *     If N is even, set K = N/2 and NISODD = .FALSE.
  255: *
  256:       IF( MOD( N, 2 ).EQ.0 ) THEN
  257:          K = N / 2
  258:          NISODD = .FALSE.
  259:       ELSE
  260:          NISODD = .TRUE.
  261:       END IF
  262: *
  263: *     Set N1 and N2 depending on LOWER
  264: *
  265:       IF( LOWER ) THEN
  266:          N2 = N / 2
  267:          N1 = N - N2
  268:       ELSE
  269:          N1 = N / 2
  270:          N2 = N - N1
  271:       END IF
  272: *
  273: *     Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
  274: *     inv(L)^C*inv(L). There are eight cases.
  275: *
  276:       IF( NISODD ) THEN
  277: *
  278: *        N is odd
  279: *
  280:          IF( NORMALTRANSR ) THEN
  281: *
  282: *           N is odd and TRANSR = 'N'
  283: *
  284:             IF( LOWER ) THEN
  285: *
  286: *              SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
  287: *              T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
  288: *              T1 -> a(0), T2 -> a(n), S -> a(N1)
  289: *
  290:                CALL DLAUUM( 'L', N1, A( 0 ), N, INFO )
  291:                CALL DSYRK( 'L', 'T', N1, N2, ONE, A( N1 ), N, ONE,
  292:      $                     A( 0 ), N )
  293:                CALL DTRMM( 'L', 'U', 'N', 'N', N2, N1, ONE, A( N ), N,
  294:      $                     A( N1 ), N )
  295:                CALL DLAUUM( 'U', N2, A( N ), N, INFO )
  296: *
  297:             ELSE
  298: *
  299: *              SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
  300: *              T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
  301: *              T1 -> a(N2), T2 -> a(N1), S -> a(0)
  302: *
  303:                CALL DLAUUM( 'L', N1, A( N2 ), N, INFO )
  304:                CALL DSYRK( 'L', 'N', N1, N2, ONE, A( 0 ), N, ONE,
  305:      $                     A( N2 ), N )
  306:                CALL DTRMM( 'R', 'U', 'T', 'N', N1, N2, ONE, A( N1 ), N,
  307:      $                     A( 0 ), N )
  308:                CALL DLAUUM( 'U', N2, A( N1 ), N, INFO )
  309: *
  310:             END IF
  311: *
  312:          ELSE
  313: *
  314: *           N is odd and TRANSR = 'T'
  315: *
  316:             IF( LOWER ) THEN
  317: *
  318: *              SRPA for LOWER, TRANSPOSE, and N is odd
  319: *              T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
  320: *
  321:                CALL DLAUUM( 'U', N1, A( 0 ), N1, INFO )
  322:                CALL DSYRK( 'U', 'N', N1, N2, ONE, A( N1*N1 ), N1, ONE,
  323:      $                     A( 0 ), N1 )
  324:                CALL DTRMM( 'R', 'L', 'N', 'N', N1, N2, ONE, A( 1 ), N1,
  325:      $                     A( N1*N1 ), N1 )
  326:                CALL DLAUUM( 'L', N2, A( 1 ), N1, INFO )
  327: *
  328:             ELSE
  329: *
  330: *              SRPA for UPPER, TRANSPOSE, and N is odd
  331: *              T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
  332: *
  333:                CALL DLAUUM( 'U', N1, A( N2*N2 ), N2, INFO )
  334:                CALL DSYRK( 'U', 'T', N1, N2, ONE, A( 0 ), N2, ONE,
  335:      $                     A( N2*N2 ), N2 )
  336:                CALL DTRMM( 'L', 'L', 'T', 'N', N2, N1, ONE, A( N1*N2 ),
  337:      $                     N2, A( 0 ), N2 )
  338:                CALL DLAUUM( 'L', N2, A( N1*N2 ), N2, INFO )
  339: *
  340:             END IF
  341: *
  342:          END IF
  343: *
  344:       ELSE
  345: *
  346: *        N is even
  347: *
  348:          IF( NORMALTRANSR ) THEN
  349: *
  350: *           N is even and TRANSR = 'N'
  351: *
  352:             IF( LOWER ) THEN
  353: *
  354: *              SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
  355: *              T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
  356: *              T1 -> a(1), T2 -> a(0), S -> a(k+1)
  357: *
  358:                CALL DLAUUM( 'L', K, A( 1 ), N+1, INFO )
  359:                CALL DSYRK( 'L', 'T', K, K, ONE, A( K+1 ), N+1, ONE,
  360:      $                     A( 1 ), N+1 )
  361:                CALL DTRMM( 'L', 'U', 'N', 'N', K, K, ONE, A( 0 ), N+1,
  362:      $                     A( K+1 ), N+1 )
  363:                CALL DLAUUM( 'U', K, A( 0 ), N+1, INFO )
  364: *
  365:             ELSE
  366: *
  367: *              SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
  368: *              T1 -> a(k+1,0) ,  T2 -> a(k,0),   S -> a(0,0)
  369: *              T1 -> a(k+1), T2 -> a(k), S -> a(0)
  370: *
  371:                CALL DLAUUM( 'L', K, A( K+1 ), N+1, INFO )
  372:                CALL DSYRK( 'L', 'N', K, K, ONE, A( 0 ), N+1, ONE,
  373:      $                     A( K+1 ), N+1 )
  374:                CALL DTRMM( 'R', 'U', 'T', 'N', K, K, ONE, A( K ), N+1,
  375:      $                     A( 0 ), N+1 )
  376:                CALL DLAUUM( 'U', K, A( K ), N+1, INFO )
  377: *
  378:             END IF
  379: *
  380:          ELSE
  381: *
  382: *           N is even and TRANSR = 'T'
  383: *
  384:             IF( LOWER ) THEN
  385: *
  386: *              SRPA for LOWER, TRANSPOSE, and N is even (see paper)
  387: *              T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
  388: *              T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
  389: *
  390:                CALL DLAUUM( 'U', K, A( K ), K, INFO )
  391:                CALL DSYRK( 'U', 'N', K, K, ONE, A( K*( K+1 ) ), K, ONE,
  392:      $                     A( K ), K )
  393:                CALL DTRMM( 'R', 'L', 'N', 'N', K, K, ONE, A( 0 ), K,
  394:      $                     A( K*( K+1 ) ), K )
  395:                CALL DLAUUM( 'L', K, A( 0 ), K, INFO )
  396: *
  397:             ELSE
  398: *
  399: *              SRPA for UPPER, TRANSPOSE, and N is even (see paper)
  400: *              T1 -> B(0,k+1),     T2 -> B(0,k),   S -> B(0,0),
  401: *              T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
  402: *
  403:                CALL DLAUUM( 'U', K, A( K*( K+1 ) ), K, INFO )
  404:                CALL DSYRK( 'U', 'T', K, K, ONE, A( 0 ), K, ONE,
  405:      $                     A( K*( K+1 ) ), K )
  406:                CALL DTRMM( 'L', 'L', 'T', 'N', K, K, ONE, A( K*K ), K,
  407:      $                     A( 0 ), K )
  408:                CALL DLAUUM( 'L', K, A( K*K ), K, INFO )
  409: *
  410:             END IF
  411: *
  412:          END IF
  413: *
  414:       END IF
  415: *
  416:       RETURN
  417: *
  418: *     End of DPFTRI
  419: *
  420:       END

CVSweb interface <joel.bertrand@systella.fr>