LAPACK  3.6.1
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
September 2012

Definition at line 144 of file dla_porcond.f.

144 *
145 * -- LAPACK computational routine (version 3.4.2) --
146 * -- LAPACK is a software package provided by Univ. of Tennessee, --
147 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148 * September 2012
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  INTEGER idamax
173  EXTERNAL lsame, idamax
174 * ..
175 * .. External Subroutines ..
176  EXTERNAL dlacn2, dpotrs, xerbla
177 * ..
178 * .. Intrinsic Functions ..
179  INTRINSIC abs, max
180 * ..
181 * .. Executable Statements ..
182 *
183  dla_porcond = 0.0d+0
184 *
185  info = 0
186  IF( n.LT.0 ) THEN
187  info = -2
188  END IF
189  IF( info.NE.0 ) THEN
190  CALL xerbla( 'DLA_PORCOND', -info )
191  RETURN
192  END IF
193 
194  IF( n.EQ.0 ) THEN
195  dla_porcond = 1.0d+0
196  RETURN
197  END IF
198  up = .false.
199  IF ( lsame( uplo, 'U' ) ) up = .true.
200 *
201 * Compute the equilibration matrix R such that
202 * inv(R)*A*C has unit 1-norm.
203 *
204  IF ( up ) THEN
205  DO i = 1, n
206  tmp = 0.0d+0
207  IF ( cmode .EQ. 1 ) THEN
208  DO j = 1, i
209  tmp = tmp + abs( a( j, i ) * c( j ) )
210  END DO
211  DO j = i+1, n
212  tmp = tmp + abs( a( i, j ) * c( j ) )
213  END DO
214  ELSE IF ( cmode .EQ. 0 ) THEN
215  DO j = 1, i
216  tmp = tmp + abs( a( j, i ) )
217  END DO
218  DO j = i+1, n
219  tmp = tmp + abs( a( i, j ) )
220  END DO
221  ELSE
222  DO j = 1, i
223  tmp = tmp + abs( a( j ,i ) / c( j ) )
224  END DO
225  DO j = i+1, n
226  tmp = tmp + abs( a( i, j ) / c( j ) )
227  END DO
228  END IF
229  work( 2*n+i ) = tmp
230  END DO
231  ELSE
232  DO i = 1, n
233  tmp = 0.0d+0
234  IF ( cmode .EQ. 1 ) THEN
235  DO j = 1, i
236  tmp = tmp + abs( a( i, j ) * c( j ) )
237  END DO
238  DO j = i+1, n
239  tmp = tmp + abs( a( j, i ) * c( j ) )
240  END DO
241  ELSE IF ( cmode .EQ. 0 ) THEN
242  DO j = 1, i
243  tmp = tmp + abs( a( i, j ) )
244  END DO
245  DO j = i+1, n
246  tmp = tmp + abs( a( j, i ) )
247  END DO
248  ELSE
249  DO j = 1, i
250  tmp = tmp + abs( a( i, j ) / c( j ) )
251  END DO
252  DO j = i+1, n
253  tmp = tmp + abs( a( j, i ) / c( j ) )
254  END DO
255  END IF
256  work( 2*n+i ) = tmp
257  END DO
258  ENDIF
259 *
260 * Estimate the norm of inv(op(A)).
261 *
262  ainvnm = 0.0d+0
263 
264  kase = 0
265  10 CONTINUE
266  CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
267  IF( kase.NE.0 ) THEN
268  IF( kase.EQ.2 ) THEN
269 *
270 * Multiply by R.
271 *
272  DO i = 1, n
273  work( i ) = work( i ) * work( 2*n+i )
274  END DO
275 
276  IF (up) THEN
277  CALL dpotrs( 'Upper', n, 1, af, ldaf, work, n, info )
278  ELSE
279  CALL dpotrs( 'Lower', n, 1, af, ldaf, work, n, info )
280  ENDIF
281 *
282 * Multiply by inv(C).
283 *
284  IF ( cmode .EQ. 1 ) THEN
285  DO i = 1, n
286  work( i ) = work( i ) / c( i )
287  END DO
288  ELSE IF ( cmode .EQ. -1 ) THEN
289  DO i = 1, n
290  work( i ) = work( i ) * c( i )
291  END DO
292  END IF
293  ELSE
294 *
295 * Multiply by inv(C**T).
296 *
297  IF ( cmode .EQ. 1 ) THEN
298  DO i = 1, n
299  work( i ) = work( i ) / c( i )
300  END DO
301  ELSE IF ( cmode .EQ. -1 ) THEN
302  DO i = 1, n
303  work( i ) = work( i ) * c( i )
304  END DO
305  END IF
306 
307  IF ( up ) THEN
308  CALL dpotrs( 'Upper', n, 1, af, ldaf, work, n, info )
309  ELSE
310  CALL dpotrs( 'Lower', n, 1, af, ldaf, work, n, info )
311  ENDIF
312 *
313 * Multiply by R.
314 *
315  DO i = 1, n
316  work( i ) = work( i ) * work( 2*n+i )
317  END DO
318  END IF
319  GO TO 10
320  END IF
321 *
322 * Compute the estimate of the reciprocal condition number.
323 *
324  IF( ainvnm .NE. 0.0d+0 )
325  $ dla_porcond = ( 1.0d+0 / ainvnm )
326 *
327  RETURN
328 *
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
idamax
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dpotrs
subroutine dpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
DPOTRS
Definition: dpotrs.f:112
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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

Here is the call graph for this function:

Here is the caller graph for this function: