Full Matrices

dim=4000;
A= RandomReal[{0,1},{dim,dim}];
A= A + Transpose[A] + DiagonalMatrix[RandomReal[dim{0.5,1.5}, dim]];
MatrixForm[A];
b=RandomReal[{0,1},dim];

Timing[
xSol=LinearSolve[A,b];
]

Timing[
ListPlot[Eigenvalues[A]]
]

Tol = (10^(-6));
x = RandomReal[{0,1},dim];
Timing[
rOld =b- A. x  ;
p = rOld;
MaxN=300;
Do[
Ap= A.p;
rOldSq=rOld.rOld;
α= (rOldSq/(p.Ap));
x = x + α p;
rNew=rOld-α Ap;
If[ Norm[rNew]<Tol, Print[i];Break[]];
β = ((rNew.rNew)/rOldSq);
p=rNew + β p;
rOld=rNew, 
{i,1,MaxN}
]
]
Norm[x-xSol]