001:       SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
002:      $                   WORK, IWORK, 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 DLACN2 in place of DLACON, 5 Feb 03, SJH.
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          NORM
013:       INTEGER            INFO, KL, KU, LDAB, N
014:       DOUBLE PRECISION   ANORM, RCOND
015: *     ..
016: *     .. Array Arguments ..
017:       INTEGER            IPIV( * ), IWORK( * )
018:       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  DGBCON estimates the reciprocal of the condition number of a real
025: *  general band matrix A, in either the 1-norm or the infinity-norm,
026: *  using the LU factorization computed by DGBTRF.
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) DOUBLE PRECISION array, dimension (LDAB,N)
051: *          Details of the LU factorization of the band matrix A, as
052: *          computed by DGBTRF.  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) DOUBLE PRECISION array, dimension (3*N)
073: *
074: *  IWORK   (workspace) INTEGER 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, T
091: *     ..
092: *     .. Local Arrays ..
093:       INTEGER            ISAVE( 3 )
094: *     ..
095: *     .. External Functions ..
096:       LOGICAL            LSAME
097:       INTEGER            IDAMAX
098:       DOUBLE PRECISION   DDOT, DLAMCH
099:       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
100: *     ..
101: *     .. External Subroutines ..
102:       EXTERNAL           DAXPY, DLACN2, DLATBS, DRSCL, XERBLA
103: *     ..
104: *     .. Intrinsic Functions ..
105:       INTRINSIC          ABS, MIN
106: *     ..
107: *     .. Executable Statements ..
108: *
109: *     Test the input parameters.
110: *
111:       INFO = 0
112:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
113:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
114:          INFO = -1
115:       ELSE IF( N.LT.0 ) THEN
116:          INFO = -2
117:       ELSE IF( KL.LT.0 ) THEN
118:          INFO = -3
119:       ELSE IF( KU.LT.0 ) THEN
120:          INFO = -4
121:       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
122:          INFO = -6
123:       ELSE IF( ANORM.LT.ZERO ) THEN
124:          INFO = -8
125:       END IF
126:       IF( INFO.NE.0 ) THEN
127:          CALL XERBLA( 'DGBCON', -INFO )
128:          RETURN
129:       END IF
130: *
131: *     Quick return if possible
132: *
133:       RCOND = ZERO
134:       IF( N.EQ.0 ) THEN
135:          RCOND = ONE
136:          RETURN
137:       ELSE IF( ANORM.EQ.ZERO ) THEN
138:          RETURN
139:       END IF
140: *
141:       SMLNUM = DLAMCH( 'Safe minimum' )
142: *
143: *     Estimate the norm of inv(A).
144: *
145:       AINVNM = ZERO
146:       NORMIN = 'N'
147:       IF( ONENRM ) THEN
148:          KASE1 = 1
149:       ELSE
150:          KASE1 = 2
151:       END IF
152:       KD = KL + KU + 1
153:       LNOTI = KL.GT.0
154:       KASE = 0
155:    10 CONTINUE
156:       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
157:       IF( KASE.NE.0 ) THEN
158:          IF( KASE.EQ.KASE1 ) THEN
159: *
160: *           Multiply by inv(L).
161: *
162:             IF( LNOTI ) THEN
163:                DO 20 J = 1, N - 1
164:                   LM = MIN( KL, N-J )
165:                   JP = IPIV( J )
166:                   T = WORK( JP )
167:                   IF( JP.NE.J ) THEN
168:                      WORK( JP ) = WORK( J )
169:                      WORK( J ) = T
170:                   END IF
171:                   CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
172:    20          CONTINUE
173:             END IF
174: *
175: *           Multiply by inv(U).
176: *
177:             CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
178:      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
179:      $                   INFO )
180:          ELSE
181: *
182: *           Multiply by inv(U').
183: *
184:             CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
185:      $                   KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
186:      $                   INFO )
187: *
188: *           Multiply by inv(L').
189: *
190:             IF( LNOTI ) THEN
191:                DO 30 J = N - 1, 1, -1
192:                   LM = MIN( KL, N-J )
193:                   WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
194:      $                        WORK( J+1 ), 1 )
195:                   JP = IPIV( J )
196:                   IF( JP.NE.J ) THEN
197:                      T = WORK( JP )
198:                      WORK( JP ) = WORK( J )
199:                      WORK( J ) = T
200:                   END IF
201:    30          CONTINUE
202:             END IF
203:          END IF
204: *
205: *        Divide X by 1/SCALE if doing so will not cause overflow.
206: *
207:          NORMIN = 'Y'
208:          IF( SCALE.NE.ONE ) THEN
209:             IX = IDAMAX( N, WORK, 1 )
210:             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
211:      $         GO TO 40
212:             CALL DRSCL( N, SCALE, WORK, 1 )
213:          END IF
214:          GO TO 10
215:       END IF
216: *
217: *     Compute the estimate of the reciprocal condition number.
218: *
219:       IF( AINVNM.NE.ZERO )
220:      $   RCOND = ( ONE / AINVNM ) / ANORM
221: *
222:    40 CONTINUE
223:       RETURN
224: *
225: *     End of DGBCON
226: *
227:       END
228: