LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 185 of file dporfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: