1: *> \brief \b DGBCON
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGBCON + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbcon.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbcon.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbcon.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
22: * WORK, IWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER NORM
26: * INTEGER INFO, KL, KU, LDAB, N
27: * DOUBLE PRECISION ANORM, RCOND
28: * ..
29: * .. Array Arguments ..
30: * INTEGER IPIV( * ), IWORK( * )
31: * DOUBLE PRECISION AB( LDAB, * ), WORK( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DGBCON estimates the reciprocal of the condition number of a real
41: *> general band matrix A, in either the 1-norm or the infinity-norm,
42: *> using the LU factorization computed by DGBTRF.
43: *>
44: *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
45: *> condition number is computed as
46: *> RCOND = 1 / ( norm(A) * norm(inv(A)) ).
47: *> \endverbatim
48: *
49: * Arguments:
50: * ==========
51: *
52: *> \param[in] NORM
53: *> \verbatim
54: *> NORM is CHARACTER*1
55: *> Specifies whether the 1-norm condition number or the
56: *> infinity-norm condition number is required:
57: *> = '1' or 'O': 1-norm;
58: *> = 'I': Infinity-norm.
59: *> \endverbatim
60: *>
61: *> \param[in] N
62: *> \verbatim
63: *> N is INTEGER
64: *> The order of the matrix A. N >= 0.
65: *> \endverbatim
66: *>
67: *> \param[in] KL
68: *> \verbatim
69: *> KL is INTEGER
70: *> The number of subdiagonals within the band of A. KL >= 0.
71: *> \endverbatim
72: *>
73: *> \param[in] KU
74: *> \verbatim
75: *> KU is INTEGER
76: *> The number of superdiagonals within the band of A. KU >= 0.
77: *> \endverbatim
78: *>
79: *> \param[in] AB
80: *> \verbatim
81: *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
82: *> Details of the LU factorization of the band matrix A, as
83: *> computed by DGBTRF. U is stored as an upper triangular band
84: *> matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
85: *> the multipliers used during the factorization are stored in
86: *> rows KL+KU+2 to 2*KL+KU+1.
87: *> \endverbatim
88: *>
89: *> \param[in] LDAB
90: *> \verbatim
91: *> LDAB is INTEGER
92: *> The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
93: *> \endverbatim
94: *>
95: *> \param[in] IPIV
96: *> \verbatim
97: *> IPIV is INTEGER array, dimension (N)
98: *> The pivot indices; for 1 <= i <= N, row i of the matrix was
99: *> interchanged with row IPIV(i).
100: *> \endverbatim
101: *>
102: *> \param[in] ANORM
103: *> \verbatim
104: *> ANORM is DOUBLE PRECISION
105: *> If NORM = '1' or 'O', the 1-norm of the original matrix A.
106: *> If NORM = 'I', the infinity-norm of the original matrix A.
107: *> \endverbatim
108: *>
109: *> \param[out] RCOND
110: *> \verbatim
111: *> RCOND is DOUBLE PRECISION
112: *> The reciprocal of the condition number of the matrix A,
113: *> computed as RCOND = 1/(norm(A) * norm(inv(A))).
114: *> \endverbatim
115: *>
116: *> \param[out] WORK
117: *> \verbatim
118: *> WORK is DOUBLE PRECISION array, dimension (3*N)
119: *> \endverbatim
120: *>
121: *> \param[out] IWORK
122: *> \verbatim
123: *> IWORK is INTEGER array, dimension (N)
124: *> \endverbatim
125: *>
126: *> \param[out] INFO
127: *> \verbatim
128: *> INFO is INTEGER
129: *> = 0: successful exit
130: *> < 0: if INFO = -i, the i-th argument had an illegal value
131: *> \endverbatim
132: *
133: * Authors:
134: * ========
135: *
136: *> \author Univ. of Tennessee
137: *> \author Univ. of California Berkeley
138: *> \author Univ. of Colorado Denver
139: *> \author NAG Ltd.
140: *
141: *> \ingroup doubleGBcomputational
142: *
143: * =====================================================================
144: SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
145: $ WORK, IWORK, INFO )
146: *
147: * -- LAPACK computational routine --
148: * -- LAPACK is a software package provided by Univ. of Tennessee, --
149: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150: *
151: * .. Scalar Arguments ..
152: CHARACTER NORM
153: INTEGER INFO, KL, KU, LDAB, N
154: DOUBLE PRECISION ANORM, RCOND
155: * ..
156: * .. Array Arguments ..
157: INTEGER IPIV( * ), IWORK( * )
158: DOUBLE PRECISION AB( LDAB, * ), WORK( * )
159: * ..
160: *
161: * =====================================================================
162: *
163: * .. Parameters ..
164: DOUBLE PRECISION ONE, ZERO
165: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
166: * ..
167: * .. Local Scalars ..
168: LOGICAL LNOTI, ONENRM
169: CHARACTER NORMIN
170: INTEGER IX, J, JP, KASE, KASE1, KD, LM
171: DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T
172: * ..
173: * .. Local Arrays ..
174: INTEGER ISAVE( 3 )
175: * ..
176: * .. External Functions ..
177: LOGICAL LSAME
178: INTEGER IDAMAX
179: DOUBLE PRECISION DDOT, DLAMCH
180: EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
181: * ..
182: * .. External Subroutines ..
183: EXTERNAL DAXPY, DLACN2, DLATBS, DRSCL, XERBLA
184: * ..
185: * .. Intrinsic Functions ..
186: INTRINSIC ABS, MIN
187: * ..
188: * .. Executable Statements ..
189: *
190: * Test the input parameters.
191: *
192: INFO = 0
193: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
194: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
195: INFO = -1
196: ELSE IF( N.LT.0 ) THEN
197: INFO = -2
198: ELSE IF( KL.LT.0 ) THEN
199: INFO = -3
200: ELSE IF( KU.LT.0 ) THEN
201: INFO = -4
202: ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
203: INFO = -6
204: ELSE IF( ANORM.LT.ZERO ) THEN
205: INFO = -8
206: END IF
207: IF( INFO.NE.0 ) THEN
208: CALL XERBLA( 'DGBCON', -INFO )
209: RETURN
210: END IF
211: *
212: * Quick return if possible
213: *
214: RCOND = ZERO
215: IF( N.EQ.0 ) THEN
216: RCOND = ONE
217: RETURN
218: ELSE IF( ANORM.EQ.ZERO ) THEN
219: RETURN
220: END IF
221: *
222: SMLNUM = DLAMCH( 'Safe minimum' )
223: *
224: * Estimate the norm of inv(A).
225: *
226: AINVNM = ZERO
227: NORMIN = 'N'
228: IF( ONENRM ) THEN
229: KASE1 = 1
230: ELSE
231: KASE1 = 2
232: END IF
233: KD = KL + KU + 1
234: LNOTI = KL.GT.0
235: KASE = 0
236: 10 CONTINUE
237: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
238: IF( KASE.NE.0 ) THEN
239: IF( KASE.EQ.KASE1 ) THEN
240: *
241: * Multiply by inv(L).
242: *
243: IF( LNOTI ) THEN
244: DO 20 J = 1, N - 1
245: LM = MIN( KL, N-J )
246: JP = IPIV( J )
247: T = WORK( JP )
248: IF( JP.NE.J ) THEN
249: WORK( JP ) = WORK( J )
250: WORK( J ) = T
251: END IF
252: CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
253: 20 CONTINUE
254: END IF
255: *
256: * Multiply by inv(U).
257: *
258: CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
259: $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
260: $ INFO )
261: ELSE
262: *
263: * Multiply by inv(U**T).
264: *
265: CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
266: $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
267: $ INFO )
268: *
269: * Multiply by inv(L**T).
270: *
271: IF( LNOTI ) THEN
272: DO 30 J = N - 1, 1, -1
273: LM = MIN( KL, N-J )
274: WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
275: $ WORK( J+1 ), 1 )
276: JP = IPIV( J )
277: IF( JP.NE.J ) THEN
278: T = WORK( JP )
279: WORK( JP ) = WORK( J )
280: WORK( J ) = T
281: END IF
282: 30 CONTINUE
283: END IF
284: END IF
285: *
286: * Divide X by 1/SCALE if doing so will not cause overflow.
287: *
288: NORMIN = 'Y'
289: IF( SCALE.NE.ONE ) THEN
290: IX = IDAMAX( N, WORK, 1 )
291: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
292: $ GO TO 40
293: CALL DRSCL( N, SCALE, WORK, 1 )
294: END IF
295: GO TO 10
296: END IF
297: *
298: * Compute the estimate of the reciprocal condition number.
299: *
300: IF( AINVNM.NE.ZERO )
301: $ RCOND = ( ONE / AINVNM ) / ANORM
302: *
303: 40 CONTINUE
304: RETURN
305: *
306: * End of DGBCON
307: *
308: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgbcon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack