LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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:142
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:136
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dpotrs
subroutine dpotrs(uplo, n, nrhs, a, lda, b, ldb, info)
DPOTRS
Definition dpotrs.f:110
Here is the call graph for this function:
Here is the caller graph for this function: