clc
clear all
meth=1; %1 for iterative method  2 for perconditioning
alphm=1%1 for old alpha    2  for   new alpha
rowtable=4;obtainrho=1;e=1e-3;
fil=fopen('optdat.txt','w');
filtable=fopen('tableinex.txt','w');

fallper={'H','T','B',' '};
sallper={'SS','SSD','DDBD'};
Diag={'AH','DAH','BDAH','A'};
Lower=   {'Z','LAH','BLAH','Z'};
Qi={'In','DAH','BDAH'};

c1={'Iterative method','preconditioning'};
c2={'oldalpha','newalpha'};
vari={'size','precon','alpha' ,'rho','CPU','IT'};
for alphm=1:2
for  meth=1
fprintf(filtable,'\\begin{table}\n');
fprintf(filtable,'\\begin{center}\n');
fprintf(filtable,'\\caption{IT,CPU,$\\rho');
fprintf(filtable,'(\\alpha)$ and $\\alpha$');
fprintf(filtable,'for the inexact EPSS ');
cap=sprintf('%s for %s.}',c1{meth},c2{alphm});
fprintf(filtable,cap);
fprintf(filtable,'\n \\scriptsize\n');
fprintf(filtable,'\\label{inex}\n');
fprintf(filtable,'\\scalebox{.84}\n');
fprintf(filtable,'{\n');
i_n=0;
for in=3:4
i_n=i_n+1;
n=2^(in-1)*100;
tolout=1e-7;
tolin=1e-3;rest=20;maxres=50;
[ B,C,D,E ] = baimatrix(n );
sprintf('for n=%g matrices are generated.',n);
A=sparse([B C; D E ]);
AH=.5*(A+A');
DAH=diag(diag(AH));
LAH=2*(tril(AH)-DAH);
BDAH=sparse(n,n);
BDAH(1:.9*n,1:.9*n)=B;
BDAH(.9*n+1:end,.9*n+1:end)=E;
BLAH=sparse(n,n);
BLAH(.9*n+1:end,1:.9*n)=D+C';
In=speye(n,n);Z=sparse(n,n);
tk=0;ic=0;fc=['|c'];f2=['&%s'];
f3=f2;f1=['\\multirow[c]{4}{*}{\\bfseries %g}&%s'];
f4=f2;row=['\\hline \n &'];
for cas=1:12
ic=ic+1;
alcas(ic)=cas;
fc=['|c|c',fc];
if cas<5
f1=[f1,sprintf('&$%%1.4f$')];
else
f1=[f1,sprintf('&$%%1.4f$')];
end;
f2=[f2,'&%1.2f'];
f3=[f3,'&%g'];
f4=[f4,'&%1.4f'];
minro1=1;minro2=1;
minit=10000;mincpu=1000;minalph=-3;
va=1;
for ia=va
tk=tk+1;
tic
iq=floor((cas-1)/4)+1;
ip=mod(cas-1,4)+1;
eval(sprintf('D=Diag{%g};',ip));
eval(sprintf('D=%s;',D));
eval(sprintf('L=Lower{%g};',ip));
eval(sprintf('L=%s;',L));
eval(sprintf('Q=Qi{%g};',iq));
eval(sprintf('Q=%s;',Q));
P=D+L;
if alphm==1
alph=sqrt(eigs(D,Q,1,'lm')*eigs(D,Q,1,'sm'));
elseif alphm==2
alph= (sqrt(eigs(P'*(Q\P),Q,1,'lm')));
end;
Sigma=alph*Q;
time_alpha=toc;time_alpha=0;
S=A-P;
if obtainrho==1
Gfun=@(v)((Sigma+P)\((Sigma-S)*((Sigma+S)\((Sigma-P)*v))));
[rho1,er,itt]=specradi(Gfun,300,e) ;
if er>e
sprintf('no') ;
Gtild=((Sigma+P)\(Sigma-S))*((Sigma+S)\(Sigma-P));
rho1=abs(eigs(Gtild, 1));
end;
else
rho1=2;
end;
if rho1<minro1
minro1=rho1;
end;
u=zeros(n,1);
b=A*ones(n,1);
r=b;k=0;
%LS=ichol(Sigma+P,struct('type','ict','droptol',.01));
if meth==1
tic
while (norm(r)/norm(b))>tolout
k=k+1;
if cas==2||cas==6||cas==10
[z,fl]=cgs(Sigma+P,r,tolin,rest);%,LS,LS');
else
[z,fl]=tgmres(Sigma+P,r,rest,tolin,maxres);
end;
u=u+z;r=b-A*u;
[z,fl]=tgmres(Sigma+S,r,rest,tolin,maxres);
u=u+z;r=b-A*u;
end;
elseif meth==2
tic
if cas==2||cas==6||cas==10
fslo=@(v)invsehss(P,S,Sigma,v,rest,tolin,maxres);
else
fslo=@(v)invsepss(P,S,Sigma,v,rest,tolin,maxres);
end;
[xsol,ff,relres,it]=lgmres(A,b,rest,tolout,maxres,[],fslo);
k=(it(1)-1)*rest+it(2);
end;
time_so=toc+time_alpha;
if k<=minit && mincpu>time_so
minit=k;
minrok1=rho1;
mincpu=time_so;
minalph=alph;
end;
end;
ALLrpre{cas}={sprintf('%s%s',fallper{ip}, sallper{iq})};
da1=sprintf('%14s%s%12.5f%11.5f%9.3f%7.0f\n',...
fallper{ip}, sallper{iq},minalph,minrok1,mincpu,minit);
if tk==1
Tl=table(n,ALLrpre{cas},minalph,...
minrok1,mincpu,minit,'VariableNames',vari)
dat=sprintf('\n size  Pricon  Optalpha   Rho1    Cpu    Iter   \n');
fprintf(fil,dat);
dat=sprintf('____   _____  ______  ______  ______  _____   \n');
fprintf(fil,dat);
else
fprintf(da1);
end;
fprintf(fil,sprintf('%g%8s%s%11.3f%11.3f%11.2f%6.0f\n',...
n,fallper{ip}, sallper{iq},minalph,minrok1,mincpu ,minit));
datatable(:,ic)=[minalph,minrok1, mincpu ,minit]';
if i_n==1 && mod(ic,4)==1
row=sprintf('%s&\\\\multicolumn{4}{c|}',row);
qh=Qi{1+floor((cas-1)/4)};
row=sprintf('%s{$\\\\Sigma=\\\\alpha %s$}',row,qh);
end;
end;
f1=sprintf('%s\\\\\\\\ \\\\cline{2-%g}\\n',f1,ic+2);
f2=sprintf('%s\\\\\\\\ \\\\cline{2-%g}\\n',f2,ic+2);
f4=sprintf('%s\\\\\\\\ \\\\cline{2-%g}\\n',f4,ic+2);
f3=[f3,' \\\\ \\hline \n'];
row=[row,'\\\\ \n'];
if i_n==1
fprintf(filtable,sprintf('\\\\begin{tabu}{%s|}\\n',fc));
fprintf(filtable,row);
fprintf(filtable, '\\hline\n$n$&$\\mathbf{ P}$');
for i_all=1:length(alcas)
qh=1+mod((i_all-1),4);
fprintf(filtable,sprintf('&%s',sprintf('${\\\\mathbf P}_%g$',qh)));
end;
fprintf(filtable,'\\\\ \n \\hline \n');
end;
for i_t=1:rowtable
if i_t==1
fprintf(filtable,f1,n,vari{i_t+2},datatable(i_t,:));
elseif i_t==2
fprintf(filtable,f4,vari{i_t+2},datatable(i_t,:));
elseif i_t==3
fprintf(filtable,f2,vari{i_t+2},datatable(i_t,:));
else
fprintf(filtable,f3,vari{i_t+2},datatable(i_t,:));
end;
end;
end;
fprintf(filtable,'\\end{tabu}\n');
fprintf(filtable,'}\n');
fprintf(filtable,'\\end{center}\n');
fprintf(filtable,'\\end{table}\n');
end;
end;
fclose(fil);fclose(filtable);
clc;
type('optdat.txt')
type('tableinex.txt')

function x= invsehss(P,S,Sigma,v0,rest,tolin,maxres)
[z1,fl]=cgs(Sigma+P,v0,tolin,rest);%,LS,LS');
[x,fl]=tgmres(Sigma+S,Sigma*z1,rest,tolin,maxres);
end
function x= invsepss(P,S,Sigma,v0,rest,tolin,maxres)
[z1,fl]=tgmres(Sigma+P,v0,rest,tolin,maxres);
[x,fl]=tgmres(Sigma+S,Sigma*z1,rest,tolin,maxres);
end

function [ B,C,D,E ] = baimatrix(n )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
q=.9*n;
W=sparse(q,q);
W(q,q)=q+1;
for iw=1:q-1
W(iw,iw)=iw+1;
W(iw,iw+1)=1;W(iw+1,iw)=1;
end;
N=sparse(n-q,n-q);
N(n-q,n-q)=n-q+1;
for in=1:n-q-1
N(in,in)=in+1;
N(in,in+1)=1; N(in+1,in)=1;
end;
F=sparse(q,n-q);
for j_f=1:n-q
F(j_f+2*q-n,j_f)=j_f;
end;

w=1./[1:n-q];
G=sparse(diag(w)) ;
B=sparse(W);
C= sparse(F*G);
D=-sparse(F');
E=sparse(N);
end




