% Initialize b, x, jac1, jac2, iter to demonstrate testfmg
% Written by J. Demmel, UC Berkeley, April 22, 1995
k=7;
kk=2^k+1;
[x,y]=meshdom(1:kk,1:kk); %
% b=sin((0:kk-1)*3*pi/kk)'; 
t=2*eye(kk-2)-diag(ones(kk-3,1),1)-diag(ones(kk-3,1),-1);
[v,d]=eig(t); %
[ds,is]=sort(diag(d)); %
vs=v(:,is); %
% b=vs*(1.02.^(0:-1:-kk+1)'.*randn(kk,1)); %
sol=zeros(kk-2,1);
sol((kk-1)/8:3*(kk-1)/8) = (1:(kk-1)/4 + 1);
sol(5*(kk-1)/8:7*(kk-1)/8) = -2*((kk-1)/4 + 1:-1:1);
xsol=[0;sol;0];
b=[0;t*sol;0];
nb = norm(b,'inf');
b=b/nb;
sol = sol/nb;
xsol = xsol/nb;
[nb,mb]=size(b);
hold off; clg
% x=b;b=zeros(size(b));
figure(1)
hold off, clg
subplot(121)
% axes('position',[.1,.1,.35,.5])
plot(b,'w')
axis([1,nb,-1,1])
title('right hand side b')
subplot(122)
% axes('position',[.55,.1,.35,.5])
plot(xsol,'w');
axis([1,nb,-1,1])
title('solution x')
figure(2)
hold off, clg
plot(vs'*t*sol,'w');
axis([1,nb,-1,1])
title('frequency components of b')
x=zeros(size(b));
jac1=1;jac2=1;iter=8;
