001:       DOUBLE PRECISION FUNCTION ZLA_GBRCOND_C( TRANS, N, KL, KU, AB, 
002:      $                             LDAB, AFB, LDAFB, IPIV, C, CAPPLY, 
003:      $                             INFO, WORK, RWORK )
004: *
005: *     -- LAPACK routine (version 3.2)                                 --
006: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
007: *     -- Jason Riedy of Univ. of California Berkeley.                 --
008: *     -- November 2008                                                --
009: *
010: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
011: *     -- Univ. of California Berkeley and NAG Ltd.                    --
012: *
013:       IMPLICIT NONE
014: *     ..
015: *     .. Scalar Arguments ..
016:       CHARACTER          TRANS
017:       LOGICAL            CAPPLY
018:       INTEGER            N, KL, KU, KD, LDAB, LDAFB, INFO
019: *     ..
020: *     .. Array Arguments ..
021:       INTEGER            IPIV( * )
022:       COMPLEX*16         AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
023:       DOUBLE PRECISION   C( * ), RWORK( * )
024: *
025: *     ZLA_GBRCOND_C Computes the infinity norm condition number of
026: *     op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
027: *     WORK is a COMPLEX*16 workspace of size 2*N, and
028: *     RWORK is a DOUBLE PRECISION workspace of size 3*N.
029: *     ..
030: *     .. Local Scalars ..
031:       LOGICAL            NOTRANS
032:       INTEGER            KASE, I, J
033:       DOUBLE PRECISION   AINVNM, ANORM, TMP
034:       COMPLEX*16         ZDUM
035: *     ..
036: *     .. Local Arrays ..
037:       INTEGER            ISAVE( 3 )
038: *     ..
039: *     .. External Functions ..
040:       LOGICAL            LSAME
041:       EXTERNAL           LSAME
042: *     ..
043: *     .. External Subroutines ..
044:       EXTERNAL           ZLACN2, ZGBTRS, XERBLA
045: *     ..
046: *     .. Intrinsic Functions ..
047:       INTRINSIC          ABS, MAX
048: *     ..
049: *     .. Statement Functions ..
050:       DOUBLE PRECISION   CABS1
051: *     ..
052: *     .. Statement Function Definitions ..
053:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
054: *     ..
055: *     .. Executable Statements ..
056:       ZLA_GBRCOND_C = 0.0D+0
057: *
058:       INFO = 0
059:       NOTRANS = LSAME( TRANS, 'N' )
060:       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
061:      $     LSAME( TRANS, 'C' ) ) THEN
062:       ELSE IF( N.LT.0 ) THEN
063:          INFO = -2
064:       END IF
065:       IF( INFO.NE.0 ) THEN
066:          CALL XERBLA( 'ZLA_GBRCOND_C', -INFO )
067:          RETURN
068:       END IF
069: *
070: *     Compute norm of op(A)*op2(C).
071: *
072:       ANORM = 0.0D+0
073:       KD = KU + 1
074:       IF ( NOTRANS ) THEN
075:          DO I = 1, N
076:             TMP = 0.0D+0
077:             IF ( CAPPLY ) THEN
078:                DO J = 1, N
079:                   IF ( I.GE.MAX( 1, J-KU )
080:      $                 .AND. I.LE.MIN( N, J+KL ) ) THEN
081:                      TMP = TMP + CABS1(AB( KD+I-J, J ) ) / C( J )
082:                   END IF
083:                END DO
084:             ELSE
085:                DO J = 1, N
086:                   IF ( I.GE.MAX( 1, J-KU )
087:      $                 .AND. I.LE.MIN( N, J+KL ) ) THEN
088:                      TMP = TMP + CABS1( AB( KD+I-J, J ) )
089:                   END IF
090:                END DO
091:             END IF
092:             RWORK( 2*N+I ) = TMP
093:             ANORM = MAX( ANORM, TMP )
094:          END DO
095:       ELSE
096:          DO I = 1, N
097:             TMP = 0.0D+0
098:             IF ( CAPPLY ) THEN
099:                DO J = 1, N
100:                   IF ( I.GE.MAX( 1, J-KU )
101:      $                 .AND. I.LE.MIN( N, J+KL ) ) THEN
102:                      TMP = TMP + CABS1( AB( J, KD+I-J ) ) / C( J )
103:                   END IF
104:                END DO
105:             ELSE
106:                DO J = 1, N
107:                   IF ( I.GE.MAX( 1, J-KU )
108:      $                 .AND. I.LE.MIN( N, J+KL ) ) THEN
109:                      TMP = TMP + CABS1( AB( J, KD+I-J ) )
110:                   END IF
111:                END DO
112:             END IF
113:             RWORK( 2*N+I ) = TMP
114:             ANORM = MAX( ANORM, TMP )
115:          END DO
116:       END IF
117: *
118: *     Quick return if possible.
119: *
120:       IF( N.EQ.0 ) THEN
121:          ZLA_GBRCOND_C = 1.0D+0
122:          RETURN
123:       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
124:          RETURN
125:       END IF
126: *
127: *     Estimate the norm of inv(op(A)).
128: *
129:       AINVNM = 0.0D+0
130: *
131:       KASE = 0
132:    10 CONTINUE
133:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
134:       IF( KASE.NE.0 ) THEN
135:          IF( KASE.EQ.2 ) THEN
136: *
137: *           Multiply by R.
138: *
139:             DO I = 1, N
140:                WORK( I ) = WORK( I ) * RWORK( 2*N+I )
141:             END DO
142: *
143:             IF ( NOTRANS ) THEN
144:                CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
145:      $              IPIV, WORK, N, INFO )
146:             ELSE
147:                CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
148:      $              LDAFB, IPIV, WORK, N, INFO )
149:             ENDIF
150: *
151: *           Multiply by inv(C).
152: *
153:             IF ( CAPPLY ) THEN
154:                DO I = 1, N
155:                   WORK( I ) = WORK( I ) * C( I )
156:                END DO
157:             END IF
158:          ELSE
159: *
160: *           Multiply by inv(C').
161: *
162:             IF ( CAPPLY ) THEN
163:                DO I = 1, N
164:                   WORK( I ) = WORK( I ) * C( I )
165:                END DO
166:             END IF
167: *
168:             IF ( NOTRANS ) THEN
169:                CALL ZGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
170:      $              LDAFB, IPIV,  WORK, N, INFO )
171:             ELSE
172:                CALL ZGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
173:      $              IPIV, WORK, N, INFO )
174:             END IF
175: *
176: *           Multiply by R.
177: *
178:             DO I = 1, N
179:                WORK( I ) = WORK( I ) * RWORK( 2*N+I )
180:             END DO
181:          END IF
182:          GO TO 10
183:       END IF
184: *
185: *     Compute the estimate of the reciprocal condition number.
186: *
187:       IF( AINVNM .NE. 0.0D+0 )
188:      $   ZLA_GBRCOND_C = 1.0D+0 / AINVNM
189: *
190:       RETURN
191: *
192:       END
193: