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

◆ cporfs()

subroutine cporfs ( character  uplo,
integer  n,
integer  nrhs,
complex, dimension( lda, * )  a,
integer  lda,
complex, dimension( ldaf, * )  af,
integer  ldaf,
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 
)

CPORFS

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

Purpose:
 CPORFS improves the computed solution to a system of linear
 equations when the coefficient matrix is Hermitian positive definite,
 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]A
          A is COMPLEX array, dimension (LDA,N)
          The Hermitian matrix A.  If UPLO = 'U', the leading N-by-N
          upper triangular part of A contains the upper triangular part
          of the matrix A, and the strictly lower triangular part of A
          is not referenced.  If UPLO = 'L', the leading N-by-N lower
          triangular part of A contains the lower triangular part of
          the matrix A, and the strictly upper triangular part of A is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX array, dimension (LDAF,N)
          The triangular factor U or L from the Cholesky factorization
          A = U**H*U or A = L*L**H, as computed by CPOTRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[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 CPOTRS.
          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 181 of file cporfs.f.

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