LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sgerfs()

subroutine sgerfs ( character  TRANS,
integer  N,
integer  NRHS,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  FERR,
real, dimension( * )  BERR,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

SGERFS

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

Purpose:
 SGERFS 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 REAL 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 REAL array, dimension (LDAF,N)
          The factors L and U from the factorization A = P*L*U
          as computed by SGETRF.
[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 SGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is REAL 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 REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SGETRS.
          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 REAL 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
December 2016

Definition at line 187 of file sgerfs.f.

187 *
188 * -- LAPACK computational routine (version 3.7.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 * December 2016
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  REAL 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  REAL zero
209  parameter( zero = 0.0e+0 )
210  REAL one
211  parameter( one = 1.0e+0 )
212  REAL two
213  parameter( two = 2.0e+0 )
214  REAL three
215  parameter( three = 3.0e+0 )
216 * ..
217 * .. Local Scalars ..
218  LOGICAL notran
219  CHARACTER transt
220  INTEGER count, i, j, k, kase, nz
221  REAL eps, lstres, s, safe1, safe2, safmin, xk
222 * ..
223 * .. Local Arrays ..
224  INTEGER isave( 3 )
225 * ..
226 * .. External Subroutines ..
227  EXTERNAL saxpy, scopy, sgemv, sgetrs, slacn2, xerbla
228 * ..
229 * .. Intrinsic Functions ..
230  INTRINSIC abs, max
231 * ..
232 * .. External Functions ..
233  LOGICAL lsame
234  REAL slamch
235  EXTERNAL lsame, slamch
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( 'SGERFS', -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 = slamch( 'Epsilon' )
284  safmin = slamch( '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 scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
302  CALL sgemv( 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 sgetrs( trans, n, 1, af, ldaf, ipiv, work( n+1 ), n,
359  $ info )
360  CALL saxpy( 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 SLACN2 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 slacn2( 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 sgetrs( 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 sgetrs( 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 SGERFS
437 *
subroutine sgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
SGETRS
Definition: sgetrs.f:123
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:91
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:138
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
Here is the call graph for this function:
Here is the caller graph for this function: