LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DSYMV
Definition: dsymv.f:152
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
Here is the call graph for this function:
Here is the caller graph for this function: