LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cpprfs()

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.

Definition at line 169 of file cpprfs.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  REAL BERR( * ), FERR( * ), RWORK( * )
182  COMPLEX AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
183  $ X( LDX, * )
184 * ..
185 *
186 * ====================================================================
187 *
188 * .. Parameters ..
189  INTEGER ITMAX
190  parameter( itmax = 5 )
191  REAL ZERO
192  parameter( zero = 0.0e+0 )
193  COMPLEX CONE
194  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
195  REAL TWO
196  parameter( two = 2.0e+0 )
197  REAL THREE
198  parameter( three = 3.0e+0 )
199 * ..
200 * .. Local Scalars ..
201  LOGICAL UPPER
202  INTEGER COUNT, I, IK, J, K, KASE, KK, NZ
203  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
204  COMPLEX ZDUM
205 * ..
206 * .. Local Arrays ..
207  INTEGER ISAVE( 3 )
208 * ..
209 * .. External Subroutines ..
210  EXTERNAL caxpy, ccopy, chpmv, clacn2, cpptrs, xerbla
211 * ..
212 * .. Intrinsic Functions ..
213  INTRINSIC abs, aimag, max, real
214 * ..
215 * .. External Functions ..
216  LOGICAL LSAME
217  REAL SLAMCH
218  EXTERNAL lsame, slamch
219 * ..
220 * .. Statement Functions ..
221  REAL CABS1
222 * ..
223 * .. Statement Function definitions ..
224  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( 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( 'CPPRFS', -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 = slamch( 'Epsilon' )
262  safmin = slamch( '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 ccopy( n, b( 1, j ), 1, work, 1 )
279  CALL chpmv( 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( real( 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( real( 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 cpptrs( uplo, n, 1, afp, work, n, info )
349  CALL caxpy( 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 CLACN2 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 clacn2( 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 cpptrs( 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 cpptrs( 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 CPPRFS
424 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine chpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
CHPMV
Definition: chpmv.f:149
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:133
subroutine cpptrs(UPLO, N, NRHS, AP, B, LDB, INFO)
CPPTRS
Definition: cpptrs.f:108
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: