Return to dpotrf2.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief \b DPOTRF2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *  Definition:
    9: *  ===========
   10: *
   11: *       RECURSIVE SUBROUTINE DPOTRF2( UPLO, N, A, LDA, INFO )
   12: *
   13: *       .. Scalar Arguments ..
   14: *       CHARACTER          UPLO
   15: *       INTEGER            INFO, LDA, N
   16: *       ..
   17: *       .. Array Arguments ..
   18: *       REAL               A( LDA, * )
   19: *       ..
   20: *
   21: *
   22: *> \par Purpose:
   23: *  =============
   24: *>
   25: *> \verbatim
   26: *>
   27: *> DPOTRF2 computes the Cholesky factorization of a real symmetric
   28: *> positive definite matrix A using the recursive algorithm.
   29: *>
   30: *> The factorization has the form
   31: *>    A = U**T * U,  if UPLO = 'U', or
   32: *>    A = L  * L**T,  if UPLO = 'L',
   33: *> where U is an upper triangular matrix and L is lower triangular.
   34: *>
   35: *> This is the recursive version of the algorithm. It divides
   36: *> the matrix into four submatrices:
   37: *>
   38: *>        [  A11 | A12  ]  where A11 is n1 by n1 and A22 is n2 by n2
   39: *>    A = [ -----|----- ]  with n1 = n/2
   40: *>        [  A21 | A22  ]       n2 = n-n1
   41: *>
   42: *> The subroutine calls itself to factor A11. Update and scale A21
   43: *> or A12, update A22 then calls itself to factor A22.
   44: *>
   45: *> \endverbatim
   46: *
   47: *  Arguments:
   48: *  ==========
   49: *
   50: *> \param[in] UPLO
   51: *> \verbatim
   52: *>          UPLO is CHARACTER*1
   53: *>          = 'U':  Upper triangle of A is stored;
   54: *>          = 'L':  Lower triangle of A is stored.
   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,out] A
   64: *> \verbatim
   65: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   66: *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
   67: *>          N-by-N upper triangular part of A contains the upper
   68: *>          triangular part of the matrix A, and the strictly lower
   69: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   70: *>          leading N-by-N lower triangular part of A contains the lower
   71: *>          triangular part of the matrix A, and the strictly upper
   72: *>          triangular part of A is not referenced.
   73: *>
   74: *>          On exit, if INFO = 0, the factor U or L from the Cholesky
   75: *>          factorization A = U**T*U or A = L*L**T.
   76: *> \endverbatim
   77: *>
   78: *> \param[in] LDA
   79: *> \verbatim
   80: *>          LDA is INTEGER
   81: *>          The leading dimension of the array A.  LDA >= max(1,N).
   82: *> \endverbatim
   83: *>
   84: *> \param[out] INFO
   85: *> \verbatim
   86: *>          INFO is INTEGER
   87: *>          = 0:  successful exit
   88: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
   89: *>          > 0:  if INFO = i, the leading minor of order i is not
   90: *>                positive definite, and the factorization could not be
   91: *>                completed.
   92: *> \endverbatim
   93: *
   94: *  Authors:
   95: *  ========
   96: *
   97: *> \author Univ. of Tennessee
   98: *> \author Univ. of California Berkeley
   99: *> \author Univ. of Colorado Denver
  100: *> \author NAG Ltd.
  101: *
  102: *> \ingroup doublePOcomputational
  103: *
  104: *  =====================================================================
  105:       RECURSIVE SUBROUTINE DPOTRF2( UPLO, N, A, LDA, INFO )
  106: *
  107: *  -- LAPACK computational routine --
  108: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  109: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  110: *
  111: *     .. Scalar Arguments ..
  112:       CHARACTER          UPLO
  113:       INTEGER            INFO, LDA, N
  114: *     ..
  115: *     .. Array Arguments ..
  116:       DOUBLE PRECISION   A( LDA, * )
  117: *     ..
  118: *
  119: *  =====================================================================
  120: *
  121: *     .. Parameters ..
  122:       DOUBLE PRECISION   ONE, ZERO
  123:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  124: *     ..
  125: *     .. Local Scalars ..
  126:       LOGICAL            UPPER
  127:       INTEGER            N1, N2, IINFO
  128: *     ..
  129: *     .. External Functions ..
  130:       LOGICAL            LSAME, DISNAN
  131:       EXTERNAL           LSAME, DISNAN
  132: *     ..
  133: *     .. External Subroutines ..
  134:       EXTERNAL           DSYRK, DTRSM, XERBLA
  135: *     ..
  136: *     .. Intrinsic Functions ..
  137:       INTRINSIC          MAX, SQRT
  138: *     ..
  139: *     .. Executable Statements ..
  140: *
  141: *     Test the input parameters
  142: *
  143:       INFO = 0
  144:       UPPER = LSAME( UPLO, 'U' )
  145:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  146:          INFO = -1
  147:       ELSE IF( N.LT.0 ) THEN
  148:          INFO = -2
  149:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  150:          INFO = -4
  151:       END IF
  152:       IF( INFO.NE.0 ) THEN
  153:          CALL XERBLA( 'DPOTRF2', -INFO )
  154:          RETURN
  155:       END IF
  156: *
  157: *     Quick return if possible
  158: *
  159:       IF( N.EQ.0 )
  160:      $   RETURN
  161: *
  162: *     N=1 case
  163: *
  164:       IF( N.EQ.1 ) THEN
  165: *
  166: *        Test for non-positive-definiteness
  167: *
  168:          IF( A( 1, 1 ).LE.ZERO.OR.DISNAN( A( 1, 1 ) ) ) THEN
  169:             INFO = 1
  170:             RETURN
  171:          END IF
  172: *
  173: *        Factor
  174: *
  175:          A( 1, 1 ) = SQRT( A( 1, 1 ) )
  176: *
  177: *     Use recursive code
  178: *
  179:       ELSE
  180:          N1 = N/2
  181:          N2 = N-N1
  182: *
  183: *        Factor A11
  184: *
  185:          CALL DPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO )
  186:          IF ( IINFO.NE.0 ) THEN
  187:             INFO = IINFO
  188:             RETURN
  189:          END IF
  190: *
  191: *        Compute the Cholesky factorization A = U**T*U
  192: *
  193:          IF( UPPER ) THEN
  194: *
  195: *           Update and scale A12
  196: *
  197:             CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE,
  198:      $                  A( 1, 1 ), LDA, A( 1, N1+1 ), LDA )
  199: *
  200: *           Update and factor A22
  201: *
  202:             CALL DSYRK( UPLO, 'T', N2, N1, -ONE, A( 1, N1+1 ), LDA,
  203:      $                  ONE, A( N1+1, N1+1 ), LDA )
  204:             CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
  205:             IF ( IINFO.NE.0 ) THEN
  206:                INFO = IINFO + N1
  207:                RETURN
  208:             END IF
  209: *
  210: *        Compute the Cholesky factorization A = L*L**T
  211: *
  212:          ELSE
  213: *
  214: *           Update and scale A21
  215: *
  216:             CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE,
  217:      $                  A( 1, 1 ), LDA, A( N1+1, 1 ), LDA )
  218: *
  219: *           Update and factor A22
  220: *
  221:             CALL DSYRK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA,
  222:      $                  ONE, A( N1+1, N1+1 ), LDA )
  223:             CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
  224:             IF ( IINFO.NE.0 ) THEN
  225:                INFO = IINFO + N1
  226:                RETURN
  227:             END IF
  228:          END IF
  229:       END IF
  230:       RETURN
  231: *
  232: *     End of DPOTRF2
  233: *
  234:       END