LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dsyrfs ( character  UPLO,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
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 
)

DSYRFS

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

Purpose: 
 DSYRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric indefinite, 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 factored form of the matrix A.  AF contains the block
          diagonal matrix D and the multipliers used to obtain the
          factor U or L from the factorization A = U*D*U**T or
          A = L*D*L**T as computed by DSYTRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by DSYTRF.
[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 DSYTRS.
          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 193 of file dsyrfs.f.

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