LAPACK  3.10.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.

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 *
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 cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:121
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
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: