BiCGSTAB convergence

Question

So I need a fast converging solver for SysLinEq as a subroutine in fortran, decided to test BiCGStab in Matlab.

Thank God I decided to test it out on first before implementing in Fortran as a subroutine.

might be inappropriate place to ask but when I run the following https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method

method from wiki for matrices it somewhat works for small systems,

but for the matrices bigger than 40x40 it does not want to converge.

might be coding/analytical issue. At first I tried regular CG, but my matrices aren't anything around symmetric positive definite, but sparse.

Matrices and vectors are generated at random for testing purpose.

If this is inappropriate place to ask please recommend the best forum - I asked on overflow, last week, no updates.
thank you in advance

    clc; clear; close all;
n = 20;
epsilon = 1e-4
c = 0;

A = rand(n);


b = rand(n);
b = b(1,:)';
fprintf('max err using linsolve is:\n %.22f',max(abs(A*linsolve(A,b)-b)))
fprintf('\n so the solution exists');

x=zeros(1,n)';

x0 = zeros(1,n)';
r0 = b - A*x0;

rhat = r0-100;

rho0 = 1;
alpha = 1;
w0 = 1;

v0 = zeros(1,n)';
p0 = zeros(1,n)';

rim1 = r0;
rhoim1=rho0;
wim1 = w0;
pim1 = p0;
vim1 = v0;
xim1 = x0;
ctr = 0;
while max(abs(A*x-b)>epsilon)

    rhoi = dot(rhat , rim1);
    beta = (rhoi/rhoim1)*(alpha/wim1);
    pi = rim1 + beta*(pim1 - wim1*vim1);
    vi = A*pi;
    alpha = rhoi / dot(rhat,vi);
    h = xim1 + alpha*pi;

    if (max(A*h-b)<epsilon)
       x = h;
       break; 
    end

    s = rim1 - alpha*vi;
    t = A*s;
    wi = dot(t,s)/dot(t,t);
    xi = h+wi*s;

    if (max(A*xi-b)<epsilon)
       x = xi;
       break; 
    end
    ri = s - wi*t;

    rim1 = ri;
    rhoim1=rhoi;
    wim1 = wi;
    pim1 = pi;
    vim1 = vi;
    xim1 = xi;
    ctr = ctr+1;
end
max(abs(A*x - b))
ctr

Matrix sparsity:

Answer 1

First of all, size 40 is pretty much microscopic for the sort of purposes that BiCGstab was invented. I assure you, this method is great for matrices of sizes in the millions and beyond.

Then: even the best iterative method greatly benefits from a preconditioner.

But your real problem is: "Matrices and vectors are generated at random for testing purpose." Nope. Nope. The whole of iterative method theory comes from operators from elliptic PDEs. And the further you get from symmetric and positive definite, the more trouble they have converging. So: compute A*e (where e is the all-ones vector) and add the result to your diagonal. Then try again.