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.

Q Code (Complex)

Clear[FeastQComplex]
FeastQComplex[{A_,Y_},{Lmin_,Lmax_}] := Module[
{Ne=8,r,c,B,n=Length[A],XW, Q,θe,Ze,Qe,Qhe},
{r,c}=0.5{Lmax-Lmin,Lmax+Lmin};
B=SparseArray[Band[{1,1}]->1.0,{n,n}];
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}}));

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];
Qhe= LinearSolve[Transpose[Ze*B-A],Y];
Q=Q-(XW[[e,2]]/4)*r*(((E^(I*θe))*Qe)-((E^(-I*θe))*Qhe));
];
Q
]

Feast Complex Code

FEASTComplex[A_,{Lmin_,Lmax_},M_]:=Module[
{n=Length[A],M0=Ceiling[1.5*M],Y,evals,evecs},
Y=RandomReal[{-1,1},{n,M0}];
Q=FeastQComplex[{A,Y},{Lmin,Lmax}];
Aq=Transpose[Q].A.Q;
Bq=Transpose[Q].Q;
{evals,evecs}=Eigensystem[{Aq,Bq}];
{evals,Q.Transpose[evecs]}
]

Test:

A=RandomComplex[{-1,1},{4,4}];
Eigenvalues[A]
Eigenvectors[A]

{Lmin,Lmax}={-1,0};
M=1;
{evals,evecs}=FEASTComplex[A,{Lmin,Lmax},M]