LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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[out] WORK
126 *> \verbatim
127 *> WORK is REAL array, dimension (3*N).
128 *> Workspace.
129 *> \endverbatim
130 *>
131 *> \param[out] 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 *> \ingroup realGEcomputational
146 *
147 * =====================================================================
148  REAL function sla_gercond( trans, n, a, lda, af, ldaf, ipiv,
149  $ cmode, c, info, work, iwork )
150 *
151 * -- LAPACK computational routine --
152 * -- LAPACK is a software package provided by Univ. of Tennessee, --
153 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
154 *
155 * .. Scalar Arguments ..
156  CHARACTER trans
157  INTEGER n, lda, ldaf, info, cmode
158 * ..
159 * .. Array Arguments ..
160  INTEGER ipiv( * ), iwork( * )
161  REAL a( lda, * ), af( ldaf, * ), work( * ),
162  $ c( * )
163 * ..
164 *
165 * =====================================================================
166 *
167 * .. Local Scalars ..
168  LOGICAL notrans
169  INTEGER kase, i, j
170  REAL ainvnm, tmp
171 * ..
172 * .. Local Arrays ..
173  INTEGER isave( 3 )
174 * ..
175 * .. External Functions ..
176  LOGICAL lsame
177  EXTERNAL lsame
178 * ..
179 * .. External Subroutines ..
180  EXTERNAL slacn2, sgetrs, xerbla
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC abs, max
184 * ..
185 * .. Executable Statements ..
186 *
187  sla_gercond = 0.0
188 *
189  info = 0
190  notrans = lsame( trans, 'N' )
191  IF ( .NOT. notrans .AND. .NOT. lsame(trans, 'T')
192  $ .AND. .NOT. lsame(trans, 'C') ) THEN
193  info = -1
194  ELSE IF( n.LT.0 ) THEN
195  info = -2
196  ELSE IF( lda.LT.max( 1, n ) ) THEN
197  info = -4
198  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
199  info = -6
200  END IF
201  IF( info.NE.0 ) THEN
202  CALL xerbla( 'SLA_GERCOND', -info )
203  RETURN
204  END IF
205  IF( n.EQ.0 ) THEN
206  sla_gercond = 1.0
207  RETURN
208  END IF
209 *
210 * Compute the equilibration matrix R such that
211 * inv(R)*A*C has unit 1-norm.
212 *
213  IF (notrans) THEN
214  DO i = 1, n
215  tmp = 0.0
216  IF ( cmode .EQ. 1 ) THEN
217  DO j = 1, n
218  tmp = tmp + abs( a( i, j ) * c( j ) )
219  END DO
220  ELSE IF ( cmode .EQ. 0 ) THEN
221  DO j = 1, n
222  tmp = tmp + abs( a( i, j ) )
223  END DO
224  ELSE
225  DO j = 1, n
226  tmp = tmp + abs( a( i, j ) / c( j ) )
227  END DO
228  END IF
229  work( 2*n+i ) = tmp
230  END DO
231  ELSE
232  DO i = 1, n
233  tmp = 0.0
234  IF ( cmode .EQ. 1 ) THEN
235  DO j = 1, n
236  tmp = tmp + abs( a( j, i ) * c( j ) )
237  END DO
238  ELSE IF ( cmode .EQ. 0 ) THEN
239  DO j = 1, n
240  tmp = tmp + abs( a( j, i ) )
241  END DO
242  ELSE
243  DO j = 1, n
244  tmp = tmp + abs( a( j, i ) / c( j ) )
245  END DO
246  END IF
247  work( 2*n+i ) = tmp
248  END DO
249  END IF
250 *
251 * Estimate the norm of inv(op(A)).
252 *
253  ainvnm = 0.0
254 
255  kase = 0
256  10 CONTINUE
257  CALL slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
258  IF( kase.NE.0 ) THEN
259  IF( kase.EQ.2 ) THEN
260 *
261 * Multiply by R.
262 *
263  DO i = 1, n
264  work(i) = work(i) * work(2*n+i)
265  END DO
266 
267  IF (notrans) THEN
268  CALL sgetrs( 'No transpose', n, 1, af, ldaf, ipiv,
269  $ work, n, info )
270  ELSE
271  CALL sgetrs( 'Transpose', n, 1, af, ldaf, ipiv,
272  $ work, n, info )
273  END IF
274 *
275 * Multiply by inv(C).
276 *
277  IF ( cmode .EQ. 1 ) THEN
278  DO i = 1, n
279  work( i ) = work( i ) / c( i )
280  END DO
281  ELSE IF ( cmode .EQ. -1 ) THEN
282  DO i = 1, n
283  work( i ) = work( i ) * c( i )
284  END DO
285  END IF
286  ELSE
287 *
288 * Multiply by inv(C**T).
289 *
290  IF ( cmode .EQ. 1 ) THEN
291  DO i = 1, n
292  work( i ) = work( i ) / c( i )
293  END DO
294  ELSE IF ( cmode .EQ. -1 ) THEN
295  DO i = 1, n
296  work( i ) = work( i ) * c( i )
297  END DO
298  END IF
299 
300  IF (notrans) THEN
301  CALL sgetrs( 'Transpose', n, 1, af, ldaf, ipiv,
302  $ work, n, info )
303  ELSE
304  CALL sgetrs( 'No transpose', n, 1, af, ldaf, ipiv,
305  $ work, n, info )
306  END IF
307 *
308 * Multiply by R.
309 *
310  DO i = 1, n
311  work( i ) = work( i ) * work( 2*n+i )
312  END DO
313  END IF
314  GO TO 10
315  END IF
316 *
317 * Compute the estimate of the reciprocal condition number.
318 *
319  IF( ainvnm .NE. 0.0 )
320  $ sla_gercond = ( 1.0 / ainvnm )
321 *
322  RETURN
323 *
324 * End of SLA_GERCOND
325 *
326  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function sla_gercond(TRANS, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
SLA_GERCOND estimates the Skeel condition number for a general matrix.
Definition: sla_gercond.f:150
subroutine sgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
SGETRS
Definition: sgetrs.f:121
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:136