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

Definition at line 183 of file sgerfs.f.

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