% Routine for complex bi-conjugate gradient method for
% solving linear systems.  Can be used for complex quantities,
% but real variables are used here.  Based on discussion by
% Jacobs in "Sparse matrices and their uses", ed. I.S. Duff.
% designed for use with square (nl x nl) systems.
% This routine works; tested on sample problems.
%
% Syntax: x = cbcg(a,b)
%
% by D. Holmgren, Mar 21 94. holmgren@phobos.astro.uwo.ca
function x = cbcg(a,b)
tol = 1.e-03;
ni = input('Number of iterations: ');
nl = length(b);
% this seems to be a good starting solution...
x0 = mean(b) .* ones(nl,1);
x = x0;
% set initial residual...
r = b - a * x;
% set initial bi-residual...
rb = r;
% set first direction vector...
p = r;
% and first bi-direction vector...
pb = rb;
for i = 1:ni
% some workspaces...
  rb1 = rb;
  r1 = r;
% one of three forms of the coefficient alpha...
  alpha = (rb' * r) / (pb' * (a * p));
% new estimate of x...
  x = x + alpha .* p;
% new residual...
  r = b - a * x;
% new bi-residual...
  rb = rb - alpha .* (a' * pb);
% one of three forms of coefficient beta (note sign)...
  beta = (rb' * r) / (rb1' * r1);
%  beta = -(rb' * (a * p)) / (pb' * (a * p));
%  beta = -((a' * pb)' * r) / (pb' * (a * p));
% set new direction vector...
  p = r + beta .* p;
% set new bi-direction vector...
  pb = rb + beta .* pb;
  sumsq = norm(r)^2;
  disp(' Sum of squares of residuals:')
  disp(sumsq)
% jump out if sum of squares is small...
  if sumsq < tol, break, end
end