load theta2
N=16;N2=N^2;b=zeros(2*N2,1);x0=b;b(1)=1;
A1=Dirac_r(-0.1,1,N,N2,theta2);
A2=Dirac_r(-0.05,1,N,N2,theta2);
A3=Dirac_r(0,1,N,N2,theta2);
e3=eig(A3);
plot(e3,'o')
tol=1e-13;nmax=2*N2;
[x_gmres1,r_gmres1,H1] = GMRES(A1,b,x0,tol,nmax);
[x_gmres2,r_gmres2,H2] = GMRES(A2,b,x0,tol,nmax);
[x_gmres3,r_gmres3,H3] = GMRES(A3,b,x0,tol,nmax);
[x_cgne1,r_cgne1] = CGNE(A1,b,x0,tol,nmax);
[x_cgne2,r_cgne2] = CGNE(A2,b,x0,tol,nmax);
[x_cgne3,r_cgne3] = CGNE(A3,b,x0,tol,nmax);
k_gmres1=max(size(r_gmres1));
k_gmres2=max(size(r_gmres2));
k_gmres3=max(size(r_gmres3));
k_cgne1=2*max(size(r_cgne1));
k_cgne2=2*max(size(r_cgne2));
k_cgne3=2*max(size(r_cgne3));
semilogy(1:k_gmres1,r_gmres1,';gmres1;')
hold
semilogy(1:k_gmres2,r_gmres2,';gmres2;')
semilogy(1:k_gmres3,r_gmres3,';gmres3;')
semilogy(1:2:k_cgne1,r_cgne1,';cgne1;')
semilogy(1:2:k_cgne2,r_cgne2,';cgne2;')
semilogy(1:2:k_cgne3,r_cgne3,';cgne3;')
xlabel('# matrix-vector')
ylabel('residual norm')
replot
gset terminal postscript
gset out 'conv_hist.ps'
replot
gset terminal x11
hold
x0=rand(2*N2,1);
[x_gmres1,r_gmres1,H1] = GMRES(A1,b,x0,tol,nmax);
[x_cgne1,r_cgne1] = CGNE(A1,b,x0,tol,nmax);
[x_bicg1,r_bicg1] = BiCGg5(A1,b,x0,tol,nmax);
[x_bicgstab1,r_bicgstab1]=BiCGstab(A1,b,x0,tol,nmax);
k_gmres1=max(size(r_gmres1));
k_cgne1=2*max(size(r_cgne1));
k_bicg1=max(size(r_bicg1));
k_bicgstab1=2*max(size(r_bicgstab1));
semilogy(1:k_gmres1,r_gmres1,';gmres1;')
hold
semilogy(1:2:k_cgne1,r_cgne1,';cgne1;')
semilogy(1:k_bicg1,r_bicg1,';bicg1;')
semilogy(1:2:k_bicgstab1,r_bicgstab1,';bicgstab1;')
xlabel('# matrix-vector')
ylabel('residual norm')
replot
gset terminal postscript
gset out 'all_solver_conv_hist.ps'
replot
gset terminal x11
hold
U1=cos(theta2)+sqrt(-1)*sin(theta2);
N=16;N2=N^2;
A=Dirac_KS(0,N,N2,U1);
norm(A+A')
ea=eig(i*A);
plot(ea,'o;ea;')
b=zeros(N2,1);x0=b;b(1)=1;
tol=1e-13;nmax=N2;
[x_gmres,r_gmres,H] = GMRES(A,b,x0,tol,nmax);
[x_cgne,r_cgne] = CGNE(A,b,x0,tol,nmax);
k_cgne=max(size(r_cgne));
semilogy(r_gmres,';gmres;')
hold
semilogy(1:2:2*k_cgne,r_cgne,';cgne;')
xlabel('# matrix-vector')
ylabel('residual norm')
replot
gset terminal postscript
gset out 'conv_hist_ks.ps'
replot
gset terminal x11
hold
b=zeros(2*N2,1);b(1)=1;
[Q,T]=Lanczos(A3'*A3,A3'*b,k_cgne3);
x_lanczos3=Q(:,1:k_cgne3)*(T\[norm(A3'*b);zeros(k_cgne3-1,1)]);
norm(x_cgne3-x_lanczos3)