Return to zpstf2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief \b ZPSTF2 computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZPSTF2 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpstf2.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpstf2.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpstf2.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       DOUBLE PRECISION   TOL
   25: *       INTEGER            INFO, LDA, N, RANK
   26: *       CHARACTER          UPLO
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       COMPLEX*16         A( LDA, * )
   30: *       DOUBLE PRECISION   WORK( 2*N )
   31: *       INTEGER            PIV( N )
   32: *       ..
   33: *
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *>
   40: *> ZPSTF2 computes the Cholesky factorization with complete
   41: *> pivoting of a complex Hermitian positive semidefinite matrix A.
   42: *>
   43: *> The factorization has the form
   44: *>    P**T * A * P = U**H * U ,  if UPLO = 'U',
   45: *>    P**T * A * P = L  * L**H,  if UPLO = 'L',
   46: *> where U is an upper triangular matrix and L is lower triangular, and
   47: *> P is stored as vector PIV.
   48: *>
   49: *> This algorithm does not attempt to check that A is positive
   50: *> semidefinite. This version of the algorithm calls level 2 BLAS.
   51: *> \endverbatim
   52: *
   53: *  Arguments:
   54: *  ==========
   55: *
   56: *> \param[in] UPLO
   57: *> \verbatim
   58: *>          UPLO is CHARACTER*1
   59: *>          Specifies whether the upper or lower triangular part of the
   60: *>          symmetric matrix A is stored.
   61: *>          = 'U':  Upper triangular
   62: *>          = 'L':  Lower triangular
   63: *> \endverbatim
   64: *>
   65: *> \param[in] N
   66: *> \verbatim
   67: *>          N is INTEGER
   68: *>          The order of the matrix A.  N >= 0.
   69: *> \endverbatim
   70: *>
   71: *> \param[in,out] A
   72: *> \verbatim
   73: *>          A is COMPLEX*16 array, dimension (LDA,N)
   74: *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
   75: *>          n by n upper triangular part of A contains the upper
   76: *>          triangular part of the matrix A, and the strictly lower
   77: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   78: *>          leading n by n lower triangular part of A contains the lower
   79: *>          triangular part of the matrix A, and the strictly upper
   80: *>          triangular part of A is not referenced.
   81: *>
   82: *>          On exit, if INFO = 0, the factor U or L from the Cholesky
   83: *>          factorization as above.
   84: *> \endverbatim
   85: *>
   86: *> \param[out] PIV
   87: *> \verbatim
   88: *>          PIV is INTEGER array, dimension (N)
   89: *>          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
   90: *> \endverbatim
   91: *>
   92: *> \param[out] RANK
   93: *> \verbatim
   94: *>          RANK is INTEGER
   95: *>          The rank of A given by the number of steps the algorithm
   96: *>          completed.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] TOL
  100: *> \verbatim
  101: *>          TOL is DOUBLE PRECISION
  102: *>          User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
  103: *>          will be used. The algorithm terminates at the (K-1)st step
  104: *>          if the pivot <= TOL.
  105: *> \endverbatim
  106: *>
  107: *> \param[in] LDA
  108: *> \verbatim
  109: *>          LDA is INTEGER
  110: *>          The leading dimension of the array A.  LDA >= max(1,N).
  111: *> \endverbatim
  112: *>
  113: *> \param[out] WORK
  114: *> \verbatim
  115: *>          WORK is DOUBLE PRECISION array, dimension (2*N)
  116: *>          Work space.
  117: *> \endverbatim
  118: *>
  119: *> \param[out] INFO
  120: *> \verbatim
  121: *>          INFO is INTEGER
  122: *>          < 0: If INFO = -K, the K-th argument had an illegal value,
  123: *>          = 0: algorithm completed successfully, and
  124: *>          > 0: the matrix A is either rank deficient with computed rank
  125: *>               as returned in RANK, or is not positive semidefinite. See
  126: *>               Section 7 of LAPACK Working Note #161 for further
  127: *>               information.
  128: *> \endverbatim
  129: *
  130: *  Authors:
  131: *  ========
  132: *
  133: *> \author Univ. of Tennessee
  134: *> \author Univ. of California Berkeley
  135: *> \author Univ. of Colorado Denver
  136: *> \author NAG Ltd.
  137: *
  138: *> \date December 2016
  139: *
  140: *> \ingroup complex16OTHERcomputational
  141: *
  142: *  =====================================================================
  143:       SUBROUTINE ZPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
  144: *
  145: *  -- LAPACK computational routine (version 3.7.0) --
  146: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  147: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  148: *     December 2016
  149: *
  150: *     .. Scalar Arguments ..
  151:       DOUBLE PRECISION   TOL
  152:       INTEGER            INFO, LDA, N, RANK
  153:       CHARACTER          UPLO
  154: *     ..
  155: *     .. Array Arguments ..
  156:       COMPLEX*16         A( LDA, * )
  157:       DOUBLE PRECISION   WORK( 2*N )
  158:       INTEGER            PIV( N )
  159: *     ..
  160: *
  161: *  =====================================================================
  162: *
  163: *     .. Parameters ..
  164:       DOUBLE PRECISION   ONE, ZERO
  165:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  166:       COMPLEX*16         CONE
  167:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
  168: *     ..
  169: *     .. Local Scalars ..
  170:       COMPLEX*16         ZTEMP
  171:       DOUBLE PRECISION   AJJ, DSTOP, DTEMP
  172:       INTEGER            I, ITEMP, J, PVT
  173:       LOGICAL            UPPER
  174: *     ..
  175: *     .. External Functions ..
  176:       DOUBLE PRECISION   DLAMCH
  177:       LOGICAL            LSAME, DISNAN
  178:       EXTERNAL           DLAMCH, LSAME, DISNAN
  179: *     ..
  180: *     .. External Subroutines ..
  181:       EXTERNAL           ZDSCAL, ZGEMV, ZLACGV, ZSWAP, XERBLA
  182: *     ..
  183: *     .. Intrinsic Functions ..
  184:       INTRINSIC          DBLE, DCONJG, MAX, SQRT
  185: *     ..
  186: *     .. Executable Statements ..
  187: *
  188: *     Test the input parameters
  189: *
  190:       INFO = 0
  191:       UPPER = LSAME( UPLO, 'U' )
  192:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  193:          INFO = -1
  194:       ELSE IF( N.LT.0 ) THEN
  195:          INFO = -2
  196:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  197:          INFO = -4
  198:       END IF
  199:       IF( INFO.NE.0 ) THEN
  200:          CALL XERBLA( 'ZPSTF2', -INFO )
  201:          RETURN
  202:       END IF
  203: *
  204: *     Quick return if possible
  205: *
  206:       IF( N.EQ.0 )
  207:      $   RETURN
  208: *
  209: *     Initialize PIV
  210: *
  211:       DO 100 I = 1, N
  212:          PIV( I ) = I
  213:   100 CONTINUE
  214: *
  215: *     Compute stopping value
  216: *
  217:       DO 110 I = 1, N
  218:          WORK( I ) = DBLE( A( I, I ) )
  219:   110 CONTINUE
  220:       PVT = MAXLOC( WORK( 1:N ), 1 )
  221:       AJJ = DBLE( A( PVT, PVT ) )
  222:       IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
  223:          RANK = 0
  224:          INFO = 1
  225:          GO TO 200
  226:       END IF
  227: *
  228: *     Compute stopping value if not supplied
  229: *
  230:       IF( TOL.LT.ZERO ) THEN
  231:          DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
  232:       ELSE
  233:          DSTOP = TOL
  234:       END IF
  235: *
  236: *     Set first half of WORK to zero, holds dot products
  237: *
  238:       DO 120 I = 1, N
  239:          WORK( I ) = 0
  240:   120 CONTINUE
  241: *
  242:       IF( UPPER ) THEN
  243: *
  244: *        Compute the Cholesky factorization P**T * A * P = U**H* U
  245: *
  246:          DO 150 J = 1, N
  247: *
  248: *        Find pivot, test for exit, else swap rows and columns
  249: *        Update dot products, compute possible pivots which are
  250: *        stored in the second half of WORK
  251: *
  252:             DO 130 I = J, N
  253: *
  254:                IF( J.GT.1 ) THEN
  255:                   WORK( I ) = WORK( I ) +
  256:      $                        DBLE( DCONJG( A( J-1, I ) )*
  257:      $                              A( J-1, I ) )
  258:                END IF
  259:                WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
  260: *
  261:   130       CONTINUE
  262: *
  263:             IF( J.GT.1 ) THEN
  264:                ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
  265:                PVT = ITEMP + J - 1
  266:                AJJ = WORK( N+PVT )
  267:                IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
  268:                   A( J, J ) = AJJ
  269:                   GO TO 190
  270:                END IF
  271:             END IF
  272: *
  273:             IF( J.NE.PVT ) THEN
  274: *
  275: *              Pivot OK, so can now swap pivot rows and columns
  276: *
  277:                A( PVT, PVT ) = A( J, J )
  278:                CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
  279:                IF( PVT.LT.N )
  280:      $            CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
  281:      $                        A( PVT, PVT+1 ), LDA )
  282:                DO 140 I = J + 1, PVT - 1
  283:                   ZTEMP = DCONJG( A( J, I ) )
  284:                   A( J, I ) = DCONJG( A( I, PVT ) )
  285:                   A( I, PVT ) = ZTEMP
  286:   140          CONTINUE
  287:                A( J, PVT ) = DCONJG( A( J, PVT ) )
  288: *
  289: *              Swap dot products and PIV
  290: *
  291:                DTEMP = WORK( J )
  292:                WORK( J ) = WORK( PVT )
  293:                WORK( PVT ) = DTEMP
  294:                ITEMP = PIV( PVT )
  295:                PIV( PVT ) = PIV( J )
  296:                PIV( J ) = ITEMP
  297:             END IF
  298: *
  299:             AJJ = SQRT( AJJ )
  300:             A( J, J ) = AJJ
  301: *
  302: *           Compute elements J+1:N of row J
  303: *
  304:             IF( J.LT.N ) THEN
  305:                CALL ZLACGV( J-1, A( 1, J ), 1 )
  306:                CALL ZGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA,
  307:      $                     A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
  308:                CALL ZLACGV( J-1, A( 1, J ), 1 )
  309:                CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
  310:             END IF
  311: *
  312:   150    CONTINUE
  313: *
  314:       ELSE
  315: *
  316: *        Compute the Cholesky factorization P**T * A * P = L * L**H
  317: *
  318:          DO 180 J = 1, N
  319: *
  320: *        Find pivot, test for exit, else swap rows and columns
  321: *        Update dot products, compute possible pivots which are
  322: *        stored in the second half of WORK
  323: *
  324:             DO 160 I = J, N
  325: *
  326:                IF( J.GT.1 ) THEN
  327:                   WORK( I ) = WORK( I ) +
  328:      $                        DBLE( DCONJG( A( I, J-1 ) )*
  329:      $                              A( I, J-1 ) )
  330:                END IF
  331:                WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
  332: *
  333:   160       CONTINUE
  334: *
  335:             IF( J.GT.1 ) THEN
  336:                ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
  337:                PVT = ITEMP + J - 1
  338:                AJJ = WORK( N+PVT )
  339:                IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
  340:                   A( J, J ) = AJJ
  341:                   GO TO 190
  342:                END IF
  343:             END IF
  344: *
  345:             IF( J.NE.PVT ) THEN
  346: *
  347: *              Pivot OK, so can now swap pivot rows and columns
  348: *
  349:                A( PVT, PVT ) = A( J, J )
  350:                CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
  351:                IF( PVT.LT.N )
  352:      $            CALL ZSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
  353:      $                        1 )
  354:                DO 170 I = J + 1, PVT - 1
  355:                   ZTEMP = DCONJG( A( I, J ) )
  356:                   A( I, J ) = DCONJG( A( PVT, I ) )
  357:                   A( PVT, I ) = ZTEMP
  358:   170          CONTINUE
  359:                A( PVT, J ) = DCONJG( A( PVT, J ) )
  360: *
  361: *              Swap dot products and PIV
  362: *
  363:                DTEMP = WORK( J )
  364:                WORK( J ) = WORK( PVT )
  365:                WORK( PVT ) = DTEMP
  366:                ITEMP = PIV( PVT )
  367:                PIV( PVT ) = PIV( J )
  368:                PIV( J ) = ITEMP
  369:             END IF
  370: *
  371:             AJJ = SQRT( AJJ )
  372:             A( J, J ) = AJJ
  373: *
  374: *           Compute elements J+1:N of column J
  375: *
  376:             IF( J.LT.N ) THEN
  377:                CALL ZLACGV( J-1, A( J, 1 ), LDA )
  378:                CALL ZGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ),
  379:      $                     LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
  380:                CALL ZLACGV( J-1, A( J, 1 ), LDA )
  381:                CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
  382:             END IF
  383: *
  384:   180    CONTINUE
  385: *
  386:       END IF
  387: *
  388: *     Ran to completion, A has full rank
  389: *
  390:       RANK = N
  391: *
  392:       GO TO 200
  393:   190 CONTINUE
  394: *
  395: *     Rank is number of steps completed.  Set INFO = 1 to signal
  396: *     that the factorization cannot be used to solve a system.
  397: *
  398:       RANK = J - 1
  399:       INFO = 1
  400: *
  401:   200 CONTINUE
  402:       RETURN
  403: *
  404: *     End of ZPSTF2
  405: *
  406:       END