Return to dsycon_3.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief \b DSYCON_3
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DSYCON_3 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsycon_3.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsycon_3.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsycon_3.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DSYCON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND,
   22: *                            WORK, IWORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          UPLO
   26: *       INTEGER            INFO, LDA, N
   27: *       DOUBLE PRECISION   ANORM, RCOND
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       INTEGER            IPIV( * ), IWORK( * )
   31: *       DOUBLE PRECISION   A( LDA, * ), E ( * ), WORK( * )
   32: *       ..
   33: *
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *> DSYCON_3 estimates the reciprocal of the condition number (in the
   40: *> 1-norm) of a real symmetric matrix A using the factorization
   41: *> computed by DSYTRF_RK or DSYTRF_BK:
   42: *>
   43: *>    A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
   44: *>
   45: *> where U (or L) is unit upper (or lower) triangular matrix,
   46: *> U**T (or L**T) is the transpose of U (or L), P is a permutation
   47: *> matrix, P**T is the transpose of P, and D is symmetric and block
   48: *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
   49: *>
   50: *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
   51: *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
   52: *> This routine uses BLAS3 solver DSYTRS_3.
   53: *> \endverbatim
   54: *
   55: *  Arguments:
   56: *  ==========
   57: *
   58: *> \param[in] UPLO
   59: *> \verbatim
   60: *>          UPLO is CHARACTER*1
   61: *>          Specifies whether the details of the factorization are
   62: *>          stored as an upper or lower triangular matrix:
   63: *>          = 'U':  Upper triangular, form is A = P*U*D*(U**T)*(P**T);
   64: *>          = 'L':  Lower triangular, form is A = P*L*D*(L**T)*(P**T).
   65: *> \endverbatim
   66: *>
   67: *> \param[in] N
   68: *> \verbatim
   69: *>          N is INTEGER
   70: *>          The order of the matrix A.  N >= 0.
   71: *> \endverbatim
   72: *>
   73: *> \param[in] A
   74: *> \verbatim
   75: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   76: *>          Diagonal of the block diagonal matrix D and factors U or L
   77: *>          as computed by DSYTRF_RK and DSYTRF_BK:
   78: *>            a) ONLY diagonal elements of the symmetric block diagonal
   79: *>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
   80: *>               (superdiagonal (or subdiagonal) elements of D
   81: *>                should be provided on entry in array E), and
   82: *>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
   83: *>               If UPLO = 'L': factor L in the subdiagonal part of A.
   84: *> \endverbatim
   85: *>
   86: *> \param[in] LDA
   87: *> \verbatim
   88: *>          LDA is INTEGER
   89: *>          The leading dimension of the array A.  LDA >= max(1,N).
   90: *> \endverbatim
   91: *>
   92: *> \param[in] E
   93: *> \verbatim
   94: *>          E is DOUBLE PRECISION array, dimension (N)
   95: *>          On entry, contains the superdiagonal (or subdiagonal)
   96: *>          elements of the symmetric block diagonal matrix D
   97: *>          with 1-by-1 or 2-by-2 diagonal blocks, where
   98: *>          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced;
   99: *>          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.
  100: *>
  101: *>          NOTE: For 1-by-1 diagonal block D(k), where
  102: *>          1 <= k <= N, the element E(k) is not referenced in both
  103: *>          UPLO = 'U' or UPLO = 'L' cases.
  104: *> \endverbatim
  105: *>
  106: *> \param[in] IPIV
  107: *> \verbatim
  108: *>          IPIV is INTEGER array, dimension (N)
  109: *>          Details of the interchanges and the block structure of D
  110: *>          as determined by DSYTRF_RK or DSYTRF_BK.
  111: *> \endverbatim
  112: *>
  113: *> \param[in] ANORM
  114: *> \verbatim
  115: *>          ANORM is DOUBLE PRECISION
  116: *>          The 1-norm of the original matrix A.
  117: *> \endverbatim
  118: *>
  119: *> \param[out] RCOND
  120: *> \verbatim
  121: *>          RCOND is DOUBLE PRECISION
  122: *>          The reciprocal of the condition number of the matrix A,
  123: *>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
  124: *>          estimate of the 1-norm of inv(A) computed in this routine.
  125: *> \endverbatim
  126: *>
  127: *> \param[out] WORK
  128: *> \verbatim
  129: *>          WORK is DOUBLE PRECISION array, dimension (2*N)
  130: *> \endverbatim
  131: *>
  132: *> \param[out] IWORK
  133: *> \verbatim
  134: *>          IWORK is INTEGER array, dimension (N)
  135: *> \endverbatim
  136: *>
  137: *> \param[out] INFO
  138: *> \verbatim
  139: *>          INFO is INTEGER
  140: *>          = 0:  successful exit
  141: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  142: *> \endverbatim
  143: *
  144: *  Authors:
  145: *  ========
  146: *
  147: *> \author Univ. of Tennessee
  148: *> \author Univ. of California Berkeley
  149: *> \author Univ. of Colorado Denver
  150: *> \author NAG Ltd.
  151: *
  152: *> \ingroup doubleSYcomputational
  153: *
  154: *> \par Contributors:
  155: *  ==================
  156: *> \verbatim
  157: *>
  158: *>  June 2017,  Igor Kozachenko,
  159: *>                  Computer Science Division,
  160: *>                  University of California, Berkeley
  161: *>
  162: *>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
  163: *>                  School of Mathematics,
  164: *>                  University of Manchester
  165: *>
  166: *> \endverbatim
  167: *
  168: *  =====================================================================
  169:       SUBROUTINE DSYCON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND,
  170:      $                     WORK, IWORK, INFO )
  171: *
  172: *  -- LAPACK computational routine --
  173: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  174: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  175: *
  176: *     .. Scalar Arguments ..
  177:       CHARACTER          UPLO
  178:       INTEGER            INFO, LDA, N
  179:       DOUBLE PRECISION   ANORM, RCOND
  180: *     ..
  181: *     .. Array Arguments ..
  182:       INTEGER            IPIV( * ), IWORK( * )
  183:       DOUBLE PRECISION   A( LDA, * ), E( * ), WORK( * )
  184: *     ..
  185: *
  186: *  =====================================================================
  187: *
  188: *     .. Parameters ..
  189:       DOUBLE PRECISION   ONE, ZERO
  190:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  191: *     ..
  192: *     .. Local Scalars ..
  193:       LOGICAL            UPPER
  194:       INTEGER            I, KASE
  195:       DOUBLE PRECISION   AINVNM
  196: *     ..
  197: *     .. Local Arrays ..
  198:       INTEGER            ISAVE( 3 )
  199: *     ..
  200: *     .. External Functions ..
  201:       LOGICAL            LSAME
  202:       EXTERNAL           LSAME
  203: *     ..
  204: *     .. External Subroutines ..
  205:       EXTERNAL           DLACN2, DSYTRS_3, XERBLA
  206: *     ..
  207: *     .. Intrinsic Functions ..
  208:       INTRINSIC          MAX
  209: *     ..
  210: *     .. Executable Statements ..
  211: *
  212: *     Test the input parameters.
  213: *
  214:       INFO = 0
  215:       UPPER = LSAME( UPLO, 'U' )
  216:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  217:          INFO = -1
  218:       ELSE IF( N.LT.0 ) THEN
  219:          INFO = -2
  220:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  221:          INFO = -4
  222:       ELSE IF( ANORM.LT.ZERO ) THEN
  223:          INFO = -7
  224:       END IF
  225:       IF( INFO.NE.0 ) THEN
  226:          CALL XERBLA( 'DSYCON_3', -INFO )
  227:          RETURN
  228:       END IF
  229: *
  230: *     Quick return if possible
  231: *
  232:       RCOND = ZERO
  233:       IF( N.EQ.0 ) THEN
  234:          RCOND = ONE
  235:          RETURN
  236:       ELSE IF( ANORM.LE.ZERO ) THEN
  237:          RETURN
  238:       END IF
  239: *
  240: *     Check that the diagonal matrix D is nonsingular.
  241: *
  242:       IF( UPPER ) THEN
  243: *
  244: *        Upper triangular storage: examine D from bottom to top
  245: *
  246:          DO I = N, 1, -1
  247:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO )
  248:      $         RETURN
  249:          END DO
  250:       ELSE
  251: *
  252: *        Lower triangular storage: examine D from top to bottom.
  253: *
  254:          DO I = 1, N
  255:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO )
  256:      $         RETURN
  257:          END DO
  258:       END IF
  259: *
  260: *     Estimate the 1-norm of the inverse.
  261: *
  262:       KASE = 0
  263:    30 CONTINUE
  264:       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
  265:       IF( KASE.NE.0 ) THEN
  266: *
  267: *        Multiply by inv(L*D*L**T) or inv(U*D*U**T).
  268: *
  269:          CALL DSYTRS_3( UPLO, N, 1, A, LDA, E, IPIV, WORK, N, INFO )
  270:          GO TO 30
  271:       END IF
  272: *
  273: *     Compute the estimate of the reciprocal condition number.
  274: *
  275:       IF( AINVNM.NE.ZERO )
  276:      $   RCOND = ( ONE / AINVNM ) / ANORM
  277: *
  278:       RETURN
  279: *
  280: *     End of DSYCON_3
  281: *
  282:       END