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