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

◆ 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*
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
chptrs
subroutine chptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
CHPTRS
Definition chptrs.f:115
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
Here is the call graph for this function:
Here is the caller graph for this function: