LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
real function sla_porcond ( character  UPLO,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldaf, * )  AF,
integer  LDAF,
integer  CMODE,
real, dimension( * )  C,
integer  INFO,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK 
)

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

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

Purpose:
    SLA_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 REAL 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 REAL 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 SPOTRF.
[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 REAL 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 REAL 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 142 of file sla_porcond.f.

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