LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
dla_gbrcond.f
Go to the documentation of this file.
1*> \brief \b DLA_GBRCOND estimates the Skeel condition number for a general banded matrix.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DLA_GBRCOND + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gbrcond.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gbrcond.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gbrcond.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* DOUBLE PRECISION FUNCTION DLA_GBRCOND( TRANS, N, KL, KU, AB, LDAB,
20* AFB, LDAFB, IPIV, CMODE, C,
21* INFO, WORK, IWORK )
22*
23* .. Scalar Arguments ..
24* CHARACTER TRANS
25* INTEGER N, LDAB, LDAFB, INFO, KL, KU, CMODE
26* ..
27* .. Array Arguments ..
28* INTEGER IWORK( * ), IPIV( * )
29* DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ),
30* $ C( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DLA_GBRCOND Estimates the Skeel condition number of op(A) * op2(C)
40*> where op2 is determined by CMODE as follows
41*> CMODE = 1 op2(C) = C
42*> CMODE = 0 op2(C) = I
43*> CMODE = -1 op2(C) = inv(C)
44*> The Skeel condition number cond(A) = norminf( |inv(A)||A| )
45*> is computed by computing scaling factors R such that
46*> diag(R)*A*op2(C) is row equilibrated and computing the standard
47*> infinity-norm condition number.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] TRANS
54*> \verbatim
55*> TRANS is CHARACTER*1
56*> Specifies the form of the system of equations:
57*> = 'N': A * X = B (No transpose)
58*> = 'T': A**T * X = B (Transpose)
59*> = 'C': A**H * X = B (Conjugate Transpose = Transpose)
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*> N is INTEGER
65*> The number of linear equations, i.e., the order of the
66*> matrix A. N >= 0.
67*> \endverbatim
68*>
69*> \param[in] KL
70*> \verbatim
71*> KL is INTEGER
72*> The number of subdiagonals within the band of A. KL >= 0.
73*> \endverbatim
74*>
75*> \param[in] KU
76*> \verbatim
77*> KU is INTEGER
78*> The number of superdiagonals within the band of A. KU >= 0.
79*> \endverbatim
80*>
81*> \param[in] AB
82*> \verbatim
83*> AB is DOUBLE PRECISION array, dimension (LDAB,N)
84*> On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
85*> The j-th column of A is stored in the j-th column of the
86*> array AB as follows:
87*> AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
88*> \endverbatim
89*>
90*> \param[in] LDAB
91*> \verbatim
92*> LDAB is INTEGER
93*> The leading dimension of the array AB. LDAB >= KL+KU+1.
94*> \endverbatim
95*>
96*> \param[in] AFB
97*> \verbatim
98*> AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
99*> Details of the LU factorization of the band matrix A, as
100*> computed by DGBTRF. U is stored as an upper triangular
101*> band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
102*> and the multipliers used during the factorization are stored
103*> in rows KL+KU+2 to 2*KL+KU+1.
104*> \endverbatim
105*>
106*> \param[in] LDAFB
107*> \verbatim
108*> LDAFB is INTEGER
109*> The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
110*> \endverbatim
111*>
112*> \param[in] IPIV
113*> \verbatim
114*> IPIV is INTEGER array, dimension (N)
115*> The pivot indices from the factorization A = P*L*U
116*> as computed by DGBTRF; row i of the matrix was interchanged
117*> with row IPIV(i).
118*> \endverbatim
119*>
120*> \param[in] CMODE
121*> \verbatim
122*> CMODE is INTEGER
123*> Determines op2(C) in the formula op(A) * op2(C) as follows:
124*> CMODE = 1 op2(C) = C
125*> CMODE = 0 op2(C) = I
126*> CMODE = -1 op2(C) = inv(C)
127*> \endverbatim
128*>
129*> \param[in] C
130*> \verbatim
131*> C is DOUBLE PRECISION array, dimension (N)
132*> The vector C in the formula op(A) * op2(C).
133*> \endverbatim
134*>
135*> \param[out] INFO
136*> \verbatim
137*> INFO is INTEGER
138*> = 0: Successful exit.
139*> i > 0: The ith argument is invalid.
140*> \endverbatim
141*>
142*> \param[out] WORK
143*> \verbatim
144*> WORK is DOUBLE PRECISION array, dimension (5*N).
145*> Workspace.
146*> \endverbatim
147*>
148*> \param[out] IWORK
149*> \verbatim
150*> IWORK is INTEGER array, dimension (N).
151*> Workspace.
152*> \endverbatim
153*
154* Authors:
155* ========
156*
157*> \author Univ. of Tennessee
158*> \author Univ. of California Berkeley
159*> \author Univ. of Colorado Denver
160*> \author NAG Ltd.
161*
162*> \ingroup la_gbrcond
163*
164* =====================================================================
165 DOUBLE PRECISION FUNCTION dla_gbrcond( TRANS, N, KL, KU, AB,
166 $ LDAB,
167 $ AFB, LDAFB, IPIV, CMODE, C,
168 $ INFO, WORK, IWORK )
169*
170* -- LAPACK computational routine --
171* -- LAPACK is a software package provided by Univ. of Tennessee, --
172* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173*
174* .. Scalar Arguments ..
175 CHARACTER trans
176 INTEGER n, ldab, ldafb, info, kl, ku, cmode
177* ..
178* .. Array Arguments ..
179 INTEGER iwork( * ), ipiv( * )
180 DOUBLE PRECISION ab( ldab, * ), afb( ldafb, * ), work( * ),
181 $ c( * )
182* ..
183*
184* =====================================================================
185*
186* .. Local Scalars ..
187 LOGICAL notrans
188 INTEGER kase, i, j, kd, ke
189 DOUBLE PRECISION ainvnm, tmp
190* ..
191* .. Local Arrays ..
192 INTEGER isave( 3 )
193* ..
194* .. External Functions ..
195 LOGICAL lsame
196 EXTERNAL lsame
197* ..
198* .. External Subroutines ..
199 EXTERNAL dlacn2, dgbtrs, xerbla
200* ..
201* .. Intrinsic Functions ..
202 INTRINSIC abs, max
203* ..
204* .. Executable Statements ..
205*
206 dla_gbrcond = 0.0d+0
207*
208 info = 0
209 notrans = lsame( trans, 'N' )
210 IF ( .NOT. notrans .AND. .NOT. lsame(trans, 'T')
211 $ .AND. .NOT. lsame(trans, 'C') ) THEN
212 info = -1
213 ELSE IF( n.LT.0 ) THEN
214 info = -2
215 ELSE IF( kl.LT.0 .OR. kl.GT.n-1 ) THEN
216 info = -3
217 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
218 info = -4
219 ELSE IF( ldab.LT.kl+ku+1 ) THEN
220 info = -6
221 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
222 info = -8
223 END IF
224 IF( info.NE.0 ) THEN
225 CALL xerbla( 'DLA_GBRCOND', -info )
226 RETURN
227 END IF
228 IF( n.EQ.0 ) THEN
229 dla_gbrcond = 1.0d+0
230 RETURN
231 END IF
232*
233* Compute the equilibration matrix R such that
234* inv(R)*A*C has unit 1-norm.
235*
236 kd = ku + 1
237 ke = kl + 1
238 IF ( notrans ) THEN
239 DO i = 1, n
240 tmp = 0.0d+0
241 IF ( cmode .EQ. 1 ) THEN
242 DO j = max( i-kl, 1 ), min( i+ku, n )
243 tmp = tmp + abs( ab( kd+i-j, j ) * c( j ) )
244 END DO
245 ELSE IF ( cmode .EQ. 0 ) THEN
246 DO j = max( i-kl, 1 ), min( i+ku, n )
247 tmp = tmp + abs( ab( kd+i-j, j ) )
248 END DO
249 ELSE
250 DO j = max( i-kl, 1 ), min( i+ku, n )
251 tmp = tmp + abs( ab( kd+i-j, j ) / c( j ) )
252 END DO
253 END IF
254 work( 2*n+i ) = tmp
255 END DO
256 ELSE
257 DO i = 1, n
258 tmp = 0.0d+0
259 IF ( cmode .EQ. 1 ) THEN
260 DO j = max( i-kl, 1 ), min( i+ku, n )
261 tmp = tmp + abs( ab( ke-i+j, i ) * c( j ) )
262 END DO
263 ELSE IF ( cmode .EQ. 0 ) THEN
264 DO j = max( i-kl, 1 ), min( i+ku, n )
265 tmp = tmp + abs( ab( ke-i+j, i ) )
266 END DO
267 ELSE
268 DO j = max( i-kl, 1 ), min( i+ku, n )
269 tmp = tmp + abs( ab( ke-i+j, i ) / c( j ) )
270 END DO
271 END IF
272 work( 2*n+i ) = tmp
273 END DO
274 END IF
275*
276* Estimate the norm of inv(op(A)).
277*
278 ainvnm = 0.0d+0
279
280 kase = 0
281 10 CONTINUE
282 CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
283 IF( kase.NE.0 ) THEN
284 IF( kase.EQ.2 ) THEN
285*
286* Multiply by R.
287*
288 DO i = 1, n
289 work( i ) = work( i ) * work( 2*n+i )
290 END DO
291
292 IF ( notrans ) THEN
293 CALL dgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
294 $ ipiv, work, n, info )
295 ELSE
296 CALL dgbtrs( 'Transpose', n, kl, ku, 1, afb, ldafb,
297 $ ipiv,
298 $ work, n, info )
299 END IF
300*
301* Multiply by inv(C).
302*
303 IF ( cmode .EQ. 1 ) THEN
304 DO i = 1, n
305 work( i ) = work( i ) / c( i )
306 END DO
307 ELSE IF ( cmode .EQ. -1 ) THEN
308 DO i = 1, n
309 work( i ) = work( i ) * c( i )
310 END DO
311 END IF
312 ELSE
313*
314* Multiply by inv(C**T).
315*
316 IF ( cmode .EQ. 1 ) THEN
317 DO i = 1, n
318 work( i ) = work( i ) / c( i )
319 END DO
320 ELSE IF ( cmode .EQ. -1 ) THEN
321 DO i = 1, n
322 work( i ) = work( i ) * c( i )
323 END DO
324 END IF
325
326 IF ( notrans ) THEN
327 CALL dgbtrs( 'Transpose', n, kl, ku, 1, afb, ldafb,
328 $ ipiv,
329 $ work, n, info )
330 ELSE
331 CALL dgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
332 $ ipiv, work, n, info )
333 END IF
334*
335* Multiply by R.
336*
337 DO i = 1, n
338 work( i ) = work( i ) * work( 2*n+i )
339 END DO
340 END IF
341 GO TO 10
342 END IF
343*
344* Compute the estimate of the reciprocal condition number.
345*
346 IF( ainvnm .NE. 0.0d+0 )
347 $ dla_gbrcond = ( 1.0d+0 / ainvnm )
348*
349 RETURN
350*
351* End of DLA_GBRCOND
352*
353 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
DGBTRS
Definition dgbtrs.f:137
double precision function dla_gbrcond(trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, cmode, c, info, work, iwork)
DLA_GBRCOND estimates the Skeel condition number for a general banded matrix.
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:134
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48