LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:121
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
Here is the call graph for this function:
Here is the caller graph for this function: