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

◆ dsyrfs()

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.

Definition at line 189 of file dsyrfs.f.

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