File:  [local] / rpl / lapack / lapack / zpotf2.f
Revision 1.19: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:34 2023 UTC (9 months 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 ZPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm).
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZPOTF2 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpotf2.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpotf2.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpotf2.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZPOTF2( UPLO, N, A, LDA, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       CHARACTER          UPLO
   25: *       INTEGER            INFO, LDA, N
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       COMPLEX*16         A( LDA, * )
   29: *       ..
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> ZPOTF2 computes the Cholesky factorization of a complex Hermitian
   38: *> positive definite matrix A.
   39: *>
   40: *> The factorization has the form
   41: *>    A = U**H * U ,  if UPLO = 'U', or
   42: *>    A = L  * L**H,  if UPLO = 'L',
   43: *> where U is an upper triangular matrix and L is lower triangular.
   44: *>
   45: *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
   46: *> \endverbatim
   47: *
   48: *  Arguments:
   49: *  ==========
   50: *
   51: *> \param[in] UPLO
   52: *> \verbatim
   53: *>          UPLO is CHARACTER*1
   54: *>          Specifies whether the upper or lower triangular part of the
   55: *>          Hermitian matrix A is stored.
   56: *>          = 'U':  Upper triangular
   57: *>          = 'L':  Lower triangular
   58: *> \endverbatim
   59: *>
   60: *> \param[in] N
   61: *> \verbatim
   62: *>          N is INTEGER
   63: *>          The order of the matrix A.  N >= 0.
   64: *> \endverbatim
   65: *>
   66: *> \param[in,out] A
   67: *> \verbatim
   68: *>          A is COMPLEX*16 array, dimension (LDA,N)
   69: *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
   70: *>          n by n upper triangular part of A contains the upper
   71: *>          triangular part of the matrix A, and the strictly lower
   72: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   73: *>          leading n by n lower triangular part of A contains the lower
   74: *>          triangular part of the matrix A, and the strictly upper
   75: *>          triangular part of A is not referenced.
   76: *>
   77: *>          On exit, if INFO = 0, the factor U or L from the Cholesky
   78: *>          factorization A = U**H *U  or A = L*L**H.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] LDA
   82: *> \verbatim
   83: *>          LDA is INTEGER
   84: *>          The leading dimension of the array A.  LDA >= max(1,N).
   85: *> \endverbatim
   86: *>
   87: *> \param[out] INFO
   88: *> \verbatim
   89: *>          INFO is INTEGER
   90: *>          = 0: successful exit
   91: *>          < 0: if INFO = -k, the k-th argument had an illegal value
   92: *>          > 0: if INFO = k, the leading minor of order k is not
   93: *>               positive definite, and the factorization could not be
   94: *>               completed.
   95: *> \endverbatim
   96: *
   97: *  Authors:
   98: *  ========
   99: *
  100: *> \author Univ. of Tennessee
  101: *> \author Univ. of California Berkeley
  102: *> \author Univ. of Colorado Denver
  103: *> \author NAG Ltd.
  104: *
  105: *> \ingroup complex16POcomputational
  106: *
  107: *  =====================================================================
  108:       SUBROUTINE ZPOTF2( UPLO, N, A, LDA, INFO )
  109: *
  110: *  -- LAPACK computational routine --
  111: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  112: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  113: *
  114: *     .. Scalar Arguments ..
  115:       CHARACTER          UPLO
  116:       INTEGER            INFO, LDA, N
  117: *     ..
  118: *     .. Array Arguments ..
  119:       COMPLEX*16         A( LDA, * )
  120: *     ..
  121: *
  122: *  =====================================================================
  123: *
  124: *     .. Parameters ..
  125:       DOUBLE PRECISION   ONE, ZERO
  126:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  127:       COMPLEX*16         CONE
  128:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
  129: *     ..
  130: *     .. Local Scalars ..
  131:       LOGICAL            UPPER
  132:       INTEGER            J
  133:       DOUBLE PRECISION   AJJ
  134: *     ..
  135: *     .. External Functions ..
  136:       LOGICAL            LSAME, DISNAN
  137:       COMPLEX*16         ZDOTC
  138:       EXTERNAL           LSAME, ZDOTC, DISNAN
  139: *     ..
  140: *     .. External Subroutines ..
  141:       EXTERNAL           XERBLA, ZDSCAL, ZGEMV, ZLACGV
  142: *     ..
  143: *     .. Intrinsic Functions ..
  144:       INTRINSIC          DBLE, MAX, SQRT
  145: *     ..
  146: *     .. Executable Statements ..
  147: *
  148: *     Test the input parameters.
  149: *
  150:       INFO = 0
  151:       UPPER = LSAME( UPLO, 'U' )
  152:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  153:          INFO = -1
  154:       ELSE IF( N.LT.0 ) THEN
  155:          INFO = -2
  156:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  157:          INFO = -4
  158:       END IF
  159:       IF( INFO.NE.0 ) THEN
  160:          CALL XERBLA( 'ZPOTF2', -INFO )
  161:          RETURN
  162:       END IF
  163: *
  164: *     Quick return if possible
  165: *
  166:       IF( N.EQ.0 )
  167:      $   RETURN
  168: *
  169:       IF( UPPER ) THEN
  170: *
  171: *        Compute the Cholesky factorization A = U**H *U.
  172: *
  173:          DO 10 J = 1, N
  174: *
  175: *           Compute U(J,J) and test for non-positive-definiteness.
  176: *
  177:             AJJ = DBLE( A( J, J ) ) - DBLE( ZDOTC( J-1, A( 1, J ), 1,
  178:      $            A( 1, J ), 1 ) )
  179:             IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
  180:                A( J, J ) = AJJ
  181:                GO TO 30
  182:             END IF
  183:             AJJ = SQRT( AJJ )
  184:             A( J, J ) = AJJ
  185: *
  186: *           Compute elements J+1:N of row J.
  187: *
  188:             IF( J.LT.N ) THEN
  189:                CALL ZLACGV( J-1, A( 1, J ), 1 )
  190:                CALL ZGEMV( 'Transpose', J-1, N-J, -CONE, A( 1, J+1 ),
  191:      $                     LDA, A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
  192:                CALL ZLACGV( J-1, A( 1, J ), 1 )
  193:                CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
  194:             END IF
  195:    10    CONTINUE
  196:       ELSE
  197: *
  198: *        Compute the Cholesky factorization A = L*L**H.
  199: *
  200:          DO 20 J = 1, N
  201: *
  202: *           Compute L(J,J) and test for non-positive-definiteness.
  203: *
  204:             AJJ = DBLE( A( J, J ) ) - DBLE( ZDOTC( J-1, A( J, 1 ), LDA,
  205:      $            A( J, 1 ), LDA ) )
  206:             IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
  207:                A( J, J ) = AJJ
  208:                GO TO 30
  209:             END IF
  210:             AJJ = SQRT( AJJ )
  211:             A( J, J ) = AJJ
  212: *
  213: *           Compute elements J+1:N of column J.
  214: *
  215:             IF( J.LT.N ) THEN
  216:                CALL ZLACGV( J-1, A( J, 1 ), LDA )
  217:                CALL ZGEMV( 'No transpose', N-J, J-1, -CONE, A( J+1, 1 ),
  218:      $                     LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
  219:                CALL ZLACGV( J-1, A( J, 1 ), LDA )
  220:                CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
  221:             END IF
  222:    20    CONTINUE
  223:       END IF
  224:       GO TO 40
  225: *
  226:    30 CONTINUE
  227:       INFO = J
  228: *
  229:    40 CONTINUE
  230:       RETURN
  231: *
  232: *     End of ZPOTF2
  233: *
  234:       END

CVSweb interface <joel.bertrand@systella.fr>