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

◆ zsyrfs()

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.

Definition at line 190 of file zsyrfs.f.

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