LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dla_porcond()

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.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (3*N).
     Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 139 of file dla_porcond.f.

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 *
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
Here is the call graph for this function:
Here is the caller graph for this function: