LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cherfs()

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

CHERFS

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

Purpose:
 CHERFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian 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 array, dimension (LDA,N)
          The Hermitian 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 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**H or
          A = L*D*L**H as computed by CHETRF.
[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 CHETRF.
[in]B
          B is COMPLEX 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 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by CHETRS.
          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 COMPLEX array, dimension (2*N)
[out]RWORK
          RWORK is REAL 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 cherfs.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  REAL BERR( * ), FERR( * ), RWORK( * )
204  COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
205  $ WORK( * ), X( LDX, * )
206 * ..
207 *
208 * =====================================================================
209 *
210 * .. Parameters ..
211  INTEGER ITMAX
212  parameter( itmax = 5 )
213  REAL ZERO
214  parameter( zero = 0.0e+0 )
215  COMPLEX ONE
216  parameter( one = ( 1.0e+0, 0.0e+0 ) )
217  REAL TWO
218  parameter( two = 2.0e+0 )
219  REAL THREE
220  parameter( three = 3.0e+0 )
221 * ..
222 * .. Local Scalars ..
223  LOGICAL UPPER
224  INTEGER COUNT, I, J, K, KASE, NZ
225  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
226  COMPLEX ZDUM
227 * ..
228 * .. Local Arrays ..
229  INTEGER ISAVE( 3 )
230 * ..
231 * .. External Subroutines ..
232  EXTERNAL caxpy, ccopy, chemv, chetrs, clacn2, xerbla
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC abs, aimag, max, real
236 * ..
237 * .. External Functions ..
238  LOGICAL LSAME
239  REAL SLAMCH
240  EXTERNAL lsame, slamch
241 * ..
242 * .. Statement Functions ..
243  REAL CABS1
244 * ..
245 * .. Statement Function definitions ..
246  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( 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( 'CHERFS', -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 = slamch( 'Epsilon' )
288  safmin = slamch( '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 ccopy( n, b( 1, j ), 1, work, 1 )
305  CALL chemv( 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 ) + abs( real( 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 ) + abs( real( 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 chetrs( uplo, n, 1, af, ldaf, ipiv, work, n, info )
367  CALL caxpy( 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 CLACN2 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 clacn2( 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**H).
411 *
412  CALL chetrs( 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 chetrs( 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 CHERFS
442 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine chemv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CHEMV
Definition: chemv.f:154
subroutine chetrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CHETRS
Definition: chetrs.f:120
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:133
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: