LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ zporfs()

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.

Definition at line 181 of file zporfs.f.

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