File:  [local] / rpl / lapack / lapack / dtbcon.f
Revision 1.14: download - view: text, annotated - select for diffs - revision graph
Sat Aug 27 15:34:41 2016 UTC (7 years, 8 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_25, HEAD
Cohérence Lapack.

    1: *> \brief \b DTBCON
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DTBCON + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtbcon.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtbcon.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtbcon.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DTBCON( NORM, UPLO, DIAG, N, KD, AB, LDAB, RCOND, WORK,
   22: *                          IWORK, INFO )
   23:    24: *       .. Scalar Arguments ..
   25: *       CHARACTER          DIAG, NORM, UPLO
   26: *       INTEGER            INFO, KD, LDAB, N
   27: *       DOUBLE PRECISION   RCOND
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       INTEGER            IWORK( * )
   31: *       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )
   32: *       ..
   33: *  
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *>
   40: *> DTBCON estimates the reciprocal of the condition number of a
   41: *> triangular band matrix A, in either the 1-norm or the infinity-norm.
   42: *>
   43: *> The norm of A is computed and an estimate is obtained for
   44: *> norm(inv(A)), then the reciprocal of the condition number is
   45: *> computed as
   46: *>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
   47: *> \endverbatim
   48: *
   49: *  Arguments:
   50: *  ==========
   51: *
   52: *> \param[in] NORM
   53: *> \verbatim
   54: *>          NORM is CHARACTER*1
   55: *>          Specifies whether the 1-norm condition number or the
   56: *>          infinity-norm condition number is required:
   57: *>          = '1' or 'O':  1-norm;
   58: *>          = 'I':         Infinity-norm.
   59: *> \endverbatim
   60: *>
   61: *> \param[in] UPLO
   62: *> \verbatim
   63: *>          UPLO is CHARACTER*1
   64: *>          = 'U':  A is upper triangular;
   65: *>          = 'L':  A is lower triangular.
   66: *> \endverbatim
   67: *>
   68: *> \param[in] DIAG
   69: *> \verbatim
   70: *>          DIAG is CHARACTER*1
   71: *>          = 'N':  A is non-unit triangular;
   72: *>          = 'U':  A is unit triangular.
   73: *> \endverbatim
   74: *>
   75: *> \param[in] N
   76: *> \verbatim
   77: *>          N is INTEGER
   78: *>          The order of the matrix A.  N >= 0.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] KD
   82: *> \verbatim
   83: *>          KD is INTEGER
   84: *>          The number of superdiagonals or subdiagonals of the
   85: *>          triangular band matrix A.  KD >= 0.
   86: *> \endverbatim
   87: *>
   88: *> \param[in] AB
   89: *> \verbatim
   90: *>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
   91: *>          The upper or lower triangular band matrix A, stored in the
   92: *>          first kd+1 rows of the array. The j-th column of A is stored
   93: *>          in the j-th column of the array AB as follows:
   94: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
   95: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
   96: *>          If DIAG = 'U', the diagonal elements of A are not referenced
   97: *>          and are assumed to be 1.
   98: *> \endverbatim
   99: *>
  100: *> \param[in] LDAB
  101: *> \verbatim
  102: *>          LDAB is INTEGER
  103: *>          The leading dimension of the array AB.  LDAB >= KD+1.
  104: *> \endverbatim
  105: *>
  106: *> \param[out] RCOND
  107: *> \verbatim
  108: *>          RCOND is DOUBLE PRECISION
  109: *>          The reciprocal of the condition number of the matrix A,
  110: *>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
  111: *> \endverbatim
  112: *>
  113: *> \param[out] WORK
  114: *> \verbatim
  115: *>          WORK is DOUBLE PRECISION array, dimension (3*N)
  116: *> \endverbatim
  117: *>
  118: *> \param[out] IWORK
  119: *> \verbatim
  120: *>          IWORK is INTEGER array, dimension (N)
  121: *> \endverbatim
  122: *>
  123: *> \param[out] INFO
  124: *> \verbatim
  125: *>          INFO is INTEGER
  126: *>          = 0:  successful exit
  127: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  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 November 2011
  139: *
  140: *> \ingroup doubleOTHERcomputational
  141: *
  142: *  =====================================================================
  143:       SUBROUTINE DTBCON( NORM, UPLO, DIAG, N, KD, AB, LDAB, RCOND, WORK,
  144:      $                   IWORK, INFO )
  145: *
  146: *  -- LAPACK computational routine (version 3.4.0) --
  147: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  148: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  149: *     November 2011
  150: *
  151: *     .. Scalar Arguments ..
  152:       CHARACTER          DIAG, NORM, UPLO
  153:       INTEGER            INFO, KD, LDAB, N
  154:       DOUBLE PRECISION   RCOND
  155: *     ..
  156: *     .. Array Arguments ..
  157:       INTEGER            IWORK( * )
  158:       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )
  159: *     ..
  160: *
  161: *  =====================================================================
  162: *
  163: *     .. Parameters ..
  164:       DOUBLE PRECISION   ONE, ZERO
  165:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  166: *     ..
  167: *     .. Local Scalars ..
  168:       LOGICAL            NOUNIT, ONENRM, UPPER
  169:       CHARACTER          NORMIN
  170:       INTEGER            IX, KASE, KASE1
  171:       DOUBLE PRECISION   AINVNM, ANORM, SCALE, SMLNUM, XNORM
  172: *     ..
  173: *     .. Local Arrays ..
  174:       INTEGER            ISAVE( 3 )
  175: *     ..
  176: *     .. External Functions ..
  177:       LOGICAL            LSAME
  178:       INTEGER            IDAMAX
  179:       DOUBLE PRECISION   DLAMCH, DLANTB
  180:       EXTERNAL           LSAME, IDAMAX, DLAMCH, DLANTB
  181: *     ..
  182: *     .. External Subroutines ..
  183:       EXTERNAL           DLACN2, DLATBS, DRSCL, XERBLA
  184: *     ..
  185: *     .. Intrinsic Functions ..
  186:       INTRINSIC          ABS, DBLE, MAX
  187: *     ..
  188: *     .. Executable Statements ..
  189: *
  190: *     Test the input parameters.
  191: *
  192:       INFO = 0
  193:       UPPER = LSAME( UPLO, 'U' )
  194:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
  195:       NOUNIT = LSAME( DIAG, 'N' )
  196: *
  197:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
  198:          INFO = -1
  199:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  200:          INFO = -2
  201:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
  202:          INFO = -3
  203:       ELSE IF( N.LT.0 ) THEN
  204:          INFO = -4
  205:       ELSE IF( KD.LT.0 ) THEN
  206:          INFO = -5
  207:       ELSE IF( LDAB.LT.KD+1 ) THEN
  208:          INFO = -7
  209:       END IF
  210:       IF( INFO.NE.0 ) THEN
  211:          CALL XERBLA( 'DTBCON', -INFO )
  212:          RETURN
  213:       END IF
  214: *
  215: *     Quick return if possible
  216: *
  217:       IF( N.EQ.0 ) THEN
  218:          RCOND = ONE
  219:          RETURN
  220:       END IF
  221: *
  222:       RCOND = ZERO
  223:       SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( 1, N ) )
  224: *
  225: *     Compute the norm of the triangular matrix A.
  226: *
  227:       ANORM = DLANTB( NORM, UPLO, DIAG, N, KD, AB, LDAB, WORK )
  228: *
  229: *     Continue only if ANORM > 0.
  230: *
  231:       IF( ANORM.GT.ZERO ) THEN
  232: *
  233: *        Estimate the norm of the inverse of A.
  234: *
  235:          AINVNM = ZERO
  236:          NORMIN = 'N'
  237:          IF( ONENRM ) THEN
  238:             KASE1 = 1
  239:          ELSE
  240:             KASE1 = 2
  241:          END IF
  242:          KASE = 0
  243:    10    CONTINUE
  244:          CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
  245:          IF( KASE.NE.0 ) THEN
  246:             IF( KASE.EQ.KASE1 ) THEN
  247: *
  248: *              Multiply by inv(A).
  249: *
  250:                CALL DLATBS( UPLO, 'No transpose', DIAG, NORMIN, N, KD,
  251:      $                      AB, LDAB, WORK, SCALE, WORK( 2*N+1 ), INFO )
  252:             ELSE
  253: *
  254: *              Multiply by inv(A**T).
  255: *
  256:                CALL DLATBS( UPLO, 'Transpose', DIAG, NORMIN, N, KD, AB,
  257:      $                      LDAB, WORK, SCALE, WORK( 2*N+1 ), INFO )
  258:             END IF
  259:             NORMIN = 'Y'
  260: *
  261: *           Multiply by 1/SCALE if doing so will not cause overflow.
  262: *
  263:             IF( SCALE.NE.ONE ) THEN
  264:                IX = IDAMAX( N, WORK, 1 )
  265:                XNORM = ABS( WORK( IX ) )
  266:                IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
  267:      $            GO TO 20
  268:                CALL DRSCL( N, SCALE, WORK, 1 )
  269:             END IF
  270:             GO TO 10
  271:          END IF
  272: *
  273: *        Compute the estimate of the reciprocal condition number.
  274: *
  275:          IF( AINVNM.NE.ZERO )
  276:      $      RCOND = ( ONE / ANORM ) / AINVNM
  277:       END IF
  278: *
  279:    20 CONTINUE
  280:       RETURN
  281: *
  282: *     End of DTBCON
  283: *
  284:       END

CVSweb interface <joel.bertrand@systella.fr>