function [X,res,R,Rt,P,Pt]=BiCG(A,b,x0,n);
% BICG Bi-Conjugate Gradient Method
%   [X,res,R,Rt,P,Pt]=BiCG(A,b,x0,n); computes an approximate solution
%   of the linear system Ax=b performing n steps of the Bi-Cojugate
%   Gradient algorithm starting with the initial vector x0. The matrix
%   X contains the approximate solutions and the vector res the norms
%   of the corresponding residuals. The matrices R, Rt and P, Pt
%   contain the residuals and directions.

x=x0;r=b-A*x0; rt=r;
p=zeros(size(b)); pt=p;
rho=1; res=norm(r);
X=x0; P=[]; Pt=[]; R=r; Rt=rt; 
for j=1:n
  rhoold=rho; rho=rt'*r;
  beta=rho/rhoold;
  p=r+beta*p; pt=rt+beta*pt;
  P=[P p]; Pt=[Pt pt];                                 % not needed
  Ap=A*p;
  alpha=rho/(pt'*Ap);
  r=r-alpha*Ap; rt=rt-alpha*A'*pt;
  R=[R r]; Rt=[Rt rt];                                 % not needed
  x=x+alpha*p;
  X=[X x];                                             % not needed
  res=[res norm(r)]; 
end;