# Structure preserving deflation of infinite eigenvalues in structured pencils.

1. Introduction. In this paper we develop numerical methods to calculate a structure preserving deflation under unitary (or in the real case real orthogonal) equivalence transformations for the infinite eigenvalue part in eigenvalue problems for structured matrix pencils

(1.1) ([lambda]N - M)x = 0,

where N = [[sigma].sub.N][N.sup.*], M = [[sigma].sub.M][M.sup.*][member of] [K.sup.n,n], [[sigma].sub.N], [[sigma].sub.M] [member of] {[+ or -]1}, and * is either T, the transpose, or *, the conjugate transpose. Here K stands for either R or C.

Eigenvalue problems of this form arise in many applications, in particular in the context of linear quadratic optimal control problems (see, e.g., [9, 12, 15, 17, 19]), [H.sub.[infinity]] control problems (see, e.g., [1, 7, 14, 20]), and other applications; see, e.g., [10, 13].

In most applications one needs information about the Kronecker structure (eigenvalues, eigenvectors, multiplicities, Jordan chains, sign characteristic) of these pencils; see [18]. Numerically, in order to obtain such information, structure preserving methods (using congruence rather than equivalence) are preferred because only these preserve all possible characteristics and they reveal the symmetries in the Kronecker structure of a structured pencil of the form (1.1). A structure preserving method may also save computational time and work since it can make use of the symmetries in the matrix structures to simplify the computation. In order to compute the Kronecker structure of (1.1) in a numerically backward stable way, a structured staircase form under unitary structure preserving transformations has been derived in [5] and implemented as production software in [4]. In the real case, the unitary transformations can be chosen to be real orthogonal, but to simplify, in the following when speaking of unitary transformations, we implicitly mean real orthogonal transformations in the case of real matrices. The following result summarizes the staircase form.

Theorem 1.1 (Structured staircase form). For a pencil [lambda]N - M with N = [[sigma].sub.N]N*, M = [[sigma].sub.M][M.sup.*] [member of] [K.sup.n,n], there exists a real orthogonal or unitary matrix U [member of] [K.sup.n,n], respectively, such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(1.2)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [q.sub.1] [greater than or equal to] [n.sub.1] [greater than or equal to] [q.sub.2] [greater than or equal to] [n.sub.2] [greater than or equal to] ... [greater than or equal to] [q.sub.m] [greater than or equal to] [n.sub.m],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the blocks [[SIGMA].sub.22], [DELTA], and [[GAMMA].sub.j], j = 1, ..., m, are nonsingular.

Proof. The proof for the real case, [[sigma].sub.N] = -1, and [[sigma].sub.M] = 1 has been given in [5]; the other cases follow analogously.

It is clear that the transformation introduced in Theorem 1.1 preserves the structure of the pencil, but note that in the real case or in the complex case with * being the complex conjugate, this transformation is a congruence transformation, while in the complex case with * being the transpose, this is just a structure preserving equivalence transformation but not a congruence transformation.

This staircase form allows to deflate the singular part and some of the infinite eigenvalue parts of the pencil in a structure preserving way so that for eigenvalue, eigenvector, and invariant subspace computations, only the central block [lambda][N.sub.m+1,m+1] - [M.sub.m+1,m+1] of the form (1.2) has to be considered, which is regular and of index at most one, i.e., it has only finite eigenvalues plus infinite eigenvalues with Kronecker blocks of size one due to the fact that [DELTA] and [[SIGMA].sub.22] are nonsingular.

The staircase form has recently been extended to parameter dependent matrix pencils arising in the control of differential-algebraic systems with variable coefficients [9], where, however, the transformation is more complex.

From the staircase form it is clear how to deflate all the parts associated with higher (than one) indices and singular parts using unitary transformations (despite the usual difficulties with rank decisions). However, for a long time it remained an open problem how to deflate the remaining index one part associated with the eigenvalue infinity in a structure preserving way using unitary structure preserving transformations, i.e., to reduce the central subpencil to a subpencil [lambda][N.sub.f] - [M.sub.f] of the same structure (with a nonsingular [N.sub.f]) that contains all the information about the finite eigenvalues. In this paper we solve this problem and develop a new technique for structure preserving deflation of infinite eigenvalues via unitary transformations. We present an error and perturbation analysis and also several numerical examples that demonstrate the properties of the new method and also compare it to non-unitary structure preserving transformations.

2. Deflation of the index one part. A simple way to achieve a structure preserving but non-unitary deflation is to use the Schur complement. Compute a factorization

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

of N with [N.sub.1] of full row rank. This can be done via the rank-revealing QR decomposition or the singular value decomposition; see [8]. The critical part in this factorization is the decision about the rank of N, which is done in the usual way by setting those part to zero that are associated to the singular values which are smaller than the machine precision times the norm of N; see [6] for a detailed discussion of this topic.

Using the structure of the pencil, we have

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. By construction, [[??].sub.11] is invertible, and since [[lambda].sub.N] - M is regular and of index at most one, [[??].sub.22] is also invertible. Then, with

(2.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

one has

(2.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with the Schur complement

(2.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Hence, [lambda]N - M is equivalent to the decoupled block diagonal pencil

(A[[??].sub.11] - S) [direct sum] ([lambda]0 - [[??].sub.22]).

This deflation procedure via a non-unitary transformation preserves the structure, but if [[??].sub.22] is ill-conditioned with respect to inversion, then the eigenvalues of the pencil [lambda][[??].sub.11]-S may be corrupted due to the ill-conditioning in the computation of the Schur complement S. We analyze the properties of the Schur complement approach in Section 4 and present some numerical examples in Section 7.

Another possibility to deflate the part associated with the infinite eigenvalues is to use equivalence transformations that are unitary but not structure preserving. Starting again with the factorization (2.1), we form [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] partitioned analogously and let V be a unitary matrix such that

(2.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [[??].sub.22] square nonsingular. Setting

(2.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

then, since [[??].sub.11] and [[??].sub.22] in (2.7) are nonsingular, we can eliminate [[??].sub.12] and [[??].sub.12] simultaneously by block Gaussian eliminations. It follows that [lambda]N-M is equivalent to the decoupled block diagonal pencil

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The computation of the subpencil [lambda] [[??].sub.11]-[[??].sub.11] can be performed in a backward stable way, but it has lost its symmetry structures, and as a consequence it is hard to obtain the resulting symmetric structure in the finite spectrum of [lambda]N-M in further computations. Thus, this method is inadequate for many of the tasks where the exact eigenvalue symmetry is needed. We present some numerical examples to demonstrate this deficiency in Section 7. Furthermore, in the resulting pencil, the sign characteristics, which are invariants of the pencil under structure preserving transformations [18], are not preserved.

3. Deflation under unitary structure preserving transformations. To derive a structure preserving deflation under unitary structure preserving transformations, consider the unitary matrices U, V obtained in (2.1) and (2.6). Form

(3.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then we have the following surprising result.

Theorem 3.1. Consider a strucaured pencil, of the form (1.1) which is regular and of index at most one. Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] be partitioned conformably to the block structure in (3.1), where U, V are obtained in (2.1) and (2.6). Then [Q.sub.11] is invertible and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where the subpencils [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are given in (3.1), (2.4), and (2.7), respectively.

Consequently, the finite eigenvalues of [lambda]N-M are exactly the finite eigenvalues of the structured subpencil [lambda][N.sub.11]-[M.sub.11].

Proof. By definition of Q, one has

(3.2) [lambda]V*NV - V*MV = Q*([lambda][??] - [??])

and

(3.3) [lambda]V*NV - V*MV = Q*([lambda][??] - [??])Q.

Based on (2.2), one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and then (2.6) becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

From this relation we obtain that

(3.4) [[sigma].sub.m] [[??].sup.*.sub.12][Q.sub.11] + [[??].sub.22][Q.sub.21] = 0

and

(3.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Suppose that [Q.sub.11] is singular, then there exists a nonzero vector x such that [Q.sub.11]x = 0. By post-multiplying x to (3.4) and using the fact that [[??].sub.22] is invertible, one has [Q.sub.21]x = 0. Then [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] contradicting the fact that Q is unitary. Therefore, [Q.sub.11] must be invertible. Comparing the (1,1) block on both sides of (3.2) leads to the first equivalence relation

[lambda][N.sub.11] - [M.sub.11] = [Q.sup.*.sub.11] ([lambda] [[??].sub.11] - [[??].sub.11]).

Since [[??].sub.22] is invertible, from the second relation in (3.5), so are [[??].sub.22] and [Q.sub.22]. From (3.4) it follows that

(3.6) [Q.sub.21][Q.sup.-1.sub.11] = -[[sigma].sub.M] [[??].sup.-1.sub.22] [[??].sup.*.sub.12].

Then the matrix [??] defined in (2.3) becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and, defining

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and using the fact that Q is unitary, it follows that

Q = [??][??].

Hence, (3.3) can be written as

(3.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and comparing the (1,1) block on both sides of (3.7) leads to the second equivalence relation

[lambda][N.sub.11] - [M.sub.11] = [Q.sup.*.sub.11]([lambda][[??].sub.11] - S)[Q.sub.11].

One can also obtain another perception of how the subpencils [lambda] [[??].sub.11] - S and [lambda][N.sub.11] - [M.sub.11] are related to the pencils [lambda]V*NV - V*MV and [lambda][??] - [??]. The relation (3.7) implies that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Setting

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we have Q = LR and also

(3.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which has the pencil ([lambda][N.sub.11] - [M.sub.11]) in the (1,1) position. The pencil [lambda]V*NV - V*MV, however, is unitarily congruent or equivalent in the complex transpose case to the original pencil, and the finite eigenvalue part can be simply extracted from the (1,1) block.

It follows from (3.8), that the columns of the matrices

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

form orthonormal bases for the right deflating subspaces corresponding to the finite and infinite eigenvalues of [lambda][??] - [??], respectively. In contrast to this, (2.4) shows that the Schur complement method simply uses the non-orthonormal basis with the columns of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

i.e., the first block column of L for a structure preserving transformation to block diagonalize [lambda][??] - [??].

For the original pencil [lambda]N - M with V = [[V.sub.1], [V.sub.2]] and U = [[U.sub.1], [U.sub.2]], the columns of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

span the right deflating subspaces corresponding to the finite and infinite eigenvalues, respectively.

The minimal angle between the two right deflating subspaces is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [parallel]*[parallel] denotes the Euclidean norm for vectors and the spectral norm for matrices. The minimal angle measures the closeness of the two deflating subspaces.

For Q, there exists a CS decomposition; see [8]. If the size [n.sub.1] of [Q.sub.11] is larger than or equal to the size [n.sub.2] of [Q.sub.22], then this has the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [W.sub.1], [W.sub.2], [Z.sub.1], [Z.sub.2] are unitary matrices, and the diagonal matrices [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] have nonnegative diagonal entries satisfying [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. If [n.sub.1] [less than or equal to] [n.sub.2], then one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Suppose that the singular values of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], then except for possibly some extra singular values at 1, the singular values of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Thus, we have

[[THETA].sub.min] = min [[theta].sub.j].

Introducing

(3.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and using (3.6), one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and hence,

[[THETA].sub.min] = [cot.sup.-1] [rho],

Thus, the smallest singular values of [Q.sub.11] and [Q.sub.22] are given by

(3.10) [[sigma].sub.min]([Q.sub.11]) = [[sigma].sub.min]([Q.sub.22]) = sin [[THETA].sub.min] = 1/[square root of 1 + [[rho].sup.2]],

and [parallel][Q.sub.12][parallel] = [parallel][Q.sub.21][parallel] = [rho]/[square root of 1 + [[rho].sup.2]]. If [rho] is large, which happens commonly when M22 is nearly singular, then the numerical computation of the Schur complement S in (2.5) may be subject to large errors. On the other hand, in this case the block [N.sub.11] = [Q.sup.*.sub.11] [[??].sub.11][Q.sup.11] may be close to a singular matrix as well, which may also lead to big numerical problems, in particular when deciding about the rank of N.

Let us summarize the methods for computing the subpencils [lambda][N.sub.11] - [M.sub.11] defined in (3.1) in the following algorithms. The first algorithm, Algorithm 1, uses the factorizations (2.1) and (2.6) to compute U and V and eventually the subpencil [lambda][N.sub.11] - [M.sub.11]. The second algorithm, Algorithm 2, is just a different version, which computes the factorization (2.2) explicitly, and the factorization (3.1) is computed using Q = U* V instead of V.

ALGORITHM 1. Consider a regular structured pencil [lambda]N - M of index at most one with N* = [[sigma].sub.n] N and [M.sup.*] = [[sigma].sub.m] M.

1. Use a rank revealing QR or singular value decomposition to compute a unitary matrix U such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [N.sub.1] is of full row rank. Partition U = [[U.sub.1] [U.sub.2]] accordingly, and set

[M.sub.2] = [U.sup.*.sub.2]M.

2. Use a rank revealing QR or singular value decomposition to determine a unitary matrix V such that

[M.sub.2]V = [0 [[??].sub.22]], det [[??].sub.22] [not equal to] 0.

3. Compute

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and return [lambda][N.sub.11] - [M.sub.11] as the pencil associated with the finite spectrum of [lambda]N-M.

Algorithm 1 requires two singular value decompositions or rank revealing QR factorizations as well as the associated matrix-matrix multiplications. The most difficult part is the rank decision for the matrix N in step 1, which will affect the number of finite eigenvalues and the size and therefore also the conditioning of [M.sub.2].

An alternative version of Algorithm 1 is given in the following algorithm.

Algorithm 2. Consider a regular structured pencil [lambda]N-M of index at most one with N* = [[sigma].sub.N] N and [M.sup.*] = [[sigma].sub.M] M.

1. Use a Schur-like decomposition of N to compute a unitary matrix U such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and compute

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

2. Use a rank revealing QR or singular value decomposition to compute a unitary matrix [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] partitioned accordingly, such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

3. Compute

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

or, if only [lambda][N.sub.11] - [M.sub.11] is needed, then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Algorithm 2 uses in contrast to Algorithm 1 a structured Schur-like decomposition in the first step and performs two unitary transformations on the pencil. The other costs are comparable to those of Algorithm 1.

We will perform the error and perturbation analysis for these algorithms as well as that of the Schur complement approach in the next section.

4. Error and perturbation analysis. In this section we present a detailed error and perturbation analysis for the different deflation procedures presented in the last sections. Our error analysis for the structure preserving method under unitary transformations will be based on Algorithm 2. The error analysis for Algorithm 1 is essentially the same.

Let u be the unit roundoff, and denote for each of the matrices defined in the previous section the computed counterpart by the corresponding calligraphic letter, e.g., let U, [[??].sub.11] be the computed U and [[??].sub.11] in step 1 of Algorithm 2. Following standard backward error analysis ([8]), there exists a unitary matrix U such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. For the computed M, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Setting X = U*[??] we obtain from the exact factorization (2.2) that

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Assume that [parallel][delta] [[??].sub.11][parallel] < [[sigma].sub.min]([[??].sub.11]). Then [[??].sub.11] + [delta] [[??].sub.11] is nonsingular. One can show that for

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Introduce the unitary matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [(I + GG*).sup.1/2] and [(I + G*G).sup.1/2] are the positive definite square roots of I + GG* and I + G*G, respectively. Then one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Comparing with (4.1), one can show that X = [??]diag([X.sub.1], [X.sub.2]), where [X.sub.1], [X.sub.2] are unitary. Without loss of generality, we set

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and introduce the condition number

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Since [parallel]G[parallel] = 0([k.sub.N]u), and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Clearly, here [k.sub.N] measures the sensitivity of the null space of N. Hence,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and it follows from [??] = X*([??] + U*[delta][M.sub.u]U)X that

(4.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Partition [??] conformable to [??], assume that

[parallel][delta] [[??].sub.22][parallel] < [[sigma].sub.min]([[??].sub.22])

so that [[??].sub.22] is invertible, and let S be the computed Schur complement based on [??]. Then, under the assumption that a backward stable linear system solver is used for computing [[??].sup.-1.sub.22] [[??].sup.*.sub.12] columnwise, one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [rho] is defined in (3.9). Some elementary calculations and (4.2) lead to

S = S + [E.sub.s] + [F.sub.s] =. S + [delta]S,

with

[parallel][F.sub.s][parallel] = O ([(1 + [rho]).sup.2] [k.sub.N] [parallel]m[parallel]u),

where [F.sub.s] results from the errors introduced in [??] by performing the transformation with [??]. Note that due to the structure of the error, we could absorb [E.sub.s] into [F.sub.s]. Altogether, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

It then follows that the columns of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

span the right deflating subspaces corresponding to the finite and infinite eigenvalues, respectively, of the pencil

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [lambda][[??].sub.11] - S, [lambda]0 - [[??].sub.22] are the corresponding computed subpencils associated to these bases.

Let [??] be unitary such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and therefore the right deflating subspace associated with the finite eigenvalues is spanned by the orthonormal columns of the matrix [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Introducing

Y = Q*[??],

it follows from (4.2) that

(4.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Similar to the matrix X, one may express

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where Gy = [([[??].sub.22] + [[epsilon].sub.2]).sup.-1] [[epsilon].sub.1]. One then has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[??].sub.M] = [parallel][[??].sup.-1.sub.22][parallel] [parallel]M[parallel] [greater than or equal to] 1.

Since [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is block upper triangular, and because of the structure of the pencil, the columns of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

are the right deflating subspaces of [lambda](N + [delta][N.sub.u]) - (M + [delta][M.sub.u]) corresponding to the finite and infinite eigenvalues, respectively. Hence, the minimal angle between these two subspaces is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The perturbation in the deflating subspace associated to the finite eigenvalues can then be measured by [[delta].sub.f] := [sin.sup.-1] [parallel][[DELTA].sub.f][parallel] with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the perturbation in the deflating subspace of the eigenvalue infinity can be measured by [[delta].sub.[infinity]]:= [sin.sup.-1] [parallel][[DELTA].sub.[infinity]][parallel] with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Inserting the obtained norm bounds, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It should be noted that in the case that N is already in block diagonal form as [??], which is for example the case when N arises from the middle block of the staircase form (1.2), then U = [??] = I, and the computed subpencil corresponding to the finite eigenvalues is given by

[lambda][[??].sub.11] - S, S = S + [E.sub.s], [parallel][E.sub.s][parallel] = O ([(1 + [rho]).sup.2] [parallel]M[parallel]u) ,

where [E.sub.s] is arising just from the computation of S.

In step 2 of Algorithm 2, for the computed version Q of Q, there exists a unitary matrix [??] such that [parallel]Q - Q[parallel] = O(u) and

(4.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [parallel]E[parallel] = O ([parallel][[??].sub.22][parallel]u) = O ([parallel][M.sub.22][parallel]u).

Introducing

Y = Q*[??],

then similar to (4.3) we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [parallel][epsilon][parallel] = 0([k.sub.N][parallel]M[parallel]u). Then one has

([0 [[??].sub.22]] + EQ + [epsilon]) Y = [0 [[??].sub.22]].

Partitioning EQ + [epsilon] = [[E.sub.11] [E.sub.12]], then [parallel][E.sub.11][parallel] = O (([parallel][M.sub.22][parallel] + [k.sub.N] [parallel]M[parallel]) u). Similar to y, we may express

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [G.sub.Y] = [([[??].sub.22] + [E.sub.12]).sup.-1] [E.sub.11]. We obtain the estimates

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Using these estimates, the computed [N.sub.11] and [M.sub.11] can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [parallel][delta][N.sub.11][parallel] = O ([parallel][Q.sub.11][[parallel].sup.2] [parallel]N[parallel]u), and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Using the relation between [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Since

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

using (3.6), one obtains

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where, by using (3.5),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Using [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], it follows that

[[DELTA].sub.M,2] = [Q.sup.*.sub.21][E.sub.11] + o(u),

and hence,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Here we have used the fact that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In the situation that only the finite eigenvalues are considered, one can reduce the effort to compute [lambda][N.sub.11] - [M.sub.11] only. So if we set

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

then the columns of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] form an orthonormal basis for the right deflating subspace of [lambda][??] - [[??].sub.1] corresponding to the finite eigenvalues, where [[??].sub.1] is a perturbed [??] of order O(M11u) because of the inexact QR factorization (4.4). The resulting subpencil is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and the pencil [lambda][N.sub.11] - [M.sub.11] is obtained by adding the extra error pencil [lambda][[DELTA].sub.N] - [[DELTA].sub.M] introduced by the numerical computation of [N.sub.11] and [M.sub.11]. Similarly, the columns of 0 span the right deflating subspace of [lambda][??] - [[??].sub.1] corresponding to the eigenvalue infinity, and the columns of the matrices

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

span the corresponding deflating subspaces of a pencil slightly perturbed from [lambda]N - M.

Setting

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we have [[member of].sub.[infinity]] = [[delta].sub.[infinity]] = 0([k.sub.N]u), and similarly,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the minimal angle between the two perturbed subspaces is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

If [lambda]N - M is already in the form [lambda][??] - [??], then U = [??] = I, and in this case,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and it follows that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We may express [N.sub.11] and [M.sub.11] as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

By using (3.10), we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Recalling from (3.10) that [parallel][Q.sup.-1.sub.11][parallel] = [square root of 1 + [[rho].sup.2]], and comparing the errors in the subpencils [lambda][[??].sub.11] - S and [lambda][N.sub.11] - [M.sub.11], it follows that [parallel][[??].sub.M][parallel] and [parallel][delta]S[parallel] have the same order, while [parallel][[??].sub.N][parallel] can be larger than [parallel][delta][[??].sub.11][parallel] by a factor [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The latter will be of equal order only when [parallel] [[??].sub.N] [parallel] can be larger than [parallel] [delta][[??].sub.11][parallel] by a factor [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The latter will be of equal order only when [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is not too large.

On the other hand, if we transform [lambda][[??].sub.11] - S to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first quantity is of the same order as [parallel][[DELTA].sub.N][parallel], while the second one can be larger than [parallel][[DELTA].sub.M][parallel] by a factor [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is not too large. So in terms of numerical stability, the error analysis does not show which method has an advantage over the other. This is an unusual circumstance in numerical analysis, where usually the methods based on unitary transformations in the worst case situation have smaller errors than the ones based on non-unitary transformations. We will demonstrate this effect in the numerical examples in Section 7.

5. An alternative point of view. Another way to understand the new unitary structure preserving method of deflating the infinite eigenvalue part is to consider the extension trick introduced in [9]. Partition the matrix U in (2.1) as U = [[U.sub.1] [U.sub.2]] with [U.sub.1] of the same size as [N.sup.T.sub.1]. For any eigenvector x of [lambda]N - M corresponding to a finite eigenvalue [lambda], it then follows from U*([lambda]N - M)x = 0 that [U.sup.*.sub.2]Mx = 0. So the original eigenvalue problem is equivalent to the non-square, non-symmetric eigenvalue problem

(5.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and we are only looking for eigenvectors in the nullspace of [U.sup.*.sub.2]M. Let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

be partitioned according to the block structure in (3.1). Then the eigenvalue problem (5.1) is equivalent to the extended eigenvalue problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Clearly then it follows that [z.sub.2] = 0 and

([lambda][N.sub.11] - [M.sub.11])[z.sub.1] = 0.

Hence, an eigenvector of [lambda]N - M corresponding to a finite eigenvalue has the form [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where [z.sub.1] is an eigenvector of the subpencil [lambda][N.sub.11] - [M.sub.11]. More generally, the columns of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] form an orthonormal basis for the deflating subspace of [lambda]N - M corresponding to the finite eigenvalues.

6. Palindromic and anti-palindromic pencils. Structured pencils closely related to those discussed before are the palindromic/anti-palindromic pencils of the form

[lambda][A.SUP.*] - [sigma]A,

where A is an arbitrary real or complex square matrix and [sigma] = [+ or -]1; see [10, 16] for a detailed discussion of the relationship. Assume further that the eigenvalue [sigma] (if it occurs) is semisimple, which corresponds to the property of being regular and of index at most one for the structured pencils discussed before. Applying a Cayley transformation (see [10, 11]), the pencil [lambda]N - M with

N := [A.sup.*] - [sigma]A = -[sigma][N.sup.*], M := [A.sup.*] + [sigma]A = [sigma][M.sup.*]

is a structured pencil of the form (1.1) with [sigma]N = -[sigma] and [sigma]M = [sigma], and the semi-simple eigenvalue [sigma] becomes the eigenvalue to and the pencil has index at most one. To this pencil we can apply the discussed transformations V, U as defined in (2.1) and (2.6). Then

(6.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Partitioning

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

conformably, one has

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and hence

[A.sup.*.sub.2] = [sigma] [B.sub.2], [M.sub.2] = 2[sigma] [B.sub.2].

Furthermore, applying V to [M.sub.2], we obtain [M.sub.2]V = [0 [[??].sub.22]], and hence [B.sub.2]V = [0 [[??].sub.22]] with [[??].sub.22] = [sigma][[??].sub.22]/2.

Setting

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and thus

[lambda][N.sub.11] - [M.sub.11] = [lambda]([A.sup.*.sub.11] - [sigma][A.sub.11]) - ([A.sup.*.sub.11] + [sigma][A.sub.11])

is the subpencil of [lambda]N - M with the finite eigenvalues only. By taking the inverse Cayley transformation, this is equivalent to the property that the palindromic subpencil

[lambda][A.sup.*.sub.11] - [sigma][A.sub.11]

contains all the eigenvalues of [lambda][A.sup.*] - [sigma]A except the eigenvalue [sigma], which has been deflated.

It should be noted though that the Cayley transformation or its inverse may lead to cancellation errors if there are eigenvalues close to [sigma] but not equal to [sigma]. Note when some eigenvalues are close [sigma] but not equal to a, since [A.sup.*] - [sigma]A has to be computed explicitly, the rank decision in computing the first factorization (6.1) is more difficult to make.

7. Numerical examples. In this section we present several numerical tests to compare the computed finite eigenvalues of the subpencils generated by the three methods: structured unitary equivalence transformation (Algorithms 1 and 2), Schur complement transformation (by using (2.4)), and the non-structured equivalence transformation (with (2.7)). All the tests were performed in MATLAB (Version 8.2.0) with double precision. All the tested pencils were real, and we use relative errors to measure the accuracy of the computed finite eigenvalues, where the "exact" eigenvalues are computed with the Matlab code eig in extended precision vpa. The tested pencils are all skew-symmetric/symmetric. The finite eigenvalues of the subpencil extracted with the two structure preserving methods are computed with the MEX code skewHamileig ([2, 3]). If the pencil associated with the finite eigenvalues is extracted via non-structured equivalence transformations, we use the MATLAB code qz. We also display the quantities [[THETA].sub.min], [[delta].sub.f], [[member of].sub.f], [[member of].sub.[infinity]], as well as p defined in (3.9). Since [[??].sub.min] and [[??].sub.min] are always very close to 0min, we do not display them. We do not show either since it is the same as [[member of].sub.[infinity]]. The following quantities were obtained from the numerical tests.

* [Eu.sub.max], [Es.sub.max], [Eh.sub.max]: the maximum relative error of the finite eigenvalues for a given pencil with the structured unitary equivalence transformation method, the Schur complement method, and the non-structured equivalence transformation method, respectively.

* [Eu.sub.min], [Es.sub.min], [Eh.sub.min]: the minimum relative error of the finite eigenvalues for a given pencil with the structured unitary equivalence transformation method, the Schur complement method, and the non-structured equivalence transformation method, respectively.

* [Eu.sub.GM], [Es.sub.GM], [Eh.sub.GM]: the geometric mean of the relative errors of the finite eigenvalues for a given pencil with the structured unitary equivalence transformation method, the Schur complement method, and the non-structured equivalence transformation method, respectively.

* [Re.sub.hat]: the maximum real part of the eigenvalues computed by the non-structured equivalence transformation method in the case when all the exact finite eigenvalues are purely imaginary.

Example 1. We tested a set of skew-symmetric/symmetric pencils [X.sup.T]([lambda]N - M)X, where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and X is a randomly generated nonsingular matrix. The finite eigenvalues of such a pencil are always [+ or -]i[square root of 6] and [+ or -]i[square root of 6]/[beta] independent of X and [alpha]. If [beta] is small in modulus, then N is close to a singular matrix, and when [alpha] is small in modulus, then it can be expected that the eigenvalue infinity part will affect the finite eigenvalues. Note that when [beta] = 1, then the pencil has two pairs of double eigenvalues.

We have performed the computations 10 times for each pair of ([alpha], [beta]) where each time a different random matrix X is generated. Tables 7.1-7.4 display the results for each pair of ([alpha], [beta]). The three methods determine the finite eigenvalues with about the same order of relative errors. The non-structured method produces slightly less accurate eigenvalues. Because it forms non-structured subpencils, the computed eigenvalues are no longer purely imaginary. The real parts of the computed finite eigenvalues grow as [absolute value of [beta]] is getting smaller. For all three methods, in this example, the errors seem insensitive to the value of [alpha], while they increase when [absolute value of [beta]] decreases. However, a does affect the deflating subspace associated with the finite eigenvalues.

Example 2. We tested a set of skew-symmetric/symmetric pencils [X.sup.T]([lambda]N - M)X, where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and X is a randomly generated matrix. The finite eigenvalues of such a pencil are always [+ or -][beta] with both algebraic and geometric multiplicities of 2.

When [absolute value of [beta]] is tiny, then all the finite eigenvalues are close to zero and the eigenvalue condition number is expected to increase. Tables 7.5-7.8 display the results for 4 different sets of ([alpha], [beta]) (with 10 pencils for each set as in Example 1), illustrating that the three methods have essentially the same accuracy. The results also show that the relative errors do increase when [absolute value of [beta]] is getting tiny, and that again, the choice of a does not affect the accuracy very much. For this example, no method yields finite eigenvalues that are exactly real.

Example 3. For skew-symmetric/symmetric pencils of the form [X.sup.T] ([lambda]N - M)X, where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and X being a randomly generated nonsingular matrix, the finite eigenvalues of such a pencil are always [+ or -][beta] with both algebraic and geometric multiplicities of 2. The numerical results are similar to those of Example 2 and are not presented here.

Example 4. For skew-symmetric/symmetric pencils of the form [X.sup.T] ([lambda]N - M)X, where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and X a randomly generated nonsingular matrix, there exist two 2 x 2 Jordan blocks, one for the eigenvalue [beta] and another for -[beta]. In this case it is expected that the computed finite eigenvalues have big errors due to the Jordan blocks, and when [absolute value of [beta]] becomes smaller, then the errors will get larger. Our numerical results presented in Tables 7.9 and 7.10 confirm these expectations for two sets of ([alpha], [beta]).

All the numerical tests confirm our expectations obtained from the error analysis and illustrate that, surprisingly, the Schur complement approach and the approach via orthogonal structure preserving transformations deliver the same accuracy despite possible ill-conditioning. The unstructured orthogonal transformation, however, shows the expected problems with purely imaginary eigenvalues.

8. Conclusions. We have presented and analyzed several methods to deflate the infinite eigenvalue part in a structured regular pencil of index at most one. We have shown via a careful error analysis that it is possible to do this deflation via a structure preserving real orthogonal or unitary transformation as well as with a Schur complement approach. Both methods yield similar results in the perturbation and error analysis and this is confirmed in the numerical tests. This is surprising and counterintuitive to the general wisdom that unitary transformations typically perform better than transformations with potentially unbounded transformations. Nonetheless, we suggest to use the unitary structure preserving transformations since they are very much in line with the staircase algorithm. Both methods are definitely preferable to the non-structured unitary transformation.

Acknowledgments. This work was started while the first author visited the second and completed while the second author was visiting the first. The authors thank both host institutions for their hospitality and generosity. The authors also thank the referees for their valuable comments and suggestions.

REFERENCES

[1] P. BENNER, R. Byers, V. Mehrmann, AND H. Xu, A robust numerical method, for the j-iteration in HTC control, Linear Algebra Appl., 425 (2007), pp. 548-570.

[2] P. BENNER, V. Sima, AND M. Voigt, FORTRAN 77 subroutines for the solutions of skew Hamiltonian/Hamiltonian eigenproblems--part I: algorithms and applications, Preprint MPIMD/13-11, Max Planck Institute, Magdeburg, 2013.

[3]--, FORTRAN 77 subroutines for the solutions of skew-Hamiltonian/Hamiltonian eigenproblems--part II: implementation and numerical results, Preprint MPIMD/13-12, Max Planck Institute, Magdeburg, 2013.

[4] T. BRULL AND V. Mehrmann, STCSSP: A FORTRAN 77 routine to compute a structured staircase form for a (skew-)symmetric/(skew-)symmetric matrix pencil, Preprint 31-2007, Institut fur Mathematik, Technical University Berlin, Berlin, 2007.

[5] R. BYERS, V. Mehrmann, AND H. Xu, A structured staircase algorithm for skew-symmetric/symmetric pencils, Electron. Trans. Numer. Anal., 26 (2007), pp. 1-33. http://etna.math.kent.edu/vol.26.2007/pp1-33.dir

[6] J. DEMMEL AND B. KAGSTROM, Stably computing the Kronecker structure and reducing subspaces of singular pencils A-[lambda]B for uncertain data, in Large Scale Eigenvalue Problems (Oberlech, 1985), J. Cullum and R. A. Willoughby, eds., vol. 127 of North-Holland Math. Stud., North-Holland, Amsterdam, 1986, pp. 283-323.

[7] P. GAHINET AND A. J. LAUB, Numerically reliable computation of optimal performance in singular HTC control, SIAM J. Control Optim., 35 (1997), pp. 1690-1710.

[8] G. H. GOLUB AND C. F. VAN LOAN, Matrix Computations, 3rd ed., Johns Hopkins University Press, Baltimore, 1996.

[9] P. KUNKEL, V. MEHRMANN, AND L. SCHOLZ, Self-adjoint differential-algebraic equations, Math. Control Signals Systems, 26 (2014), pp. 47-76.

[10] D. S. MACKEY, N. MACKEY, C. MEHL, AND V. MEHRMANN, Structured polynomial eigenvalue problems: good vibrations from good linearizations, SIAM J. Matrix Anal. Appl., 28 (2006), pp. 1029-1051.

[11]--, Mobius transformations of matrix polynomials, Linear Algebra Appl., in press, doi: 10.1016/j.laa.2014.05.013

[12] V. MEHRMANN, The Autonomous Linear Quadratic Control Problem. Theory and Numerical Solution, Springer, Heidelberg, 1991.

[13] V. MEHRMANN AND D. WATKINS, Structure-preserving methods for computing eigenpairs of large sparse skew-Hamiltonian/Hamiltonian pencils, SIAM J. Sci. Comput., 22 (2000), pp. 1905-1925.

[14] V. MEHRMANN AND H. XU, Numerical methods in control, J. Comput. Appl. Math., 123 (2000), pp. 371-394.

[15] P. H. PETKOV, N. D. CHRISTOV, AND M. M. Konstantinov, Computational Methods for Linear Control Systems, Prentice-Hall Int., Hertfordshire, 1991.

[16] C. SCHRODER, Palindromic and Even Eigenvalue Problems--Analysis and Numerical Methods, Ph.D. Thesis, Department of Mathematics, Technical University Berlin, Berlin, 2008.

[17] V. SIMA, Algorithms for Linear-Quadratic Optimization, Marcel Dekker, New York, 1996.

[18] R. C. THOMPSON, Pencils of complex and real symmetric and skew matrices, Linear Algebra Appl., 147 (1991), pp. 323-371.

[19] P. VAN DOOREN, A generalized eigenvalue approach for solving Riccati equations, SIAM J. Sci. Statist. Comput., 2 (1981), pp. 121-135.

[20] K. ZHOU, J. C. DOYLE, AND K. GLOVER, Robust and Optimal Control, Prentice-Hall, Englewood Cliffs, 1996.

VOLKER MEHRMANN ([dagger]) AND HONGGUO XU ([double dagger])

* Received May 28, 2014. Accepted October 13, 2014. Published online on January, 23, 2015. Recommended by D. Szyld. Supported by Deutsche Forschungsgemeinschaft through the DFG Research Center Matheon Mathematics for Key Technologies in Berlin and by the University of Kansas, Dept. of Mathematics, during a research visit in 2013. The second author is partially supported by the Alexander von Humboldt Foundation and the Deutsche Forschungsgemeinschaft, through the DFG Research Center Matheon Mathematics for Key Technologies in Berlin.

([dagger]) Institut fur Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, Germany (mehrmann@math.tu-berlin.de).

([double dagger]) Department of Mathematics, University of Kansas, Lawrence, KS 66045, USA (xu@math.ku.edu).
```Table 7.1

([alpha], [beta]) = ([10.sup.-3], 1): eigenvalue errors,
errors in deflating subspaces, [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]          4e-15    9e-15    2e-14    4e-14    6e-14
[Es.sub.max]          7e-15    4e-14    1e-14    6e-14    6e-14
[Eh.sub.max]          1e-14    8e-14    2e-14    1e-13    6e-14
[Eu.sub.min]          2e-15    9e-15    3e-15    1e-14    5e-15
[Es.sub.min]            0      1e-14    5e-15    9e-15    4e-15
[Eh.sub.min]          1e-15    4e-14    1e-14    1e-14    6e-15
[Eu.sub.GM]           2e-15    9e-15    8e-15    2e-14    2e-14
[Es.sub.GM]             0      2e-14    8e-15    2e-14    2e-14
[Eh.sub.GM]           4e-15    6e-14    2e-14    4e-14    2e-14

[Re.sub.hat]          2e-14    2e-13    2e-14    2e-13    1e-14

[[member of].sub.f]   3e-11    3e-11    2e-12    2e-11    4e-11
[[delta].sub.f]       3e-11    2e-11    2e-12    7e-12    2e-11
[[member of].sub.     6e-16    6e-16    5e-16    7e-16    5e-16
[infinity]]

[THETA]               6e-01    2e-01    2e-01    2e-01    8e-02
[rho]                  1.4      4.7      4.0      6.4       13

6        7        8        9        10

[Eu.sub.max]          6e-14    7e-14    8e-14    2e-13    4e-13
[Es.sub.max]          6e-14    7e-14    1e-13    6e-14    5e-13
[Eh.sub.max]          8e-14    6e-14    2e-13    1e-12    5e-13
[Eu.sub.min]          6e-15    2e-14    1e-14    9e-16    3e-15
[Es.sub.min]          4e-15    2e-14    6e-15    4e-15    4e-15
[Eh.sub.min]          3e-15    5e-15    9e-14    6e-15    2e-15
[Eu.sub.GM]           2e-14    3e-14    3e-14    1e-14    3e-14
[Es.sub.GM]           2e-14    4e-14    3e-14    2e-14    4e-14
[Eh.sub.GM]           2e-14    2e-14    1e-13    8e-14    4e-14

[Re.sub.hat]          4e-14    9e-15    2e-13    4e-13    1e-13

[[member of].sub.f]   4e-11    1e-11    2e-11    1e-10    7e-11
[[delta].sub.f]       2e-11    1e-11    2e-11    1e-10    6e-11
[[member of].sub.     5e-16    3e-16    3e-16    2e-16    5e-16
[infinity]]

[THETA]               3e-01    2e-01    8e-02    8e-02    1e-01
[rho]                  3.5      6.5       13       13      7.6

Table 7.2

([alpha], [beta]) = ([10.sup.-3], [10.sup.-5]): eigenvalue errors,
errors of deflating subspaces, and [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]           8e-13    3e-12    6e-12    7e-12    9e-12
[Es.sub.max]           3e-12    3e-12    7e-12    9e-12    4e-11
[Eh.sub.max]           8e-12    7e-11    3e-10    8e-12    1e-10
[Eu.sub.min]           6e-14    4e-16    8e-14    4e-13    2e-14
[Es.sub.min]           6e-14    7e-16    3e-14    2e-13    2e-15
[Eh.sub.min]           5e-14    4e-14    2e-12    4e-13    4e-14
[Eu.sub.GM]            2e-13    3e-14    7e-13    2e-12    4e-13
[Es.sub.GM]            4e-13    5e-14    4e-13    1e-12    3e-13
[Eh.sub.GM]            6e-13    2e-12    2e-11    2e-12    2e-12

[Re.sub.hat]           2e-06    3e-06    7e-05    2e-06    2e-05

[[member of].sub.f]    1e-08    8e-08    5e-08    8e-10    2e-07
[[delta].sub.f]        1e-08    8e-08    5e-08    8e-10    2e-07
[[member of].sub.      2e-12    2e-11    6e-12    6e-12    4e-11
[infinity]]

[THETA]                2e-01    3e-01    2e-01    1e-01    1e-01
[rho]                   4.2      3.1      6.1      9.8      8.6

6        7        8        9        10

[Eu.sub.max]           2e-11    8e-11    8e-11    3e-10    2e-09
[Es.sub.max]           3e-11    4e-11    4e-11    6e-10    2e-09
[Eh.sub.max]           3e-11    1e-10    4e-11    1e-09    3e-09
[Eu.sub.min]           9e-15    9e-13    3e-15    5e-16    1e-11
[Es.sub.min]           9e-15    1e-12    6e-15    7e-16    1e-11
[Eh.sub.min]           8e-15    9e-13    9e-15    2e-15    1e-11
[Eu.sub.GM]            4e-13    9e-12    5e-13    4e-13    2e-10
[Es.sub.GM]            5e-13    8e-12    5e-13    7e-13    1e-10
[Eh.sub.GM]            5e-13    1e-11    6e-13    2e-12    2e-10

[Re.sub.hat]           7e-06    1e-05    1e-05    3e-04    6e-04

[[member of].sub.f]    1e-07    1e-07    1e-07    4e-07    6e-06
[[delta].sub.f]        1e-07    1e-07    1e-07    4e-07    6e-06
[[member of].sub.      1e-11    7e-12    6e-12    1e-10    3e-11
[infinity]]

[THETA]                1e-01    5e-02    1e-01    8e-02    6e-03
[rho]                   7.2       19      8.3       12      165

Table 7.3

([alpha], [beta]) = ([10.sup.-7], 1): eigenvalue errors,
errors of deflating subspaces, and [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]          3e-15    4e-15    4e-15    6e-15    7e-15
[Es.sub.max]          1e-15    4e-15    8e-15    5e-15    4e-15
[Eh.sub.max]          3e-14    8e-15    3e-14    4e-15    5e-15
[Eu.sub.min]          7e-16    9e-16    2e-16    7e-16      0
[Es.sub.min]          2e-16    4e-16    1e-15    9e-16    1e-15
[Eh.sub.min]          4e-16    4e-15    3e-15    2e-15    4e-15
[Eu.sub.GM]           2e-15    2e-15    9e-16    2e-15      0
[Es.sub.GM]           5e-16    1e-15    3e-15    2e-15    2e-15
[Eh.sub.GM]           3e-15    5e-15    1e-14    3e-15    4e-15

[Re.sub.hat]          3e-14    2e-14    3e-14    3e-15    1e-14

[[member of].sub.f]   5e-08    4e-08    1e-07    2e-07    2e-08
[[delta].sub.f]       4e-08    3e-08    1e-07    9e-08    2e-08
[[member of].sub.     2e-16    5e-16    3e-16    3e-16    3e-16
[infinity]]

[THETA]               2e-01    7e-01    3e-01    2e-01    4e-01
[rho]                  4.6      1.2      3.7      5.8      2.4

6        7        8        9        10

[Eu.sub.max]          9e-15    9e-15    9e-15    4e-14    6e-14
[Es.sub.max]          8e-15    4e-15    2e-14    2e-14    4e-14
[Eh.sub.max]          2e-14    3e-14    3e-14    4e-14    3e-14
[Eu.sub.min]          2e-15    5e-15    4e-15    5e-15    2e-15
[Es.sub.min]          1e-15    2e-15    3e-15    2e-14    3e-15
[Eh.sub.min]          4e-15    4e-15    1e-15    8e-15    1e-14
[Eu.sub.GM]           5e-15    6e-15    6e-15    1e-14    1e-14
[Es.sub.GM]           3e-15    3e-15    8e-15    2e-14    1e-14
[Eh.sub.GM]           9e-15    1e-14    7e-15    2e-14    2e-14

[Re.sub.hat]          1e-14    3e-14    7e-15    3e-14    3e-14

[[member of].sub.f]   2e-07    4e-07    1e-07    1e-06    4e-07
[[delta].sub.f]       1e-07    3e-07    1e-07    6e-07    3e-07
[[member of].sub.     7e-16    3e-16    2e-16    2e-15    5e-16
[infinity]]

[THETA]               3e-01    2e-01    7e-02    2e-01    2e-01
[rho]                  3.8      4.9       14      4.6      5.0

Table 7.4

([alpha], [beta]) = ([10.sup.-7], [10.sup.-5]): eigenvalue
errors, errors of deflating subspaces, and [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]           1e-12    2e-12    9e-12    3e-11    3e-11
[Es.sub.max]           2e-11    2e-12    7e-12    5e-12    2e-11
[Eh.sub.max]           2e-11    3e-11    3e-11    7e-11    3e-10
[Eu.sub.min]           5e-16    1e-15    2e-15    1e-13    4e-15
[Es.sub.min]           2e-15    1e-14    9e-16    2e-13    5e-15
[Eh.sub.min]           6e-15    9e-14    7e-15    3e-13    3e-14
[Eu.sub.GM]            3e-14    4e-14    1e-13    2e-12    4e-13
[Es.sub.GM]            2e-13    1e-13    8e-14    9e-13    3e-13
[Eh.sub.GM]            3e-13    2e-12    8e-13    5e-12    2e-12

[Re.sub.hat]           4e-06    6e-06    2e-05    1e-05    6e-05

[[member of].sub.f]    2e-04    4e-04    1e-03    4e-04    2e-03
[[delta].sub.f]        2e-04    4e-04    1e-03    4e-04    2e-03
[[member of].sub.      6e-12    2e-12    2e-11    9e-12    3e-11
[infinity]]

[THETA]                3e-01    2e-01    5e-01    5e-02    1e-01
[rho]                   3.8      5.6      1.7       21       10

6        7        8        9        10

[Eu.sub.max]           3e-11    52e-11   9e-11    1e-10    2e-10
[Es.sub.max]           3e-11    7e-11    1e-10    8e-11    1e-10
[Eh.sub.max]           1e-10    6e-11    1e-10    9e-11    1e-10
[Eu.sub.min]           2e-15    9e-14    1e-14    3e-13    1e-15
[Es.sub.min]           3e-15    4e-14    4e-15    2e-13    6e-15
[Eh.sub.min]           1e-14    7e-14    5e-15    7e-13    1e-14
[Eu.sub.GM]            3e-13    2e-12    1e-12    5e-12    4e-13
[Es.sub.GM]            3e-13    2e-12    7e-13    4e-12    8e-13
[Eh.sub.GM]            1e-12    2e-12    8e-13    8e-12    1e-12

[Re.sub.hat]           2e-05    4e-06    2e-05    2e-05    2e-05

[[member of].sub.f]    5e-04    2e-03    2e-04    1e-04    2e-04
[[delta].sub.f]        5e-04    2e-03    2e-04    1e-04    2e-04
[[member of].sub.      2e-11    1e-11    3e-11    4e-11    3e-11
[infinity]]

[THETA]                4e-01    7e-02    5e-02    2e-01    7e-02
[rho]                   2.5       14       20      4.0       14

Table 7.5

([alpha], [beta]) = ([10.sup.-3], 1): eigenvalue errors,
errors of deflating subspaces, and [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]          2e-15    2e-14    2e-14    2e-14    2e-14
[Es.sub.max]          2e-15    3e-14    3e-14    9e-16    2e-14
[Eh.sub.max]          2e-14    6e-14    6e-14    5e-14    2e-13
[Eu.sub.min]          2e-15    2e-14    2e-14    2e-14    2e-14
[Es.sub.min]          9e-16    3e-14    6e-15    9e-16    2e-14
[Eh.sub.min]          2e-15    8e-15    2e-14    5e-16    1e-14
[Eu.sub.GM]           2e-15    2e-14    2e-14    2e-14    2e-14
[Es.sub.GM]           2e-15    3e-14    2e-14    9e-16    2e-14
[Eh.sub.GM]           5e-15    3e-14    3e-14    8e-15    4e-14

[[member of].sub.f]   8e-12    9e-11    9e-12    2e-11    4e-11
[[delta].sub.f]       4e-12    4e-11    9e-12    1e-11    2e-11
[[member of].sub.     4e-16    5e-16    6e-16    3e-16    4e-16
[infinity]]

[THETA]               3e-01    3e-01    5e-01    2e-01    3e-01
[rho]                   4        4        3        8        5

6        7        8        9        10

[Eu.sub.max]          2e-14    3e-14    5e-14    6e-14    7e-12
[Es.sub.max]          2e-14    4e-14    9e-14    7e-14    3e-12
[Eh.sub.max]          2e-14    9e-14    6e-14    2e-13    6e-12
[Eu.sub.min]          2e-15    9e-15    2e-14    2e-14    2e-14
[Es.sub.min]          7e-15    2e-14    2e-14    2e-15    5e-14
[Eh.sub.min]          5e-16    7e-15    6e-14    1e-14      0
[Eu.sub.GM]           6e-15    2e-14    3e-14    3e-14    4e-13
[Es.sub.GM]           9e-15    3e-14    4e-14    41e-14   4e-13
[Eh.sub.GM]           5e-15    3e-14    6e-14    4e-14      0

[[member of].sub.f]   4e-12    5e-12    3e-11    9e-12    2e-10
[[delta].sub.f]       4e-12    4e-12    2e-11    6e-12    2e-10
[[member of].sub.     4e-16    2e-16    5e-16    3e-16    5e-16
[infinity]]

[THETA]               3e-01    3e-01    5e-01    2e-01    3e-02
[rho]                   5        4        3        6        40

Table 7.6

([alpha], [beta]) = ([10.sup.-3], [10.sup.-5]): eigenvalue
errors, errors of deflating subspaces, and [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]          9e-05    2e-04    2e-04    3e-04    4e-04
[Es.sub.max]          5e-05    5e-04    5e-05    2e-04    3e-04
[Eh.sub.max]          2e-04    3e-04    2e-04    6e-03    4e-04
[Eu.sub.min]          9e-05    2e-04    5e-06    7e-05    4e-04
[Es.sub.min]          5e-05    4e-05    5e-05    2e-04    3e-04
[Eh.sub.min]          2e-04    3e-04    9e-06    2e-04    4e-04
[Eu.sub.GM]           9e-05    2e-04    3e-05    2e-04    4e-04
[Es.sub.GM]           5e-05    2e-04    5e-05    2e-04    3e-04
[Eh.sub.GM]           2e-04    3e-04    5e-05    9e-04    4e-04

[[member of].sub.f]   2e-12    3e-11    5e-12    4e-11    2e-11
[[delta].sub.f]       2e-13    2e-11    4e-12    2e-11    2e-11
[[member of].sub.     3e-16    4e-16    3e-16    4e-16    4e-16
[infinity]]

[THETA]               5e-01    5e-01    3e-01    4e-01    2e-01
[rho]                   3        3        5        4        9

6        7        8        9        10

[Eu.sub.max]          4e-04    5e-04    7e-04    3e-03    2e-02
[Es.sub.max]          3e-04    6e-04    2e-04    2e-03    6e-03
[Eh.sub.max]          2e-03    2e-03    2e-03    4e-03    3e-02
[Eu.sub.min]          5e-05    5e-04    5e-05    2e-05    2e-04
[Es.sub.min]          5e-05    6e-04    2e-04    5e-05    2e-04
[Eh.sub.min]          2e-04    3e-05    4e-05    7e-04    9e-04
[Eu.sub.GM]           2e-04    5e-04    2e-04    2e-04    2e-03
[Es.sub.GM]           2e-04    6e-04    2e-04    3e-04    8e-04
[Eh.sub.GM]           4e-04    3e-04    2e-04    2e-03    5e-03

[[member of].sub.f]   2e-11    5e-11    6e-12    2e-11    2e-10
[[delta].sub.f]       2e-11    4e-11    5e-12    1e-11    7e-11
[[member of].sub.     3e-16    3e-16    5e-16    4e-16    7e-16
[infinity]]

[THETA]               3e-01    7e-02    3e-01    2e-01    2e-02
[rho]                   5        20       4        6        80

Table 7.7

([alpha], [beta]) = ([10.sup.-7], 1): eigenvalue errors,
errors of deflating subspaces, and [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]           9e-15    1e-14    1e-14    2e-14    3e-14
[Es.sub.max]           2e-14    9e-15    3e-14    2e-14    3e-14
[Eh.sub.max]           2e-14    2e-14    3e-13    3e-13    2e-13
[Eu.sub.min]           6e-15    1e-14    7e-15    2e-14    3e-14
[Es.sub.min]           2e-15    3e-15    6e-15    2e-14    3e-14
[Eh.sub.min]           9e-15    7e-15    3e-14    6e-15    3e-14
[Eu.sub.GM]            7e-15    1e-14    9e-15    2e-14    3e-14
[Es.sub.GM]            5e-15    5e-15    2e-14    2e-14    3e-14
[Eh.sub.GM]            2e-14    2e-14    6e-14    5e-14    7e-14

[[member of].sub.f]    9e-08    5e-07    9e-08    2e-07    9e-07
[[delta].sub.f]        7e-08    2e-07    9e-08    9e-08    3e-07
[[member of].sub.      5e-16    3e-16    4e-16    7e-16    4e-16
[infinity]]

[THETA]                3e-01    5e-01    2e-01    7e-01    3e-01
[rho]                    5e       3        9        2        4

6        7        8        9        10

[Eu.sub.max]           3e-14    5e-14    5e-14    2e-13    4e-13
[Es.sub.max]           5e-14    3e-14    2e-13    2e-14    4e-13
[Eh.sub.max]           1e-14    6e-14    2e-13    3e-13    4e-13
[Eu.sub.min]           3e-14    4e-15    5e-14    4e-15    8e-15
[Es.sub.min]           3e-15    2e-15    4e-15    2e-14    7e-15
[Eh.sub.min]           7e-15    3e-15    2e-13    2e-15    7e-15
[Eu.sub.GM]            3e-14    2e-14    5e-14    3e-14    6e-14
[Es.sub.GM]            2e-14    8e-15    3e-14    2e-14    5e-14
[Eh.sub.GM]            8e-15    2e-14    2e-13    3e-14    5e-14

[[member of].sub.f]    2e-07    5e-07    2e-08    2e-07    7e-07
[[delta].sub.f]        1e-07    5e-07    8e-09    5e-08    7e-07
[[member of].sub.      3e-16    3e-16    4e-16    3e-16    6e-16
[infinity]]

[THETA]                7e-02    3e-01    7e-02    5e-01    2e-01
[rho]                    20       5        20       2        8

Table 7.8

([alpha], [beta]) = ([10.sup.-7], 10-5): eigenvalue errors,
errors of deflating subspaces, and [THETA], [rho].

1        2        3        4        5

[Eu.sub.max]           3e-05    3e-05    5e-05    2e-04    3e-04
[Es.sub.max]           8e-06    3e-05    6e-05    2e-04    6e-04
[Eh.sub.max]           5e-04    4e-05    4e-03    8e-04    2e-03
[Eu.sub.min]           2e-06    3e-05    5e-05    2e-04    1e-04
[Es.sub.min]           8e-06    3e-05    6e-05    2e-04    7e-05
[Eh.sub.min]           2e-05    4e-05    2e-04    2e-04    3e-04
[Eu.sub.GM]            6e-06    3e-05    5e-05    2e-04    2e-04
[Es.sub.GM]            8e-06    3e-05    6e-05    2e-04    2e-04
[Eh.sub.GM]            9e-05    4e-05    8e-04    4e-04    7e-04

[[member of].sub.f]    4e-08    2e-08    5e-07    5e-07    2e-07
[[delta].sub.f]        3e-08    3e-08    4e-07    5e-07    3e-07
[[member of].sub.      2e-16    1e-15    3e-16    7e-16    4e-16
[infinity]]

[THETA]                2e-01    5e-01    6e-01    2e-01    1e-01
[rho]                    6        3        2        5        20

6        7        8        9        10

[Eu.sub.max]           4e-04    7e-04    3e-03    1e-01    2e-01
[Es.sub.max]           4e-04    8e-04    2e-03    3e-02    3e-01
[Eh.sub.max]           4e-04    3e-03    2e-03    8e-02    2e+00
[Eu.sub.min]           4e-04    7e-04    3e-05    5e-05    9e-05
[Es.sub.min]           4e-04    8e-04    6e-05    4e-05    2e-03
[Eh.sub.min]           4e-04    7e-04    4e-04    9e-05    2e-04
[Eu.sub.GM]            4e-04    7e-04    3e-04    3e-03    4e-03
[Es.sub.GM]            4e-04    8e-04    3e-04    2e-03    2e-02
[Eh.sub.GM]            4e-04    2e-03    9e-04    3e-03    2e-02

[[member of].sub.f]    7e-07    2e-07    6e-07    9e-07    2e-06
[[delta].sub.f]        8e-07    2e-07    3e-07    4e-07    2e-06
[[member of].sub.      7e-16    3e-16    5e-16    6e-16    3e-16
[infinity]]

[THETA]                8e-02    7e-02    6e-02    2e-02    5e-03
[rho]                    20       20       20       90      300

Table 7.9

([alpha], [beta]) = ([10.sup.-5], 1): eigenvalue errors,
errors  of deflating subspaces, and [THETA], [rho]

1        2        3        4        5

[Eu.sub.max]           8e-08    8e-08    9e-08    1e-07    2e-07
[Es.sub.max]           2e-07    4e-08    7e-08    6e-08    1e-07
[Eh.sub.max]           2e-07    8e-08    2e-07    2e-07    2e-07
[Eu.sub.min]           8e-08    8e-08    9e-08    1e-07    2e-07
[Es.sub.min]           2e-07    4e-08    7e-08    6e-08    1e-07
[Eh.sub.min]           8e-08    7e-08    2e-07    4e-08    2e-07
[Eu.sub.GM]            8e-08    8e-08    9e-08    1e-07    2e-07
[Es.sub.GM]            2e-07    4e-08    7e-08    6e-08    1e-07
[Eh.sub.GM]            2e-07    7e-08    2e-07    7e-08    2e-07

[[member of].sub.f]    2e-09    2e-09    3e-09    8e-10    3e-10
[[delta].sub.f]        1e-09    3e-09    3e-09    9e-10    4e-10
[[member of].sub.      5e-16    4e-16    4e-16    2e-15    1e-15
[infinity]]

[THETA]                3e-01    4e-01    4e-01    6e-01    2e-01
[rho]                    5        3        3        2        7

6        7        8        9        10

[Eu.sub.max]           2e-07    2e-07    2e-07    2e-07    1e-06
[Es.sub.max]           2e-07    5e-08    2e-07    2e-07    8e-07
[Eh.sub.max]           4e-07    3e-07    4e-07    6e-07    2e-06
[Eu.sub.min]           2e-07    2e-07    2e-07    2e-07    1e-06
[Es.sub.min]           2e-07    5e-08    2e-07    2e-07    8e-07
[Eh.sub.min]           3e-07    3e-07    2e-07    5e-07    1e-06
[Eu.sub.GM]            2e-07    2e-07    2e-07    2e-07    1e-06
[Es.sub.GM]            2e-07    5e-08    2e-07    2e-07    8e-07
[Eh.sub.GM]            3e-07    3e-07    3e-07    5e-07    2e-06

[[member of].sub.f]    8e-11    7e-09    3e-10    4e-09    3e-09
[[delta].sub.f]        2e-10    2e-09    3e-10    2e-09    3e-09
[[member of].sub.      3e-16    5e-16    6e-16    3e-16    8e-16
[infinity]]

[THETA]                3e-01    2e-01    5e-01    1e-01    2e-01
[rho]                    5        7        2        20       10

Table 7.10

([alpha], [beta]) = ([10.sup.-5], [10.sup.-6]): eigenvalue
errors, errors of deflating subspaces, and [THETA], [rho]

1        2        3        4        5

[Eu.sub.max]          5e-02    6e-02    7e-02    8e-02    1e-01
[Es.sub.max]          5e-02    5e-02    8e-02    8e-02    8e-02
[Eh.sub.max]          5e-02    6e-02    3e-01    3e-01    2e-01
[Eu.sub.min]          5e-02    6e-02    7e-02    8e-02    1e-01
[Es.sub.min]          5e-02    5e-02    8e-02    8e-02    8e-02
[Eh.sub.min]          2e-02    3e-02    2e-01    9e-02    5e-02
[Eu.sub.GM]           5e-02    6e-02    7e-02    8e-02    1e-01
[Es.sub.GM]           5e-02    5e-02    8e-02    8e-02    8e-02
[Eh.sub.GM]           3e-02    5e-02    3e-01    2e-01    8e-02

[[member of].sub.f]   3e-10    9e-10    1e-09    2e-09    3e-09
[[delta].sub.f]       3e-10    8e-10    4e-10    2e-09    2e-09
[[member of].sub.     4e-16    7e-16    5e-16    6e-16    3e-16
[infinity]]

[THETA]               4e-01    5e-01    2e-01    2e-01    4e-01
[rho]                   3        3        6        6        4

6        7        8        9        10

[Eu.sub.max]          1e-01    2e-01    3e-01    4e-01    5e-01
[Es.sub.max]          2e-01    7e-02    3e-01    2e-01    6e-01
[Eh.sub.max]          9e-02    2e-01    6e-01    3e-01    1e+00
[Eu.sub.min]          1e-01    2e-01    3e-01    4e-01    5e-01
[Es.sub.min]          2e-01    7e-02    3e-01    2e-01    5e-01
[Eh.sub.min]          8e-02    2e-01    3e-01    3e-02    2e-01
[Eu.sub.GM]           1e-01    2e-01    3e-01    4e-01    5e-01
[Es.sub.GM]           2e-01    7e-02    3e-01    2e-01    6e-01
[Eh.sub.GM]           9e-02    2e-01    5e-01    9e-02    4e-01

[[member of].sub.f]   2e-10    1e-09    3e-09    5e-09    3e-08
[[delta].sub.f]       3e-10    6e-10    3e-09    2e-09    3e-08
[[member of].sub.     9e-16    2e-16    4e-16    7e-16    5e-16
[infinity]]

[THETA]               3e-01    6e-01    1e-01    6e-02    3e-02
[rho]                   5        2        20       20       50
```
COPYRIGHT 2015 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.