Real Fast Structure-Preserving Algorithm for Eigenproblem of Complex Hermitian Matrices

Abstract

It is well known that the flops for complex operations are usually 4 times of real cases. In the paper, using real operations instead of complex, a real fast structure-preserving algorithm for eigenproblem of complex Hermitian matrices is given. We make use of the real symmetric and skew-Hamiltonian structure transformed by Wilkinson's way, focus on symplectic orthogonal similarity transformations and their structure-preserving property, and then reduce it into a two-by-two block tridiagonal symmetric matrix. Finally a real algorithm can be quickly obtained for eigenvalue problems of the original Hermitian matrix. Numerical experiments show that the fast algorithm can solve real complex Hermitian matrix efficiently, stably, and with high precision.

1. Introduction

We consider eigenvalue problems of the complex Hermitian matrix of the form 𝐻𝑥=𝜆𝑥,𝐻∈ℂ^𝑛×𝑛,𝑥∈ℂ^𝑛, (1) where 𝜆 is a real scalar and 𝐻∈ℂ^𝑛×𝑛 is a complex Hermitian matrix of the form 𝐻=𝐴+𝑖𝐵, (2) where 𝐴 is a real symmetric matrix; that is, 𝐴^𝑇=𝐴∈ℝ^𝑛 and 𝐵 is a real skew-symmetric nonzero matrix, that is, 𝐵=−𝐵^𝑇≠0∈ℝ^𝑛, given by 𝐻𝐴=^𝑇+𝐻2,𝐵=𝐻−𝐻^𝑇,2𝑖 (3) and 𝑥∈ℂ^𝑛 can be written as 𝑥=𝑢+𝑖𝑣, where 𝑢, 𝑣∈ℝ^𝑛.

We know that the property of the complex Hermitian matrix is also possessed by the symmetry matrix; therefore little work has been done on complex Hermitian matrix. However, from the view of numerical computation, the biggest difference between these two types of matrices is a complex Hermitian matrix that involves complex operations, while the latter does not. It is well known that a complex Hermitian matrix is unitary similar to a real diagonal matrix; that is, it has real eigenvalue but complex eigenvector. Because the flops for the complex operation are usually 4 times of the real case, so how to reduce the computation cost of the eigenvalue problem of complex Hermitian matrices is very meaningful and worthy of research.

In 1965, the famous numerical analysis expert Wilkinson [1] gave a way to transform a complex Hermitian matrix into real case. To be more precise, if (𝜆,𝑢+𝑖𝑣) is the eigenpair of 𝐻, then 𝐻𝑥=(𝐴+𝑖𝐵)(𝑢+𝑖𝑣)=𝜆(𝑢+𝑖𝑣). (4)

Since the eigenvalues of 𝐻 are real, the previous equation can be written as 𝑢𝑣𝑢𝑣𝐴−𝐵𝐵𝐴=𝜆. (5)

So the eigenproblem (1) is transformed into the eigenproblem of the real symmetry matrix 𝑆=𝐴−𝐵𝐵𝐴∈ℝ^2𝑛×2𝑛. (6)

Note that the matrix 𝑆 in (6) is not only real symmetry but also skew-Hamiltonian. Therefore, in the following discussion, inspired by the algorithms for the eigenproblem of Hamiltonian (skew-Hamiltonian) matrix [2–11], we devise an algorithm by taking full advantage of the special structure of (6).

In this paper we first prove that orthogonal symplectic similarity transformations not only preserve the symmetry and skew-Hamiltonian structure of (6), but also eigenvectors of (6), since the eigenvectors of (6) have special structure; that is, it is also orthogonal symplectic. Therefore, by these transformations, a real stable, accurate, and fast method is devised to reduce (6) into a block symmetry tridiagonal matrix, and then we get the eigenproblem of the original Hermitian matrix. When implementing the method, although its dimension doubles the one of 𝐻, we can make full use of the special structure of 𝑆 to avoid superfluous calculation and storage.

The remainder of the paper is organized as follows. In Section 2, we show the basic orthogonal symplectic matrix and its algorithm. In Section 3, we present the overall real algorithm for eigenproblem of complex Hermitian matrix. Some implementation aspects are discussed in Section 4. Finally, in Section 5, the results of numerical experiments are discussed to end the paper.

2. Basic Orthogonal Symplectic Matrix and Its Algorithm

It is well known that by choosing a proper unitary vector 𝑤, Householder transformation can zero some entries of a given vector; likewise, we defined Householder symplectic transformation which has the following form.

Definition 1. A real 2𝑛×2𝑛 matrix 𝑆 is a symplectic matrix if 𝑆𝐽𝑆^𝑇=𝐽, where 𝑆^𝑇 is the transpose of 𝑆 and 𝐽 is the orthogonal matrix: 𝐽=0𝐼_𝑛−𝐼_𝑛0∈ℝ^2𝑛×2𝑛. (7) Here 𝐼_𝑛 is the identity 𝑛×𝑛 matrix and 0 is the zero 𝑛×𝑛 matrix.

Definition 2. Letting 𝑘∈{1,2,…,𝑛}, a matrix 𝐻(𝑘,𝑤)∈ℝ^2𝑛×2𝑛 is called Householder symplectic matrix, if it has the form 𝐻()=𝐼𝑘,𝑤diag_𝑘−10𝐼,𝑃0diag_𝑘−1,𝑃∈ℝ^2𝑛×2𝑛, (8) where 𝑃=𝐼−2𝑤𝑤^𝑇𝑤^𝑇𝑤,𝑤∈ℝ^{𝑛−𝑘+1}. (9) Obviously, 𝐻=diag(𝐼_𝑘−1,𝑃) is Householder transformation of order 𝑛, when 𝑤=0, 𝐻(𝑘,𝑤)=𝐼_2𝑛.

Note that Householder symplectic matrix 𝐻(𝑘,𝑤) is just a direct sum of two “ordinary” 𝑛-by-𝑛 Householder matrices, specially when 𝑤=0, 𝐻(𝑘,𝑤)=𝐼_2𝑛. So based on the algorithm of “ordinary” Householder transformation, we obtain the following algorithm for Householder symplectic transformation 𝐻(𝑘,𝑤) of order 2𝑛.

Algorithm 3 (Householder symplectic transformation [8]). Given 𝑘(1≤𝑘<𝑛) and 𝑦,𝑧∈ℝ^𝑛, the following algorithm determines 𝑤=(𝑤_𝑘,…,𝑤_𝑛)^𝑇, such that if 𝐻()𝑦𝑧=𝑣𝑥𝑘,𝑤, (10) then 𝑥_𝑖=0,𝑖=𝑘+1,…,𝑛: 𝜎∶=𝑧²_𝑘+𝑧²_𝑘+1+⋯+𝑧²_𝑛^1/2,𝑤_𝑘∶=𝑧_𝑘+sign𝑧_𝑘𝑤𝜎,_𝑖∶=𝑧_𝑖 (11) end.

𝐻(𝑘,𝑤) is completely determined by 𝑘th to 𝑛th components of vector 𝑧; it makes 𝑥 satisfy 𝑥_𝑖=0 for 𝑖=𝑘+1,…,𝑛. Note that by interchanging the roles of 𝑦 and 𝑧, the previous algorithm can be used to determine 𝐻(𝑘,𝑤) such that 𝑣_𝑖=0 for 𝑖=𝑘+1,…,𝑛.

Householder symplectic transformation can be used to zero large portions of a vector as Householder transformation; likewise, we can define Givens symplectic transformation, which can be used to zero single entries.

Definition 4. Given 𝑘∈{1,2,…,𝑛}, a matrix 𝐽(𝑘,𝜃)∈ℝ^2𝑛×2𝑛 is called Givens symplectic matrix, if it has the form 𝐽()==𝐼𝑘,𝜃𝐶𝑆−𝑆𝐶_𝑘−10_𝑘−1𝐼cos𝜃sin𝜃_𝑛−𝑘0_𝑛−𝑘0_𝑘−1𝐼_𝑘−10−sin𝜃cos𝜃_𝑛−𝑘𝐼_𝑛−𝑘, (12) where 𝐶=diag𝐼_𝑘−1,cos𝜃,𝐼_𝑛−𝑘,𝑆=diag0_𝑘−1,sin𝜃,0_𝑛−𝑘. (13)
𝐽(𝑘,𝜃) is an “ordinary” 2𝑛-by-2𝑛 Givens rotation that rotates in planes 𝑘 and 𝑛+𝑘. The corresponding algorithm is as follows.

Algorithm 5 (Givens symplectic transformation [8]). Given 𝑘(1≤𝑘<𝑛) and 𝑦,𝑧∈ℝ^𝑛, the following algorithm determines 𝑐=cos𝜃, 𝑠=sin𝜃 such that if𝐽()𝑦𝑧=𝑣𝑥𝑘,𝑤, (14) then 𝑥_𝑘=0: 𝜎∶=𝑦²_𝑘+𝑧²_𝑘^1/2,if𝜎=0then𝑐∶=1and𝑠∶=0else𝑐∶=𝑦_𝑘/𝜎and𝑠∶=𝑧_𝑘/𝜎end. (15)

3. The Reduction of the Symmetry and Skew-Hamiltonian Matrix 𝑆

Now we consider some properties of the symmetry and skew-Hamiltonian matrix 𝑆 of (6).

Lemma 6. If (𝜆,[^𝑢_𝑣]) is the eigenpair of 𝑆=𝐴−𝐵𝐵𝐴, (16) then (𝜆,[_𝑢^−𝑣]) is also its eigenpair.

Proof. If 𝑢𝑣𝑢𝑣𝐴−𝐵𝐵𝐴=𝜆, (17) then we have 𝑢𝑢𝐴−𝐵𝐵𝐴−𝑣=𝜆−𝑣. (18)

According to Lemma 6, there are two different eigenvectors corresponding to the same eigenvalue, so it can be easy to prove that 𝐴−𝐵𝐵𝐴𝑈−𝑉𝑉𝑈=Λ𝑈−𝑉𝑉𝑈, (19) where Λ={𝜆₁,…,𝜆_𝑛} is the set of eigenvalues of (6) and 𝑈−𝑉𝑉𝑈∈ℝ^2𝑛×2𝑛 (20) is the eigenvectors of (6). Note that the matrix in (20) is symplectic orthogonal. The following lemma tells us that the product of symplectic orthogonal matrices is also symplectic orthogonal.

Lemma 7. Supposing that 2𝑛-by-2𝑛 real matrices of the form 𝑄𝑄=¹𝑄₂−𝑄₂𝑄₁,𝑄𝑄=¹𝑄₂−𝑄₂𝑄₁ (21) are both orthogonal and symplectic, then their product 𝑄𝑄 is both orthogonal and symplectic.

Lemma 7 shows that the special structure of eigenvectors of 𝑆 can be preserved by the symplectic orthogonal similarity transformations. The following theorem is the foundation of the fast and stable structure-preserving algorithm given in the following discussion.

Theorem 8. Suppose 𝐴^𝑇=𝐴∈𝑅^𝑛×𝑛,𝐵^𝑇=−𝐵≠0∈𝑅^𝑛×𝑛, and 𝑄𝑄=¹𝑄₂−𝑄₂𝑄₁ (22) is an orthogonal symplectic matrix. If 𝐴₁₁𝐴₁₂𝐴₂₁𝐴₂₂ (23) satisfy 𝐴₁₁𝐴₁₂𝐴₂₁𝐴₂₂=𝑄₁𝑄₂−𝑄₂𝑄₁^𝑇𝐴−𝐵𝑄𝐵𝐴¹𝑄₂−𝑄₂𝑄₁, (24) then 𝐴₁₁=𝐴₂₂,𝐴₂₁=−𝐴₁₂,𝐴^𝑇₁₁=𝐴₁₁,𝐴^𝑇₁₂=−𝐴₁₂. (25)

Proof. It is easy to obtain 𝐴₁₁=𝑄^𝑇₁𝐴𝑄₁+𝑄^𝑇₁𝐵𝑄₂−𝑄^𝑇₂𝐵𝑄₁+𝑄^𝑇₂𝐴𝑄₂,𝐴₂₂=𝑄^𝑇₂𝐴𝑄₂−𝑄^𝑇₂𝐵𝑄₁+𝑄^𝑇₁𝐴𝑄₁+𝑄^𝑇₁𝐵𝑄₂,𝐴₁₂=𝑄^𝑇₁𝐴𝑄₂−𝑄^𝑇₁𝐵𝑄₁−𝑄₂𝐵𝑄₂−𝑄^𝑇₂𝐴𝑄₁,𝐴₂₁=𝑄^𝑇₂𝐴𝑄₁+𝑄^𝑇₂𝐵𝑄₂+𝑄^𝑇₁𝐵𝑄₁−𝑄^𝑇₁𝐴𝑄₂. (26) Then 𝐴₁₁=𝐴₂₂,𝐴₂₁=−𝐴₁₂. (27)
Moreover 𝐴^𝑇=𝐴,𝐵^𝑇=−𝐵, hence we have 𝐴^𝑇₁₁=𝐴₁₁,𝐴^𝑇₁₂=−𝐴₁₂. (28)

The previous theorem shows that the symplectic similarity transformation preserves the symmetry and skew-Hamiltonian structure of 𝑆=,𝐴𝐴−𝐵𝐵𝐴^𝑇=𝐴,𝐵^𝑇=−𝐵. (29)

This structure-preserving property provides a way to transform the matrix into a special skew-Hamiltonian matrix: 𝑇00𝑇∈ℝ^2𝑛×2𝑛, (30) where 𝑇∈ℝ^𝑛×𝑛 is a real symmetry tridiagonal matrix. So we have the following results, which provide a theoretical basis for the real fast structure-preserving algorithm for eigenproblems of a Hermitian matrix.

Theorem 9. Letting 𝑆=^𝐴−𝐵_𝐵𝐴∈ℝ^2𝑛×2𝑛, where 𝐴^𝑇=𝐴,𝐵^𝑇=−𝐵, then there exists an orthogonal symplectic matrix 𝑄=^𝑄1𝑄2_{−𝑄2𝑄1}∈ℝ^2𝑛×2𝑛 such that 𝑄^𝑇𝑆𝑄=𝑇00𝑇∈ℝ^2𝑛×2𝑛, (31) where 𝑇 is a symmetric tridiagonal matrix. Suppose that 𝑇 has the eigenvalue decomposition 𝑄₁^𝑇𝑇𝑄₁=Λ, (32) that is 𝑄₁00𝑄₁^𝑇𝑇0𝑄0𝑇¹00𝑄₁=.Λ00Λ (33) Then eigenvalues of 𝑆 are 𝜆(𝑆)=𝜆(𝑇)⊔𝜆(𝑇), (34) and eigenvectors are the columns of 𝑄𝑄𝑄=¹𝑄₂−𝑄₂𝑄₁𝑄₁00𝑄₁=𝑄₁𝑄₁𝑄₂𝑄₁−𝑄₂𝑄₁𝑄₁𝑄₁. (35)

We now consider how an orthogonal symplectic 𝑄∈ℝ^2𝑛×2𝑛 can be determined to reduce 𝑆=^𝐴−𝐵_𝐵𝐴 (𝐴^𝑇=𝐴,𝐵^𝑇=−𝐵) to 𝑄^𝑇𝑆𝑄=𝑇00𝑇, (36) where 𝑇 is a symmetry tridiagonal matrix. It is worth illustrating a few steps of the algorithm before presenting it in full generality. We depict an 8-by-8 symmetry and skew-Hamiltonian matrix (𝑛=4) as follows: 𝑆==𝐴−𝐵𝐵𝐴𝑥𝑥𝑥𝑥0𝑥𝑥𝑥𝑥𝑥𝑥𝑥𝑥0𝑥𝑥𝑥𝑥𝑥𝑥𝑥𝑥0𝑥𝑥𝑥𝑥𝑥𝑥𝑥𝑥00𝑥𝑥𝑥𝑥𝑥𝑥𝑥𝑥0𝑥𝑥𝑥𝑥𝑥𝑥𝑥𝑥0𝑥𝑥𝑥𝑥𝑥𝑥𝑥𝑥0().𝑥𝑥𝑥𝑥𝑛=4 (37) The zero diagonals in 𝐵 follow from skew symmetry. It is easy to verify that the (2, 1) and (1, 2) block of this matrix can be zeroed by applying a sequence of Householder and Givens symplectic similarity transformations.

The first step is to zero 𝑏₃₁ and 𝑏₄₁ using a Householder symplectic 𝐻₁=𝐻(2,𝑤). This matrix can be determined by executing Algorithm 3 with 𝑘=2, 𝑦=𝐴𝑒₁, and 𝑧=𝐵𝑒₁. Note that after performing, only the rows and columns 2, 3, and 4 of 𝐴, 𝐵, and −𝐵 are affected, the result being 𝑆=𝐴−𝐵𝐵𝐴∶=𝐻₁𝑆𝐻^𝑇₁=𝑥𝑥𝑥𝑥0𝑥00𝑥𝑥𝑥𝑥𝑥0𝑥𝑥𝑥𝑥𝑥𝑥0𝑥0𝑥𝑥𝑥𝑥𝑥0𝑥𝑥00𝑥00𝑥𝑥𝑥𝑥𝑥0𝑥𝑥𝑥𝑥𝑥𝑥0𝑥0𝑥𝑥𝑥𝑥𝑥0𝑥𝑥0().𝑥𝑥𝑥𝑥𝑛=4 (38) Then (1,3) and (1,4) positions of 𝐵 are zeroed; meanwhile the same positions of −𝐵 are zeroed because Householder symplectic similarity transformation preserves the structure of 𝑆.

The next step is to zero 𝑏₂₁ using Givens similarity 𝐽₁=𝐽(2,𝜃), which can be achieved by applying Algorithm 5 with 𝑘=2, 𝑦=𝐴𝑒₁, and 𝑧=𝐵𝑒₁. And notice that only the second row and column of 𝐴, 𝐵 and −𝐵 be affected. This implies that 𝑆=𝐴−𝐵𝐵𝐴∶=𝐽₁𝑆𝐽^𝑇₁=𝑥𝑥𝑥𝑥0000𝑥𝑥𝑥𝑥00𝑥𝑥𝑥𝑥𝑥𝑥0𝑥0𝑥𝑥𝑥𝑥𝑥0𝑥𝑥00000𝑥𝑥𝑥𝑥00𝑥𝑥𝑥𝑥𝑥𝑥0𝑥0𝑥𝑥𝑥𝑥𝑥0𝑥𝑥0().𝑥𝑥𝑥𝑥𝑛=4 (39)

Next we compute a Householder symplectic 𝐻(2,𝑤) to zero 𝑎₃₁ and 𝑎₄₁, which can be achieved by applying Algorithm 3 with 𝑘=2, 𝑦=𝐵𝑒₁, and 𝑧=𝐴𝑒₁. Denote 𝐺₁=𝐻(2,𝑤), the result being 𝑆=𝐴−𝐵𝐵𝐴∶=𝐺₁𝑆𝐺^𝑇₁=𝑥𝑥000000𝑥𝑥𝑥𝑥00𝑥𝑥0𝑥𝑥𝑥0𝑥0𝑥0𝑥𝑥𝑥0𝑥𝑥00000𝑥𝑥0000𝑥𝑥𝑥𝑥𝑥𝑥0𝑥0𝑥0𝑥𝑥𝑥0𝑥𝑥0().0𝑥𝑥𝑥𝑛=4 (40)

Note that the (1,3) and (1,4) positions of 𝐴 are zeroed. This completes the zeroing in the first column of 𝑆. Likewise we can zero the 2th of 𝑆 with the property that 𝑆==𝐴−𝐵𝐵𝐴𝑥𝑥000000𝑥𝑥𝑥000000𝑥𝑥𝑥000𝑥00𝑥𝑥00𝑥00000𝑥𝑥000000𝑥𝑥𝑥0000𝑥0𝑥𝑥𝑥00𝑥0().00𝑥𝑥𝑛=4 (41)

Notice that only the (3,4) and (4,3) positions of 𝐵 are nonzero; they are zeroed only by applying Givens symplectic similarity transformation 𝐺₃=𝐽(4,𝜃) to 𝑆 and we get 𝑆=𝐺₃𝑆𝐺^𝑇₃==𝑇00𝑇𝑥𝑥000000𝑥𝑥𝑥000000𝑥𝑥𝑥000000𝑥𝑥00000000𝑥𝑥000000𝑥𝑥𝑥000000𝑥𝑥𝑥0000().00𝑥𝑥𝑛=4 (42)

This completes the zeros of 𝑆 with 𝑛=4. The similar method can be applied to zero the general matrix of order 2𝑛. Overall we have the following procedure “Reduction to block symmetry tridiagonal matrix.” Given matrix 𝑆=𝐴−𝐵𝐵𝐴with𝐴^𝑇=𝐴,𝐵^𝑇=−𝐵. (43)

The following algorithm updates 𝐴 to the symmetry tridiagonal matrix 𝑇; meanwhile 𝐵 and −𝐵 are reduced into an 𝑛×𝑛 zero matrix; moreover we have 𝜆(𝑆)=𝜆(𝑇)⊔𝜆(𝑇).

Algorithm 10 (reduction to block symmetry tridiagonal matrix). For 𝑘=1,…,𝑛−2 do:
 if 𝑘≤𝑛−2 do:
(1)  Apply Algorithm 3 with 𝑦=𝐴𝑒_𝑘 and 𝑧=𝐵𝑒_𝑘 to determine 𝐻_𝑘=𝐻(𝑘+1,𝑤). Update is as follows: 𝐴−𝐵𝐵𝐴∶=𝐻_𝑘𝐻𝐴−𝐵𝐵𝐴^𝑇𝑘. (44)
(2)  Apply Algorithm 5 with 𝑦=𝐴𝑒_𝑘 and 𝑧=𝐵𝑒_𝑘 to determine 𝐽_𝑘=𝐽(𝑘+1,𝜃).  Updates are as follows: 𝐴−𝐵𝐵𝐴∶=𝐽_𝑘𝐽𝐴−𝐵𝐵𝐴^𝑇𝑘. (45)
(3)  Apply Algorithm 3 with 𝑦=𝐵𝑒_𝑘,𝑧=𝐴𝑒_𝑘 to determine 𝐺_𝑘=𝐻(𝑘+1,𝑤). Updates are as follows: 𝐴−𝐵𝐵𝐴∶=𝐺_𝑘𝐺𝐴−𝐵𝐵𝐴^𝑇𝑘. (46)
If 𝑘=𝑛−1 do the folowing.
(4) Apply Algorithm 5 with 𝑦=𝐴𝑒_𝑛−1 and 𝑧=𝐵𝑒_𝑛−1 to determine 𝐽_𝑛−1=𝐽(𝑛,𝜃). Update is as follows: 𝐴−𝐵𝐵𝐴∶=𝐽_𝑛−1𝐽𝐴−𝐵𝐵𝐴^𝑇𝑛−1. (47)

After completion, we get the block symmetry tridiagonal matrix 𝑇00𝑇, where 𝑇 is a symmetry tridiagonal matrix. The orthogonal matrix 𝑄^𝑇 is given by 𝑄^𝑇=𝐽_𝑛−1𝐺_𝑛−2𝐽_𝑛−2𝐻_𝑛−2⋯𝐺₂𝐽₂𝐻₂𝐺₁𝐽₁𝐻₁. (48)

The method of reduction 𝑆 into block symmetry tridiagonal matrix is given by Algorithm 10; in the former 𝑛−2 steps, each step we premultiply by one Givens symplectic similarity and two Householder symplectic similarity and by postmultiply their transpose. Finally in the 𝑛−1 step, only the Givens symplectic similarity is done. According to (48), it is clear that 𝑄^𝑇 is the product of 3(𝑛−2)+1 orthogonal matrices, which is backward stability.

Algorithm 11 (the overall algorithm). (1) Transform Hermitian matrix 𝐻 into the symmetry and skew-Hamiltonian matrix: 𝑆=𝐴−𝐵𝐵𝐴where𝐴^𝑇=𝐴,𝐵^𝑇=−𝐵. (49)
(2)  Apply Algorithm 10 to reduce 𝑆 into a block symmetry tridiagonal matrix: 𝑄^𝑇𝐻𝑄=𝑇00𝑇. (50)
(3)  Compute eigenpairs (Λ, 𝑋) of 𝑇 via 𝑄𝑅.
(4)  Setting 𝑄₁₁𝑄₂₁=𝑄₁𝑄₂−𝑄₂𝑄₁𝑋0=𝑄₁𝑋−𝑄₂𝑋, (51) then Λ is eigenvalue of 𝐻, and 𝑄₁₁+𝑖𝑄₂₁ is the eigenvector of 𝐻.

4. Implementation Aspects

When implementing Algorithm 10, the matrix 𝑆 must be taken full advantage of in order to avoid superfluous calculation and storage. In this section there are two parts; the first is to consider the computational cost and storage for the reduction to this block symmetry tridiagonal matrix, and the second is the computational cost and storage of the eigenvector 𝑄^𝑇.

Due to the special structure 𝐴=𝐴^𝑇 and 𝐵^𝑇=−𝐵, shared by 𝐴 and 𝐵, only 𝑛² workspace is required for matrix 𝑆=𝐴−𝐵𝐵𝐴where𝐴^𝑇=𝐴,𝐵^𝑇=−𝐵. (52)

Now we consider the computation cost for reducing 𝑆 to the block symmetry tridiagonal form, for Householder symplectic similarity, let 𝐻_𝑘=𝐻(𝑘+1,𝑤)=diag𝐼_𝑘,𝑃,𝐼_𝑘,𝐴,𝑃𝑆=¹¹𝐴₂₁𝐵₁₁𝐵₂₁𝐴^𝑇₂₁𝐴₂₂−𝐵^𝑇₂₁𝐵₂₂𝐵^𝑇₁₁−𝐵^𝑇₂₁𝐴₁₁𝐴^𝑇₁₂𝐵^𝑇₂₁𝐵^𝑇₂₂𝐴₁₂𝐴₂₂𝑘𝑘𝑛−𝑘𝑛−𝑘𝑘𝑛−𝑘𝑘𝑛−𝑘, (53)(54) then 𝐻_𝑘𝑆𝐻_𝑘=𝐴₁₁𝑃𝐴^𝑇₁₂𝐵₁₁𝑃𝐵₂₁𝐴₁₂𝑃𝑃𝐴₂₂𝑃−𝐵^𝑇₂₁𝑃𝑃𝐵₂₂𝑃𝐵^𝑇₁₁−𝑃𝐵^𝑇₂₁𝐴₁₁𝑃𝐴^𝑇₁₂𝐵^𝑇₂₁𝑃𝑃𝐵^𝑇₂₂𝐴₁₂𝑃𝐴₂₂𝑃𝑘𝑘𝑛−𝑘𝑛−𝑘𝑘𝑛−𝑘𝑘𝑛−𝑘. (55)

In 𝑘th step, updating 𝑆 to 𝐻_𝑘𝑆𝐻_𝑘 is equal to two ordinary Householder similarities for 𝑛×𝑛 matrix, one for symmetry matrix 𝐴 and one for skew-symmetry matrix 𝐵, so we only need to compute 𝑃𝐴^𝑇₁₂, 𝑃𝐴₂₂𝑃, 𝑃𝐵₂₁, and 𝑃𝐵₂₂𝑃. Likewise, we can update 𝑆 to 𝐺_𝑘𝑆𝐺_𝑘 by the same way.

Algorithm 12 (Householder symplectic similarity). For𝑘=1∶𝑛−2(𝐵(,()𝑣,𝛽=house𝑘+1∶𝑛,𝑘))𝑤=𝛽𝐵𝑘+1∶𝑛,𝑘∶𝑛^𝑇𝐵()()𝑣,𝑘+1∶𝑛,𝑘∶𝑛=𝐵𝑘+1∶𝑛,𝑘∶𝑛−𝑣𝑤^𝑇,𝑤=𝛽𝐵(𝑘∶𝑛,𝑘+1∶𝑛)𝑣^𝑇,𝐵(𝑘∶𝑛,𝑘+1∶𝑛)=𝐵(𝑘∶𝑛,𝑘+1∶𝑛)−𝑤𝑣^𝑇,𝑤=𝛽𝐴(𝑘+1∶𝑛,𝑘∶𝑛)^𝑇𝐴()𝑣,𝑝+1∶𝑛,𝑝∶𝑛=𝐴𝑘+1∶𝑛,𝑘∶𝑛−𝑣𝑤^𝑇,()𝐴()()𝑤=𝛽𝐴𝑘∶𝑛,𝑘+1∶𝑛𝑣,𝑘∶𝑛,𝑘+1∶𝑛=𝐴𝑘∶𝑛,𝑘+1∶𝑛−𝑤𝑣^𝑇 (56) end.

Updating by Givens symplectic similarity is equally simple; setting 𝑐=cos(𝜃), 𝑠=sin(𝜃)  then 𝐽_𝑘()=𝐼=𝐽𝑘+1,𝜃^𝑘0_𝑘𝐼𝑐𝑠_{𝑛−𝑘−1}0_{𝑛−𝑘−1}0_𝑘𝐼_𝑘0−𝑠𝑐_{𝑛−𝑘−1}𝐼_{𝑛−𝑘−1},𝐴𝑆=¹¹𝐴₂₁𝐴₃₁𝐵₁₁𝐵₂₁𝐵₃₁𝐴₁₂𝐴₂₂𝐴₃₂𝐵₁₂𝐵₂₂𝐵₃₂𝐴₁₃𝐴₂₃𝐴₃₃𝐵₁₃𝐵₂₃𝐵₃₃−𝐵₁₁−𝐵₂₁−𝐵₃₁𝐴₁₁𝐴₂₁𝐴₃₁−𝐵₁₂−𝐵₂₂−𝐵₃₂𝐴₁₂𝐴₂₂𝐴₃₂−𝐵₁₃−𝐵₂₃−𝐵₃₃𝐴₁₃𝐴₂₃𝐴₃₃11.𝑘+1𝑛−𝑘−2𝑘+1𝑛−𝑘−2𝑘+11𝑛−𝑘−2𝑘1𝑛−𝑘−2 (57)

Denoting 𝑆_𝐽=(𝐽_𝑘𝑆𝐽^𝑇_𝑘), then we have 𝑆_𝐽()=𝐴1∶2𝑛,1∶𝑛₁₁𝑐𝐴₂₁+𝑠𝐵₂₁𝐴₃₁𝐵₁₁𝑐𝐵₂₁−𝑠𝐴₂₁𝐵₃₁𝑐𝐴₁₂−𝑠𝐵₁₂𝐴₂₂𝑐𝐴₃₂−𝑠𝐵₃₂𝑐𝐵₁₂+𝑠𝐴₁₂𝐵₂₂𝑐𝐵₃₂+𝑠𝐴₃₂𝐴^𝑇₃₁𝐴^𝑇₃₂𝐴₃₃𝐵₁₃𝑐𝐵₂₃−𝑆𝐴₂₃𝐵₃₃11.𝑘+1𝑛−𝑘−2𝑘+1𝑛−𝑘−2𝑘+11𝑛−𝑘−2 (58)

Note that updating 𝑆 into 𝐽_𝑘𝑆𝐽^𝑇_𝑘 such that 𝐵(𝑘+1,𝑘)=0 affects the rows  (resp. the columns) 𝑘+1, 𝑛+𝑘+1 of 𝑆, that is, the rows(columns) 𝑘+1, 𝑛+𝑘+1 of 𝐴 and 𝐵. According to (58), we only need to compute 𝑐𝐴₂₁+𝑠𝐵₂₁, 𝑐𝐴₃₂−𝑠𝐵₃₂, 𝑐𝐵₂₁−𝑠𝐴₂₁ and 𝑐𝐵₃₂+𝑠𝐴₃₂ because 𝑆_𝐽 is also symmetry and skew-Hamiltonian matrix. Therefore, only 2(𝑛−𝑘) flops are required in 𝑘th step, altogether ∑^𝑛−1_𝑘=12(𝑛−𝑘)=𝑛²−𝑛 flops. Bellow is the algorithm for 𝐽_𝑘𝑆𝐽^𝑇_𝑘 in 𝑘th step.

Algorithm 13 (Givens symplectic similarity). Consider ()(),[](),𝐴()()(),𝐴()(),𝑆()()(),𝐵()()(),𝐵()(),𝐴()(),𝐴()()𝐵()𝐵()𝑎=𝐴𝑘+1,𝑘,𝑏=𝐵𝑘+1,𝑘𝑐,𝑠=givens𝑎,𝑏𝑘+1,𝑘=𝑐𝐴𝑘+1,𝑘+𝑠𝐵𝑘+1,𝑘𝑘,𝑘+1=𝐴𝑘+1,𝑘for𝑖=𝑘+2∶𝑛𝑖,𝑘+1=𝑐𝐴𝑖,𝑘+1−𝑠𝐵𝑖,𝑘+1𝑖,𝑘+1=𝑐𝐵𝑖,𝑘+1+𝑠𝐴𝑖,𝑘+1𝑘+1,𝑖=−𝐵𝑖,𝑘+1𝑖,𝑘+1=𝑆𝑖,𝑘+1𝑘+1,𝑖=𝐴𝑖,𝑘+1end𝑘+1,𝑘=0,𝑘,𝑘+1=0. (59)

Overall (4/3)𝑛³×4+𝒪(𝑛²)=(16/3)𝑛³+𝒪(𝑛²) flops are required for reducing 𝑆 to 𝑇00𝑇.

This reduction includes 4 times computational cost of Householder tridiagonalized for 𝑛×𝑛 symmetry matrix and 𝑛²−𝑛 flops for Givens similarity transformation. Since one complex operation is four times of real operation, (16/3)𝑛³+𝒪(𝑛²) flops are required for complex Hermitian tridiagonalization.

Now we consider the computation and storage problem for eigenvectors 𝑄^𝑇=𝐽_𝑛−1𝐺_𝑛−2𝐽_𝑛−2𝐻_𝑛−2⋯𝐺₂𝐽₂𝐻₂𝐺₁𝐽₁𝐻₁. (60)

To obtain 𝑄^𝑇, we can premultiply only by 𝑛×𝑛 matrices by using the special structure of 𝐽_𝑘, 𝐺_𝑘, and 𝐻_𝑘 though they are 2𝑛-by-2𝑛 matrices. When premultiplying by Householder symplectics 𝐻_𝑘, let a 2𝑛×2𝑛 matrix 𝑄^𝑇 be of the form 𝑄^𝑇=𝑄_𝑁×𝑁−𝑈_𝑁×𝑁𝑈_𝑁×𝑁𝑄_𝑁×𝑁=𝑄₁₁𝑄₂₁𝑈₁₁𝑈₂₁𝑄₁₂𝑄₂₂𝑈₁₂𝑈₂₂−𝑈₁₁−𝑈₂₁𝑄₁₁𝑄₂₁−𝑈₁₂−𝑈₂₂𝑄₁₂𝑄₂₂𝑘𝑘.𝑛−𝑘𝑛−𝑘𝑘𝑛−𝑘𝑘𝑛−𝑘 (61)

So two 𝑛² workspace is required for 𝑄^𝑇, one for 𝑄_𝑁×𝑁 and one for 𝑈_𝑁×𝑁. Premultiplying by 𝐻_𝑘=diag(𝐼_𝑘+1,𝑃_𝐴,𝐼_𝑘+1,𝑃_𝐴) in 𝑘th step, we have 𝐻_𝑘𝑄^𝑇=𝑄₁₁𝑃_𝐴𝑄₂₁𝑈₁₁𝑃_𝐴𝑈₂₁𝑄₁₂𝑃_𝐴𝑄₂₂𝑈₁₂𝑃_𝐴𝑈₂₂−𝑈₁₁−𝑃_𝐴𝑈₂₁𝑄₁₁𝑃_𝐴𝑄₂₁−𝑈₁₂−𝑃_𝐴𝑈₂₂𝑄₁₂𝑃_𝐴𝑄₂₂𝑘𝑘𝑛−𝑘𝑛−𝑘𝑘𝑛−𝑘𝑘𝑛−𝑘. (62) Therefore, only 𝑄₂₁, 𝑄₂₂ and the real part of 𝑄^𝑇 are updated to 𝑃_𝐴𝑄₂₁, 𝑃_𝐴𝑄₂₂ and also imaginary part 𝑈₂₁, 𝑈₂₂ updated to 𝑃_𝐴𝑈₂₁, 𝑃_𝐴𝑈₂₂. The main work is computing the product of Householder matrix 𝐻=𝐼−𝛽𝑣𝑣^𝑇∈ℝ^{(𝑛−𝑘)×(𝑛−𝑘)} and [𝑄₂₁ 𝑄₂₂] ∈ℝ^{(𝑛−𝑘)×(𝑛)}, and the product of 𝐻 and [𝑈₂₁ 𝑈₂₂] ∈ℝ^{(𝑛−𝑘)×(𝑛)}, where the first columns of 𝑄₂₁ and 𝑈₂₁ are zero. So 4(𝑛−𝑘)(𝑛−𝑘−1) flops are required each step, altogether ∑^𝑛−2_𝑘=18(𝑛−𝑘)(𝑛−𝑘−1)=(8/3)𝑛³+𝒪(𝑛²) flops. Likewise premultiplying by 𝐺_𝑘 requires other (8/3)𝑛³+𝒪(𝑛²) flops.

Algorithm 14 (computation of 𝐻_𝑘𝑄^𝑇).  For 𝑘=1∶𝑛−2(𝐵(,()𝑣,𝛽=house𝑘+1∶𝑛,𝑘))𝑤=𝛽𝑄1𝑘+1∶𝑛,2∶𝑛^𝑇()()𝑣,𝑄1𝑘+1∶𝑛,2∶𝑛=𝑄1𝑘+1∶𝑛,2∶𝑛−𝑣𝑤^𝑇,𝑤=𝛽𝑄2(𝑘+1∶𝑛,2∶𝑛)^𝑇()()𝑣,𝑄2𝑘+1∶𝑛,2∶𝑛=𝑄2𝑘+1∶𝑛,2∶𝑛−𝑣𝑤^𝑇 (63) end.

For Givens symplectic matrix 𝐽_𝑘, set𝑄^𝑇=𝑄₁₁𝑄₂₁𝑄₃₁𝑈₁₁𝑈₂₁𝑈₃₁𝑄₁₂𝑄₂₂𝑄₃₂𝑈₁₂𝑈₂₂𝑈₃₂𝑄₁₃𝑄₂₃𝑄₃₃𝑈₁₃𝑈₂₃𝑈₃₃−𝑈₁₁−𝑈₂₁−𝑈₃₁𝑄₁₁𝑄₂₁𝑄₃₁−𝑈₁₂−𝑈₂₂−𝑈₃₂𝑄₁₂𝑄₂₂𝑄₃₂−𝑈₁₃−𝑈₂₃−𝑈₃₃𝑄₁₃𝑄₂₃𝑄₃₃𝑘1𝑘1𝑛−𝑘−1𝑛−𝑘−1𝑘1𝑛−𝑘−1𝑘1𝑛−𝑘−1. (64)

Premultiplying by 𝐽_𝑘, denoting 𝑄^𝑇_𝐽=𝐽_𝑘𝑄^𝑇, we have 𝑄^𝑇_𝐽()=𝑄1∶2𝑛,1∶𝑛₁₁𝑐𝑄₂₁+𝑠𝑈₂₁𝑄₃₁0𝑐𝑈₂₁−𝑠𝑄₂₁0𝑄₁₂𝑐𝑄₂₂+𝑠𝑈₂₂𝑄₃₂0𝑐𝑈₂₂−𝑠𝑄₂₂0𝑄₁₃𝑐𝑄₂₃+𝑠𝑈₂₃𝑄₃₃0𝑐𝑈₂₃−𝑠𝑄₂₃0𝑘1𝑘1.𝑛−𝑘−1𝑛−𝑘−1𝑘1𝑛−𝑘−1 (65) Therefore, each step updating the 𝑘+1th and 𝑛+𝑘+1th row of 𝑄^𝑇 requires 4(𝑛−1) flops, altogether 4(𝑛−1)² flops. Below is the algorithm.

Algorithm 15 (computation of 𝐽_𝑘𝑄^𝑇). For 𝑘=1∶𝑛−2()(),()(),()()()(),()(),()()(),()()()𝑇1=𝑄1,𝑇2=𝑄2,𝑄1𝑘+1,2∶𝑛=𝑐𝑄1𝑘+1,2∶𝑛𝑄2𝑘+1,2∶𝑛=𝑠𝑄2𝑘+1,2∶𝑛𝑄2𝑘+2∶𝑛,2∶𝑛=0,𝑄1𝑘+1∶𝑛,2∶𝑛=𝑄1𝑘+1∶𝑛,2∶𝑛+𝑄2𝑘+1∶𝑛,2∶𝑛𝑇1𝑘+1,2∶𝑛=−𝑠𝑇1𝑘+1,2∶𝑛𝑇1𝑘+2∶𝑛,2∶𝑛=0,𝑇2𝑘+1,2∶𝑛=𝑐𝑇2𝑘+1,2∶𝑛𝑄2𝑘+1∶𝑛,2∶𝑛=𝑇1𝑘+1∶𝑛,2∶𝑛+𝑇2𝑘+1∶𝑛,2∶𝑛 (66) end.
Overall in each step we only need to premultiply two 𝑛×𝑛 Householder matrix and one 𝑛×𝑛 Givens matrix in the procedure of computing the eigenvector.

5. Numerical Experiments

All codes were run in MATLAB 7.04 in double precision and all experiments were performed on a personal computer with 1.41 Ghz central processing unit (AMD Sempron(tm) Processor 2500+), 520 MB memory, and Windows XP system. Machine epsilon 𝑒𝑝𝑠≈2.22×10⁻¹⁶.

Example 16. To illustrate Algorithm 11 more clearly, we give a simple example and show the processes and results of Algorithm 11 of 𝐻₁: 𝐻₁=1−3𝑖−4𝑖3𝑖14𝑖4𝑖4𝑖1. (67)
We first transform 𝐻₁ into symmetry and skew-Hamiltonian matrix 𝑆: 𝐻→𝑆=100034010−304001−4−400−3−410030−4010440001. (68) Then 𝑆→𝐻^𝑇𝑆𝐻=100050010−50−40010400−501005040100−40001→𝐽^𝑇𝐻^𝑇==,𝑆𝐻𝐽150000514000041000000150000514000041𝑇00𝑇 (69) where ,.𝐻=10000000.60.800000.8−0.600000010000000.60.800000.8−0.6𝐽=1000000000−10001000000100010000000001 (70)
Denoting that (𝐷,𝑉) is the eigenpair of 𝑇,,𝑉=0.5522−0.62470.5522−0.70710.00000.70710.44170.78090.4417𝐷=−5.40310001.00000007.4031 (71) then 𝑄𝑉0=𝐻₁𝐽𝑉0=0.5522−0.62470.55220.35340.62470.3534−0.2650−0.4685−0.2650000−0.42430.00000.4243−0.56570.00000.5657 (72) is the first half of eigenvectors of 𝑆. The computing eigenproblem of 𝐻₁ is showed in Table 1.
According to Tables 1 and 2, it can be seen that the eigenvectors are different except the one corresponding to 1. In fact, by Table 1, the eigenvector corresponding to 7.4031 is 𝑣⁽¹⁾=(0.5522,0.3534+0.4243𝑖,−0.2650+0.5657𝑖)^𝑇, and 𝑣⁽²⁾=(−0.2343−0.5000𝑖,0.2343−0.5000𝑖,0.6247)^𝑇 by Table 2. It is easy to prove that rank(𝑣⁽¹⁾,𝑣⁽²⁾)=1, that is, 𝑣⁽¹⁾=𝑒^𝑖𝜃𝑣⁽²⁾. So the complex eigenvector after normalizing is not unique, which is the biggest difference from real eigenvector.

Table 1 

Example 16—real fast algorithm for 𝐻₁ (Algorithm 11).

Table 2 

Example 16—complex algorithm (eig) for eigenproblem 𝐻₁.

Example 17. Now we give a little bigger matrix 𝐻₂; we depict its eigenproblem in Tables 3 and 4, respectively: 𝐻₂=12+𝑖3+𝑖4+𝑖5+𝑖2−𝑖23−1𝑖463−𝑖3+1𝑖3154−𝑖4131+𝑖5−𝑖651−𝑖3. (73)
According to tables, we can see that the eigenvectors are complex and it is easy to prove that the two eigenvectors corresponding to the same eigenvalue are linear dependent, so the eigenvectors, which are computed by the real fast algorithm (Algorithm 11) span the same space as complex algorithm (eig).

Table 3 

Example 17—the real fast algorithm for eigenproblem of 𝐻₂.

Table 4 

Example 17—the real fast algorithm for eigenproblem of 𝐻₂.

Example 18. In this example we construct complex Hermitian matrices of the form 𝐻₃()+1=diag1,2,…,2∗𝑛0−2−1210−2−121⋱⋱⋱⋱−2−120_2𝑛×2𝑛+𝑖∗,0𝐵−𝐵0 (74) where 𝐵=(1/8)diag(eye(𝑛)). In tables we only give a few eigenvalues with largest modulus and compare it with complex algorithm (eig). Let 𝑛=10, 100; that is, the order of matrix 𝐻₃ is 20, 200, respectively. The eigenvalues and errors, which are computed respectively by real fast algorithm and complex algorithm (eig), are showed in Tables 5 and 6. All the eigenvalues of matrix of order 20 are given, while only the 20th largest modulus eigenvalues of matrix order 200 are present. It is easy to prove that eigenvectors corresponding to the same eigenvalue are linear dependent. According to the tables, we have the same result as the complex algorithm (eig), but only the real operation is used.

Table 5 

Example 17—the real fast algorithm for eigenproblem of 𝐻₂.

Table 6 

Example 17—the real fast algorithm for eigenproblem of 𝐻₂.

Example 19. Complex Hermitian matrices are constructed by Hilbert matrix, and the MATLAB code is as follows: (𝑛);𝐻𝑛=10;𝐻=hilb₄=𝐻+tril(𝐻)−diagdiag(𝐻)^′𝑖−tril(𝐻)−diagdiag(𝐻). (75)(76)(77) When 𝑛=10, the 6 largest modulus eigenvalues are present and 3 eigenvectors corresponding to the 3 largest modulus eigenvalues are given here: 𝜆𝐻_4eig=,𝜆𝐻2.137981825207990.51461933549659−0.49873888113430−0.188658336083330.17747332069964−0.08500068573387_4Algorithm11=,𝑉2.137981825207990.51461933549659−0.49873888113430−0.188658336083330.17747332069964−0.08500068573387_(eig)=,𝑄=.−0.1267−0.0857𝑖−0.0156−0.1526𝑖0.1174+0.0100𝑖0.0529+0.2493𝑖−0.0069+0.2192𝑖−0.1508−0.1187𝑖0.0647−0.2646𝑖−0.2177−0.2291𝑖0.1777+0.1521𝑖−0.3229+0.1315𝑖0.3057−0.0144𝑖−0.0986−0.2811𝑖0.2230+0.2709𝑖−0.1554+0.3121𝑖−0.0529+0.2981𝑖0.2089−0.2644𝑖−0.2960−0.2141𝑖0.2578−0.2379𝑖−0.2801−0.2083𝑖0.2142−0.2740𝑖−0.3615−0.0351𝑖−0.2635+0.2411𝑖0.2822+0.1932𝑖0.2113+0.3073𝑖0.1389+0.3183𝑖−0.1273+0.3242𝑖0.1640−0.3632𝑖0.3304−0.3474−0.4016−0.6616−0.4938−0.3499−0.3850−0.2226𝑖0.1121−0.3048𝑖−0.0667+0.5036𝑖−0.2449−0.2266𝑖0.2977+0.0431𝑖0.2830+0.3297𝑖−0.1682−0.2102𝑖0.1960+0.2453𝑖0.3451+0.0933𝑖−0.1203−0.1921𝑖0.0573+0.3070𝑖0.2918−0.0624𝑖−0.0880−0.1754𝑖−0.0550+0.2958𝑖0.2086−0.1463𝑖−0.0649−0.1608𝑖−0.1329+0.2524𝑖0.1273−0.1825𝑖−0.0477−0.1480𝑖−0.1816+0.1980𝑖0.0581−0.1897𝑖−0.0346−0.1368𝑖−0.2082+0.1428𝑖0.0028−0.1802𝑖−0.0243−0.1270𝑖−0.2190+0.0916𝑖−0.0394−0.1614𝑖 (78)(79)(80)
We prove that rank([𝑉(∶,𝑘)𝑄(∶,𝑘)])=1 for 𝑘=1,2,3. It is very interesting that the eigenvectors computed from our algorithm are unique but not the complex algorithm (eig).

Example 20. At last we construct some rand Hermitian matrices of order 𝑁; its MATLAB code is as follows. 𝐴=rand(𝑛),𝑆𝐴=𝐴^′(𝑛)+𝐴,𝐵=rand,𝑆𝐵=𝐵−𝐵^′,𝐻=𝑆𝐴+𝑖^∗𝑆𝐵. (81)(82)(83)
Letting 2∗𝑁=200∶100∶1000, denote max(relerr)=max_𝑗𝜆^eig_𝑗−𝜆^Sympl_𝑗𝜆^Sympl_𝑗(),minrelerr=min_𝑗𝜆^eig_𝑗−𝜆^Sympl_𝑗𝜆^Sympl_𝑗. (84)(85)
All the relative errors we compute are present by a factor of order 𝒪(10^𝑑). The average relative errors for all examples of each dimension we computed for each example in Table 7 are denoted by average (relerr). The max, min, and average errors of 𝑟=𝐻𝑣_𝑖−𝑣_𝑖𝜆 are also considered in Table 8. According to two tables, all the “relerr” and “residual” are order 𝒪(10⁻¹³) and very small. Hence it can be concluded that eigenproblems computed by the real fast algorithm have high precision, and the method always converges.

Table 7 

Example 20—relative error.

Table 8 

Example 20—residual 𝑟=𝐻𝑣_𝑖−𝑣_𝑖𝜆.

Acknowledgments

This work is supported by National Natural Science Foundation of China nos. 11101087 and 10261012 and supported by Research Start Foundation of Fuzhou University nos. 022427 and 600629.