Return to dpotrf2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack (3.6.0) et du numéro de version du RPL/2.

    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: *> \date November 2015
  103: *
  104: *> \ingroup doublePOcomputational
  105: *
  106: *  =====================================================================
  107:       RECURSIVE SUBROUTINE DPOTRF2( UPLO, N, A, LDA, INFO )
  108: *
  109: *  -- LAPACK computational routine (version 3.6.0) --
  110: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  111: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  112: *     November 2015
  113: *
  114: *     .. Scalar Arguments ..
  115:       CHARACTER          UPLO
  116:       INTEGER            INFO, LDA, N
  117: *     ..
  118: *     .. Array Arguments ..
  119:       DOUBLE PRECISION   A( LDA, * )
  120: *     ..
  121: *
  122: *  =====================================================================
  123: *
  124: *     .. Parameters ..
  125:       DOUBLE PRECISION   ONE, ZERO
  126:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  127: *     ..
  128: *     .. Local Scalars ..
  129:       LOGICAL            UPPER            
  130:       INTEGER            N1, N2, IINFO
  131: *     ..
  132: *     .. External Functions ..
  133:       LOGICAL            LSAME, DISNAN
  134:       EXTERNAL           LSAME, DISNAN
  135: *     ..
  136: *     .. External Subroutines ..
  137:       EXTERNAL           DSYRK, DTRSM, XERBLA
  138: *     ..
  139: *     .. Intrinsic Functions ..
  140:       INTRINSIC          MAX, SQRT
  141: *     ..
  142: *     .. Executable Statements ..
  143: *
  144: *     Test the input parameters
  145: *
  146:       INFO = 0
  147:       UPPER = LSAME( UPLO, 'U' )
  148:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  149:          INFO = -1
  150:       ELSE IF( N.LT.0 ) THEN
  151:          INFO = -2
  152:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  153:          INFO = -4
  154:       END IF
  155:       IF( INFO.NE.0 ) THEN
  156:          CALL XERBLA( 'DPOTRF2', -INFO )
  157:          RETURN
  158:       END IF
  159: *
  160: *     Quick return if possible
  161: *
  162:       IF( N.EQ.0 )
  163:      $   RETURN
  164: *
  165: *     N=1 case
  166: *
  167:       IF( N.EQ.1 ) THEN
  168: *
  169: *        Test for non-positive-definiteness
  170: *
  171:          IF( A( 1, 1 ).LE.ZERO.OR.DISNAN( A( 1, 1 ) ) ) THEN
  172:             INFO = 1
  173:             RETURN
  174:          END IF
  175: *
  176: *        Factor
  177: *
  178:          A( 1, 1 ) = SQRT( A( 1, 1 ) )
  179: *
  180: *     Use recursive code
  181: *
  182:       ELSE
  183:          N1 = N/2
  184:          N2 = N-N1
  185: *
  186: *        Factor A11
  187: *
  188:          CALL DPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO )
  189:          IF ( IINFO.NE.0 ) THEN
  190:             INFO = IINFO
  191:             RETURN
  192:          END IF    
  193: *
  194: *        Compute the Cholesky factorization A = U**T*U
  195: *
  196:          IF( UPPER ) THEN
  197: *
  198: *           Update and scale A12
  199: *
  200:             CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE,
  201:      $                  A( 1, 1 ), LDA, A( 1, N1+1 ), LDA )            
  202: *
  203: *           Update and factor A22
  204: *          
  205:             CALL DSYRK( UPLO, 'T', N2, N1, -ONE, A( 1, N1+1 ), LDA,
  206:      $                  ONE, A( N1+1, N1+1 ), LDA )
  207:             CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
  208:             IF ( IINFO.NE.0 ) THEN
  209:                INFO = IINFO + N1
  210:                RETURN
  211:             END IF
  212: *
  213: *        Compute the Cholesky factorization A = L*L**T
  214: *
  215:          ELSE
  216: *
  217: *           Update and scale A21
  218: *
  219:             CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE, 
  220:      $                  A( 1, 1 ), LDA, A( N1+1, 1 ), LDA )
  221: *
  222: *           Update and factor A22
  223: *
  224:             CALL DSYRK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA,
  225:      $                  ONE, A( N1+1, N1+1 ), LDA )
  226:             CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
  227:             IF ( IINFO.NE.0 ) THEN
  228:                INFO = IINFO + N1
  229:                RETURN
  230:             END IF
  231:          END IF
  232:       END IF
  233:       RETURN
  234: *
  235: *     End of DPOTRF2
  236: *
  237:       END