LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dgerfs ( character  TRANS,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DGERFS

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

Purpose: 
 DGERFS 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 = 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDAF,N)
          The factors L and U from the factorization A = P*L*U
          as computed by DGETRF.
[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 DGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DGETRS.
          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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (3*N)
[out]IWORK
          IWORK is INTEGER 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
November 2011

Definition at line 187 of file dgerfs.f.

187 *
188 * -- LAPACK computational routine (version 3.4.0) --
189 * -- LAPACK is a software package provided by Univ. of Tennessee, --
190 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191 * November 2011
192 *
193 * .. Scalar Arguments ..
194  CHARACTER trans
195  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs
196 * ..
197 * .. Array Arguments ..
198  INTEGER ipiv( * ), iwork( * )
199  DOUBLE PRECISION a( lda, * ), af( ldaf, * ), b( ldb, * ),
200  $ berr( * ), ferr( * ), work( * ), x( ldx, * )
201 * ..
202 *
203 * =====================================================================
204 *
205 * .. Parameters ..
206  INTEGER itmax
207  parameter ( itmax = 5 )
208  DOUBLE PRECISION zero
209  parameter ( zero = 0.0d+0 )
210  DOUBLE PRECISION one
211  parameter ( one = 1.0d+0 )
212  DOUBLE PRECISION two
213  parameter ( two = 2.0d+0 )
214  DOUBLE PRECISION three
215  parameter ( three = 3.0d+0 )
216 * ..
217 * .. Local Scalars ..
218  LOGICAL notran
219  CHARACTER transt
220  INTEGER count, i, j, k, kase, nz
221  DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
222 * ..
223 * .. Local Arrays ..
224  INTEGER isave( 3 )
225 * ..
226 * .. External Subroutines ..
227  EXTERNAL daxpy, dcopy, dgemv, dgetrs, dlacn2, xerbla
228 * ..
229 * .. Intrinsic Functions ..
230  INTRINSIC abs, max
231 * ..
232 * .. External Functions ..
233  LOGICAL lsame
234  DOUBLE PRECISION dlamch
235  EXTERNAL lsame, dlamch
236 * ..
237 * .. Executable Statements ..
238 *
239 * Test the input parameters.
240 *
241  info = 0
242  notran = lsame( trans, 'N' )
243  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
244  $ lsame( trans, 'C' ) ) 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 = -10
256  ELSE IF( ldx.LT.max( 1, n ) ) THEN
257  info = -12
258  END IF
259  IF( info.NE.0 ) THEN
260  CALL xerbla( 'DGERFS', -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  IF( notran ) THEN
275  transt = 'T'
276  ELSE
277  transt = 'N'
278  END IF
279 *
280 * NZ = maximum number of nonzero elements in each row of A, plus 1
281 *
282  nz = n + 1
283  eps = dlamch( 'Epsilon' )
284  safmin = dlamch( 'Safe minimum' )
285  safe1 = nz*safmin
286  safe2 = safe1 / eps
287 *
288 * Do for each right hand side
289 *
290  DO 140 j = 1, nrhs
291 *
292  count = 1
293  lstres = three
294  20 CONTINUE
295 *
296 * Loop until stopping criterion is satisfied.
297 *
298 * Compute residual R = B - op(A) * X,
299 * where op(A) = A, A**T, or A**H, depending on TRANS.
300 *
301  CALL dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
302  CALL dgemv( trans, n, n, -one, a, lda, x( 1, j ), 1, one,
303  $ work( n+1 ), 1 )
304 *
305 * Compute componentwise relative backward error from formula
306 *
307 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
308 *
309 * where abs(Z) is the componentwise absolute value of the matrix
310 * or vector Z. If the i-th component of the denominator is less
311 * than SAFE2, then SAFE1 is added to the i-th components of the
312 * numerator and denominator before dividing.
313 *
314  DO 30 i = 1, n
315  work( i ) = abs( b( i, j ) )
316  30 CONTINUE
317 *
318 * Compute abs(op(A))*abs(X) + abs(B).
319 *
320  IF( notran ) THEN
321  DO 50 k = 1, n
322  xk = abs( x( k, j ) )
323  DO 40 i = 1, n
324  work( i ) = work( i ) + abs( a( i, k ) )*xk
325  40 CONTINUE
326  50 CONTINUE
327  ELSE
328  DO 70 k = 1, n
329  s = zero
330  DO 60 i = 1, n
331  s = s + abs( a( i, k ) )*abs( x( i, j ) )
332  60 CONTINUE
333  work( k ) = work( k ) + s
334  70 CONTINUE
335  END IF
336  s = zero
337  DO 80 i = 1, n
338  IF( work( i ).GT.safe2 ) THEN
339  s = max( s, abs( work( n+i ) ) / work( i ) )
340  ELSE
341  s = max( s, ( abs( work( n+i ) )+safe1 ) /
342  $ ( work( 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 dgetrs( trans, n, 1, af, ldaf, ipiv, work( n+1 ), n,
359  $ info )
360  CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
361  lstres = berr( j )
362  count = count + 1
363  GO TO 20
364  END IF
365 *
366 * Bound error from formula
367 *
368 * norm(X - XTRUE) / norm(X) .le. FERR =
369 * norm( abs(inv(op(A)))*
370 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
371 *
372 * where
373 * norm(Z) is the magnitude of the largest component of Z
374 * inv(op(A)) is the inverse of op(A)
375 * abs(Z) is the componentwise absolute value of the matrix or
376 * vector Z
377 * NZ is the maximum number of nonzeros in any row of A, plus 1
378 * EPS is machine epsilon
379 *
380 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
381 * is incremented by SAFE1 if the i-th component of
382 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
383 *
384 * Use DLACN2 to estimate the infinity-norm of the matrix
385 * inv(op(A)) * diag(W),
386 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
387 *
388  DO 90 i = 1, n
389  IF( work( i ).GT.safe2 ) THEN
390  work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
391  ELSE
392  work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
393  END IF
394  90 CONTINUE
395 *
396  kase = 0
397  100 CONTINUE
398  CALL dlacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
399  $ kase, isave )
400  IF( kase.NE.0 ) THEN
401  IF( kase.EQ.1 ) THEN
402 *
403 * Multiply by diag(W)*inv(op(A)**T).
404 *
405  CALL dgetrs( transt, n, 1, af, ldaf, ipiv, work( n+1 ),
406  $ n, info )
407  DO 110 i = 1, n
408  work( n+i ) = work( i )*work( n+i )
409  110 CONTINUE
410  ELSE
411 *
412 * Multiply by inv(op(A))*diag(W).
413 *
414  DO 120 i = 1, n
415  work( n+i ) = work( i )*work( n+i )
416  120 CONTINUE
417  CALL dgetrs( trans, n, 1, af, ldaf, ipiv, work( n+1 ), n,
418  $ info )
419  END IF
420  GO TO 100
421  END IF
422 *
423 * Normalize error.
424 *
425  lstres = zero
426  DO 130 i = 1, n
427  lstres = max( lstres, abs( x( i, j ) ) )
428  130 CONTINUE
429  IF( lstres.NE.zero )
430  $ ferr( j ) = ferr( j ) / lstres
431 *
432  140 CONTINUE
433 *
434  RETURN
435 *
436 * End of DGERFS
437 *
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
dgetrs
subroutine dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
DGETRS
Definition: dgetrs.f:123
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dlacn2
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:138

Here is the call graph for this function:

Here is the caller graph for this function: