next up previous contents index
Next: Error Bounds for Linear Up: Accuracy and Stability Previous: Improved Error Bounds

Error Bounds for Linear Equation Solving

 

Let Ax=b be the system to be solved, and tex2html_wrap_inline17482 the computed solution. Let n be the dimension of A. An approximate error bound  for tex2html_wrap_inline17482 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:

  1. Suppose that Ax=b is solved using the simple driver PSGESV    (section 3.2.1). Then the approximate error boundgif


    displaymath18202
    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

    tex2html_wrap_inline18242,
    displaymath18203
    Then (to 4 decimal places)
    displaymath18204
    tex2html_wrap_inline18244, tex2html_wrap_inline18246, the true reciprocal condition number tex2html_wrap_inline18248, tex2html_wrap_inline18250, and the true error tex2html_wrap_inline18252.  

  2. 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:  
    displaymath18205
    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, tex2html_wrap_inline18258, tex2html_wrap_inline18260, and the actual error is tex2html_wrap_inline18262.

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 tex2html_wrap_inline17482 be the solution computed by ScaLAPACK (or LAPACK) using any of their linear equation solvers. Let r be the residual tex2html_wrap_inline18270. In the absence of rounding error, r would be zero and tex2html_wrap_inline17482 would equal x; with rounding error, one can only say the following:

The normwise backward error of the computed solution tex2html_wrap_inline17482,     with respect to the infinity norm, is the pair E,f, which minimizes
displaymath18206
subject to the constraint tex2html_wrap_inline18075. The minimal value of tex2html_wrap_inline18284 is given by
displaymath18207
One can show that the computed solution tex2html_wrap_inline17482 satisfies tex2html_wrap_inline18288, where p(n) is a modestly growing function of n. The corresponding condition number is tex2html_wrap_inline18294.   The error tex2html_wrap_inline18296 is bounded by
displaymath18208
In the first code fragment in the preceding section, tex2html_wrap_inline18300, which is tex2html_wrap_inline18302 in the numerical example, is approximated by tex2html_wrap_inline18304. Approximations  of tex2html_wrap_inline18121 -- 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 tex2html_wrap_inline18308 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 tex2html_wrap_inline18314, resulting in tex2html_wrap_inline18316.

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 tex2html_wrap_inline17482 with the following properties:

The componentwise backward error of the computed solution tex2html_wrap_inline17482 is the pair E,f which minimizes    
displaymath18209
(where we interpret 0/0 as 0) subject to the constraint tex2html_wrap_inline18075. The minimal value of tex2html_wrap_inline18330 is given by
displaymath18210
One can show that for most problems the tex2html_wrap_inline17482 computed by PxyySVX satisfies tex2html_wrap_inline18334, where p(n) is a modestly growing function of n. In other words, tex2html_wrap_inline17482 is the exact solution of the perturbed problem tex2html_wrap_inline18075, where E and f are small relative perturbations in each entry of A and b, respectively. The corresponding condition number is tex2html_wrap_inline18352.   The error tex2html_wrap_inline18296 is bounded by
displaymath18211

The routines PxyyRFS and PxyySVX return tex2html_wrap_inline18356,                              which is called BERR  (for Backward ERRor), and a bound on the the actual error tex2html_wrap_inline18358, 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 tex2html_wrap_inline18360 given above:
displaymath18212
Here, tex2html_wrap_inline18362 is the computed value of the residual tex2html_wrap_inline18364, 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 tex2html_wrap_inline18366.

Even in the rare cases where PxyyRFS fails to make BERR close to its minimum tex2html_wrap_inline17202, the error bound FERR may remain small. See [9] for details.


next up previous contents index
Next: Error Bounds for Linear Up: Accuracy and Stability Previous: Improved Error Bounds

Susan Blackford
Tue May 13 09:21:01 EDT 1997