Return to dpotrf.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    1:       SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       CHARACTER          UPLO
   10:       INTEGER            INFO, LDA, N
   11: *     ..
   12: *     .. Array Arguments ..
   13:       DOUBLE PRECISION   A( LDA, * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  DPOTRF computes the Cholesky factorization of a real symmetric
   20: *  positive definite matrix A.
   21: *
   22: *  The factorization has the form
   23: *     A = U**T * U,  if UPLO = 'U', or
   24: *     A = L  * L**T,  if UPLO = 'L',
   25: *  where U is an upper triangular matrix and L is lower triangular.
   26: *
   27: *  This is the block version of the algorithm, calling Level 3 BLAS.
   28: *
   29: *  Arguments
   30: *  =========
   31: *
   32: *  UPLO    (input) CHARACTER*1
   33: *          = 'U':  Upper triangle of A is stored;
   34: *          = 'L':  Lower triangle of A is stored.
   35: *
   36: *  N       (input) INTEGER
   37: *          The order of the matrix A.  N >= 0.
   38: *
   39: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
   40: *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
   41: *          N-by-N upper triangular part of A contains the upper
   42: *          triangular part of the matrix A, and the strictly lower
   43: *          triangular part of A is not referenced.  If UPLO = 'L', the
   44: *          leading N-by-N lower triangular part of A contains the lower
   45: *          triangular part of the matrix A, and the strictly upper
   46: *          triangular part of A is not referenced.
   47: *
   48: *          On exit, if INFO = 0, the factor U or L from the Cholesky
   49: *          factorization A = U**T*U or A = L*L**T.
   50: *
   51: *  LDA     (input) INTEGER
   52: *          The leading dimension of the array A.  LDA >= max(1,N).
   53: *
   54: *  INFO    (output) INTEGER
   55: *          = 0:  successful exit
   56: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   57: *          > 0:  if INFO = i, the leading minor of order i is not
   58: *                positive definite, and the factorization could not be
   59: *                completed.
   60: *
   61: *  =====================================================================
   62: *
   63: *     .. Parameters ..
   64:       DOUBLE PRECISION   ONE
   65:       PARAMETER          ( ONE = 1.0D+0 )
   66: *     ..
   67: *     .. Local Scalars ..
   68:       LOGICAL            UPPER
   69:       INTEGER            J, JB, NB
   70: *     ..
   71: *     .. External Functions ..
   72:       LOGICAL            LSAME
   73:       INTEGER            ILAENV
   74:       EXTERNAL           LSAME, ILAENV
   75: *     ..
   76: *     .. External Subroutines ..
   77:       EXTERNAL           DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA
   78: *     ..
   79: *     .. Intrinsic Functions ..
   80:       INTRINSIC          MAX, MIN
   81: *     ..
   82: *     .. Executable Statements ..
   83: *
   84: *     Test the input parameters.
   85: *
   86:       INFO = 0
   87:       UPPER = LSAME( UPLO, 'U' )
   88:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
   89:          INFO = -1
   90:       ELSE IF( N.LT.0 ) THEN
   91:          INFO = -2
   92:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
   93:          INFO = -4
   94:       END IF
   95:       IF( INFO.NE.0 ) THEN
   96:          CALL XERBLA( 'DPOTRF', -INFO )
   97:          RETURN
   98:       END IF
   99: *
  100: *     Quick return if possible
  101: *
  102:       IF( N.EQ.0 )
  103:      $   RETURN
  104: *
  105: *     Determine the block size for this environment.
  106: *
  107:       NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 )
  108:       IF( NB.LE.1 .OR. NB.GE.N ) THEN
  109: *
  110: *        Use unblocked code.
  111: *
  112:          CALL DPOTF2( UPLO, N, A, LDA, INFO )
  113:       ELSE
  114: *
  115: *        Use blocked code.
  116: *
  117:          IF( UPPER ) THEN
  118: *
  119: *           Compute the Cholesky factorization A = U'*U.
  120: *
  121:             DO 10 J = 1, N, NB
  122: *
  123: *              Update and factorize the current diagonal block and test
  124: *              for non-positive-definiteness.
  125: *
  126:                JB = MIN( NB, N-J+1 )
  127:                CALL DSYRK( 'Upper', 'Transpose', JB, J-1, -ONE,
  128:      $                     A( 1, J ), LDA, ONE, A( J, J ), LDA )
  129:                CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
  130:                IF( INFO.NE.0 )
  131:      $            GO TO 30
  132:                IF( J+JB.LE.N ) THEN
  133: *
  134: *                 Compute the current block row.
  135: *
  136:                   CALL DGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1,
  137:      $                        J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ),
  138:      $                        LDA, ONE, A( J, J+JB ), LDA )
  139:                   CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit',
  140:      $                        JB, N-J-JB+1, ONE, A( J, J ), LDA,
  141:      $                        A( J, J+JB ), LDA )
  142:                END IF
  143:    10       CONTINUE
  144: *
  145:          ELSE
  146: *
  147: *           Compute the Cholesky factorization A = L*L'.
  148: *
  149:             DO 20 J = 1, N, NB
  150: *
  151: *              Update and factorize the current diagonal block and test
  152: *              for non-positive-definiteness.
  153: *
  154:                JB = MIN( NB, N-J+1 )
  155:                CALL DSYRK( 'Lower', 'No transpose', JB, J-1, -ONE,
  156:      $                     A( J, 1 ), LDA, ONE, A( J, J ), LDA )
  157:                CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
  158:                IF( INFO.NE.0 )
  159:      $            GO TO 30
  160:                IF( J+JB.LE.N ) THEN
  161: *
  162: *                 Compute the current block column.
  163: *
  164:                   CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
  165:      $                        J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ),
  166:      $                        LDA, ONE, A( J+JB, J ), LDA )
  167:                   CALL DTRSM( 'Right', 'Lower', 'Transpose', 'Non-unit',
  168:      $                        N-J-JB+1, JB, ONE, A( J, J ), LDA,
  169:      $                        A( J+JB, J ), LDA )
  170:                END IF
  171:    20       CONTINUE
  172:          END IF
  173:       END IF
  174:       GO TO 40
  175: *
  176:    30 CONTINUE
  177:       INFO = INFO + J - 1
  178: *
  179:    40 CONTINUE
  180:       RETURN
  181: *
  182: *     End of DPOTRF
  183: *
  184:       END