LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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