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