|
|
Next: Error Bounds for Linear
Up: Accuracy and Stability
Previous: Improved Error Bounds
Let Ax=b be the system to be solved, and the computed
solution. Let n be the dimension of A.
An approximate error bound
for may be obtained in one of the following two ways,
depending on whether the solution is computed by a simple driver or
an expert driver:
- Suppose that Ax=b is solved using the simple driver PSGESV
(section 3.2.1).
Then the approximate error bound
can be computed by the following code fragment.
EPSMCH = PSLAMCH( ICTXT, 'E' )
* Get infinity-norm of A
ANORM = PSLANGE( 'I', N, N, A, IA, JA, DESCA, WORK )
* Solve system; The solution X overwrites B
CALL PSGESV( N, 1, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO )
IF( INFO.GT.0 ) THEN
PRINT *,'Singular Matrix'
ELSE IF( N.GT.0 ) THEN
* Get reciprocal condition number RCOND of A
CALL PSGECON( 'I', N, A, IA, JA, DESCA, ANORM, RCOND, WORK,
$ LWORK, IWORK, LIWORK, INFO )
RCOND = MAX( RCOND, EPSMCH )
ERRBD = EPSMCH / RCOND
END IF
For example, suppose
,
Then (to 4 decimal places)
,
,
the true reciprocal condition number ,
, and the true error
.
- Suppose that Ax=b is solved using the expert driver PSGESVX
(section 3.2.1).
This routine provides an explicit error bound FERR, measured
with the infinity-norm:
For example, the following code fragment solves
Ax=b and computes an approximate error bound FERR:
CALL PSGESVX( 'E', 'N', N, 1, A, IA, JA, DESCA, AF, IAF, JAF,
$ DESCAF, IPIV, EQUED, R, C, B, IB, JB, DESCB, X, IX,
$ JX, DESCX, RCOND, FERR, BERR, WORK, LWORK, IWORK,
$ LIWORK, INFO )
IF( INFO.GT.0 ) PRINT *,'(Nearly) Singular Matrix'
For the same A and b as above,
,
,
and the actual error is .
This example illustrates
that the expert driver provides an error bound with less programming
effort than the simple driver, and also that it may produce a significantly
more accurate answer.
Similar code fragments, with obvious adaptations,
may be used with all the driver routines PxPOSV and PxPOSVX
in Table 3.2.
For example, if a symmetric positive definite or Hermitian positive
definite system is solved by using the simple driver PxPOSV,
then PxLANSY
or PxLANHE, respectively, must
be used to compute ANORM, and
PxPOCON must
be used to compute RCOND.
The drivers
PxGBSV
(for solving general band matrices with partial pivoting),
PxPBSV
(for solving positive definite band matrices) and
PxPTSV
(for solving positive definite tridiagonal matrices),
do not yet have the corresponding routines needed to compute
error bounds, namely,
PxLAnHE to compute ANORM and PxyyCON to compute RCOND.
The drivers
PxDBSV
(for solving general band matrices) and
PxDTSV
(for solving general tridiagonal matrices) do not pivot for numerical
stability, and so may be faster but less accurate than their
pivoting counterparts above. These routines may be used safely when
any diagonal pivot sequence leads to a stable factorization;
diagonally dominant matrices and symmetric positive definite matrices
[71] have this property, for example.
Further Details: Error Bounds for Linear Equation Solving
The conventional error analysis of linear
equation solving goes as follows.
Let Ax=b be the system to be solved. Let be the solution
computed by ScaLAPACK (or LAPACK) using any of their linear equation solvers.
Let r be
the residual . In the absence of rounding error, r
would be zero and would equal x; with rounding error, one can
only say the following:
The normwise backward error of the computed solution ,
with respect to the infinity norm,
is the pair E,f, which minimizes
subject to the constraint .
The minimal value of
is given by
One can show that the computed solution
satisfies ,
where p(n) is a modestly growing function of n.
The corresponding condition number is
.
The error is bounded by
In the first code fragment in the preceding section, ,
which is in the numerical example,
is approximated by .
Approximations
of -- or, strictly speaking, its reciprocal RCOND --
are returned by computational routines
PxyyCON (section 3.3.1) or driver routines
PxyySVX (section 3.2.1). The code fragment
makes sure RCOND is at least EPSMCH to
avoid overflow in computing
ERRBD.
This limits
ERRBD to a maximum of 1, which is no loss of generality because
a relative error of 1 or more indicates the same thing:
a complete loss of accuracy.
Note that the
value of RCOND returned by PxyySVX may apply to a linear
system obtained from Ax=b by equilibration, namely,
scaling the rows and columns of A in order to make the
condition number smaller. This is the case in the second
code fragment in the preceding section, where the program
chose to scale the rows by the factors returned in
,
resulting in .
As stated in section 6.4.2,
this approach does not respect the presence
of zero or tiny entries in A. In contrast,
the ScaLAPACK computational routines
PxyyRFS (section 3.3.1) or driver routines PxyySVX
(section 3.2.1) will (except in rare cases)
compute a solution with the following properties:
The componentwise backward error
of the computed solution is the pair E,f which minimizes
(where we interpret 0/0 as 0)
subject to the constraint .
The minimal value of
is given by
One can show that for most problems the computed by PxyySVX
satisfies ,
where p(n) is a modestly growing function of n.
In other words, is the exact solution of the
perturbed problem ,
where E and f are small relative perturbations in each entry of A and
b, respectively.
The corresponding condition number is
.
The error is bounded by
The routines PxyyRFS and PxyySVX return ,
which is called BERR
(for Backward ERRor),
and a bound on the the actual error
, called FERR
(for Forward ERRor), as
in the second code fragment in the last section.
FERR is actually calculated by the following formula, which can
be smaller than the bound given above:
Here, is the computed value of the residual , and
the norm in the numerator is estimated by using the same estimation
subroutine used for RCOND.
The value of
BERR for the example in the preceding section is .
Even in the rare cases where PxyyRFS fails to make
BERR close to its minimum , the error bound FERR
may remain small. See [9]
for details.
Next: Error Bounds for Linear
Up: Accuracy and Stability
Previous: Improved Error Bounds
Susan Blackford
Tue May 13 09:21:01 EDT 1997
|