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

CPPRFS

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

Purpose: 
 CPPRFS 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 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 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 SPPTRF/CPPTRF,
          packed columnwise in a linear array in the same format as A
          (see AP).
[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 CPPTRS.
          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 173 of file cpprfs.f.

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