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