LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
chpmv
subroutine chpmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
CHPMV
Definition chpmv.f:149
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:133
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cpptrs
subroutine cpptrs(uplo, n, nrhs, ap, b, ldb, info)
CPPTRS
Definition cpptrs.f:108
Here is the call graph for this function:
Here is the caller graph for this function: