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