clear;
%  syms sx sy sa sb;
%  su=sin(2*pi*sx)*sin(2*pi*sy);
%  sf = - diff(su,sx,2) - diff(su,sy,2) + sa*diff(su,sx) + sb*diff(su,sy);
%  ufunc=inline(char(vectorize(subs(su,{sx,sy,sa,sb},{'x','y','a','b'}))),'x','y','a','b');
%  ffunc=inline(char(vectorize(subs(sf,{sx,sy,sa,sb},{'x','y','a','b'}))),'x','y','a','b');
%  
n=100;
nx=n;ny=n;nz=1;
%  
%  hx=1/(nx+1);
%  hy=1/(ny+1);
%  if(nz>1)
%    dim=3
%    hz=1/(nz+1);
%  else
%    dim=2;
%    hz=0;
%  end
%  
%  h=(hx^2 + hy^2 + hz^2)/dim;

a=.1;
b=.1;
A=genA(nx,ny,nz,a,b,0,0);

%  xx=linspace(0,1,n+2);
%  [X,Y]=meshgrid(xx(2:end-1));
%  Uex = ufunc(X,Y,a,b);
%  F = h*ffunc(X,Y,a,b);
%  f=reshape(F,n*n,1);
%  
%  u=A\f;
%  
%  U=reshape(u,n,n);
%  
%  [Ux,Uy]=gradient(U,hx,hy);
%  E=Uex-U;
%  
%  norm(E);
%  figure;
%  subplot(1,2,1),
%  surf(X,Y,U);
%  subplot(1,2,2),
%  surf(X,Y,Uex);
%  
%  figure;
%  subplot(1,2,1),
%  contour(X,Y,U);hold on;quiver(X,Y,Ux,Uy);
%  subplot(1,2,2),
%  surf(X,Y,E);

condA=condest(A);
tic;
[uA,flagA,relresA,iterA,resvecA]=bicgstab(A,ones(n*n,1),1e-8,n*n,[],[],rand(n*n,1));
tA=toc;
disp(sprintf(' A: cond = %g, iter=%d, nnz(A)=%d, time=%g',condA,ceil(iterA),nnz(A),tA));

for myeps=[.6 .5 .4 .3 .2 .12];
  tic;
  M=spai(A,myeps,10,10,1,100000,0);
  myp = @(x,M) M*x;
  [uB,flagB,relresB,iterB,resvecB]=bicgstab(A,ones(n*n,1),1e-8,n*n,@(x,M)myp(x,M),[],rand(n*n,1),M);
  tB=toc;
  B=M*A;
  condB=condest(B);
  disp(sprintf('MA: cond = %g, iter=%d, nnz(MA)=%d, time=%g, flag?  %d',condB,ceil(iterB),nnz(B),tB,flagB));
end

figure;
subplot(1,2,1),
spy(A);
subplot(1,2,2),
spy(B);