jdqz.m

数据降维工具箱

      [MZ,QMZ]=FormPM(q,z); 
      %%% solve preconditioned system
      [t,xtol] = feval(lsolver,theta,Q,Z,MZ,QMZ,r,spar(par,nit));

  end
    
return
%------------------------------------------------------------------------
function [MZ,QMZ]=FormPM(q,z)
% compute vectors and matrices for skew projection

global Qschur MinvZ QastMinvZ

  Minv_z=SolvePrecond(z);
  QMZ=[QastMinvZ,Qschur'*Minv_z;q'*MinvZ,q'*Minv_z];
  MZ=[MinvZ,Minv_z];
    
return
%%%======================================================================
%%%======== LINEAR SOLVERS ==============================================
%%%======================================================================
function [x,xtol] = exact(theta,Q,Z,r)
% produces the exact solution if matrices are given
% Is only feasible for low dimensional matrices
% Only of interest for experimental purposes
%
global Operator_A Operator_B

n=size(r,1); 

if ischar(Operator_A)
   [MZ,QMZ]=FormPM(Q(:,end),Z(:,end)); 
   if n>200
     [x,xtol]=SolvePCE(theta,Q,Z,MZ,QMZ,r,'cgstab',[1.0e-10,500,4]); 
   else
     [x,xtol]=SolvePCE(theta,Q,Z,MZ,QMZ,r,'gmres',[1.0e-10,100]); 
   end
   return
end

k=size(Q,2);
Aug=[theta(2)*Operator_A-theta(1)*Operator_B,Z;Q',zeros(k,k)];
x=Aug\[r;zeros(k,1)]; x([n+1:n+k],:)=[]; xtol=1;

% L=eig(full(Aug)); plot(real(L),imag(L),'*'), pause
%%% [At,Bt]=MV(x); At=theta(2)*At-theta(1)*Bt; 
%%% xtol=norm(r-At+Z*(Z'*At))/norm(r); 

return

%%%===== Iterative methods ==============================================
function [r,xtol] = olsen(theta,Q,Z,MZ,M,r,par)
% returns the preconditioned residual as approximate solution
% May be sufficient in case of an excellent preconditioner

  r=SkewProj(Q,MZ,M,SolvePrecond(r)); xtol=0;

return
%------------------------------------------------------------------------
function [x,rnrm] = cgstab(theta,Q,Z,MZ,M,r,par)
% BiCGstab(ell) with preconditioning
% [x,rnrm] = cgstab(theta,Q,Z,MZ,M,r,par)
% Computes iteratively an approximation to the solution 
% of the linear system Q'*x = 0 and Atilde*x=r 
% where Atilde=(I-Z*Z)*(A-theta*B)*(I-Q*Q').
% using (I-MZ*(M\Q'))*inv(K) as preconditioner
%
% This function is specialized for use in JDQZ.
% integer nmv: number of matrix multiplications
% rnrm: relative residual norm
%
%  par=[tol,mxmv,ell] where 
%    integer m: max number of iteration steps
%    real tol: residual reduction
%
% rnrm: obtained residual reduction
%
% -- References: ETNA

% Gerard Sleijpen ([e-mail])
% Copyright (c) 1998, Gerard Sleijpen

% -- Initialization --
%

global Precond_Type

tol=par(1); max_it=par(2); l=par(3); n=size(r,1);
rnrm=1; nmv=0;
 
if max_it < 2 | tol>=1, x=r; return, end
%%% 0 step of bicgstab eq. 1 step of bicgstab
%%% Then x is a multiple of b

TP=Precond_Type;
if TP==0, r=SkewProj(Q,MZ,M,SolvePrecond(r)); tr=r;
else, tr=RepGS(Z,r); end
rnrm=norm(r); snrm=rnrm; tol=tol*snrm;

sigma=1; omega=1; 
x=zeros(n,1); u=zeros(n,1);
J1=2:l+1; 
   
%%% HIST=[0,1];

if TP <2 %% explicit preconditioning
% -- Iteration loop
while (nmv < max_it)

   sigma=-omega*sigma;
   for j = 1:l,
      rho=tr'*r(:,j);  bet=rho/sigma;
      u=r-bet*u;
      u(:,j+1)=PreMV(theta,Q,MZ,M,u(:,j));
      sigma=tr'*u(:,j+1);  alp=rho/sigma;
      r=r-alp*u(:,2:j+1);
      r(:,j+1)=PreMV(theta,Q,MZ,M,r(:,j));
      x=x+alp*u(:,1);
      G(1,1)=r(:,1)'*r(:,1); rnrm=sqrt(G(1,1));
      if rnrm<tol, l=j; J1=2:l+1; r=r(:,1:l+1); break, end
   end
   nmv = nmv+2*l;

   for i=2:l+1 
     G(i,1:i)=r(:,i)'*r(:,1:i); G(1:i,i)=G(i,1:i)'; 
   end
   if TP, g=Z'*r; G=G-g'*g; end
   d=G(J1,1); gamma=G(J1,J1)\d;  
   rnrm=sqrt(real(G(1,1)-d'*gamma));   %%% compute norm in l-space
   %%% HIST=[HIST;[nmv,rnrm/snrm]];

   x=x+r(:,1:l)*gamma;
   if rnrm < tol, break, end     %%% sufficient accuracy. No need to update r,u
   omega=gamma(l,1); gamma=[1;-gamma];
   u=u*gamma; r=r*gamma; 
   if TP, g=g*gamma; r=r-Z*g; end

   % rnrm = norm(r); 
end

else %% implicit preconditioning

I=eye(2*l); v0=I(:,1:l); s0=I(:,l+1:2*l);
y0=zeros(2*l,1); V=zeros(n,2*l); 

while (nmv < max_it)

   sigma=-omega*sigma;
   y=y0; v=v0; s=s0;
   for j = 1:l,
      rho=tr'*r(:,j);  bet=rho/sigma;
      u=r-bet*u;
      if j>1,                 %%% collect the updates for x in l-space
         v(:,1:j-1)=s(:,1:j-1)-bet*v(:,1:j-1); 
      end
      [u(:,j+1),V(:,j)]=PreMV(theta,Q,MZ,M,u(:,j));
      sigma=tr'*u(:,j+1);  alp=rho/sigma;
      r=r-alp*u(:,2:j+1);
      if j>1, 
         s(:,1:j-1)=s(:,1:j-1)-alp*v(:,2:j); 
      end
      [r(:,j+1),V(:,l+j)]=PreMV(theta,Q,MZ,M,r(:,j));
      y=y+alp*v(:,1);  
      G(1,1)=r(:,1)'*r(:,1); rnrm=sqrt(G(1,1));
      if rnrm<tol, l=j; J1=2:l+1; s=s(:,1:l); break, end
   end
   nmv = nmv+2*l;

   for i=2:l+1 
     G(i,1:i)=r(:,i)'*r(:,1:i); G(1:i,i)=G(i,1:i)'; 
   end
   g=Z'*r; G=G-g'*g;         %%% but, do the orth to Z implicitly
   d=G(J1,1); gamma=G(J1,J1)\d;  
   rnrm=sqrt(real(G(1,1)-d'*gamma)); %%% compute norm in l-space
   x=x+V*(y+s*gamma);

   %%% HIST=[HIST;[nmv,rnrm/snrm]];

   if rnrm < tol, break, end  %%% sufficient accuracy. No need to update r,u
   omega=gamma(l,1); gamma=[1;-gamma];
   u=u*gamma; r=r*gamma; 
   g=g*gamma; r=r-Z*g;        %%% Do the orth to Z explicitly
                              %%% In exact arithmetic not needed, but
                              %%% appears to be more stable.

end
end

if TP==1, x=SkewProj(Q,MZ,M,SolvePrecond(x)); end
rnrm = rnrm/snrm;
%%% plot(HIST(:,1),log10(HIST(:,2)+eps),'*'), drawnow
return
%----------------------------------------------------------------------
function [v,rnrm] = gmres0(theta,Q,Z,MZ,M,v,par)
% GMRES
% [x,rnrm] = gmres(theta,Q,Z,MZ,M,v,par)
% Computes iteratively an approximation to the solution 
% of the linear system Q'*x = 0 and Atilde*x=b 
% where Atilde=(I-Z*Z)*(A-theta*B)*(I-Q*Q').
% using (I-MZ*(M\Q'))*inv(K) as preconditioner
%
% If used as implicit preconditioner then FGMRES.
%
% par=[tol,m] where
%  integer m: degree of the minimal residual polynomial
%  real tol: residual reduction
%
% rnrm: obtained residual reduction
%
% -- References: Saad & Schultz SISC 1986

% Gerard Sleijpen ([e-mail])
% Copyright (c) 1998, Gerard Sleijpen

% -- Initialization

global Precond_Type

tol=par(1); max_it=par(2); n = size(v,1);
rnrm = 1; j=0;

if max_it < 2 | tol>=1, return, end 
%%% 0 step of gmres eq. 1 step of gmres
%%% Then x is a multiple of b
 
H = zeros(max_it +1,max_it); Rot=[ones(1,max_it);zeros(1,max_it)];

TP=Precond_Type;

TP=Precond_Type; 
if TP==0
  v=SkewProj(Q,MZ,M,SolvePrecond(v)); rho0 = norm(v); v = v/rho0;
else
  v=RepGS(Z,v); 
end

V = [v];
tol = tol * rnrm; 
y = [ rnrm ; zeros(max_it,1) ];

while (j < max_it) & (rnrm > tol),
  j=j+1;
  [v,w]=PreMV(theta,Q,MZ,M,v); 
  if TP 
    if TP == 2, W=[W,w]; end 
    v=RepGS(Z,v,0); 
  end
  [v,h] = RepGS(V,v); H(1:size(h,1),j) = h;
  V = [V, v]; 
  for i = 1:j-1,
    a = Rot(:,i);
    H(i:i+1,j) = [a'; -a(2) a(1)]*H(i:i+1,j);
  end
  J=[j, j+1];
  a=H(J,j);
  if a(2) ~= 0
     cs = norm(a); 
     a = a/cs; Rot(:,j) = a;
     H(J,j) = [cs; 0];
     y(J) = [a'; -a(2) a(1)]*y(J);
  end 
  rnrm = abs(y(j+1));
end

J=[1:j];  
if TP == 2
  v = W(:,J)*(H(J,J)\y(J));
else
  v = V(:,J)*(H(J,J)\y(J));
end

if TP==1, v=SkewProj(Q,MZ,M,SolvePrecond(v)); end

return
%%%======================================================================
function [v,rnrm] = gmres(theta,Q,Z,MZ,M,v,par)
% GMRES 
% [x,nmv,rnrm] = gmres(theta,Q,Z,MZ,M,v,par)
% Computes iteratively an approximation to the solution 
% of the linear system Q'*x = 0 and Atilde*x=r 
% where Atilde=(I-Z*Z)*(A-theta*B)*(I-Q*Q').
% using (I-MZ*(M\Q'))*inv(K) as preconditioner.
%
% If used as implicit preconditioner, then FGMRES.
%
% par=[tol,m] where
%  integer m: degree of the minimal residual polynomial
%  real tol: residual reduction
%
% nmv:  number of MV with Atilde
% rnrm: obtained residual reduction
%
% -- References: Saad
% Same as gmres0. However this variant uses MATLAB built-in functions
% slightly more efficient (see Sleijpen and van den Eshof).
%
%
% Gerard Sleijpen ([e-mail])
% Copyright (c) 2002, Gerard Sleijpen

global Precond_Type

% -- Initialization
tol=par(1); max_it=par(2); n = size(v,1);
j=0;

if max_it < 2 | tol>=1, rnrm=1; return, end 
%%% 0 step of gmres eq. 1 step of gmres
%%% Then x is a multiple of b
 
H = zeros(max_it +1,max_it); Gamma=1; rho=1;

TP=Precond_Type; 
if TP==0
  v=SkewProj(Q,MZ,M,SolvePrecond(v)); rho0 = norm(v); v = v/rho0;
else
  v=RepGS(Z,v); rho0=1;
end

V = zeros(n,0); W=zeros(n,0);
tol0 = 1/(tol*tol); 
%% HIST=1;
while (j < max_it) & (rho < tol0) 

  V=[V,v]; j=j+1;
  [v,w]=PreMV(theta,Q,MZ,M,v);
  if TP 
    if TP == 2, W=[W,w]; end 
    v=RepGS(Z,v,0); 
  end
  [v,h] = RepGS(V,v); 
  H(1:size(h,1),j)=h; gamma=H(j+1,j);
  
  if gamma==0, break %%% Lucky break-down
  else
    gamma= -Gamma*h(1:j)/gamma; 
    Gamma=[Gamma,gamma];
    rho=rho+gamma'*gamma;
  end     
          
  %% HIST=[HIST;(gamma~=0)/sqrt(rho)]; 
    
end

if gamma==0; %%% Lucky break-down
   e1=zeros(j,1); e1(1)=rho0; rnrm=0; 
   if TP == 2
     v=W*(H(1:j,1:j)\e1); 
   else
     v=V*(H(1:j,1:j)\e1); 
   end 
else %%% solve in least square sense 
   e1=zeros(j+1,1); e1(1)=rho0; rnrm=1/sqrt(rho);
   if TP == 2
     v=W*(H(1:j+1,1:j)\e1); 
   else
     v=V*(H(1:j+1,1:j)\e1); 
   end 
end

if TP==1, v=SkewProj(Q,MZ,M,SolvePrecond(v)); end
%% HIST=log10(HIST+eps); J=[0:size(HIST,1)-1]';
%% plot(J,HIST(:,1),'*'); drawnow
return
%%%======== END SOLVE CORRECTION EQUATION ===============================       

%%%======================================================================
%%%======== BASIC OPERATIONS ============================================
%%%======================================================================
function [Av,Bv]=MV(v)
% [y,z]=MV(x)
%  y=A*x, z=B*x

%  y=MV(x,theta)
%  y=(A-theta*B)*x
%

global Operator_Form Operator_MVs Operator_A Operator_B Operator_Params

  Bv=v;
  switch Operator_Form
     case 1 % both Operator_A and B are strings 
        Operator_Params{1}=v;
        Av=feval(Operator_A,Operator_Params{:}); 
        Bv=feval(Operator_B,Operator_Params{:});
     case 2
        Operator_Params{1}=v;
        [Av,Bv]=feval(Operator_A,Operator_Params{:});
     case 3
        Operator_Params{1}=v;
        Operator_Params{2}='A';
        Av=feval(Operator_A,Operator_Params{:});
        Operator_Params{2}='B';
        Bv=feval(Operator_A,Operator_Params{:});
     case 4
        Operator_Params{1}=v;
        Av=feval(Operator_A,Operator_Params{:}); 
        Bv=Operator_B*v;
     case 5
        Operator_Params{1}=v;
        Av=feval(Operator_A,Operator_Params{:});
     case 6
        Av=Operator_A*v; 
        Operator_Params{1}=v;
        Bv=feval(Operator_B,Operator_Params{:});
     case 7
        Av=Operator_A*v; 
        Bv=Operator_B*v;
     case 8
        Av=Operator_A*v;
  end

  Operator_MVs = Operator_MVs +size(v,2);

% [Av(1:5,1),Bv(1:5,1)], pause
return
%------------------------------------------------------------------------
function y=SolvePrecond(y);