LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ chprfs()

subroutine chprfs ( character  UPLO,
integer  N,
integer  NRHS,
complex, dimension( * )  AP,
complex, dimension( * )  AFP,
integer, dimension( * )  IPIV,
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 
)

CHPRFS

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

Purpose:
 CHPRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian indefinite
 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
[in]AFP
          AFP is COMPLEX array, dimension (N*(N+1)/2)
          The factored form of the matrix A.  AFP contains the block
          diagonal matrix D and the multipliers used to obtain the
          factor U or L from the factorization A = U*D*U**H or
          A = L*D*L**H as computed by CHPTRF, stored as a packed
          triangular matrix.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by CHPTRF.
[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 CHPTRS.
          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 178 of file chprfs.f.

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