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