LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ssyrfs()

subroutine ssyrfs ( character  UPLO,
integer  N,
integer  NRHS,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

SSYRFS

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

Purpose:
 SSYRFS 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 REAL 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 REAL 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 SSYTRF.
[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 SSYTRF.
[in]B
          B is REAL 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 REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SSYTRS.
          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 REAL 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 REAL 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 REAL 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 ssyrfs.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  REAL 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  REAL ZERO
212  parameter( zero = 0.0e+0 )
213  REAL ONE
214  parameter( one = 1.0e+0 )
215  REAL TWO
216  parameter( two = 2.0e+0 )
217  REAL THREE
218  parameter( three = 3.0e+0 )
219 * ..
220 * .. Local Scalars ..
221  LOGICAL UPPER
222  INTEGER COUNT, I, J, K, KASE, NZ
223  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
224 * ..
225 * .. Local Arrays ..
226  INTEGER ISAVE( 3 )
227 * ..
228 * .. External Subroutines ..
229  EXTERNAL saxpy, scopy, slacn2, ssymv, ssytrs, xerbla
230 * ..
231 * .. Intrinsic Functions ..
232  INTRINSIC abs, max
233 * ..
234 * .. External Functions ..
235  LOGICAL LSAME
236  REAL SLAMCH
237  EXTERNAL lsame, slamch
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( 'SSYRFS', -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 = slamch( 'Epsilon' )
279  safmin = slamch( '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 scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
296  CALL ssymv( 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 ssytrs( uplo, n, 1, af, ldaf, ipiv, work( n+1 ), n,
359  $ info )
360  CALL saxpy( 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 SLACN2 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 slacn2( 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 ssytrs( 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 ssytrs( 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 SSYRFS
437 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:136
subroutine ssytrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
SSYTRS
Definition: ssytrs.f:120
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine ssymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SSYMV
Definition: ssymv.f:152
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: