LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
December 2016

Definition at line 188 of file cgerfs.f.

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