LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dla_gercond.f
Go to the documentation of this file.
1 *> \brief \b DLA_GERCOND estimates the Skeel condition number for a general 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_GERCOND + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gercond.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gercond.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gercond.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * DOUBLE PRECISION FUNCTION DLA_GERCOND ( TRANS, N, A, LDA, AF,
22 * LDAF, IPIV, CMODE, C,
23 * INFO, WORK, IWORK )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER TRANS
27 * INTEGER N, LDA, LDAF, INFO, CMODE
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IPIV( * ), IWORK( * )
31 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
32 * $ C( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DLA_GERCOND 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] A
72 *> \verbatim
73 *> A is DOUBLE PRECISION array, dimension (LDA,N)
74 *> On entry, the N-by-N matrix A.
75 *> \endverbatim
76 *>
77 *> \param[in] LDA
78 *> \verbatim
79 *> LDA is INTEGER
80 *> The leading dimension of the array A. LDA >= max(1,N).
81 *> \endverbatim
82 *>
83 *> \param[in] AF
84 *> \verbatim
85 *> AF is DOUBLE PRECISION array, dimension (LDAF,N)
86 *> The factors L and U from the factorization
87 *> A = P*L*U as computed by DGETRF.
88 *> \endverbatim
89 *>
90 *> \param[in] LDAF
91 *> \verbatim
92 *> LDAF is INTEGER
93 *> The leading dimension of the array AF. LDAF >= max(1,N).
94 *> \endverbatim
95 *>
96 *> \param[in] IPIV
97 *> \verbatim
98 *> IPIV is INTEGER array, dimension (N)
99 *> The pivot indices from the factorization A = P*L*U
100 *> as computed by DGETRF; row i of the matrix was interchanged
101 *> with row IPIV(i).
102 *> \endverbatim
103 *>
104 *> \param[in] CMODE
105 *> \verbatim
106 *> CMODE is INTEGER
107 *> Determines op2(C) in the formula op(A) * op2(C) as follows:
108 *> CMODE = 1 op2(C) = C
109 *> CMODE = 0 op2(C) = I
110 *> CMODE = -1 op2(C) = inv(C)
111 *> \endverbatim
112 *>
113 *> \param[in] C
114 *> \verbatim
115 *> C is DOUBLE PRECISION array, dimension (N)
116 *> The vector C in the formula op(A) * op2(C).
117 *> \endverbatim
118 *>
119 *> \param[out] INFO
120 *> \verbatim
121 *> INFO is INTEGER
122 *> = 0: Successful exit.
123 *> i > 0: The ith argument is invalid.
124 *> \endverbatim
125 *>
126 *> \param[in] WORK
127 *> \verbatim
128 *> WORK is DOUBLE PRECISION array, dimension (3*N).
129 *> Workspace.
130 *> \endverbatim
131 *>
132 *> \param[in] IWORK
133 *> \verbatim
134 *> IWORK is INTEGER array, dimension (N).
135 *> Workspace.
136 *> \endverbatim
137 *
138 * Authors:
139 * ========
140 *
141 *> \author Univ. of Tennessee
142 *> \author Univ. of California Berkeley
143 *> \author Univ. of Colorado Denver
144 *> \author NAG Ltd.
145 *
146 *> \date September 2012
147 *
148 *> \ingroup doubleGEcomputational
149 *
150 * =====================================================================
151  DOUBLE PRECISION FUNCTION dla_gercond ( TRANS, N, A, LDA, AF,
152  $ ldaf, ipiv, cmode, c,
153  $ info, work, iwork )
154 *
155 * -- LAPACK computational routine (version 3.4.2) --
156 * -- LAPACK is a software package provided by Univ. of Tennessee, --
157 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
158 * September 2012
159 *
160 * .. Scalar Arguments ..
161  CHARACTER trans
162  INTEGER n, lda, ldaf, info, cmode
163 * ..
164 * .. Array Arguments ..
165  INTEGER ipiv( * ), iwork( * )
166  DOUBLE PRECISION a( lda, * ), af( ldaf, * ), work( * ),
167  $ c( * )
168 * ..
169 *
170 * =====================================================================
171 *
172 * .. Local Scalars ..
173  LOGICAL notrans
174  INTEGER kase, i, j
175  DOUBLE PRECISION ainvnm, tmp
176 * ..
177 * .. Local Arrays ..
178  INTEGER isave( 3 )
179 * ..
180 * .. External Functions ..
181  LOGICAL lsame
182  EXTERNAL lsame
183 * ..
184 * .. External Subroutines ..
185  EXTERNAL dlacn2, dgetrs, xerbla
186 * ..
187 * .. Intrinsic Functions ..
188  INTRINSIC abs, max
189 * ..
190 * .. Executable Statements ..
191 *
192  dla_gercond = 0.0d+0
193 *
194  info = 0
195  notrans = lsame( trans, 'N' )
196  IF ( .NOT. notrans .AND. .NOT. lsame(trans, 'T')
197  $ .AND. .NOT. lsame(trans, 'C') ) THEN
198  info = -1
199  ELSE IF( n.LT.0 ) THEN
200  info = -2
201  ELSE IF( lda.LT.max( 1, n ) ) THEN
202  info = -4
203  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
204  info = -6
205  END IF
206  IF( info.NE.0 ) THEN
207  CALL xerbla( 'DLA_GERCOND', -info )
208  return
209  END IF
210  IF( n.EQ.0 ) THEN
211  dla_gercond = 1.0d+0
212  return
213  END IF
214 *
215 * Compute the equilibration matrix R such that
216 * inv(R)*A*C has unit 1-norm.
217 *
218  IF (notrans) THEN
219  DO i = 1, n
220  tmp = 0.0d+0
221  IF ( cmode .EQ. 1 ) THEN
222  DO j = 1, n
223  tmp = tmp + abs( a( i, j ) * c( j ) )
224  END DO
225  ELSE IF ( cmode .EQ. 0 ) THEN
226  DO j = 1, n
227  tmp = tmp + abs( a( i, j ) )
228  END DO
229  ELSE
230  DO j = 1, n
231  tmp = tmp + abs( a( i, j ) / c( j ) )
232  END DO
233  END IF
234  work( 2*n+i ) = tmp
235  END DO
236  ELSE
237  DO i = 1, n
238  tmp = 0.0d+0
239  IF ( cmode .EQ. 1 ) THEN
240  DO j = 1, n
241  tmp = tmp + abs( a( j, i ) * c( j ) )
242  END DO
243  ELSE IF ( cmode .EQ. 0 ) THEN
244  DO j = 1, n
245  tmp = tmp + abs( a( j, i ) )
246  END DO
247  ELSE
248  DO j = 1, n
249  tmp = tmp + abs( a( j, i ) / c( j ) )
250  END DO
251  END IF
252  work( 2*n+i ) = tmp
253  END DO
254  END IF
255 *
256 * Estimate the norm of inv(op(A)).
257 *
258  ainvnm = 0.0d+0
259 
260  kase = 0
261  10 continue
262  CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
263  IF( kase.NE.0 ) THEN
264  IF( kase.EQ.2 ) THEN
265 *
266 * Multiply by R.
267 *
268  DO i = 1, n
269  work(i) = work(i) * work(2*n+i)
270  END DO
271 
272  IF (notrans) THEN
273  CALL dgetrs( 'No transpose', n, 1, af, ldaf, ipiv,
274  $ work, n, info )
275  ELSE
276  CALL dgetrs( 'Transpose', n, 1, af, ldaf, ipiv,
277  $ work, n, info )
278  END IF
279 *
280 * Multiply by inv(C).
281 *
282  IF ( cmode .EQ. 1 ) THEN
283  DO i = 1, n
284  work( i ) = work( i ) / c( i )
285  END DO
286  ELSE IF ( cmode .EQ. -1 ) THEN
287  DO i = 1, n
288  work( i ) = work( i ) * c( i )
289  END DO
290  END IF
291  ELSE
292 *
293 * Multiply by inv(C**T).
294 *
295  IF ( cmode .EQ. 1 ) THEN
296  DO i = 1, n
297  work( i ) = work( i ) / c( i )
298  END DO
299  ELSE IF ( cmode .EQ. -1 ) THEN
300  DO i = 1, n
301  work( i ) = work( i ) * c( i )
302  END DO
303  END IF
304 
305  IF (notrans) THEN
306  CALL dgetrs( 'Transpose', n, 1, af, ldaf, ipiv,
307  $ work, n, info )
308  ELSE
309  CALL dgetrs( 'No transpose', n, 1, af, ldaf, ipiv,
310  $ work, n, info )
311  END IF
312 *
313 * Multiply by R.
314 *
315  DO i = 1, n
316  work( i ) = work( i ) * work( 2*n+i )
317  END DO
318  END IF
319  go to 10
320  END IF
321 *
322 * Compute the estimate of the reciprocal condition number.
323 *
324  IF( ainvnm .NE. 0.0d+0 )
325  $ dla_gercond = ( 1.0d+0 / ainvnm )
326 *
327  return
328 *
329  END