LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
double precision function dla_porcond ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
integer  CMODE,
double precision, dimension( * )  C,
integer  INFO,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK 
)

DLA_PORCOND estimates the Skeel condition number for a symmetric positive-definite matrix.

Download DLA_PORCOND + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
    DLA_PORCOND Estimates the Skeel condition number of  op(A) * op2(C)
    where op2 is determined by CMODE as follows
    CMODE =  1    op2(C) = C
    CMODE =  0    op2(C) = I
    CMODE = -1    op2(C) = inv(C)
    The Skeel condition number  cond(A) = norminf( |inv(A)||A| )
    is computed by computing scaling factors R such that
    diag(R)*A*op2(C) is row equilibrated and computing the standard
    infinity-norm condition number.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
       = 'U':  Upper triangle of A is stored;
       = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
     The number of linear equations, i.e., the order of the
     matrix A.  N >= 0.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
     On entry, the N-by-N matrix A.
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is DOUBLE PRECISION array, dimension (LDAF,N)
     The triangular factor U or L from the Cholesky factorization
     A = U**T*U or A = L*L**T, as computed by DPOTRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= max(1,N).
[in]CMODE
          CMODE is INTEGER
     Determines op2(C) in the formula op(A) * op2(C) as follows:
     CMODE =  1    op2(C) = C
     CMODE =  0    op2(C) = I
     CMODE = -1    op2(C) = inv(C)
[in]C
          C is DOUBLE PRECISION array, dimension (N)
     The vector C in the formula op(A) * op2(C).
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit.
     i > 0:  The ith argument is invalid.
[in]WORK
          WORK is DOUBLE PRECISION array, dimension (3*N).
     Workspace.
[in]IWORK
          IWORK is INTEGER array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 144 of file dla_porcond.f.

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 *
dla_porcond
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
dlacn2
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
dpotrs
subroutine dpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
DPOTRS
Definition: dpotrs.f:112
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: