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