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

◆ cgerfs()

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

CGERFS

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

Purpose:
 CGERFS improves the computed solution to a system of linear
 equations and provides error bounds and backward error estimates for
 the solution.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations:
          = 'N':  A * X = B     (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose)
[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 original N-by-N matrix A.
[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 factors L and U from the factorization A = P*L*U
          as computed by CGETRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices from CGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[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 CGETRS.
          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 184 of file cgerfs.f.

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