001:       SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
002:      $                   WORK, RWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          NORM
012:       INTEGER            INFO, KL, KU, LDAB, N
013:       DOUBLE PRECISION   ANORM, RCOND
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            IPIV( * )
017:       DOUBLE PRECISION   RWORK( * )
018:       COMPLEX*16         AB( LDAB, * ), WORK( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  ZGBCON estimates the reciprocal of the condition number of a complex
025: *  general band matrix A, in either the 1-norm or the infinity-norm,
026: *  using the LU factorization computed by ZGBTRF.
027: *
028: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
029: *  condition number is computed as
030: *     RCOND = 1 / ( norm(A) * norm(inv(A)) ).
031: *
032: *  Arguments
033: *  =========
034: *
035: *  NORM    (input) CHARACTER*1
036: *          Specifies whether the 1-norm condition number or the
037: *          infinity-norm condition number is required:
038: *          = '1' or 'O':  1-norm;
039: *          = 'I':         Infinity-norm.
040: *
041: *  N       (input) INTEGER
042: *          The order of the matrix A.  N >= 0.
043: *
044: *  KL      (input) INTEGER
045: *          The number of subdiagonals within the band of A.  KL >= 0.
046: *
047: *  KU      (input) INTEGER
048: *          The number of superdiagonals within the band of A.  KU >= 0.
049: *
050: *  AB      (input) COMPLEX*16 array, dimension (LDAB,N)
051: *          Details of the LU factorization of the band matrix A, as
052: *          computed by ZGBTRF.  U is stored as an upper triangular band
053: *          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
054: *          the multipliers used during the factorization are stored in
055: *          rows KL+KU+2 to 2*KL+KU+1.
056: *
057: *  LDAB    (input) INTEGER
058: *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
059: *
060: *  IPIV    (input) INTEGER array, dimension (N)
061: *          The pivot indices; for 1 <= i <= N, row i of the matrix was
062: *          interchanged with row IPIV(i).
063: *
064: *  ANORM   (input) DOUBLE PRECISION
065: *          If NORM = '1' or 'O', the 1-norm of the original matrix A.
066: *          If NORM = 'I', the infinity-norm of the original matrix A.
067: *
068: *  RCOND   (output) DOUBLE PRECISION
069: *          The reciprocal of the condition number of the matrix A,
070: *          computed as RCOND = 1/(norm(A) * norm(inv(A))).
071: *
072: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
073: *
074: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
075: *
076: *  INFO    (output) INTEGER
077: *          = 0:  successful exit
078: *          < 0: if INFO = -i, the i-th argument had an illegal value
079: *
080: *  =====================================================================
081: *
082: *     .. Parameters ..
083:       DOUBLE PRECISION   ONE, ZERO
084:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
085: *     ..
086: *     .. Local Scalars ..
087:       LOGICAL            LNOTI, ONENRM
088:       CHARACTER          NORMIN
089:       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
090:       DOUBLE PRECISION   AINVNM, SCALE, SMLNUM
091:       COMPLEX*16         T, ZDUM
092: *     ..
093: *     .. Local Arrays ..
094:       INTEGER            ISAVE( 3 )
095: *     ..
096: *     .. External Functions ..
097:       LOGICAL            LSAME
098:       INTEGER            IZAMAX
099:       DOUBLE PRECISION   DLAMCH
100:       COMPLEX*16         ZDOTC
101:       EXTERNAL           LSAME, IZAMAX, DLAMCH, ZDOTC
102: *     ..
103: *     .. External Subroutines ..
104:       EXTERNAL           XERBLA, ZAXPY, ZDRSCL, ZLACN2, ZLATBS
105: *     ..
106: *     .. Intrinsic Functions ..
107:       INTRINSIC          ABS, DBLE, DIMAG, MIN
108: *     ..
109: *     .. Statement Functions ..
110:       DOUBLE PRECISION   CABS1
111: *     ..
112: *     .. Statement Function definitions ..
113:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
114: *     ..
115: *     .. Executable Statements ..
116: *
117: *     Test the input parameters.
118: *
119:       INFO = 0
120:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
121:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
122:          INFO = -1
123:       ELSE IF( N.LT.0 ) THEN
124:          INFO = -2
125:       ELSE IF( KL.LT.0 ) THEN
126:          INFO = -3
127:       ELSE IF( KU.LT.0 ) THEN
128:          INFO = -4
129:       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
130:          INFO = -6
131:       ELSE IF( ANORM.LT.ZERO ) THEN
132:          INFO = -8
133:       END IF
134:       IF( INFO.NE.0 ) THEN
135:          CALL XERBLA( 'ZGBCON', -INFO )
136:          RETURN
137:       END IF
138: *
139: *     Quick return if possible
140: *
141:       RCOND = ZERO
142:       IF( N.EQ.0 ) THEN
143:          RCOND = ONE
144:          RETURN
145:       ELSE IF( ANORM.EQ.ZERO ) THEN
146:          RETURN
147:       END IF
148: *
149:       SMLNUM = DLAMCH( 'Safe minimum' )
150: *
151: *     Estimate the norm of inv(A).
152: *
153:       AINVNM = ZERO
154:       NORMIN = 'N'
155:       IF( ONENRM ) THEN
156:          KASE1 = 1
157:       ELSE
158:          KASE1 = 2
159:       END IF
160:       KD = KL + KU + 1
161:       LNOTI = KL.GT.0
162:       KASE = 0
163:    10 CONTINUE
164:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
165:       IF( KASE.NE.0 ) THEN
166:          IF( KASE.EQ.KASE1 ) THEN
167: *
168: *           Multiply by inv(L).
169: *
170:             IF( LNOTI ) THEN
171:                DO 20 J = 1, N - 1
172:                   LM = MIN( KL, N-J )
173:                   JP = IPIV( J )
174:                   T = WORK( JP )
175:                   IF( JP.NE.J ) THEN
176:                      WORK( JP ) = WORK( J )
177:                      WORK( J ) = T
178:                   END IF
179:                   CALL ZAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
180:    20          CONTINUE
181:             END IF
182: *
183: *           Multiply by inv(U).
184: *
185:             CALL ZLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
186:      $                   KL+KU, AB, LDAB, WORK, SCALE, RWORK, INFO )
187:          ELSE
188: *
189: *           Multiply by inv(U').
190: *
191:             CALL ZLATBS( 'Upper', 'Conjugate transpose', 'Non-unit',
192:      $                   NORMIN, N, KL+KU, AB, LDAB, WORK, SCALE, RWORK,
193:      $                   INFO )
194: *
195: *           Multiply by inv(L').
196: *
197:             IF( LNOTI ) THEN
198:                DO 30 J = N - 1, 1, -1
199:                   LM = MIN( KL, N-J )
200:                   WORK( J ) = WORK( J ) - ZDOTC( LM, AB( KD+1, J ), 1,
201:      $                        WORK( J+1 ), 1 )
202:                   JP = IPIV( J )
203:                   IF( JP.NE.J ) THEN
204:                      T = WORK( JP )
205:                      WORK( JP ) = WORK( J )
206:                      WORK( J ) = T
207:                   END IF
208:    30          CONTINUE
209:             END IF
210:          END IF
211: *
212: *        Divide X by 1/SCALE if doing so will not cause overflow.
213: *
214:          NORMIN = 'Y'
215:          IF( SCALE.NE.ONE ) THEN
216:             IX = IZAMAX( N, WORK, 1 )
217:             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
218:      $         GO TO 40
219:             CALL ZDRSCL( N, SCALE, WORK, 1 )
220:          END IF
221:          GO TO 10
222:       END IF
223: *
224: *     Compute the estimate of the reciprocal condition number.
225: *
226:       IF( AINVNM.NE.ZERO )
227:      $   RCOND = ( ONE / AINVNM ) / ANORM
228: *
229:    40 CONTINUE
230:       RETURN
231: *
232: *     End of ZGBCON
233: *
234:       END
235: