LAPACK  3.7.0 LAPACK: Linear Algebra PACKage
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
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[in] WORK
117 *> \verbatim
118 *> WORK is DOUBLE PRECISION array, dimension (3*N).
119 *> Workspace.
120 *> \endverbatim
121 *>
122 *> \param[in] 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 *> \date December 2016
137 *
138 *> \ingroup doublePOcomputational
139 *
140 * =====================================================================
141  DOUBLE PRECISION FUNCTION dla_porcond( UPLO, N, A, LDA, AF, LDAF,
142  \$ cmode, c, info, work,
143  \$ iwork )
144 *
145 * -- LAPACK computational routine (version 3.7.0) --
146 * -- LAPACK is a software package provided by Univ. of Tennessee, --
147 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148 * December 2016
149 *
150 * .. Scalar Arguments ..
151  CHARACTER UPLO
152  INTEGER N, LDA, LDAF, INFO, CMODE
153  DOUBLE PRECISION A( lda, * ), AF( ldaf, * ), WORK( * ),
154  \$ c( * )
155 * ..
156 * .. Array Arguments ..
157  INTEGER IWORK( * )
158 * ..
159 *
160 * =====================================================================
161 *
162 * .. Local Scalars ..
163  INTEGER KASE, I, J
164  DOUBLE PRECISION AINVNM, TMP
165  LOGICAL UP
166 * ..
167 * .. Array Arguments ..
168  INTEGER ISAVE( 3 )
169 * ..
170 * .. External Functions ..
171  LOGICAL LSAME
172  EXTERNAL lsame
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL dlacn2, dpotrs, xerbla
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC abs, max
179 * ..
180 * .. Executable Statements ..
181 *
182  dla_porcond = 0.0d+0
183 *
184  info = 0
185  IF( n.LT.0 ) THEN
186  info = -2
187  END IF
188  IF( info.NE.0 ) THEN
189  CALL xerbla( 'DLA_PORCOND', -info )
190  RETURN
191  END IF
192
193  IF( n.EQ.0 ) THEN
194  dla_porcond = 1.0d+0
195  RETURN
196  END IF
197  up = .false.
198  IF ( lsame( uplo, 'U' ) ) up = .true.
199 *
200 * Compute the equilibration matrix R such that
201 * inv(R)*A*C has unit 1-norm.
202 *
203  IF ( up ) THEN
204  DO i = 1, n
205  tmp = 0.0d+0
206  IF ( cmode .EQ. 1 ) THEN
207  DO j = 1, i
208  tmp = tmp + abs( a( j, i ) * c( j ) )
209  END DO
210  DO j = i+1, n
211  tmp = tmp + abs( a( i, j ) * c( j ) )
212  END DO
213  ELSE IF ( cmode .EQ. 0 ) THEN
214  DO j = 1, i
215  tmp = tmp + abs( a( j, i ) )
216  END DO
217  DO j = i+1, n
218  tmp = tmp + abs( a( i, j ) )
219  END DO
220  ELSE
221  DO j = 1, i
222  tmp = tmp + abs( a( j ,i ) / c( j ) )
223  END DO
224  DO j = i+1, n
225  tmp = tmp + abs( a( i, j ) / c( j ) )
226  END DO
227  END IF
228  work( 2*n+i ) = tmp
229  END DO
230  ELSE
231  DO i = 1, n
232  tmp = 0.0d+0
233  IF ( cmode .EQ. 1 ) THEN
234  DO j = 1, i
235  tmp = tmp + abs( a( i, j ) * c( j ) )
236  END DO
237  DO j = i+1, n
238  tmp = tmp + abs( a( j, i ) * c( j ) )
239  END DO
240  ELSE IF ( cmode .EQ. 0 ) THEN
241  DO j = 1, i
242  tmp = tmp + abs( a( i, j ) )
243  END DO
244  DO j = i+1, n
245  tmp = tmp + abs( a( j, i ) )
246  END DO
247  ELSE
248  DO j = 1, i
249  tmp = tmp + abs( a( i, j ) / c( j ) )
250  END DO
251  DO j = i+1, n
252  tmp = tmp + abs( a( j, i ) / c( j ) )
253  END DO
254  END IF
255  work( 2*n+i ) = tmp
256  END DO
257  ENDIF
258 *
259 * Estimate the norm of inv(op(A)).
260 *
261  ainvnm = 0.0d+0
262
263  kase = 0
264  10 CONTINUE
265  CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
266  IF( kase.NE.0 ) THEN
267  IF( kase.EQ.2 ) THEN
268 *
269 * Multiply by R.
270 *
271  DO i = 1, n
272  work( i ) = work( i ) * work( 2*n+i )
273  END DO
274
275  IF (up) THEN
276  CALL dpotrs( 'Upper', n, 1, af, ldaf, work, n, info )
277  ELSE
278  CALL dpotrs( 'Lower', n, 1, af, ldaf, work, n, info )
279  ENDIF
280 *
281 * Multiply by inv(C).
282 *
283  IF ( cmode .EQ. 1 ) THEN
284  DO i = 1, n
285  work( i ) = work( i ) / c( i )
286  END DO
287  ELSE IF ( cmode .EQ. -1 ) THEN
288  DO i = 1, n
289  work( i ) = work( i ) * c( i )
290  END DO
291  END IF
292  ELSE
293 *
294 * Multiply by inv(C**T).
295 *
296  IF ( cmode .EQ. 1 ) THEN
297  DO i = 1, n
298  work( i ) = work( i ) / c( i )
299  END DO
300  ELSE IF ( cmode .EQ. -1 ) THEN
301  DO i = 1, n
302  work( i ) = work( i ) * c( i )
303  END DO
304  END IF
305
306  IF ( up ) THEN
307  CALL dpotrs( 'Upper', n, 1, af, ldaf, work, n, info )
308  ELSE
309  CALL dpotrs( 'Lower', n, 1, af, ldaf, work, n, info )
310  ENDIF
311 *
312 * Multiply by R.
313 *
314  DO i = 1, n
315  work( i ) = work( i ) * work( 2*n+i )
316  END DO
317  END IF
318  GO TO 10
319  END IF
320 *
321 * Compute the estimate of the reciprocal condition number.
322 *
323  IF( ainvnm .NE. 0.0d+0 )
324  \$ dla_porcond = ( 1.0d+0 / ainvnm )
325 *
326  RETURN
327 *
328  END
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:144
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 dpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
DPOTRS
Definition: dpotrs.f:112
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62