qemsvd_mod

function dq=qemsvd_mod(a,b,neig,plot)
 given response matrix a and vector b to be corrected the function finds 
 the vector dq so that b-a*dq=0

 if the input plot is = 1 the plot of correction effect and correctors
 strengths versus number of eigenvectors is also computed.

 svd part is carbon copy of qemsvd by L.Farvacque in qempanel

0001 function dq=qemsvd_mod(a,b,neig,plot)
0002 %function dq=qemsvd_mod(a,b,neig,plot)
0003 % given response matrix a and vector b to be corrected the function finds
0004 % the vector dq so that b-a*dq=0
0005 %
0006 % if the input plot is = 1 the plot of correction effect and correctors
0007 % strengths versus number of eigenvectors is also computed.
0008 %
0009 % svd part is carbon copy of qemsvd by L.Farvacque in qempanel
0010 
0011 [u,s,v]=svd(a,0);
0012 lambda=diag(s);
0013 nmax=length(lambda);
0014 eigsorb=u'*b;
0015 if neig > nmax
0016     neig=nmax;
0017     warning('Svd:maxsize',['number of vectors limited to ' num2str(nmax)]);
0018 end
0019 eigscor=eigsorb(1:neig)./lambda(1:neig);
0020 dq=v(:,1:neig)*eigscor;
0021 
0022 
0023 % plot correction effect to set appropriate number of eigenvectors.
0024 if nargin<4
0025     plot=0;
0026 end
0027 
0028 if plot
0029     numeigen0=neig;
0030     neig=sort([1:5:nmax,numeigen0]);
0031     dqstd=zeros(size(neig));
0032     ostd=zeros(size(neig));
0033     for ineig=1:length(neig)% loop number of eigenvectors used in the correction
0034         eigscor=eigsorb(1:neig(ineig))./lambda(1:neig(ineig));
0035         dqe=v(:,1:neig(ineig))*eigscor;
0036         o=b-a*dqe;
0037         dqstd(ineig)=std(dqe); 
0038         ostd(ineig)=std(o);
0039     end
0040     
0041     f=figure('visible','on');
0042     AX=plotyy(neig,ostd,neig,dqstd);
0043     xlabel('number of eigenvectors');
0044     set(get(AX(1),'Ylabel'),'String','rms vector to correct')
0045     set(get(AX(2),'Ylabel'),'String','rms correctors')
0046     %set(AX(1),'YScale','log');
0047     set(AX(2),'YScale','log');
0048     title('Correctors strength and residual vs number of eigenvectors')
0049     saveas(gca,[datestr(now,'yyyymmddTHHMMSS') '_eigplot.fig']);
0050     export_fig([datestr(now,'yyyymmddTHHMMSS') '_eigplot.jpg']);
0051     close(f);
0052     
0053 end