Return to zsycon.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE ZSYCON( UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK,
    2:      $                   INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
   10: *
   11: *     .. Scalar Arguments ..
   12:       CHARACTER          UPLO
   13:       INTEGER            INFO, LDA, N
   14:       DOUBLE PRECISION   ANORM, RCOND
   15: *     ..
   16: *     .. Array Arguments ..
   17:       INTEGER            IPIV( * )
   18:       COMPLEX*16         A( LDA, * ), WORK( * )
   19: *     ..
   20: *
   21: *  Purpose
   22: *  =======
   23: *
   24: *  ZSYCON estimates the reciprocal of the condition number (in the
   25: *  1-norm) of a complex symmetric matrix A using the factorization
   26: *  A = U*D*U**T or A = L*D*L**T computed by ZSYTRF.
   27: *
   28: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
   29: *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
   30: *
   31: *  Arguments
   32: *  =========
   33: *
   34: *  UPLO    (input) CHARACTER*1
   35: *          Specifies whether the details of the factorization are stored
   36: *          as an upper or lower triangular matrix.
   37: *          = 'U':  Upper triangular, form is A = U*D*U**T;
   38: *          = 'L':  Lower triangular, form is A = L*D*L**T.
   39: *
   40: *  N       (input) INTEGER
   41: *          The order of the matrix A.  N >= 0.
   42: *
   43: *  A       (input) COMPLEX*16 array, dimension (LDA,N)
   44: *          The block diagonal matrix D and the multipliers used to
   45: *          obtain the factor U or L as computed by ZSYTRF.
   46: *
   47: *  LDA     (input) INTEGER
   48: *          The leading dimension of the array A.  LDA >= max(1,N).
   49: *
   50: *  IPIV    (input) INTEGER array, dimension (N)
   51: *          Details of the interchanges and the block structure of D
   52: *          as determined by ZSYTRF.
   53: *
   54: *  ANORM   (input) DOUBLE PRECISION
   55: *          The 1-norm of the original matrix A.
   56: *
   57: *  RCOND   (output) DOUBLE PRECISION
   58: *          The reciprocal of the condition number of the matrix A,
   59: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
   60: *          estimate of the 1-norm of inv(A) computed in this routine.
   61: *
   62: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
   63: *
   64: *  INFO    (output) INTEGER
   65: *          = 0:  successful exit
   66: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   67: *
   68: *  =====================================================================
   69: *
   70: *     .. Parameters ..
   71:       DOUBLE PRECISION   ONE, ZERO
   72:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   73: *     ..
   74: *     .. Local Scalars ..
   75:       LOGICAL            UPPER
   76:       INTEGER            I, KASE
   77:       DOUBLE PRECISION   AINVNM
   78: *     ..
   79: *     .. Local Arrays ..
   80:       INTEGER            ISAVE( 3 )
   81: *     ..
   82: *     .. External Functions ..
   83:       LOGICAL            LSAME
   84:       EXTERNAL           LSAME
   85: *     ..
   86: *     .. External Subroutines ..
   87:       EXTERNAL           XERBLA, ZLACN2, ZSYTRS
   88: *     ..
   89: *     .. Intrinsic Functions ..
   90:       INTRINSIC          MAX
   91: *     ..
   92: *     .. Executable Statements ..
   93: *
   94: *     Test the input parameters.
   95: *
   96:       INFO = 0
   97:       UPPER = LSAME( UPLO, 'U' )
   98:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
   99:          INFO = -1
  100:       ELSE IF( N.LT.0 ) THEN
  101:          INFO = -2
  102:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  103:          INFO = -4
  104:       ELSE IF( ANORM.LT.ZERO ) THEN
  105:          INFO = -6
  106:       END IF
  107:       IF( INFO.NE.0 ) THEN
  108:          CALL XERBLA( 'ZSYCON', -INFO )
  109:          RETURN
  110:       END IF
  111: *
  112: *     Quick return if possible
  113: *
  114:       RCOND = ZERO
  115:       IF( N.EQ.0 ) THEN
  116:          RCOND = ONE
  117:          RETURN
  118:       ELSE IF( ANORM.LE.ZERO ) THEN
  119:          RETURN
  120:       END IF
  121: *
  122: *     Check that the diagonal matrix D is nonsingular.
  123: *
  124:       IF( UPPER ) THEN
  125: *
  126: *        Upper triangular storage: examine D from bottom to top
  127: *
  128:          DO 10 I = N, 1, -1
  129:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO )
  130:      $         RETURN
  131:    10    CONTINUE
  132:       ELSE
  133: *
  134: *        Lower triangular storage: examine D from top to bottom.
  135: *
  136:          DO 20 I = 1, N
  137:             IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO )
  138:      $         RETURN
  139:    20    CONTINUE
  140:       END IF
  141: *
  142: *     Estimate the 1-norm of the inverse.
  143: *
  144:       KASE = 0
  145:    30 CONTINUE
  146:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  147:       IF( KASE.NE.0 ) THEN
  148: *
  149: *        Multiply by inv(L*D*L') or inv(U*D*U').
  150: *
  151:          CALL ZSYTRS( UPLO, N, 1, A, LDA, IPIV, WORK, N, INFO )
  152:          GO TO 30
  153:       END IF
  154: *
  155: *     Compute the estimate of the reciprocal condition number.
  156: *
  157:       IF( AINVNM.NE.ZERO )
  158:      $   RCOND = ( ONE / AINVNM ) / ANORM
  159: *
  160:       RETURN
  161: *
  162: *     End of ZSYCON
  163: *
  164:       END