Annotation of rpl/lapack/lapack/zla_gbrcond_x.f, revision 1.17

1.9       bertrand    1: *> \brief \b ZLA_GBRCOND_X computes the infinity norm condition number of op(A)*diag(x) for general banded matrices.
1.6       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.13      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.6       bertrand    7: *
                      8: *> \htmlonly
1.13      bertrand    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">
1.6       bertrand   15: *> [TXT]</a>
1.13      bertrand   16: *> \endhtmlonly
1.6       bertrand   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 )
1.13      bertrand   24: *
1.6       bertrand   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( * )
1.13      bertrand   34: *
                     35: *
1.6       bertrand   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: *>
1.16      bertrand  129: *> \param[out] WORK
1.6       bertrand  130: *> \verbatim
                    131: *>          WORK is COMPLEX*16 array, dimension (2*N).
                    132: *>     Workspace.
                    133: *> \endverbatim
                    134: *>
1.16      bertrand  135: *> \param[out] RWORK
1.6       bertrand  136: *> \verbatim
                    137: *>          RWORK is DOUBLE PRECISION array, dimension (N).
                    138: *>     Workspace.
                    139: *> \endverbatim
                    140: *
                    141: *  Authors:
                    142: *  ========
                    143: *
1.13      bertrand  144: *> \author Univ. of Tennessee
                    145: *> \author Univ. of California Berkeley
                    146: *> \author Univ. of Colorado Denver
                    147: *> \author NAG Ltd.
1.6       bertrand  148: *
                    149: *> \ingroup complex16GBcomputational
                    150: *
                    151: *  =====================================================================
1.1       bertrand  152:       DOUBLE PRECISION FUNCTION ZLA_GBRCOND_X( TRANS, N, KL, KU, AB,
                    153:      $                                         LDAB, AFB, LDAFB, IPIV,
                    154:      $                                         X, INFO, WORK, RWORK )
                    155: *
1.17    ! bertrand  156: *  -- LAPACK computational routine --
1.6       bertrand  157: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    158: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.1       bertrand  159: *
                    160: *     .. Scalar Arguments ..
                    161:       CHARACTER          TRANS
                    162:       INTEGER            N, KL, KU, KD, KE, LDAB, LDAFB, INFO
                    163: *     ..
                    164: *     .. Array Arguments ..
                    165:       INTEGER            IPIV( * )
                    166:       COMPLEX*16         AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ),
                    167:      $                   X( * )
                    168:       DOUBLE PRECISION   RWORK( * )
                    169: *
                    170: *
                    171: *  =====================================================================
                    172: *
                    173: *     .. Local Scalars ..
                    174:       LOGICAL            NOTRANS
                    175:       INTEGER            KASE, I, J
                    176:       DOUBLE PRECISION   AINVNM, ANORM, TMP
                    177:       COMPLEX*16         ZDUM
                    178: *     ..
                    179: *     .. Local Arrays ..
                    180:       INTEGER            ISAVE( 3 )
                    181: *     ..
                    182: *     .. External Functions ..
                    183:       LOGICAL            LSAME
                    184:       EXTERNAL           LSAME
                    185: *     ..
                    186: *     .. External Subroutines ..
                    187:       EXTERNAL           ZLACN2, ZGBTRS, XERBLA
                    188: *     ..
                    189: *     .. Intrinsic Functions ..
                    190:       INTRINSIC          ABS, MAX
                    191: *     ..
                    192: *     .. Statement Functions ..
                    193:       DOUBLE PRECISION   CABS1
                    194: *     ..
                    195: *     .. Statement Function Definitions ..
                    196:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
                    197: *     ..
                    198: *     .. Executable Statements ..
                    199: *
                    200:       ZLA_GBRCOND_X = 0.0D+0
                    201: *
                    202:       INFO = 0
                    203:       NOTRANS = LSAME( TRANS, 'N' )
                    204:       IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T') .AND. .NOT.
                    205:      $     LSAME( TRANS, 'C' ) ) THEN
                    206:          INFO = -1
                    207:       ELSE IF( N.LT.0 ) THEN
                    208:          INFO = -2
                    209:       ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
                    210:          INFO = -3
                    211:       ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
                    212:          INFO = -4
                    213:       ELSE IF( LDAB.LT.KL+KU+1 ) THEN
                    214:          INFO = -6
                    215:       ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
                    216:          INFO = -8
                    217:       END IF
                    218:       IF( INFO.NE.0 ) THEN
                    219:          CALL XERBLA( 'ZLA_GBRCOND_X', -INFO )
                    220:          RETURN
                    221:       END IF
                    222: *
                    223: *     Compute norm of op(A)*op2(C).
                    224: *
                    225:       KD = KU + 1
                    226:       KE = KL + 1
                    227:       ANORM = 0.0D+0
                    228:       IF ( NOTRANS ) THEN
                    229:          DO I = 1, N
                    230:             TMP = 0.0D+0
                    231:             DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
                    232:                TMP = TMP + CABS1( AB( KD+I-J, J) * X( J ) )
                    233:             END DO
                    234:             RWORK( I ) = TMP
                    235:             ANORM = MAX( ANORM, TMP )
                    236:          END DO
                    237:       ELSE
                    238:          DO I = 1, N
                    239:             TMP = 0.0D+0
                    240:             DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
                    241:                TMP = TMP + CABS1( AB( KE-I+J, I ) * X( J ) )
                    242:             END DO
                    243:             RWORK( I ) = TMP
                    244:             ANORM = MAX( ANORM, TMP )
                    245:          END DO
                    246:       END IF
                    247: *
                    248: *     Quick return if possible.
                    249: *
                    250:       IF( N.EQ.0 ) THEN
                    251:          ZLA_GBRCOND_X = 1.0D+0
                    252:          RETURN
                    253:       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
                    254:          RETURN
                    255:       END IF
                    256: *
                    257: *     Estimate the norm of inv(op(A)).
                    258: *
                    259:       AINVNM = 0.0D+0
                    260: *
                    261:       KASE = 0
                    262:    10 CONTINUE
                    263:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
                    264:       IF( KASE.NE.0 ) THEN
                    265:          IF( KASE.EQ.2 ) THEN
                    266: *
                    267: *           Multiply by R.
                    268: *
                    269:             DO I = 1, N
                    270:                WORK( I ) = WORK( I ) * RWORK( I )
                    271:             END DO
                    272: *
                    273:             IF ( NOTRANS ) THEN
                    274:                CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
                    275:      $              IPIV, WORK, N, INFO )
                    276:             ELSE
                    277:                CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
                    278:      $              LDAFB, IPIV, WORK, N, INFO )
                    279:             ENDIF
                    280: *
                    281: *           Multiply by inv(X).
                    282: *
                    283:             DO I = 1, N
                    284:                WORK( I ) = WORK( I ) / X( I )
                    285:             END DO
                    286:          ELSE
                    287: *
1.5       bertrand  288: *           Multiply by inv(X**H).
1.1       bertrand  289: *
                    290:             DO I = 1, N
                    291:                WORK( I ) = WORK( I ) / X( I )
                    292:             END DO
                    293: *
                    294:             IF ( NOTRANS ) THEN
                    295:                CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
                    296:      $              LDAFB, IPIV, WORK, N, INFO )
                    297:             ELSE
                    298:                CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
                    299:      $              IPIV, WORK, N, INFO )
                    300:             END IF
                    301: *
                    302: *           Multiply by R.
                    303: *
                    304:             DO I = 1, N
                    305:                WORK( I ) = WORK( I ) * RWORK( I )
                    306:             END DO
                    307:          END IF
                    308:          GO TO 10
                    309:       END IF
                    310: *
                    311: *     Compute the estimate of the reciprocal condition number.
                    312: *
                    313:       IF( AINVNM .NE. 0.0D+0 )
                    314:      $   ZLA_GBRCOND_X = 1.0D+0 / AINVNM
                    315: *
                    316:       RETURN
                    317: *
1.17    ! bertrand  318: *     End of ZLA_GBRCOND_X
        !           319: *
1.1       bertrand  320:       END

CVSweb interface <joel.bertrand@systella.fr>