Feast

N = matrix size
Lmin = minimum guess for eigenvalue
Lmax = maximum guess for eigenvalue
Search Interval is [Lmin, Lmax]
Ne = number of points for Gauss - Legendre quadrature
M = number of eigenvalues in search interval
M0 = size of working subspace (greater than M)

Polizzi uses M0 ≥ 1.5M and Ne=8.

Notes AAS:
1) I made all the variables internal.
2) I added some comments.
3) I made a thing that made a symmetric matrix with a set of desired eigenvalues.
4) I broke the Feast code up into little pieces.

Making test Matrices

Clear[SymMatWithEvals]
SymMatWithEvals[evals_List]:= Module[
{n=Length[evals],Q,A},
(* Form a random rotation  *)
Q=QRDecomposition[RandomVariate[NormalDistribution[0,1],{n,n}]][[1]];
(* Similarity transforms the input eigenvalues to be unrecognizable *)
A=Q.DiagonalMatrix[evals].Transpose[Q];
(* A is already symmetric in exact arithmetic *)
(* but for some purposes I need to fix any floating point asymmetry *)
0.5 (A + Transpose[A])
]

Testing my code.  Note the evals always come out of the Eigenvalue commands sorted in decreasing magnitude.

A=SymMatWithEvals[{1,50,400}]
Eigenvalues[A]

Look as as though it works.

n=34;
evals=Sort[RandomReal[{-12,12},n],Abs[#1]>Abs[#2]&]
A=SymMatWithEvals[evals];
Norm[Eigenvalues[A]-evals]
Norm[Transpose[A]-A]

Feast Code

Q Code

The first thing I did was spit out the Q code.  It is complicated and messy but it looks better on its own.

Clear[FeastQ]
FeastQ[{A_,Y_},{Lmin_,Lmax_}] := Module[
{Ne=8,r,c,B,n=Length[A],XW, Q,θe,Ze,Qe},
{r,c}=0.5{Lmax-Lmin,Lmax+Lmin};
(* Note doing the special case with B=Id *)
B=SparseArray[Band[{1,1}]->1.0,{n,n}];
(* XW are the Gaussian integration data from Polizzi *)
XW=(({{0.183434642495649, 0.362683783378361}, {-0.183434642495649, 0.362683783378361}, {0.525532409916328, 0.313706645877887}, {-0.525532409916328, 0.313706645877887}, {0.796666477413626, 0.222381034453374}, {-0.796666477413626, 0.222381034453374}, {0.960289856497536, 0.101228536290376}, {-0.960289856497536, 0.101228536290376}}));
(* Initializing then forming Q*)
Q=0*Y;
For[e=1, e≤Ne,e++,
θe=-(π/2)(XW[[e,1]]-1);
Ze=c+r*E^(I*θe);
Qe=LinearSolve[Ze*B-A,Y];
Q=Q-(XW[[e,2]]/2)*Re[r*E^(I*θe)*Qe];
];
Q
]

All I can really test is that it runs and produces an appropriate Q.

{n,M0}={14,6};
evals=Sort[RandomReal[{-12,12},n],Abs[#1]>Abs[#2]&];
A=SymMatWithEvals[evals];
Y=RandomReal[{-1,1},{n,M0}];
{Lmin, Lmax}={0.1, 2.3};
Q=FeastQ[{A,Y},{Lmin, Lmax}];
Dimensions[Q]

Main Feast Code

I tidied it up quite a bit and let it return all the evals it finds.

FEAST[A_,{Lmin_,Lmax_},M_]:=Module[
{n=Length[A],M0=Ceiling[1.5*M],Y,evals,evecs},
(* Initial Random Directions *)
Y=RandomReal[{-1,1},{n,M0}];
(* Initializing then forming Q*)
Q=FeastQ[{A,Y},{Lmin,Lmax}];
(* forming Aq and Bq *)
(* Note doing special case with B=Id *)
Aq=Transpose[Q].A.Q;
Bq=Transpose[Q].Q;
{evals,evecs}=Eigensystem[{Aq,Bq}];
{evals,Q.Transpose[evecs]}
]

Testing.  Note that it always finds the three it should but it sometimes finds ones that do not exist!

evals={1,50,400};
n=Length[evals];
A=SymMatWithEvals[evals];
{Lmin,Lmax}={0.99999999999999999999999999, 400.000000000000000000000000000001};
M=3;
{evals,evecs}=FEAST[A,{Lmin,Lmax},M];
evals

x=RandomReal[{-1,1},{n,1}]
Norm[A.x-evals[[1]]*x,1]/Norm[A.x,1]

evecs

Testing that the algorithm works

evals={1,25,50,400,1000};
n=Length[evals];
A=SymMatWithEvals[evals];
{Lmin,Lmax}={0,1200};
M=5;
{evals,evecs}=FEAST[A,{Lmin,Lmax},M];
evals

{4 9.999999999999964`,25.000000000002128`,24.690183067044295`,1.0000000000000198`,0.`}

The following computes the residual norm between the known eigenvalues and the given eigenvalues

Norm[50-evals[[1]]]/Norm[50]
Norm[25-evals[[2]]]/Norm[25]
Norm[1-evals[[4]]]/Norm[1]

The following computes Polizzi’s norm.  We’re doing this for comparison.

Norm[A.(evecs[[All,1]]/Norm[evecs[[All,1]]])-evals[[1]]*(evecs[[All,1]]/Norm[evecs[[All,1]]]),1]/Norm[A.(evecs[[All,1]]/Norm[evecs[[All,1]]]),1]
Norm[A.(evecs[[All,2]]/Norm[evecs[[All,2]]])-evals[[2]]*(evecs[[All,2]]/Norm[evecs[[All,2]]]),1]/Norm[A.(evecs[[All,2]]/Norm[evecs[[All,2]]]),1]
Norm[A.(evecs[[All,4]]/Norm[evecs[[All,4]]])-evals[[4]]*(evecs[[All,4]]/Norm[evecs[[All,4]]]),1]/Norm[A.(evecs[[All,4]]/Norm[evecs[[All,4]]]),1]

0.011914984061113687`

0.002148514867156731`