LAPACK 3.12.1
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 137 of file dla_porcond.f.

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