[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 projectionglobal 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_Bn=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   returnendk=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 (sleijpen@math.uu.nl)% Copyright (c) 1998, Gerard Sleijpen% -- Initialization --%global Precond_Typetol=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 bTP=Precond_Type;if TP==0, r=SkewProj(Q,MZ,M,SolvePrecond(r)); tr=r;else, tr=RepGS(Z,r); endrnrm=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 loopwhile (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); endelse %% implicit preconditioningI=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.endendif TP==1, x=SkewProj(Q,MZ,M,SolvePrecond(x)); endrnrm = rnrm/snrm;%%% plot(HIST(:,1),log10(HIST(:,2)+eps),'*'), drawnowreturn%----------------------------------------------------------------------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 (sleijpen@math.uu.nl)% Copyright (c) 1998, Gerard Sleijpen% -- Initializationglobal Precond_Typetol=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); endV = [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));endJ=[1:j];  if TP == 2  v = W(:,J)*(H(J,J)\y(J));else  v = V(:,J)*(H(J,J)\y(J));endif TP==1, v=SkewProj(Q,MZ,M,SolvePrecond(v)); endreturn%%%======================================================================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 (sleijpen@math.uu.nl)% Copyright (c) 2002, Gerard Sleijpenglobal Precond_Type% -- Initializationtol=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;endV = 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)];     endif 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 endif TP==1, v=SkewProj(Q,MZ,M,SolvePrecond(v)); end%% HIST=log10(HIST+eps); J=[0:size(HIST,1)-1]';%% plot(J,HIST(:,1),'*'); drawnowreturn%%%======== 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)], pausereturn%------------------------------------------------------------------------function y=SolvePrecond(y);