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

◆ zgerfs()

subroutine zgerfs ( character  trans,
integer  n,
integer  nrhs,
complex*16, dimension( lda, * )  a,
integer  lda,
complex*16, dimension( ldaf, * )  af,
integer  ldaf,
integer, dimension( * )  ipiv,
complex*16, dimension( ldb, * )  b,
integer  ldb,
complex*16, dimension( ldx, * )  x,
integer  ldx,
double precision, dimension( * )  ferr,
double precision, dimension( * )  berr,
complex*16, dimension( * )  work,
double precision, dimension( * )  rwork,
integer  info 
)

ZGERFS

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

Purpose:
 ZGERFS 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*16 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*16 array, dimension (LDAF,N)
          The factors L and U from the factorization A = P*L*U
          as computed by ZGETRF.
[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 ZGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is COMPLEX*16 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*16 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by ZGETRS.
          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 COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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 zgerfs.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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
198 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
199 $ WORK( * ), X( LDX, * )
200* ..
201*
202* =====================================================================
203*
204* .. Parameters ..
205 INTEGER ITMAX
206 parameter( itmax = 5 )
207 DOUBLE PRECISION ZERO
208 parameter( zero = 0.0d+0 )
209 COMPLEX*16 ONE
210 parameter( one = ( 1.0d+0, 0.0d+0 ) )
211 DOUBLE PRECISION TWO
212 parameter( two = 2.0d+0 )
213 DOUBLE PRECISION THREE
214 parameter( three = 3.0d+0 )
215* ..
216* .. Local Scalars ..
217 LOGICAL NOTRAN
218 CHARACTER TRANSN, TRANST
219 INTEGER COUNT, I, J, K, KASE, NZ
220 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
221 COMPLEX*16 ZDUM
222* ..
223* .. Local Arrays ..
224 INTEGER ISAVE( 3 )
225* ..
226* .. External Functions ..
227 LOGICAL LSAME
228 DOUBLE PRECISION DLAMCH
229 EXTERNAL lsame, dlamch
230* ..
231* .. External Subroutines ..
232 EXTERNAL xerbla, zaxpy, zcopy, zgemv, zgetrs, zlacn2
233* ..
234* .. Intrinsic Functions ..
235 INTRINSIC abs, dble, dimag, max
236* ..
237* .. Statement Functions ..
238 DOUBLE PRECISION CABS1
239* ..
240* .. Statement Function definitions ..
241 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( 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( 'ZGERFS', -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 = dlamch( 'Epsilon' )
292 safmin = dlamch( '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 zcopy( n, b( 1, j ), 1, work, 1 )
310 CALL zgemv( 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 zgetrs( trans, n, 1, af, ldaf, ipiv, work, n, info )
367 CALL zaxpy( 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 ZLACN2 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 zlacn2( 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 zgetrs( 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 zgetrs( 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 ZGERFS
444*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
zgetrs
subroutine zgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
ZGETRS
Definition zgetrs.f:121
zlacn2
subroutine zlacn2(n, v, x, est, kase, isave)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition zlacn2.f:133
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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: