Annotation of rpl/lapack/lapack/dgbcon.f, revision 1.11
1.9 bertrand 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: *> \date November 2011
142: *
143: *> \ingroup doubleGBcomputational
144: *
145: * =====================================================================
1.1 bertrand 146: SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
147: $ WORK, IWORK, INFO )
148: *
1.9 bertrand 149: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 150: * -- LAPACK is a software package provided by Univ. of Tennessee, --
151: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 bertrand 152: * November 2011
1.1 bertrand 153: *
154: * .. Scalar Arguments ..
155: CHARACTER NORM
156: INTEGER INFO, KL, KU, LDAB, N
157: DOUBLE PRECISION ANORM, RCOND
158: * ..
159: * .. Array Arguments ..
160: INTEGER IPIV( * ), IWORK( * )
161: DOUBLE PRECISION AB( LDAB, * ), WORK( * )
162: * ..
163: *
164: * =====================================================================
165: *
166: * .. Parameters ..
167: DOUBLE PRECISION ONE, ZERO
168: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
169: * ..
170: * .. Local Scalars ..
171: LOGICAL LNOTI, ONENRM
172: CHARACTER NORMIN
173: INTEGER IX, J, JP, KASE, KASE1, KD, LM
174: DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T
175: * ..
176: * .. Local Arrays ..
177: INTEGER ISAVE( 3 )
178: * ..
179: * .. External Functions ..
180: LOGICAL LSAME
181: INTEGER IDAMAX
182: DOUBLE PRECISION DDOT, DLAMCH
183: EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
184: * ..
185: * .. External Subroutines ..
186: EXTERNAL DAXPY, DLACN2, DLATBS, DRSCL, XERBLA
187: * ..
188: * .. Intrinsic Functions ..
189: INTRINSIC ABS, MIN
190: * ..
191: * .. Executable Statements ..
192: *
193: * Test the input parameters.
194: *
195: INFO = 0
196: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
197: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
198: INFO = -1
199: ELSE IF( N.LT.0 ) THEN
200: INFO = -2
201: ELSE IF( KL.LT.0 ) THEN
202: INFO = -3
203: ELSE IF( KU.LT.0 ) THEN
204: INFO = -4
205: ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
206: INFO = -6
207: ELSE IF( ANORM.LT.ZERO ) THEN
208: INFO = -8
209: END IF
210: IF( INFO.NE.0 ) THEN
211: CALL XERBLA( 'DGBCON', -INFO )
212: RETURN
213: END IF
214: *
215: * Quick return if possible
216: *
217: RCOND = ZERO
218: IF( N.EQ.0 ) THEN
219: RCOND = ONE
220: RETURN
221: ELSE IF( ANORM.EQ.ZERO ) THEN
222: RETURN
223: END IF
224: *
225: SMLNUM = DLAMCH( 'Safe minimum' )
226: *
227: * Estimate the norm of inv(A).
228: *
229: AINVNM = ZERO
230: NORMIN = 'N'
231: IF( ONENRM ) THEN
232: KASE1 = 1
233: ELSE
234: KASE1 = 2
235: END IF
236: KD = KL + KU + 1
237: LNOTI = KL.GT.0
238: KASE = 0
239: 10 CONTINUE
240: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
241: IF( KASE.NE.0 ) THEN
242: IF( KASE.EQ.KASE1 ) THEN
243: *
244: * Multiply by inv(L).
245: *
246: IF( LNOTI ) THEN
247: DO 20 J = 1, N - 1
248: LM = MIN( KL, N-J )
249: JP = IPIV( J )
250: T = WORK( JP )
251: IF( JP.NE.J ) THEN
252: WORK( JP ) = WORK( J )
253: WORK( J ) = T
254: END IF
255: CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
256: 20 CONTINUE
257: END IF
258: *
259: * Multiply by inv(U).
260: *
261: CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
262: $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
263: $ INFO )
264: ELSE
265: *
1.8 bertrand 266: * Multiply by inv(U**T).
1.1 bertrand 267: *
268: CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
269: $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
270: $ INFO )
271: *
1.8 bertrand 272: * Multiply by inv(L**T).
1.1 bertrand 273: *
274: IF( LNOTI ) THEN
275: DO 30 J = N - 1, 1, -1
276: LM = MIN( KL, N-J )
277: WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
278: $ WORK( J+1 ), 1 )
279: JP = IPIV( J )
280: IF( JP.NE.J ) THEN
281: T = WORK( JP )
282: WORK( JP ) = WORK( J )
283: WORK( J ) = T
284: END IF
285: 30 CONTINUE
286: END IF
287: END IF
288: *
289: * Divide X by 1/SCALE if doing so will not cause overflow.
290: *
291: NORMIN = 'Y'
292: IF( SCALE.NE.ONE ) THEN
293: IX = IDAMAX( N, WORK, 1 )
294: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
295: $ GO TO 40
296: CALL DRSCL( N, SCALE, WORK, 1 )
297: END IF
298: GO TO 10
299: END IF
300: *
301: * Compute the estimate of the reciprocal condition number.
302: *
303: IF( AINVNM.NE.ZERO )
304: $ RCOND = ( ONE / AINVNM ) / ANORM
305: *
306: 40 CONTINUE
307: RETURN
308: *
309: * End of DGBCON
310: *
311: 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