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 pairE,f, which minimizes

subject to the constraint . The minimal value of is given by

One can show that the computed solution satisfies , wherep(n) is a modestly growing function ofn. 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 reciprocalRCOND-- are returned by computational routines PxyyCON (section 3.3.1) or driver routines PxyySVX (section 3.2.1). The code fragment makes sureRCONDis at leastEPSMCHto avoid overflow in computingERRBD. This limitsERRBDto 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 ofRCONDreturned by PxyySVX may apply to a linear system obtained fromAx=bbyequilibration, namely, scaling the rows and columns ofAin 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 pairE,fwhich 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 , wherep(n) is a modestly growing function ofn. In other words, is the exact solution of the perturbed problem , whereEandfare small relative perturbations in each entry ofAandb, 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 , calledFERR(for Forward ERRor), as in the second code fragment in the last section.FERRis 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 forRCOND.The value of

BERRfor the example in the preceding section is .Even in the rare cases where PxyyRFS fails to make

BERRclose to its minimum , the error boundFERRmay remain small. See [9] for details.

Tue May 13 09:21:01 EDT 1997