LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zporfs ( character  UPLO,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
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 
)

ZPORFS

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

Purpose: 
 ZPORFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian positive definite,
 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 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*16 array, dimension (LDAF,N)
          The triangular factor U or L from the Cholesky factorization
          A = U**H*U or A = L*L**H, as computed by ZPOTRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[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 ZPOTRS.
          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 185 of file zporfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: