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

◆ dporfs()

subroutine dporfs ( character  uplo,
integer  n,
integer  nrhs,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( ldaf, * )  af,
integer  ldaf,
double precision, dimension( ldb, * )  b,
integer  ldb,
double precision, dimension( ldx, * )  x,
integer  ldx,
double precision, dimension( * )  ferr,
double precision, dimension( * )  berr,
double precision, dimension( * )  work,
integer, dimension( * )  iwork,
integer  info 
)

DPORFS

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

Purpose:
 DPORFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric positive definite,
 and provides error bounds and backward error estimates for the
 solution.
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 order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
          upper triangular part of A contains the upper triangular part
          of the matrix A, and the strictly lower triangular part of A
          is not referenced.  If UPLO = 'L', the leading N-by-N lower
          triangular part of A contains the lower triangular part of
          the matrix A, and the strictly upper triangular part of A is
          not referenced.
[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]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
          The right hand side matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DPOTRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is DOUBLE PRECISION array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (3*N)
[out]IWORK
          IWORK is INTEGER array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Internal Parameters:
  ITMAX is the maximum number of steps of iterative refinement.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 181 of file dporfs.f.

183*
184* -- LAPACK computational routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER UPLO
190 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
191* ..
192* .. Array Arguments ..
193 INTEGER IWORK( * )
194 DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
195 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 INTEGER ITMAX
202 parameter( itmax = 5 )
203 DOUBLE PRECISION ZERO
204 parameter( zero = 0.0d+0 )
205 DOUBLE PRECISION ONE
206 parameter( one = 1.0d+0 )
207 DOUBLE PRECISION TWO
208 parameter( two = 2.0d+0 )
209 DOUBLE PRECISION THREE
210 parameter( three = 3.0d+0 )
211* ..
212* .. Local Scalars ..
213 LOGICAL UPPER
214 INTEGER COUNT, I, J, K, KASE, NZ
215 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
216* ..
217* .. Local Arrays ..
218 INTEGER ISAVE( 3 )
219* ..
220* .. External Subroutines ..
221 EXTERNAL daxpy, dcopy, dlacn2, dpotrs, dsymv, xerbla
222* ..
223* .. Intrinsic Functions ..
224 INTRINSIC abs, max
225* ..
226* .. External Functions ..
227 LOGICAL LSAME
228 DOUBLE PRECISION DLAMCH
229 EXTERNAL lsame, dlamch
230* ..
231* .. Executable Statements ..
232*
233* Test the input parameters.
234*
235 info = 0
236 upper = lsame( uplo, 'U' )
237 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
238 info = -1
239 ELSE IF( n.LT.0 ) THEN
240 info = -2
241 ELSE IF( nrhs.LT.0 ) THEN
242 info = -3
243 ELSE IF( lda.LT.max( 1, n ) ) THEN
244 info = -5
245 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
246 info = -7
247 ELSE IF( ldb.LT.max( 1, n ) ) THEN
248 info = -9
249 ELSE IF( ldx.LT.max( 1, n ) ) THEN
250 info = -11
251 END IF
252 IF( info.NE.0 ) THEN
253 CALL xerbla( 'DPORFS', -info )
254 RETURN
255 END IF
256*
257* Quick return if possible
258*
259 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
260 DO 10 j = 1, nrhs
261 ferr( j ) = zero
262 berr( j ) = zero
263 10 CONTINUE
264 RETURN
265 END IF
266*
267* NZ = maximum number of nonzero elements in each row of A, plus 1
268*
269 nz = n + 1
270 eps = dlamch( 'Epsilon' )
271 safmin = dlamch( 'Safe minimum' )
272 safe1 = nz*safmin
273 safe2 = safe1 / eps
274*
275* Do for each right hand side
276*
277 DO 140 j = 1, nrhs
278*
279 count = 1
280 lstres = three
281 20 CONTINUE
282*
283* Loop until stopping criterion is satisfied.
284*
285* Compute residual R = B - A * X
286*
287 CALL dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
288 CALL dsymv( uplo, n, -one, a, lda, x( 1, j ), 1, one,
289 $ work( n+1 ), 1 )
290*
291* Compute componentwise relative backward error from formula
292*
293* max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
294*
295* where abs(Z) is the componentwise absolute value of the matrix
296* or vector Z. If the i-th component of the denominator is less
297* than SAFE2, then SAFE1 is added to the i-th components of the
298* numerator and denominator before dividing.
299*
300 DO 30 i = 1, n
301 work( i ) = abs( b( i, j ) )
302 30 CONTINUE
303*
304* Compute abs(A)*abs(X) + abs(B).
305*
306 IF( upper ) THEN
307 DO 50 k = 1, n
308 s = zero
309 xk = abs( x( k, j ) )
310 DO 40 i = 1, k - 1
311 work( i ) = work( i ) + abs( a( i, k ) )*xk
312 s = s + abs( a( i, k ) )*abs( x( i, j ) )
313 40 CONTINUE
314 work( k ) = work( k ) + abs( a( k, k ) )*xk + s
315 50 CONTINUE
316 ELSE
317 DO 70 k = 1, n
318 s = zero
319 xk = abs( x( k, j ) )
320 work( k ) = work( k ) + abs( a( k, k ) )*xk
321 DO 60 i = k + 1, n
322 work( i ) = work( i ) + abs( a( i, k ) )*xk
323 s = s + abs( a( i, k ) )*abs( x( i, j ) )
324 60 CONTINUE
325 work( k ) = work( k ) + s
326 70 CONTINUE
327 END IF
328 s = zero
329 DO 80 i = 1, n
330 IF( work( i ).GT.safe2 ) THEN
331 s = max( s, abs( work( n+i ) ) / work( i ) )
332 ELSE
333 s = max( s, ( abs( work( n+i ) )+safe1 ) /
334 $ ( work( i )+safe1 ) )
335 END IF
336 80 CONTINUE
337 berr( j ) = s
338*
339* Test stopping criterion. Continue iterating if
340* 1) The residual BERR(J) is larger than machine epsilon, and
341* 2) BERR(J) decreased by at least a factor of 2 during the
342* last iteration, and
343* 3) At most ITMAX iterations tried.
344*
345 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
346 $ count.LE.itmax ) THEN
347*
348* Update solution and try again.
349*
350 CALL dpotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
351 CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
352 lstres = berr( j )
353 count = count + 1
354 GO TO 20
355 END IF
356*
357* Bound error from formula
358*
359* norm(X - XTRUE) / norm(X) .le. FERR =
360* norm( abs(inv(A))*
361* ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
362*
363* where
364* norm(Z) is the magnitude of the largest component of Z
365* inv(A) is the inverse of A
366* abs(Z) is the componentwise absolute value of the matrix or
367* vector Z
368* NZ is the maximum number of nonzeros in any row of A, plus 1
369* EPS is machine epsilon
370*
371* The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
372* is incremented by SAFE1 if the i-th component of
373* abs(A)*abs(X) + abs(B) is less than SAFE2.
374*
375* Use DLACN2 to estimate the infinity-norm of the matrix
376* inv(A) * diag(W),
377* where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
378*
379 DO 90 i = 1, n
380 IF( work( i ).GT.safe2 ) THEN
381 work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
382 ELSE
383 work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
384 END IF
385 90 CONTINUE
386*
387 kase = 0
388 100 CONTINUE
389 CALL dlacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
390 $ kase, isave )
391 IF( kase.NE.0 ) THEN
392 IF( kase.EQ.1 ) THEN
393*
394* Multiply by diag(W)*inv(A**T).
395*
396 CALL dpotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
397 DO 110 i = 1, n
398 work( n+i ) = work( i )*work( n+i )
399 110 CONTINUE
400 ELSE IF( kase.EQ.2 ) THEN
401*
402* Multiply by inv(A)*diag(W).
403*
404 DO 120 i = 1, n
405 work( n+i ) = work( i )*work( n+i )
406 120 CONTINUE
407 CALL dpotrs( uplo, n, 1, af, ldaf, work( n+1 ), n, info )
408 END IF
409 GO TO 100
410 END IF
411*
412* Normalize error.
413*
414 lstres = zero
415 DO 130 i = 1, n
416 lstres = max( lstres, abs( x( i, j ) ) )
417 130 CONTINUE
418 IF( lstres.NE.zero )
419 $ ferr( j ) = ferr( j ) / lstres
420*
421 140 CONTINUE
422*
423 RETURN
424*
425* End of DPORFS
426*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dsymv
subroutine dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DSYMV
Definition dsymv.f:152
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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: