Tschur=[[Tschur;Zero],a1\(Zschur'*Bv-[Tschur*a;0])]; Zero=[Zero,0];      k=k+1;      if ischar(target), Target(k,:)=[nt,0,0];      else, Target(k,:)=[0,target]; end      if DISP, ShowEig(theta,target,k); end      if (k>=nselect), break; end;        UpdateMinvZ;      J=[2:j]; j=j-1; rKNOWN=0; DETECTED=1;       Ul=Ul(:,J);       V=V*Ul; AV=AV*Ul; BV=BV*Ul;       WAV=Ul'*WAV*Ul; WBV=Ul'*WBV*Ul;       Ul=eye(j); Ur=Ul;      %%% check for conjugate pair      if PAIRS & (abs(imag(theta(2)/theta(1)))>tol)         t=ImagVector(q);          if norm(t)>tol,              %%% t perp Zschur, t in span(Q0,imag(q))             t=t-Q0*(ZastQ\(Zschur'*t));            if norm(t)>100*tol               target=ScaleEig(conj(theta));               EXPAND=1; DETECTED=0; USE_OLD=0;                if DISP, fprintf('--- Checking for conjugate pair ---\n'), end            end         end      end      INITIATE = ( j==0 & DETECTED);   elseif DETECTED %%% To detect whether another eigenpair is accurate enough      INITIATE=1;   end % if (nr<tol)   %%% restart if dim(V)> jmax   if j==jmax      j=jmin; J=[1:j];      Ur=Ur(:,J);       V=V*Ur; AV=AV*Ur; BV=BV*Ur;       WAV=Ur'*WAV*Ur; WBV=Ur'*WBV*Ur;       Ur=eye(j);    end % if jmaxend % while kQschur=Q0;endtime_needed=etime(clock,time);if JDV0 & extra>0 & DISP  fprintf('\n\n# j-dim. proj.: %2i\n\n',extra)endI=CheckSortSchur(Sigma,kappa); Target(1:length(I),:)=Target(I,:);XKNOWN=0;if nargout == 0    if ~DISP   eigenvalues=diag(Sschur)./diag(Tschur)   % Result(eigenvalues)   return, endelse  Jordan=[]; X=zeros(n,0);  if SCHUR ~= 1    if k>0      [Z,D,Jor]=FindJordan(Sschur,Tschur,SCHUR);       DT=abs(diag(D)); DS=abs(diag(Jor));      JT=find(DT<=tol & DS>tol); JS=find(DS<=tol & DT<=tol);      msg=''; DT=~isempty(JT); DS=~isempty(JS);       if DT          msg1='The eigenvalues'; msg2=sprintf(', %i',JT);          msg=[msg1,msg2,' are numerically ''Inf'''];      end,       if DS          msg1='The pencil is numerically degenerated in the directions';          msg2=sprintf(', %i',JS);           if DT, msg=[msg,sprintf('\n\n')]; end, msg=[msg,msg1,msg2,'.'];      end,       if (DT | DS), warndlg(msg,'Unreliable directions'), end      Jordan=Jor/D; X=Qschur*Z; XKNOWN=1;    end   end  [varargout{1:nargout}]=output(history,SCHUR,X,Jordan);end%-------------- display results -----------------------------------------if DISP & size(history,1)>0  rs=history(:,1); mrs=max(rs);  if mrs>0, rs=rs+0.1*eps*mrs;    subplot(2,1,1); t=history(:,2);     plot(t,log10(rs),'*-',t,log10(tol)+0*t,':')    legend('log_{10} || r_{#it} ||_2')    String=sprintf('The test subspace is computed as %s.',testspace);    title(String)    subplot(2,1,2); t=history(:,3);    plot(t,log10(rs),'-*',t,log10(tol)+0*t,':')    legend('log_{10} || r_{#MV} ||_2')      String=sprintf('JDQZ with jmin=%g, jmax=%g, residual tolerance %g.',...            jmin,jmax,tol);     title(String)     String=sprintf('Correction equation solved with %s.',lsolver);    xlabel(String),     date=fix(clock);    String=sprintf('%2i-%2i-%2i, %2i:%2i:%2i',date(3:-1:1),date(4:6));    ax=axis; text(0.2*ax(1)+0.8*ax(2),1.2*ax(3)-0.2*ax(4),String)    drawnow  end  Result(Sigma,Target,diag(Sschur),diag(Tschur),tol)end%------------------------ TEST ACCURACY ---------------------------------if k>nselect & DISP   fprintf('\n%i additional eigenpairs have been detected.\n',k-nselect)endif k<nselect & DISP   fprintf('\nFailed to detect %i eigenpairs.\n',nselect-k)endif (k>0) & DISP   Str='time_needed';                      texttest(Str,eval(Str))   fprintf('\n%39s: %9i','Number of Operator actions',Operator_MVs)   if Precond_Solves   fprintf('\n%39s: %9i','Number of preconditioner solves',Precond_Solves)   end   if 1   if SCHUR ~= 1 & XKNOWN     % Str='norm(Sschur*Z-Tschur*Z*Jordan)'; texttest(Str,eval(Str),tol0)     ok=1; eval('[AX,BX]=MV(X);','ok=0;')     if ~ok, for j=1:size(X,2), [AX(:,j),BX(:,j)]=MV(X(:,j)); end, end     Str='norm(AX*D-BX*Jor)';              texttest(Str,eval(Str),tol0)   end   ok=1; eval('[AQ,BQ]=MV(Qschur);','ok=0;')   if ~ok, for j=1:size(Qschur,2), [AQ(:,j),BQ(:,j)]=MV(Qschur(:,j)); end, end   if kappa == 1     Str='norm(AQ-Zschur*Sschur)';         texttest(Str,eval(Str),tol0)   else     Str='norm(AQ-Zschur*Sschur)/kappa';   texttest(Str,eval(Str),tol0)   end   Str='norm(BQ-Zschur*Tschur)';           texttest(Str,eval(Str),tol0)   I=eye(k);   Str='norm(Qschur''*Qschur-I)';          texttest(Str,eval(Str))     Str='norm(Zschur''*Zschur-I)';          texttest(Str,eval(Str))   nrmSschur=max(norm(Sschur),1.e-8);   nrmTschur=max(norm(Tschur),1.e-8);   Str='norm(tril(Sschur,-1))/nrmSschur';  texttest(Str,eval(Str))   Str='norm(tril(Tschur,-1))/nrmTschur';  texttest(Str,eval(Str))   end   fprintf('\n==================================================\n')endif k==0  disp('no eigenvalue could be detected with the required precision')endreturn%%%======== END JDQZ ====================================================%%%======================================================================%%%======== PREPROCESSING ===============================================%%%======================================================================%%%======== ARNOLDI (for initial spaces) ================================        function [V,AV,BV]=Arnoldi(v,Av,Bv,sigma,jmin,nselect,tol)% Apply Arnoldi with M\(A*sigma(1)'+B*sigma(2)'), to construct an % initial search subspace%global Qschurif ischar(sigma), sigma=[0,1]; end[n,j]=size(v); k=size(Qschur,2); jmin=min(jmin,n-k);if j==0 & k>0  v=RepGS(Qschur,rand(n,1)); [Av,Bv]=MV(v); j=1; endV=v; AV=Av; BV=Bv;while j<jmin;   v=[Av,Bv]*sigma';   v0=SolvePrecond(v);   if sigma(1)==0 & norm(v0-v)<tol,    %%%% then precond=I and target = 0: apply Arnoldi with A      sigma=[1,0]; v0=Av;   end   v=RepGS([Qschur,V],v0); V=[V,v];    [Av,Bv]=MV(v); AV=[AV,Av]; BV=[BV,Bv]; j=j+1; end % whilereturn%%%======== END ARNOLDI =================================================%%%======================================================================%%%======== POSTPROCESSING ==============================================%%%======================================================================%%%======== SORT QZ DECOMPOSITION INTERACTION MATRICES ==================function I=CheckSortSchur(Sigma,kappa)% I=CheckSortSchur(Sigma)%   Scales Qschur, Sschur, and Tschur such that diag(Tschur) in [0,1]%   Reorders the Partial Schur decomposition such that the `eigenvalues'%   (diag(S),diag(T)) appear in increasing chordal distance w.r.t. to%   Sigma.  %   If diag(T) is non-singular then Lambda=diag(S)./diag(T) are the%   eigenvalues.global Qschur Zschur Sschur Tschurk=size(Sschur,1); if k==0, I=[]; return, end% [AQ,BQ]=MV(Qschur);%   Str='norm(AQ-Zschur*Sschur)';                texttest(Str,eval(Str))%   Str='norm(BQ-Zschur*Tschur)';                texttest(Str,eval(Str))%--- scale such that diag(Tschur) in [0,1] ----[Tschur,D]=ScaleT(Tschur); Sschur=D\Sschur;% kappa=max(norm(Sschur,inf)/norm(Tschur,inf),1);s=diag(Sschur); t=diag(Tschur);I=(1:k)'; l=size(Sigma,1);for j=1:k   J0=(j:k)';  J=SortEig(s(I(J0)),t(I(J0)),Sigma(min(j,l),:),kappa);  I(J0)=I(J0(J));endif ~min((1:k)'==I)   [Q,Z,Sschur,Tschur]=SwapQZ(eye(k),eye(k),Sschur,Tschur,I);    [Tschur,D2]=ScaleT(Tschur); Sschur=D2\Sschur;   Qschur=Qschur*Q; Zschur=Zschur*(D*Z*D2);else   Zschur=Zschur*D;endreturn%========================================================================function [T,D]=ScaleT(T)% scale such that diag(T) in [0,1] ----   n=sign(diag(T)); n=n+(n==0); D=diag(n);   T=D\T; IT=imag(T); RT=real(T);   T=RT+IT.*(abs(IT)>eps*abs(RT))*sqrt(-1);return%%%======== COMPUTE SORTED JORDAN FORM ==================================function [X,D,Jor]=FindJordan(S,T,SCHUR)% [X,D,J]=FINDJORDAN(S,T)%   For S and T k by k upper triangular matrices%   FINDJORDAN computes the Jordan decomposition.%   X is a k by k matrix of eigenvectors and principal vectors%   D and J are k by matrices, D is diagonal, J is Jordan%   such that S*X*D=T*X*J. (diag(D),diag(J)) are the eigenvalues.%   If D is non-singular then Lambda=diag(J)./diag(D)%   are the eigenvalues.% coded by Gerard Sleijpen, May, 2002k=size(S,1);s=diag(S); t=diag(T); n=sign(t); n=n+(n==0); D=sqrt(conj(s).*s+conj(t).*t).*n;S=diag(D)\S; T=diag(D)\T; D=diag(diag(T));if k<1,   if k==0, X=[]; D=[]; Jor=[]; end  if k==1, X=1; Jor=s; end  returnendtol=k*(norm(S,1)+norm(T,1))*eps;[X,Jor,I]=PseudoJordan(S,T,tol);if SCHUR == 0for l=1:length(I)-1  if I(l)<I(l+1)-1,     J=[I(l):I(l+1)-1];      [U,JJor]=JordanBlock(Jor(J,J),tol);    X(:,J)=X(:,J)*U; Jor(J,J)=JJor;  endendendJor=Jor+diag(diag(S)); Jor=Jor.*(abs(Jor)>tol);return%==================================================function [X,Jor,J]=PseudoJordan(S,T,delta)% Computes a pseudo-Jordan decomposition for the upper triangular % matrices S and T with ordered diagonal elements. % S*X*(diag(diag(T)))=T*X*(diag(diag(S))+Jor) % with X(:,i:j) orthonormal if its % columns span an invariant subspace of (S,T).k=size(S,1); s=diag(S); t=diag(T); Jor=zeros(k); X=eye(k); J=1;for i=2:k  I=[1:i];   C=t(i,1)*S(I,I)-s(i,1)*T(I,I); C(i,i)=norm(C,inf);  if C(i,i)>0    tol=delta*C(i,i);    for j=i:-1:1       if j==1 | abs(C(j-1,j-1))>tol, break; end    end    e=zeros(i,1); e(i,1)=1;    if j==i      J=[J,i]; q=C\e; X(I,i)=q/norm(q);    else      q=X(I,j:i-1);       q=[C,T(I,I)*q;q',zeros(i-j)]\[e;zeros(i-j,1)];      q=q/norm(q(I,1)); X(I,i)=q(I,1);      Jor(j:i-1,i)=-q(i+1:2*i-j,1);    end  endendJ=[J,k+1];return%==================================================function [X,Jor,U]=JordanBlock(A,tol)%  If A is nilpotent, then A*X=X*Jor with%  Jor a Jordan block%k=size(A,1); Id=eye(k); U=Id; aa=A; j=k; jj=[]; J=1:k;while j>0  [u,s,v]=svd(aa); U(:,J)=U(:,J)*v;  sigma=diag(s); delta=tol;  J=find(sigma<delta);   if isempty(J),j=0; else, j=min(J)-1; end  jj=[jj,j]; if j==0, break, end  aa=v'*u*s; J=1:j; aa=aa(J,J);end Jor=U'*A*U; Jor=Jor.*(abs(Jor)>tol);l=length(jj); jj=[jj(l:-1:1),k];l2=jj(2)-jj(1); J=jj(1)+(1:l2); JX=Id(:,J); X=Id;for j=2:l  l1=l2+1; l2=jj(j+1)-jj(j);   J2=l1:l2; J=jj(j)+(1:l2);  JX=Jor*JX; D=diag(sqrt(diag(JX'*JX))); JX=JX/D;  [Q,S,V]=svd(JX(J,:));  JX=[JX,Id(:,J)*Q(:,J2)]; X(:,J)=JX;endJ=[];for i=1:l2  for k=l:-1:1    j=jj(k)+i; if j<=jj(k+1), J=[J,j]; end  endendX=X(:,J); Jor=X\(Jor*X); X=U*X;Jor=Jor.*(abs(Jor)>100*tol);return%%%======== END JORDAN FORM =============================================%%%======== OUTPUT ======================================================function varargout=output(history,SCHUR,X,Lambda)global Qschur Zschur Sschur Tschurif nargout == 1, varargout{1}=diag(Sschur)./diag(Tschur); return, endif nargout > 2,  varargout{nargout}=history;        endif nargout < 6 & SCHUR == 1  if nargout >1, varargout{1}=Qschur; varargout{2}=Zschur; end  if nargout >2, varargout{3}=Sschur; end  if nargout >3, varargout{4}=Tschur; endend%-------------- compute eigenpairs --------------------------------------if SCHUR ~= 1   varargout{1}=X; varargout{2}=Lambda;    if nargout >3, varargout{3}=Qschur; varargout{4}=Zschur; end   if nargout >4, varargout{5}=Sschur; end   if nargout >5, varargout{6}=Tschur; endendreturn%%%======================================================================%%%======== UPDATE PRECONDITIONED SCHUR VECTORS =========================%%%======================================================================function   UpdateMinvZglobal Qschur Zschur MinvZ QastMinvZ  [n,k]=size(Qschur);  if k==1, MinvZ=zeros(n,0); QastMinvZ = []; end  Minv_z=SolvePrecond(Zschur(:,k));   QastMinvZ=[[QastMinvZ;Qschur(:,k)'*MinvZ],Qschur'*Minv_z];  MinvZ=[MinvZ,Minv_z]; return%%%======================================================================%%%======== SOLVE CORRECTION EQUATION ===================================%%%======================================================================function [t,xtol]=SolvePCE(theta,q,z,r,lsolver,par,nit)global Qschur Zschur  Q=[Qschur,q]; Z=[Zschur,z];  switch lsolver    case 'exact'         [t,xtol] = exact(theta,Q,Z,r);    case {'gmres','cgstab','olsen'}