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