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

1.6     ! bertrand    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: *  =====================================================================
1.1       bertrand  154:       DOUBLE PRECISION FUNCTION ZLA_GBRCOND_X( TRANS, N, KL, KU, AB,
                    155:      $                                         LDAB, AFB, LDAFB, IPIV,
                    156:      $                                         X, INFO, WORK, RWORK )
                    157: *
1.6     ! bertrand  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
1.1       bertrand  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: *
1.5       bertrand  291: *           Multiply by inv(X**H).
1.1       bertrand  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>