Clear[n,x]
n=8;
f[x_]=LegendreP[n,x];
MyPrecision=32;
Sub=NSolve[f[x]==0,x,
WorkingPrecision->MyPrecision];
thing=Transpose[{x,(2/((1-(x^2)) ((f'[x])^2)))}/.Sub]

Note (Subscript[x, 1])parameter in FEAST paper is wrong in the last digot.

Derivative[0,1][LegendreP][n,x]

D[f[n,x],x]//FullForm

D[LegendreP[n,x],x]

symboliclegendre[n_,x_]:=Solve[LegendreP[n,x]==0];

LegendreP

legendreprime[n_,a_]:=D[LegendreP[n,x],x]/.x->a;
weights[n_,x_]:=2/((1-x^2) legendreprime[n,x]^2);

(*how many terms should be generated*)
h=64;

(*what numerical precision is desired?*)
precision=256;

str=OpenWrite["lgvalues.txt"];
Write[str,"abcissae"];
Do[Write[str];Write[str,n];Write[str];
nlist=symboliclegendre[n,x];
xnlist=x/.nlist;
Do[Write[str,N[Part[xnlist,i],precision]],{i,Length[xnlist]}];,{n,2,h}];
Write[str,"weights"];
Do[Write[str];Write[str,n];Write[str];
slist:=symboliclegendre[n,x];
xslist=x/.slist;
Do[Write[str,N[weights[n,Part[xslist,i]],precision]],{i,Length[xslist]}];,{n,2,h}];
Close[str];

Implementing the FEAST Algorithm

Density matrix based algorithm for solving eigenvalue problems. Phys Rev B v79, 115112, 2009, Eric Polizzi.

FEAST[A_,{λMin_,λMax_}]:= Module[
{c=(λMin+λMax)/2.0, r=(λMax-λMin)/2.0,Q},
Table[
z=c+r (E^(i θ))

Table[
precision=8;
{absc,weights,errweights}=NIntegrate`GaussBerntsenEspelidRuleData[n,precision];
{n,Sort[2*absc-1],Length[absc]},
{n,3,7,1}]

precision=8;
n=3
{absc,weights,errweights}=NIntegrate`GaussRuleData[2n+1,precision];
ListPlot[absc]
ListPlot[2*absc -1]

Gauss

RandomSymDefThm144[n_]:= Module[{M},
M= RandomReal[{-1,1},{n,n}];
M=Transpose[M].M;
M=0.5*(M+Transpose[M])
]

MyCholesky[A_]:=Module[
{AA,s,R,n=Length[A]},
(* This is the 2nd outerprodduct version from before *)
AA=A;R=0*A;
Do[
R[[i,i]]=(Sqrt[AA[[i,i]]]);
s=AA[[i,i+1;;-1]]/R[[i,i]];
R[[i,i+1;;-1]]=s;
AA[[i+1;;-1,i+1;;-1]]-=Outer[Times, s,s],
{i,1,n}];
R
]

A little bit of testing and timing.

n=12;
A=RandomSymDefThm144[n];
MyR=MyCholesky[A];
R=CholeskyDecomposition[A];
Norm[R-MyR]

MyTimes=Table[
n=Floor[ 12 (1.5^p)];
A=RandomSymDefThm144[n];
{n, Timing[MyCholesky[A];][[1]]},
{p, 3, 10}];
BuiltTimes=Table[
n=Floor[ 12 (1.5^p)];
A=RandomSymDefThm144[n];
{n, Timing[CholeskyDecomposition[A];][[1]]},
{p, 5, 16}];

CubeFn[n_]:=(0.001/(50^3))(n^3);
ListLogLogPlot[{MyTimes,BuiltTimes, {{50,CubeFn[50]},{5000,CubeFn[5000]}}},
PlotStyle->{Red, Blue,Green},
Mesh->All,
Joined->True,
Mesh->All,
GridLines->Automatic,
Frame->True]

A= Import["http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/bcsstruc1/bcsstk04.mtx.gz"]
Timing[MyR=MyCholesky[A];]
Timing[R=CholeskyDecomposition[A];]
{MatrixPlot[A],MatrixPlot[MyR]}

A= Import["http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/bcsstruc2/bcsstk14.mtx.gz"]
Timing[MyR=MyCholesky[A];]
Timing[R=CholeskyDecomposition[A];]
{MatrixPlot[A],MatrixPlot[MyR]}