File:  [local] / rpl / lapack / lapack / zla_gbrcond_x.f
Revision 1.6: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 20:43:13 2011 UTC (12 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Mise à jour de Lapack.

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

CVSweb interface <joel.bertrand@systella.fr>