LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zpprfs()

subroutine zpprfs ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( * )  AP,
complex*16, dimension( * )  AFP,
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 
)

ZPPRFS

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

Purpose:
 ZPPRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian positive definite
 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)*(2n-j)/2) = A(i,j) for j<=i<=n.
[in]AFP
          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)
          The triangular factor U or L from the Cholesky factorization
          A = U**H*U or A = L*L**H, as computed by DPPTRF/ZPPTRF,
          packed columnwise in a linear array in the same format as A
          (see AP).
[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 ZPPTRS.
          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 169 of file zpprfs.f.

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