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