The GSVD algorithm used in LAPACK ([12][10][62]) is backward stable:

Let the computed GSVD ofAandBbe and . This is nearly the exact GSVD ofA+EandB+Fin the following sense.EandFare small:there exist small , , and such that , , and are exactly orthogonal (or unitary):

and

is the exact GSVD of

A+EandB+F. Herep(n) is a modestly growing function ofn, and we takep(n) = 1 in the above code fragment.Let and be the square roots of the diagonal entries of the exact and , and let and the square roots of the diagonal entries of the computed and . Let

Then provided

Gand have full rankn, one can show [61][74] thatIn the code fragment we approximate the numerator of the last expression by and approximate the denominator by in order to compute

SERRBD; STRCON returns an approximationRCONDto .

We assume that the rank `r` of `G` equals `n`, because otherwise the
s and s are not well determined. For example, if

then `A` and `B` have
and , whereas
and have
and , which
are completely different, even though and
. In this case, ,
so `G` is nearly rank-deficient.

The reason the code fragment assumes `m` > = `n` is that in this case is
stored overwritten on `A`, and can be passed to STRCON in order to compute
`RCOND`. If `m` < = `n`, then the
first `m` rows of are
stored in `A`, and the last `m` - `n` rows of are stored in `B`. This
complicates the computation of `RCOND`: either must be copied to
a single array before calling STRCON, or else the lower level subroutine SLACON
must be used with code capable of solving linear equations with
and as coefficient matrices.

Tue Nov 29 14:03:33 EST 1994