LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zhprfs ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( * )  AP,
complex*16, dimension( * )  AFP,
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 
)

ZHPRFS

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

Purpose: 
 ZHPRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian indefinite
 and packed, 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]AP
          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
          The upper or lower triangle of the Hermitian matrix A, packed
          columnwise in a linear array.  The j-th column of A is stored
          in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
[in]AFP
          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)
          The factored form of the matrix A.  AFP 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 ZHPTRF, stored as a packed
          triangular matrix.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by ZHPTRF.
[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 ZHPTRS.
          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 182 of file zhprfs.f.

182 *
183 * -- LAPACK computational routine (version 3.4.0) --
184 * -- LAPACK is a software package provided by Univ. of Tennessee, --
185 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186 * November 2011
187 *
188 * .. Scalar Arguments ..
189  CHARACTER uplo
190  INTEGER info, ldb, ldx, n, nrhs
191 * ..
192 * .. Array Arguments ..
193  INTEGER ipiv( * )
194  DOUBLE PRECISION berr( * ), ferr( * ), rwork( * )
195  COMPLEX*16 afp( * ), ap( * ), b( ldb, * ), work( * ),
196  $ x( ldx, * )
197 * ..
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202  INTEGER itmax
203  parameter ( itmax = 5 )
204  DOUBLE PRECISION zero
205  parameter ( zero = 0.0d+0 )
206  COMPLEX*16 one
207  parameter ( one = ( 1.0d+0, 0.0d+0 ) )
208  DOUBLE PRECISION two
209  parameter ( two = 2.0d+0 )
210  DOUBLE PRECISION three
211  parameter ( three = 3.0d+0 )
212 * ..
213 * .. Local Scalars ..
214  LOGICAL upper
215  INTEGER count, i, ik, j, k, kase, kk, nz
216  DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
217  COMPLEX*16 zdum
218 * ..
219 * .. Local Arrays ..
220  INTEGER isave( 3 )
221 * ..
222 * .. External Subroutines ..
223  EXTERNAL xerbla, zaxpy, zcopy, zhpmv, zhptrs, zlacn2
224 * ..
225 * .. Intrinsic Functions ..
226  INTRINSIC abs, dble, dimag, max
227 * ..
228 * .. External Functions ..
229  LOGICAL lsame
230  DOUBLE PRECISION dlamch
231  EXTERNAL lsame, dlamch
232 * ..
233 * .. Statement Functions ..
234  DOUBLE PRECISION cabs1
235 * ..
236 * .. Statement Function definitions ..
237  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
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( ldb.LT.max( 1, n ) ) THEN
252  info = -8
253  ELSE IF( ldx.LT.max( 1, n ) ) THEN
254  info = -10
255  END IF
256  IF( info.NE.0 ) THEN
257  CALL xerbla( 'ZHPRFS', -info )
258  RETURN
259  END IF
260 *
261 * Quick return if possible
262 *
263  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
264  DO 10 j = 1, nrhs
265  ferr( j ) = zero
266  berr( j ) = zero
267  10 CONTINUE
268  RETURN
269  END IF
270 *
271 * NZ = maximum number of nonzero elements in each row of A, plus 1
272 *
273  nz = n + 1
274  eps = dlamch( 'Epsilon' )
275  safmin = dlamch( 'Safe minimum' )
276  safe1 = nz*safmin
277  safe2 = safe1 / eps
278 *
279 * Do for each right hand side
280 *
281  DO 140 j = 1, nrhs
282 *
283  count = 1
284  lstres = three
285  20 CONTINUE
286 *
287 * Loop until stopping criterion is satisfied.
288 *
289 * Compute residual R = B - A * X
290 *
291  CALL zcopy( n, b( 1, j ), 1, work, 1 )
292  CALL zhpmv( uplo, n, -one, ap, x( 1, j ), 1, one, work, 1 )
293 *
294 * Compute componentwise relative backward error from formula
295 *
296 * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
297 *
298 * where abs(Z) is the componentwise absolute value of the matrix
299 * or vector Z. If the i-th component of the denominator is less
300 * than SAFE2, then SAFE1 is added to the i-th components of the
301 * numerator and denominator before dividing.
302 *
303  DO 30 i = 1, n
304  rwork( i ) = cabs1( b( i, j ) )
305  30 CONTINUE
306 *
307 * Compute abs(A)*abs(X) + abs(B).
308 *
309  kk = 1
310  IF( upper ) THEN
311  DO 50 k = 1, n
312  s = zero
313  xk = cabs1( x( k, j ) )
314  ik = kk
315  DO 40 i = 1, k - 1
316  rwork( i ) = rwork( i ) + cabs1( ap( ik ) )*xk
317  s = s + cabs1( ap( ik ) )*cabs1( x( i, j ) )
318  ik = ik + 1
319  40 CONTINUE
320  rwork( k ) = rwork( k ) + abs( dble( ap( kk+k-1 ) ) )*
321  $ xk + s
322  kk = kk + k
323  50 CONTINUE
324  ELSE
325  DO 70 k = 1, n
326  s = zero
327  xk = cabs1( x( k, j ) )
328  rwork( k ) = rwork( k ) + abs( dble( ap( kk ) ) )*xk
329  ik = kk + 1
330  DO 60 i = k + 1, n
331  rwork( i ) = rwork( i ) + cabs1( ap( ik ) )*xk
332  s = s + cabs1( ap( ik ) )*cabs1( x( i, j ) )
333  ik = ik + 1
334  60 CONTINUE
335  rwork( k ) = rwork( k ) + s
336  kk = kk + ( n-k+1 )
337  70 CONTINUE
338  END IF
339  s = zero
340  DO 80 i = 1, n
341  IF( rwork( i ).GT.safe2 ) THEN
342  s = max( s, cabs1( work( i ) ) / rwork( i ) )
343  ELSE
344  s = max( s, ( cabs1( work( i ) )+safe1 ) /
345  $ ( rwork( i )+safe1 ) )
346  END IF
347  80 CONTINUE
348  berr( j ) = s
349 *
350 * Test stopping criterion. Continue iterating if
351 * 1) The residual BERR(J) is larger than machine epsilon, and
352 * 2) BERR(J) decreased by at least a factor of 2 during the
353 * last iteration, and
354 * 3) At most ITMAX iterations tried.
355 *
356  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
357  $ count.LE.itmax ) THEN
358 *
359 * Update solution and try again.
360 *
361  CALL zhptrs( uplo, n, 1, afp, ipiv, work, n, info )
362  CALL zaxpy( n, one, work, 1, x( 1, j ), 1 )
363  lstres = berr( j )
364  count = count + 1
365  GO TO 20
366  END IF
367 *
368 * Bound error from formula
369 *
370 * norm(X - XTRUE) / norm(X) .le. FERR =
371 * norm( abs(inv(A))*
372 * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
373 *
374 * where
375 * norm(Z) is the magnitude of the largest component of Z
376 * inv(A) is the inverse of A
377 * abs(Z) is the componentwise absolute value of the matrix or
378 * vector Z
379 * NZ is the maximum number of nonzeros in any row of A, plus 1
380 * EPS is machine epsilon
381 *
382 * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
383 * is incremented by SAFE1 if the i-th component of
384 * abs(A)*abs(X) + abs(B) is less than SAFE2.
385 *
386 * Use ZLACN2 to estimate the infinity-norm of the matrix
387 * inv(A) * diag(W),
388 * where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
389 *
390  DO 90 i = 1, n
391  IF( rwork( i ).GT.safe2 ) THEN
392  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
393  ELSE
394  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
395  $ safe1
396  END IF
397  90 CONTINUE
398 *
399  kase = 0
400  100 CONTINUE
401  CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
402  IF( kase.NE.0 ) THEN
403  IF( kase.EQ.1 ) THEN
404 *
405 * Multiply by diag(W)*inv(A**H).
406 *
407  CALL zhptrs( uplo, n, 1, afp, ipiv, work, n, info )
408  DO 110 i = 1, n
409  work( i ) = rwork( i )*work( i )
410  110 CONTINUE
411  ELSE IF( kase.EQ.2 ) THEN
412 *
413 * Multiply by inv(A)*diag(W).
414 *
415  DO 120 i = 1, n
416  work( i ) = rwork( i )*work( i )
417  120 CONTINUE
418  CALL zhptrs( uplo, n, 1, afp, ipiv, work, n, 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, cabs1( 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 ZHPRFS
437 *
zhptrs
subroutine zhptrs(UPLO, N, NRHS, AP, IPIV, B, LDB, INFO)
ZHPTRS
Definition: zhptrs.f:117
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
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
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
zhpmv
subroutine zhpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
ZHPMV
Definition: zhpmv.f:151
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: