Return to zpotrf.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZPOTRF( 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:       COMPLEX*16         A( LDA, * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  ZPOTRF computes the Cholesky factorization of a complex Hermitian
   20: *  positive definite matrix A.
   21: *
   22: *  The factorization has the form
   23: *     A = U**H * U,  if UPLO = 'U', or
   24: *     A = L  * L**H,  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) COMPLEX*16 array, dimension (LDA,N)
   40: *          On entry, the Hermitian 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**H*U or A = L*L**H.
   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:       COMPLEX*16         CONE
   66:       PARAMETER          ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ) )
   67: *     ..
   68: *     .. Local Scalars ..
   69:       LOGICAL            UPPER
   70:       INTEGER            J, JB, NB
   71: *     ..
   72: *     .. External Functions ..
   73:       LOGICAL            LSAME
   74:       INTEGER            ILAENV
   75:       EXTERNAL           LSAME, ILAENV
   76: *     ..
   77: *     .. External Subroutines ..
   78:       EXTERNAL           XERBLA, ZGEMM, ZHERK, ZPOTF2, ZTRSM
   79: *     ..
   80: *     .. Intrinsic Functions ..
   81:       INTRINSIC          MAX, MIN
   82: *     ..
   83: *     .. Executable Statements ..
   84: *
   85: *     Test the input parameters.
   86: *
   87:       INFO = 0
   88:       UPPER = LSAME( UPLO, 'U' )
   89:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
   90:          INFO = -1
   91:       ELSE IF( N.LT.0 ) THEN
   92:          INFO = -2
   93:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
   94:          INFO = -4
   95:       END IF
   96:       IF( INFO.NE.0 ) THEN
   97:          CALL XERBLA( 'ZPOTRF', -INFO )
   98:          RETURN
   99:       END IF
  100: *
  101: *     Quick return if possible
  102: *
  103:       IF( N.EQ.0 )
  104:      $   RETURN
  105: *
  106: *     Determine the block size for this environment.
  107: *
  108:       NB = ILAENV( 1, 'ZPOTRF', UPLO, N, -1, -1, -1 )
  109:       IF( NB.LE.1 .OR. NB.GE.N ) THEN
  110: *
  111: *        Use unblocked code.
  112: *
  113:          CALL ZPOTF2( UPLO, N, A, LDA, INFO )
  114:       ELSE
  115: *
  116: *        Use blocked code.
  117: *
  118:          IF( UPPER ) THEN
  119: *
  120: *           Compute the Cholesky factorization A = U'*U.
  121: *
  122:             DO 10 J = 1, N, NB
  123: *
  124: *              Update and factorize the current diagonal block and test
  125: *              for non-positive-definiteness.
  126: *
  127:                JB = MIN( NB, N-J+1 )
  128:                CALL ZHERK( 'Upper', 'Conjugate transpose', JB, J-1,
  129:      $                     -ONE, A( 1, J ), LDA, ONE, A( J, J ), LDA )
  130:                CALL ZPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
  131:                IF( INFO.NE.0 )
  132:      $            GO TO 30
  133:                IF( J+JB.LE.N ) THEN
  134: *
  135: *                 Compute the current block row.
  136: *
  137:                   CALL ZGEMM( 'Conjugate transpose', 'No transpose', JB,
  138:      $                        N-J-JB+1, J-1, -CONE, A( 1, J ), LDA,
  139:      $                        A( 1, J+JB ), LDA, CONE, A( J, J+JB ),
  140:      $                        LDA )
  141:                   CALL ZTRSM( 'Left', 'Upper', 'Conjugate transpose',
  142:      $                        'Non-unit', JB, N-J-JB+1, CONE, A( J, J ),
  143:      $                        LDA, A( J, J+JB ), LDA )
  144:                END IF
  145:    10       CONTINUE
  146: *
  147:          ELSE
  148: *
  149: *           Compute the Cholesky factorization A = L*L'.
  150: *
  151:             DO 20 J = 1, N, NB
  152: *
  153: *              Update and factorize the current diagonal block and test
  154: *              for non-positive-definiteness.
  155: *
  156:                JB = MIN( NB, N-J+1 )
  157:                CALL ZHERK( 'Lower', 'No transpose', JB, J-1, -ONE,
  158:      $                     A( J, 1 ), LDA, ONE, A( J, J ), LDA )
  159:                CALL ZPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
  160:                IF( INFO.NE.0 )
  161:      $            GO TO 30
  162:                IF( J+JB.LE.N ) THEN
  163: *
  164: *                 Compute the current block column.
  165: *
  166:                   CALL ZGEMM( 'No transpose', 'Conjugate transpose',
  167:      $                        N-J-JB+1, JB, J-1, -CONE, A( J+JB, 1 ),
  168:      $                        LDA, A( J, 1 ), LDA, CONE, A( J+JB, J ),
  169:      $                        LDA )
  170:                   CALL ZTRSM( 'Right', 'Lower', 'Conjugate transpose',
  171:      $                        'Non-unit', N-J-JB+1, JB, CONE, A( J, J ),
  172:      $                        LDA, A( J+JB, J ), LDA )
  173:                END IF
  174:    20       CONTINUE
  175:          END IF
  176:       END IF
  177:       GO TO 40
  178: *
  179:    30 CONTINUE
  180:       INFO = INFO + J - 1
  181: *
  182:    40 CONTINUE
  183:       RETURN
  184: *
  185: *     End of ZPOTRF
  186: *
  187:       END