# BiCGSTAB convergence

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.

    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:

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

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.

• 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 Jul 28, 2022 at 2:35
• @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. Jul 28, 2022 at 4:03
• 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. Jul 28, 2022 at 23:25