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');


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;
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;

    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;
    ri = s - wi*t;

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

Matrix sparsity: Matrix sparsity

  • 2
    $\begingroup$ Clearly your matrix is not random, or it would not have this kind of sparsity pattern. $\endgroup$ Jul 28, 2022 at 3:35

1 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.

  • $\begingroup$ the article I tried to use the method from said must end up with SPD matrix, but as you can see the sparsity, I either mess up something or article lied :( Thanks for suggest! but can't really be confident ending up with SPD matrix system $\endgroup$
    – 2Napasa
    Jul 28, 2022 at 2:35
  • 1
    $\begingroup$ @2Napasa You don't need exactly SPD. BiCGstab can handle non-symmetric matrices, and they can be indefinite, but they shouldn't be too far from SPD. Unfortunately there is little theory that says how far. But random matrices definitely don't work. $\endgroup$ Jul 28, 2022 at 4:03
  • $\begingroup$ thanks, I recently managed to make matrix SPD and used regular CG method work for the system, but convergence time is a disaster for 20x20 grid took 64 seconds over 500 different systems with tolerance around $1E^{-6}$. left 50x50 grid to run over 100it, will wait and see. $\endgroup$
    – 2Napasa
    Jul 28, 2022 at 23:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.