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

CPORFS

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

Purpose: 
 CPORFS 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 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 triangular factor U or L from the Cholesky factorization
          A = U**H*U or A = L*L**H, as computed by CPOTRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[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 CPOTRS.
          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.
Date
November 2011

Definition at line 185 of file cporfs.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  REAL berr( * ), ferr( * ), rwork( * )
197  COMPLEX a( lda, * ), af( ldaf, * ), b( ldb, * ),
198  $ work( * ), x( ldx, * )
199 * ..
200 *
201 * ====================================================================
202 *
203 * .. Parameters ..
204  INTEGER itmax
205  parameter ( itmax = 5 )
206  REAL zero
207  parameter ( zero = 0.0e+0 )
208  COMPLEX one
209  parameter ( one = ( 1.0e+0, 0.0e+0 ) )
210  REAL two
211  parameter ( two = 2.0e+0 )
212  REAL three
213  parameter ( three = 3.0e+0 )
214 * ..
215 * .. Local Scalars ..
216  LOGICAL upper
217  INTEGER count, i, j, k, kase, nz
218  REAL eps, lstres, s, safe1, safe2, safmin, xk
219  COMPLEX zdum
220 * ..
221 * .. Local Arrays ..
222  INTEGER isave( 3 )
223 * ..
224 * .. External Subroutines ..
225  EXTERNAL caxpy, ccopy, chemv, clacn2, cpotrs, xerbla
226 * ..
227 * .. Intrinsic Functions ..
228  INTRINSIC abs, aimag, max, real
229 * ..
230 * .. External Functions ..
231  LOGICAL lsame
232  REAL slamch
233  EXTERNAL lsame, slamch
234 * ..
235 * .. Statement Functions ..
236  REAL cabs1
237 * ..
238 * .. Statement Function definitions ..
239  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( 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( 'CPORFS', -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 = slamch( 'Epsilon' )
281  safmin = slamch( '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 ccopy( n, b( 1, j ), 1, work, 1 )
298  CALL chemv( 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( REAL( 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( REAL( 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 cpotrs( uplo, n, 1, af, ldaf, work, n, info )
360  CALL caxpy( 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 CLACN2 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 clacn2( 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 cpotrs( 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 cpotrs( 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 CPORFS
435 *
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cpotrs
subroutine cpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
CPOTRS
Definition: cpotrs.f:112
chemv
subroutine chemv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CHEMV
Definition: chemv.f:156
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
clacn2
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:135

Here is the call graph for this function:

Here is the caller graph for this function: