LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zsyrfs ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZSYRFS

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

Purpose: 
 ZSYRFS 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 COMPLEX*16 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 COMPLEX*16 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 ZSYTRF.
[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 ZSYTRF.
[in]B
          B is COMPLEX*16 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 COMPLEX*16 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by ZSYTRS.
          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 COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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 194 of file zsyrfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: