LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
sla_syrcond.f
Go to the documentation of this file.
1 *> \brief \b SLA_SYRCOND estimates the Skeel condition number for a symmetric indefinite 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_SYRCOND + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sla_syrcond.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sla_syrcond.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sla_syrcond.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * REAL FUNCTION SLA_SYRCOND( UPLO, N, A, LDA, AF, LDAF, IPIV, CMODE,
22 * C, INFO, WORK, IWORK )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER N, LDA, LDAF, INFO, CMODE
27 * ..
28 * .. Array Arguments
29 * INTEGER IWORK( * ), IPIV( * )
30 * REAL A( LDA, * ), AF( LDAF, * ), WORK( * ), C( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLA_SYRCOND estimates the Skeel condition number of op(A) * op2(C)
40 *> where op2 is determined by CMODE as follows
41 *> CMODE = 1 op2(C) = C
42 *> CMODE = 0 op2(C) = I
43 *> CMODE = -1 op2(C) = inv(C)
44 *> The Skeel condition number cond(A) = norminf( |inv(A)||A| )
45 *> is computed by computing scaling factors R such that
46 *> diag(R)*A*op2(C) is row equilibrated and computing the standard
47 *> infinity-norm condition number.
48 *> \endverbatim
49 *
50 * Arguments:
51 * ==========
52 *
53 *> \param[in] UPLO
54 *> \verbatim
55 *> UPLO is CHARACTER*1
56 *> = 'U': Upper triangle of A is stored;
57 *> = 'L': Lower triangle of A is stored.
58 *> \endverbatim
59 *>
60 *> \param[in] N
61 *> \verbatim
62 *> N is INTEGER
63 *> The number of linear equations, i.e., the order of the
64 *> matrix A. N >= 0.
65 *> \endverbatim
66 *>
67 *> \param[in] A
68 *> \verbatim
69 *> A is REAL array, dimension (LDA,N)
70 *> On entry, the N-by-N matrix A.
71 *> \endverbatim
72 *>
73 *> \param[in] LDA
74 *> \verbatim
75 *> LDA is INTEGER
76 *> The leading dimension of the array A. LDA >= max(1,N).
77 *> \endverbatim
78 *>
79 *> \param[in] AF
80 *> \verbatim
81 *> AF is REAL array, dimension (LDAF,N)
82 *> The block diagonal matrix D and the multipliers used to
83 *> obtain the factor U or L as computed by SSYTRF.
84 *> \endverbatim
85 *>
86 *> \param[in] LDAF
87 *> \verbatim
88 *> LDAF is INTEGER
89 *> The leading dimension of the array AF. LDAF >= max(1,N).
90 *> \endverbatim
91 *>
92 *> \param[in] IPIV
93 *> \verbatim
94 *> IPIV is INTEGER array, dimension (N)
95 *> Details of the interchanges and the block structure of D
96 *> as determined by SSYTRF.
97 *> \endverbatim
98 *>
99 *> \param[in] CMODE
100 *> \verbatim
101 *> CMODE is INTEGER
102 *> Determines op2(C) in the formula op(A) * op2(C) as follows:
103 *> CMODE = 1 op2(C) = C
104 *> CMODE = 0 op2(C) = I
105 *> CMODE = -1 op2(C) = inv(C)
106 *> \endverbatim
107 *>
108 *> \param[in] C
109 *> \verbatim
110 *> C is REAL array, dimension (N)
111 *> The vector C in the formula op(A) * op2(C).
112 *> \endverbatim
113 *>
114 *> \param[out] INFO
115 *> \verbatim
116 *> INFO is INTEGER
117 *> = 0: Successful exit.
118 *> i > 0: The ith argument is invalid.
119 *> \endverbatim
120 *>
121 *> \param[in] WORK
122 *> \verbatim
123 *> WORK is REAL array, dimension (3*N).
124 *> Workspace.
125 *> \endverbatim
126 *>
127 *> \param[in] IWORK
128 *> \verbatim
129 *> IWORK is INTEGER array, dimension (N).
130 *> Workspace.
131 *> \endverbatim
132 *
133 * Authors:
134 * ========
135 *
136 *> \author Univ. of Tennessee
137 *> \author Univ. of California Berkeley
138 *> \author Univ. of Colorado Denver
139 *> \author NAG Ltd.
140 *
141 *> \date December 2016
142 *
143 *> \ingroup realSYcomputational
144 *
145 * =====================================================================
146  REAL FUNCTION sla_syrcond( UPLO, N, A, LDA, AF, LDAF, IPIV, CMODE,
147  $ C, INFO, WORK, IWORK )
148 *
149 * -- LAPACK computational routine (version 3.7.0) --
150 * -- LAPACK is a software package provided by Univ. of Tennessee, --
151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
152 * December 2016
153 *
154 * .. Scalar Arguments ..
155  CHARACTER UPLO
156  INTEGER N, LDA, LDAF, INFO, CMODE
157 * ..
158 * .. Array Arguments
159  INTEGER IWORK( * ), IPIV( * )
160  REAL A( lda, * ), AF( ldaf, * ), WORK( * ), C( * )
161 * ..
162 *
163 * =====================================================================
164 *
165 * .. Local Scalars ..
166  CHARACTER NORMIN
167  INTEGER KASE, I, J
168  REAL AINVNM, SMLNUM, TMP
169  LOGICAL UP
170 * ..
171 * .. Local Arrays ..
172  INTEGER ISAVE( 3 )
173 * ..
174 * .. External Functions ..
175  LOGICAL LSAME
176  REAL SLAMCH
177  EXTERNAL lsame, slamch
178 * ..
179 * .. External Subroutines ..
180  EXTERNAL slacn2, xerbla, ssytrs
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC abs, max
184 * ..
185 * .. Executable Statements ..
186 *
187  sla_syrcond = 0.0
188 *
189  info = 0
190  IF( n.LT.0 ) THEN
191  info = -2
192  ELSE IF( lda.LT.max( 1, n ) ) THEN
193  info = -4
194  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
195  info = -6
196  END IF
197  IF( info.NE.0 ) THEN
198  CALL xerbla( 'SLA_SYRCOND', -info )
199  RETURN
200  END IF
201  IF( n.EQ.0 ) THEN
202  sla_syrcond = 1.0
203  RETURN
204  END IF
205  up = .false.
206  IF ( lsame( uplo, 'U' ) ) up = .true.
207 *
208 * Compute the equilibration matrix R such that
209 * inv(R)*A*C has unit 1-norm.
210 *
211  IF ( up ) THEN
212  DO i = 1, n
213  tmp = 0.0
214  IF ( cmode .EQ. 1 ) THEN
215  DO j = 1, i
216  tmp = tmp + abs( a( j, i ) * c( j ) )
217  END DO
218  DO j = i+1, n
219  tmp = tmp + abs( a( i, j ) * c( j ) )
220  END DO
221  ELSE IF ( cmode .EQ. 0 ) THEN
222  DO j = 1, i
223  tmp = tmp + abs( a( j, i ) )
224  END DO
225  DO j = i+1, n
226  tmp = tmp + abs( a( i, j ) )
227  END DO
228  ELSE
229  DO j = 1, i
230  tmp = tmp + abs( a( j, i ) / c( j ) )
231  END DO
232  DO j = i+1, n
233  tmp = tmp + abs( a( i, j ) / c( j ) )
234  END DO
235  END IF
236  work( 2*n+i ) = tmp
237  END DO
238  ELSE
239  DO i = 1, n
240  tmp = 0.0
241  IF ( cmode .EQ. 1 ) THEN
242  DO j = 1, i
243  tmp = tmp + abs( a( i, j ) * c( j ) )
244  END DO
245  DO j = i+1, n
246  tmp = tmp + abs( a( j, i ) * c( j ) )
247  END DO
248  ELSE IF ( cmode .EQ. 0 ) THEN
249  DO j = 1, i
250  tmp = tmp + abs( a( i, j ) )
251  END DO
252  DO j = i+1, n
253  tmp = tmp + abs( a( j, i ) )
254  END DO
255  ELSE
256  DO j = 1, i
257  tmp = tmp + abs( a( i, j) / c( j ) )
258  END DO
259  DO j = i+1, n
260  tmp = tmp + abs( a( j, i) / c( j ) )
261  END DO
262  END IF
263  work( 2*n+i ) = tmp
264  END DO
265  ENDIF
266 *
267 * Estimate the norm of inv(op(A)).
268 *
269  smlnum = slamch( 'Safe minimum' )
270  ainvnm = 0.0
271  normin = 'N'
272 
273  kase = 0
274  10 CONTINUE
275  CALL slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
276  IF( kase.NE.0 ) THEN
277  IF( kase.EQ.2 ) THEN
278 *
279 * Multiply by R.
280 *
281  DO i = 1, n
282  work( i ) = work( i ) * work( 2*n+i )
283  END DO
284 
285  IF ( up ) THEN
286  CALL ssytrs( 'U', n, 1, af, ldaf, ipiv, work, n, info )
287  ELSE
288  CALL ssytrs( 'L', n, 1, af, ldaf, ipiv, work, n, info )
289  ENDIF
290 *
291 * Multiply by inv(C).
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  ELSE
303 *
304 * Multiply by inv(C**T).
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 
316  IF ( up ) THEN
317  CALL ssytrs( 'U', n, 1, af, ldaf, ipiv, work, n, info )
318  ELSE
319  CALL ssytrs( 'L', n, 1, af, ldaf, ipiv, work, n, info )
320  ENDIF
321 *
322 * Multiply by R.
323 *
324  DO i = 1, n
325  work( i ) = work( i ) * work( 2*n+i )
326  END DO
327  END IF
328 *
329  GO TO 10
330  END IF
331 *
332 * Compute the estimate of the reciprocal condition number.
333 *
334  IF( ainvnm .NE. 0.0 )
335  $ sla_syrcond = ( 1.0 / ainvnm )
336 *
337  RETURN
338 *
339  END
subroutine ssytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
SSYTRS
Definition: ssytrs.f:122
real function sla_syrcond(UPLO, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
SLA_SYRCOND estimates the Skeel condition number for a symmetric indefinite matrix.
Definition: sla_syrcond.f:148
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:138