Return to zgbcon.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b ZGBCON
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZGBCON + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbcon.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbcon.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbcon.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
   22: *                          WORK, RWORK, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          NORM
   26: *       INTEGER            INFO, KL, KU, LDAB, N
   27: *       DOUBLE PRECISION   ANORM, RCOND
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       INTEGER            IPIV( * )
   31: *       DOUBLE PRECISION   RWORK( * )
   32: *       COMPLEX*16         AB( LDAB, * ), WORK( * )
   33: *       ..
   34: *  
   35: *
   36: *> \par Purpose:
   37: *  =============
   38: *>
   39: *> \verbatim
   40: *>
   41: *> ZGBCON estimates the reciprocal of the condition number of a complex
   42: *> general band matrix A, in either the 1-norm or the infinity-norm,
   43: *> using the LU factorization computed by ZGBTRF.
   44: *>
   45: *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
   46: *> condition number is computed as
   47: *>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
   48: *> \endverbatim
   49: *
   50: *  Arguments:
   51: *  ==========
   52: *
   53: *> \param[in] NORM
   54: *> \verbatim
   55: *>          NORM is CHARACTER*1
   56: *>          Specifies whether the 1-norm condition number or the
   57: *>          infinity-norm condition number is required:
   58: *>          = '1' or 'O':  1-norm;
   59: *>          = 'I':         Infinity-norm.
   60: *> \endverbatim
   61: *>
   62: *> \param[in] N
   63: *> \verbatim
   64: *>          N is INTEGER
   65: *>          The order of the matrix A.  N >= 0.
   66: *> \endverbatim
   67: *>
   68: *> \param[in] KL
   69: *> \verbatim
   70: *>          KL is INTEGER
   71: *>          The number of subdiagonals within the band of A.  KL >= 0.
   72: *> \endverbatim
   73: *>
   74: *> \param[in] KU
   75: *> \verbatim
   76: *>          KU is INTEGER
   77: *>          The number of superdiagonals within the band of A.  KU >= 0.
   78: *> \endverbatim
   79: *>
   80: *> \param[in] AB
   81: *> \verbatim
   82: *>          AB is COMPLEX*16 array, dimension (LDAB,N)
   83: *>          Details of the LU factorization of the band matrix A, as
   84: *>          computed by ZGBTRF.  U is stored as an upper triangular band
   85: *>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
   86: *>          the multipliers used during the factorization are stored in
   87: *>          rows KL+KU+2 to 2*KL+KU+1.
   88: *> \endverbatim
   89: *>
   90: *> \param[in] LDAB
   91: *> \verbatim
   92: *>          LDAB is INTEGER
   93: *>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
   94: *> \endverbatim
   95: *>
   96: *> \param[in] IPIV
   97: *> \verbatim
   98: *>          IPIV is INTEGER array, dimension (N)
   99: *>          The pivot indices; for 1 <= i <= N, row i of the matrix was
  100: *>          interchanged with row IPIV(i).
  101: *> \endverbatim
  102: *>
  103: *> \param[in] ANORM
  104: *> \verbatim
  105: *>          ANORM is DOUBLE PRECISION
  106: *>          If NORM = '1' or 'O', the 1-norm of the original matrix A.
  107: *>          If NORM = 'I', the infinity-norm of the original matrix A.
  108: *> \endverbatim
  109: *>
  110: *> \param[out] RCOND
  111: *> \verbatim
  112: *>          RCOND is DOUBLE PRECISION
  113: *>          The reciprocal of the condition number of the matrix A,
  114: *>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
  115: *> \endverbatim
  116: *>
  117: *> \param[out] WORK
  118: *> \verbatim
  119: *>          WORK is COMPLEX*16 array, dimension (2*N)
  120: *> \endverbatim
  121: *>
  122: *> \param[out] RWORK
  123: *> \verbatim
  124: *>          RWORK is DOUBLE PRECISION array, dimension (N)
  125: *> \endverbatim
  126: *>
  127: *> \param[out] INFO
  128: *> \verbatim
  129: *>          INFO is INTEGER
  130: *>          = 0:  successful exit
  131: *>          < 0: if INFO = -i, the i-th argument had an illegal value
  132: *> \endverbatim
  133: *
  134: *  Authors:
  135: *  ========
  136: *
  137: *> \author Univ. of Tennessee 
  138: *> \author Univ. of California Berkeley 
  139: *> \author Univ. of Colorado Denver 
  140: *> \author NAG Ltd. 
  141: *
  142: *> \date November 2011
  143: *
  144: *> \ingroup complex16GBcomputational
  145: *
  146: *  =====================================================================
  147:       SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
  148:      $                   WORK, RWORK, INFO )
  149: *
  150: *  -- LAPACK computational routine (version 3.4.0) --
  151: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  152: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  153: *     November 2011
  154: *
  155: *     .. Scalar Arguments ..
  156:       CHARACTER          NORM
  157:       INTEGER            INFO, KL, KU, LDAB, N
  158:       DOUBLE PRECISION   ANORM, RCOND
  159: *     ..
  160: *     .. Array Arguments ..
  161:       INTEGER            IPIV( * )
  162:       DOUBLE PRECISION   RWORK( * )
  163:       COMPLEX*16         AB( LDAB, * ), WORK( * )
  164: *     ..
  165: *
  166: *  =====================================================================
  167: *
  168: *     .. Parameters ..
  169:       DOUBLE PRECISION   ONE, ZERO
  170:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  171: *     ..
  172: *     .. Local Scalars ..
  173:       LOGICAL            LNOTI, ONENRM
  174:       CHARACTER          NORMIN
  175:       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
  176:       DOUBLE PRECISION   AINVNM, SCALE, SMLNUM
  177:       COMPLEX*16         T, ZDUM
  178: *     ..
  179: *     .. Local Arrays ..
  180:       INTEGER            ISAVE( 3 )
  181: *     ..
  182: *     .. External Functions ..
  183:       LOGICAL            LSAME
  184:       INTEGER            IZAMAX
  185:       DOUBLE PRECISION   DLAMCH
  186:       COMPLEX*16         ZDOTC
  187:       EXTERNAL           LSAME, IZAMAX, DLAMCH, ZDOTC
  188: *     ..
  189: *     .. External Subroutines ..
  190:       EXTERNAL           XERBLA, ZAXPY, ZDRSCL, ZLACN2, ZLATBS
  191: *     ..
  192: *     .. Intrinsic Functions ..
  193:       INTRINSIC          ABS, DBLE, DIMAG, MIN
  194: *     ..
  195: *     .. Statement Functions ..
  196:       DOUBLE PRECISION   CABS1
  197: *     ..
  198: *     .. Statement Function definitions ..
  199:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  200: *     ..
  201: *     .. Executable Statements ..
  202: *
  203: *     Test the input parameters.
  204: *
  205:       INFO = 0
  206:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
  207:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
  208:          INFO = -1
  209:       ELSE IF( N.LT.0 ) THEN
  210:          INFO = -2
  211:       ELSE IF( KL.LT.0 ) THEN
  212:          INFO = -3
  213:       ELSE IF( KU.LT.0 ) THEN
  214:          INFO = -4
  215:       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
  216:          INFO = -6
  217:       ELSE IF( ANORM.LT.ZERO ) THEN
  218:          INFO = -8
  219:       END IF
  220:       IF( INFO.NE.0 ) THEN
  221:          CALL XERBLA( 'ZGBCON', -INFO )
  222:          RETURN
  223:       END IF
  224: *
  225: *     Quick return if possible
  226: *
  227:       RCOND = ZERO
  228:       IF( N.EQ.0 ) THEN
  229:          RCOND = ONE
  230:          RETURN
  231:       ELSE IF( ANORM.EQ.ZERO ) THEN
  232:          RETURN
  233:       END IF
  234: *
  235:       SMLNUM = DLAMCH( 'Safe minimum' )
  236: *
  237: *     Estimate the norm of inv(A).
  238: *
  239:       AINVNM = ZERO
  240:       NORMIN = 'N'
  241:       IF( ONENRM ) THEN
  242:          KASE1 = 1
  243:       ELSE
  244:          KASE1 = 2
  245:       END IF
  246:       KD = KL + KU + 1
  247:       LNOTI = KL.GT.0
  248:       KASE = 0
  249:    10 CONTINUE
  250:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  251:       IF( KASE.NE.0 ) THEN
  252:          IF( KASE.EQ.KASE1 ) THEN
  253: *
  254: *           Multiply by inv(L).
  255: *
  256:             IF( LNOTI ) THEN
  257:                DO 20 J = 1, N - 1
  258:                   LM = MIN( KL, N-J )
  259:                   JP = IPIV( J )
  260:                   T = WORK( JP )
  261:                   IF( JP.NE.J ) THEN
  262:                      WORK( JP ) = WORK( J )
  263:                      WORK( J ) = T
  264:                   END IF
  265:                   CALL ZAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
  266:    20          CONTINUE
  267:             END IF
  268: *
  269: *           Multiply by inv(U).
  270: *
  271:             CALL ZLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
  272:      $                   KL+KU, AB, LDAB, WORK, SCALE, RWORK, INFO )
  273:          ELSE
  274: *
  275: *           Multiply by inv(U**H).
  276: *
  277:             CALL ZLATBS( 'Upper', 'Conjugate transpose', 'Non-unit',
  278:      $                   NORMIN, N, KL+KU, AB, LDAB, WORK, SCALE, RWORK,
  279:      $                   INFO )
  280: *
  281: *           Multiply by inv(L**H).
  282: *
  283:             IF( LNOTI ) THEN
  284:                DO 30 J = N - 1, 1, -1
  285:                   LM = MIN( KL, N-J )
  286:                   WORK( J ) = WORK( J ) - ZDOTC( LM, AB( KD+1, J ), 1,
  287:      $                        WORK( J+1 ), 1 )
  288:                   JP = IPIV( J )
  289:                   IF( JP.NE.J ) THEN
  290:                      T = WORK( JP )
  291:                      WORK( JP ) = WORK( J )
  292:                      WORK( J ) = T
  293:                   END IF
  294:    30          CONTINUE
  295:             END IF
  296:          END IF
  297: *
  298: *        Divide X by 1/SCALE if doing so will not cause overflow.
  299: *
  300:          NORMIN = 'Y'
  301:          IF( SCALE.NE.ONE ) THEN
  302:             IX = IZAMAX( N, WORK, 1 )
  303:             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
  304:      $         GO TO 40
  305:             CALL ZDRSCL( N, SCALE, WORK, 1 )
  306:          END IF
  307:          GO TO 10
  308:       END IF
  309: *
  310: *     Compute the estimate of the reciprocal condition number.
  311: *
  312:       IF( AINVNM.NE.ZERO )
  313:      $   RCOND = ( ONE / AINVNM ) / ANORM
  314: *
  315:    40 CONTINUE
  316:       RETURN
  317: *
  318: *     End of ZGBCON
  319: *
  320:       END