LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
sgetrs
subroutine sgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
SGETRS
Definition sgetrs.f:121
slacn2
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
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: