LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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[in] WORK
144 *> \verbatim
145 *> WORK is REAL array, dimension (5*N).
146 *> Workspace.
147 *> \endverbatim
148 *>
149 *> \param[in] 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 *> \date September 2012
164 *
165 *> \ingroup realGBcomputational
166 *
167 * =====================================================================
168  REAL FUNCTION sla_gbrcond( TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB,
169  $ ipiv, cmode, c, info, work, iwork )
170 *
171 * -- LAPACK computational routine (version 3.4.2) --
172 * -- LAPACK is a software package provided by Univ. of Tennessee, --
173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174 * September 2012
175 *
176 * .. Scalar Arguments ..
177  CHARACTER trans
178  INTEGER n, ldab, ldafb, info, kl, ku, cmode
179 * ..
180 * .. Array Arguments ..
181  INTEGER iwork( * ), ipiv( * )
182  REAL ab( ldab, * ), afb( ldafb, * ), work( * ),
183  $ c( * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * .. Local Scalars ..
189  LOGICAL notrans
190  INTEGER kase, i, j, kd, ke
191  REAL ainvnm, tmp
192 * ..
193 * .. Local Arrays ..
194  INTEGER isave( 3 )
195 * ..
196 * .. External Functions ..
197  LOGICAL lsame
198  EXTERNAL lsame
199 * ..
200 * .. External Subroutines ..
201  EXTERNAL slacn2, sgbtrs, xerbla
202 * ..
203 * .. Intrinsic Functions ..
204  INTRINSIC abs, max
205 * ..
206 * .. Executable Statements ..
207 *
208  sla_gbrcond = 0.0
209 *
210  info = 0
211  notrans = lsame( trans, 'N' )
212  IF ( .NOT. notrans .AND. .NOT. lsame(trans, 'T')
213  $ .AND. .NOT. lsame(trans, 'C') ) THEN
214  info = -1
215  ELSE IF( n.LT.0 ) THEN
216  info = -2
217  ELSE IF( kl.LT.0 .OR. kl.GT.n-1 ) THEN
218  info = -3
219  ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
220  info = -4
221  ELSE IF( ldab.LT.kl+ku+1 ) THEN
222  info = -6
223  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
224  info = -8
225  END IF
226  IF( info.NE.0 ) THEN
227  CALL xerbla( 'SLA_GBRCOND', -info )
228  RETURN
229  END IF
230  IF( n.EQ.0 ) THEN
231  sla_gbrcond = 1.0
232  RETURN
233  END IF
234 *
235 * Compute the equilibration matrix R such that
236 * inv(R)*A*C has unit 1-norm.
237 *
238  kd = ku + 1
239  ke = kl + 1
240  IF ( notrans ) THEN
241  DO i = 1, n
242  tmp = 0.0
243  IF ( cmode .EQ. 1 ) THEN
244  DO j = max( i-kl, 1 ), min( i+ku, n )
245  tmp = tmp + abs( ab( kd+i-j, j ) * c( j ) )
246  END DO
247  ELSE IF ( cmode .EQ. 0 ) THEN
248  DO j = max( i-kl, 1 ), min( i+ku, n )
249  tmp = tmp + abs( ab( kd+i-j, j ) )
250  END DO
251  ELSE
252  DO j = max( i-kl, 1 ), min( i+ku, n )
253  tmp = tmp + abs( ab( kd+i-j, j ) / c( j ) )
254  END DO
255  END IF
256  work( 2*n+i ) = tmp
257  END DO
258  ELSE
259  DO i = 1, n
260  tmp = 0.0
261  IF ( cmode .EQ. 1 ) THEN
262  DO j = max( i-kl, 1 ), min( i+ku, n )
263  tmp = tmp + abs( ab( ke-i+j, i ) * c( j ) )
264  END DO
265  ELSE IF ( cmode .EQ. 0 ) THEN
266  DO j = max( i-kl, 1 ), min( i+ku, n )
267  tmp = tmp + abs( ab( ke-i+j, i ) )
268  END DO
269  ELSE
270  DO j = max( i-kl, 1 ), min( i+ku, n )
271  tmp = tmp + abs( ab( ke-i+j, i ) / c( j ) )
272  END DO
273  END IF
274  work( 2*n+i ) = tmp
275  END DO
276  END IF
277 *
278 * Estimate the norm of inv(op(A)).
279 *
280  ainvnm = 0.0
281 
282  kase = 0
283  10 CONTINUE
284  CALL slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
285  IF( kase.NE.0 ) THEN
286  IF( kase.EQ.2 ) THEN
287 *
288 * Multiply by R.
289 *
290  DO i = 1, n
291  work( i ) = work( i ) * work( 2*n+i )
292  END DO
293 
294  IF ( notrans ) THEN
295  CALL sgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
296  $ ipiv, work, n, info )
297  ELSE
298  CALL sgbtrs( 'Transpose', n, kl, ku, 1, afb, ldafb, ipiv,
299  $ work, n, info )
300  END IF
301 *
302 * Multiply by inv(C).
303 *
304  IF ( cmode .EQ. 1 ) THEN
305  DO i = 1, n
306  work( i ) = work( i ) / c( i )
307  END DO
308  ELSE IF ( cmode .EQ. -1 ) THEN
309  DO i = 1, n
310  work( i ) = work( i ) * c( i )
311  END DO
312  END IF
313  ELSE
314 *
315 * Multiply by inv(C**T).
316 *
317  IF ( cmode .EQ. 1 ) THEN
318  DO i = 1, n
319  work( i ) = work( i ) / c( i )
320  END DO
321  ELSE IF ( cmode .EQ. -1 ) THEN
322  DO i = 1, n
323  work( i ) = work( i ) * c( i )
324  END DO
325  END IF
326 
327  IF ( notrans ) THEN
328  CALL sgbtrs( 'Transpose', n, kl, ku, 1, afb, ldafb, ipiv,
329  $ work, n, info )
330  ELSE
331  CALL sgbtrs( '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.0 )
347  $ sla_gbrcond = ( 1.0 / ainvnm )
348 *
349  RETURN
350 *
351  END