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