Return to zsycon_3.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief \b ZSYCON_3
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZSYCON_3 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsycon_3.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsycon_3.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsycon_3.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZSYCON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND,
   22: *                            WORK, 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( * )
   31: *       COMPLEX*16         A( LDA, * ), E ( * ), WORK( * )
   32: *       ..
   33: *
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *> ZSYCON_3 estimates the reciprocal of the condition number (in the
   40: *> 1-norm) of a complex symmetric matrix A using the factorization
   41: *> computed by ZSYTRF_RK or ZSYTRF_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 ZSYTRS_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 COMPLEX*16 array, dimension (LDA,N)
   76: *>          Diagonal of the block diagonal matrix D and factors U or L
   77: *>          as computed by ZSYTRF_RK and ZSYTRF_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 COMPLEX*16 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 ZSYTRF_RK or ZSYTRF_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 COMPLEX*16 array, dimension (2*N)
  130: *> \endverbatim
  131: *>
  132: *> \param[out] INFO
  133: *> \verbatim
  134: *>          INFO is INTEGER
  135: *>          = 0:  successful exit
  136: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  137: *> \endverbatim
  138: *
  139: *  Authors:
  140: *  ========
  141: *
  142: *> \author Univ. of Tennessee
  143: *> \author Univ. of California Berkeley
  144: *> \author Univ. of Colorado Denver
  145: *> \author NAG Ltd.
  146: *
  147: *> \ingroup complex16SYcomputational
  148: *
  149: *> \par Contributors:
  150: *  ==================
  151: *> \verbatim
  152: *>
  153: *>  June 2017,  Igor Kozachenko,
  154: *>                  Computer Science Division,
  155: *>                  University of California, Berkeley
  156: *>
  157: *>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
  158: *>                  School of Mathematics,
  159: *>                  University of Manchester
  160: *>
  161: *> \endverbatim
  162: *
  163: *  =====================================================================
  164:       SUBROUTINE ZSYCON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND,
  165:      $                     WORK, INFO )
  166: *
  167: *  -- LAPACK computational routine --
  168: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  169: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  170: *
  171: *     .. Scalar Arguments ..
  172:       CHARACTER          UPLO
  173:       INTEGER            INFO, LDA, N
  174:       DOUBLE PRECISION   ANORM, RCOND
  175: *     ..
  176: *     .. Array Arguments ..
  177:       INTEGER            IPIV( * )
  178:       COMPLEX*16         A( LDA, * ), E( * ), WORK( * )
  179: *     ..
  180: *
  181: *  =====================================================================
  182: *
  183: *     .. Parameters ..
  184:       DOUBLE PRECISION   ONE, ZERO
  185:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  186:       COMPLEX*16         CZERO
  187:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
  188: *     ..
  189: *     .. Local Scalars ..
  190:       LOGICAL            UPPER
  191:       INTEGER            I, KASE
  192:       DOUBLE PRECISION   AINVNM
  193: *     ..
  194: *     .. Local Arrays ..
  195:       INTEGER            ISAVE( 3 )
  196: *     ..
  197: *     .. External Functions ..
  198:       LOGICAL            LSAME
  199:       EXTERNAL           LSAME
  200: *     ..
  201: *     .. External Subroutines ..
  202:       EXTERNAL           ZLACN2, ZSYTRS_3, XERBLA
  203: *     ..
  204: *     .. Intrinsic Functions ..
  205:       INTRINSIC          MAX
  206: *     ..
  207: *     .. Executable Statements ..
  208: *
  209: *     Test the input parameters.
  210: *
  211:       INFO = 0
  212:       UPPER = LSAME( UPLO, 'U' )
  213:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  214:          INFO = -1
  215:       ELSE IF( N.LT.0 ) THEN
  216:          INFO = -2
  217:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  218:          INFO = -4
  219:       ELSE IF( ANORM.LT.ZERO ) THEN
  220:          INFO = -7
  221:       END IF
  222:       IF( INFO.NE.0 ) THEN
  223:          CALL XERBLA( 'ZSYCON_3', -INFO )
  224:          RETURN
  225:       END IF
  226: *
  227: *     Quick return if possible
  228: *
  229:       RCOND = ZERO
  230:       IF( N.EQ.0 ) THEN
  231:          RCOND = ONE
  232:          RETURN
  233:       ELSE IF( ANORM.LE.ZERO ) THEN
  234:          RETURN
  235:       END IF
  236: *
  237: *     Check that the diagonal matrix D is nonsingular.
  238: *
  239:       IF( UPPER ) THEN
  240: *
  241: *        Upper triangular storage: examine D from bottom to top
  242: *
  243:          DO I = N, 1, -1
  244:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO )
  245:      $         RETURN
  246:          END DO
  247:       ELSE
  248: *
  249: *        Lower triangular storage: examine D from top to bottom.
  250: *
  251:          DO I = 1, N
  252:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO )
  253:      $         RETURN
  254:          END DO
  255:       END IF
  256: *
  257: *     Estimate the 1-norm of the inverse.
  258: *
  259:       KASE = 0
  260:    30 CONTINUE
  261:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  262:       IF( KASE.NE.0 ) THEN
  263: *
  264: *        Multiply by inv(L*D*L**T) or inv(U*D*U**T).
  265: *
  266:          CALL ZSYTRS_3( UPLO, N, 1, A, LDA, E, IPIV, WORK, N, INFO )
  267:          GO TO 30
  268:       END IF
  269: *
  270: *     Compute the estimate of the reciprocal condition number.
  271: *
  272:       IF( AINVNM.NE.ZERO )
  273:      $   RCOND = ( ONE / AINVNM ) / ANORM
  274: *
  275:       RETURN
  276: *
  277: *     End of ZSYCON_3
  278: *
  279:       END