# A FEAST Algorithm for the Linear Response Eigenvalue Problem
## 4 Numerical Examples  
In this section, we present some numerical examples to illustrate the effectiveness of Algorithm 1 and the upper bounds in Theorem 1. **Example 1**. : _We first examine the upper bounds of Theorem 1. For simplicity, we consider diagonal matrices for \(K\) and \(M\). Take \(M=K=\mathrm{diag}(\lambda_{1},\dots,\lambda_{N})\) with \(N=100\) where:_  
\[\lambda_{1}=1+\eta,\;\lambda_{2}=1,\;\lambda_{3}=1-\eta,\quad\lambda_{j}= \frac{N+4-j}{2N},\;\text{for}\;j=4,\dots,N.\]  
_In such a case, there are two eigenvalue clusters: \(\{\pm\lambda_{1},\pm\lambda_{2},\pm\lambda_{3}\}\) and \(\{\pm\lambda_{4},\dots,\pm\lambda_{N}\}\), and \(\Psi=K^{-1/2}\).__We seek the approximations associated with the first cluster \(\{\pm\lambda_{1},\pm\lambda_{2},\pm\lambda_{3}\}\) and vary the parameter \(\eta>0\) to control the tightness among eigenvalues. To make the numerical example repeatable, the initial block \(Y_{{}_{(0)}}\) is chosen to be:_  
\[Y_{{}_{(0)}}= \begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&1\\ \frac{1}{N}&\sin 1&\cos 1\\ \vdots&\vdots&\vdots\\ \frac{N-3}{N}&\sin(N-3)&\cos(N-3)\end{bmatrix} .\]  
_In such a way, \(Y_{{}_{(0)}}\) satisfies the condition that \(Y_{{}_{(0)}}^{\mathrm{T}}M\Psi_{{}_{(1;3)}}\) is nonsingular. We run Algorithm 1 and check the bound for \(\|\tan\Theta_{{}_{\mathrm{M}}}(\Psi_{{}_{(1;13)}},V_{{}_{(0)}})\|_{\mathrm{F}}\) given by (11) with \(k=1\) and \(k=2\), respectively. In Algorithm 1, the trapezoidal rule is applied to calculate \(\tilde{f}(KM)Y_{{}_{(0)}}\) with the parameters \(c=1\) and \(r=0.2\). We compute the following factors:_  
\[\epsilon_{1} =\|\tan\Theta_{{}_{\mathrm{M}}}(\Psi_{{}_{(1;13)}},V_{{}_{(0)}}) \|_{\mathrm{F}}, \tag{21a}\] \[\epsilon_{2} =\left(\frac{|\tilde{f}(\lambda_{4}^{2})|}{|\tilde{f}(\lambda_{3 }^{2})|}\right)^{k}\times\|\tan\Theta_{{}_{\mathrm{M}}}(\Psi_{{}_{(c,1;3)}}, \widehat{Y})\|_{\mathrm{F}}, \tag{21b}\]  
_in Table 1. From Table 1, we can see that, as \(\eta\) goes to zero, the bounds \(\epsilon_{2}\) for \(k=1,2\) are sharp and insensitive to \(\eta\)._  
**Example 2**. : _To test the effectiveness of Algorithm 1, we chose three test problems, i.e., test 1 to test 3 in Table 2, which come from the linear response analysis for Na\({}_{2}\), Na\({}_{4}\), and silane (SiH4) compounds, respectively [9]. The matrices \(K\) and \(M\) of test 1, test 2, and test 3 are both symmetric positive definite with order \(N=1862\), 2834, and 5660, respectively. We compute the eigenvalues whose square lies in the interval \((\underline{\omega}_{l},\,\overline{\omega}_{l})\) for \(t=1,2,3\) and the associated eigenvectors by using the MATLAB parallel pool. The interval \((\underline{\omega}_{l},\,\overline{\omega}_{l})\) and their corresponding eigenvalue counts \(s_{l}\) for \(t=1,2,3\) are detailed in Table 2. All our experiments were performed on a Windows 10 (64 bit) PC-Intel(R) Core(TM) i7-6700 CPU 3.40 GHz, 16 GB of RAM using MATLAB Version 8.5 (R2015a) with machine epsilon \(2.22\times 10^{-16}\) in double precision floating point arithmetic. In demonstrating the quality of computed approximations, we calculated the normalized residual norms \(r(\rho_{j})\) defined in (10) and relative eigenvalues errors:_  
\[e(\rho_{j})=\frac{|\rho_{j}-\lambda_{j}|}{\lambda_{j}}\]  
_for the jth approximation eigenvalue \((\rho_{j},\mathbf{z}_{j})\) where \(\lambda_{j}\) were computed by MATLAB's function eig on \(H\), and considered to be the \"exact\" eigenvalues for testing purposes._  
_In this example, we first investigated how the maximal relative eigenvalue errors and normalized residual norms, i.e., \(\max e(\rho_{j})\) and \(\max r(\rho_{j})\) in Table 3, respectively, changed for an increase in the number of quadrature points \(q\) in (6). We used the direct methods, such as Gaussian elimination, to solve the shifted linear systems (7)in this example.