LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgbcon.f
Go to the documentation of this file.
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 gbcon
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
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dgbcon(norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, iwork, info)
DGBCON
Definition dgbcon.f:146
subroutine dlacn2(n, v, x, isgn, est, kase, isave)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition dlacn2.f:136
subroutine dlatbs(uplo, trans, diag, normin, n, kd, ab, ldab, x, scale, cnorm, info)
DLATBS solves a triangular banded system of equations.
Definition dlatbs.f:242
subroutine drscl(n, sa, sx, incx)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition drscl.f:84