function varargout=jdqz(varargin)%JDQZ computes a partial generalized Schur decomposition (or QZ%  decomposition) of a pair of square matrices or operators.%  %  LAMBDA=JDQZ(A,B) and JDQZ(A,B) return K eigenvalues of the matrix pair%  (A,B), where K=min(5,N) and N=size(A,1) if K has not been specified.%  %  [X,JORDAN]=JDQZ(A,B) returns the eigenvectors X and the Jordan%  structure JORDAN:  A*X=B*X*JORDAN. The diagonal of JORDAN contains the%  eigenvalues: LAMBDA=DIAG(JORDAN). JORDAN is an K by K matrix with the%  eigenvalues on the diagonal and zero or one on the first upper diagonal%  elements. The other entries are zero.%  %  [X,JORDAN,HISTORY]=JDQZ(A,B) returns also the convergence history.%  %  [X,JORDAN,Q,Z,S,T,HISTORY]=JDQZ(A,B) %  If between four and seven output arguments are required, then Q and Z%  are N by K orthonormal, S and T are K by K upper triangular such that%  they form a partial generalized Schur decomposition: A*Q=Z*S and%  B*Q=Z*T. Then LAMBDA=DIAG(S)./DIAG(T) and X=Q*Y with Y the eigenvectors%  of the pair (S,T): S*Y=T*Y*JORDAN (see also OPTIONS.Schur).%  %  JDQZ(A,B) %  JDQZ('Afun','Bfun')%  The first input argument is either a square matrix (which can be full%  or sparse, symmetric or nonsymmetric, real or complex), or a string%  containing the name of an M-file which applies a linear operator to the%  columns of a given matrix. In the latter case, the M-file, say Afun.m,%  must return the dimension N of the problem with N = Afun([],'dimension').%  For example, JDQZ('fft',...) is much faster than JDQZ(F,...), where F is%  the explicit FFT matrix.%  If another input argument is a square N by N matrix or the name of an%  M-file, then B is this argument (regardless whether A is an M-file or a%  matrix). If B has not been specified, then B is assumed to be the%  identity unless A is an M-file with two output vectors of dimension N%  with [AV,BV]=Afun(V), or with AV=Afun(V,'A') and BV=Afun(V,'B').%  %  The remaining input arguments are optional and can be given in%  practically any order:%  %  [X,JORDAN,Q,Z,S,T,HISTORY] = JDQZ(A,B,K,SIGMA,OPTIONS)%  [X,JORDAN,Q,Z,S,T,HISTORY] = JDQZ('Afun','Bfun',K,SIGMA,OPTIONS)%  %  where%  %      K         an integer, the number of desired eigenvalues.%      SIGMA     a scalar shift or a two letter string.%      OPTIONS   a structure containing additional parameters.%  %  If K is not specified, then K = MIN(N,5) eigenvalues are computed.%  %  If SIGMA is not specified, then the Kth eigenvalues largest in%  magnitude are computed. If SIGMA is a real or complex scalar, then the%  Kth eigenvalues nearest SIGMA are computed. If SIGMA is column vector%  of size (L,1), then the Jth eigenvalue nearest to SIGMA(MIN(J,L))%  is computed for J=1:K. SIGMA is the "target" for the desired eigenvalues.%  If SIGMA is one of the following strings, then it specifies the desired %  eigenvalues.%  %    SIGMA            Specified eigenvalues%  %    'LM'             Largest Magnitude  %    'SM'             Smallest Magnitude (same as SIGMA = 0)%    'LR'             Largest Real part%    'SR'             Smallest Real part%    'BE'             Both Ends. Computes K/2 eigenvalues%                     from each end of the spectrum (one more%                     from the high end if K is odd.)%  %  If 'TestSpace' is 'Harmonic' (see OPTIONS), then SIGMA = 0 is the%  default, otherwise SIGMA = 'LM' is the default.%  %  %  The OPTIONS structure specifies certain parameters in the algorithm.%  %   Field name            Parameter                             Default%  %   OPTIONS.Tol           Convergence tolerance:                1e-8 %                           norm(r) <= Tol/SQRT(K)   %   OPTIONS.jmin          Minimum dimension search subspace V   K+5%   OPTIONS.jmax          Maximum dimension search subspace V   jmin+5%   OPTIONS.MaxIt         Maximum number of iterations.         100%   OPTIONS.v0            Starting space                        ones+0.1*rand%   OPTIONS.Schur         Gives schur decomposition             'no'%                           If 'yes', then X and JORDAN are%                           not computed and [Q,Z,S,T,HISTORY]%                           is the list of output arguments.%   OPTIONS.TestSpace     Defines the test subspace W           'Harmonic'%                           'Standard':    W=sigma*A*V+B*V%                           'Harmonic':    W=A*V-sigma*B*V%                           'SearchSpace': W=V%                            W=V is justified if B is positive%                            definite.%   OPTIONS.Disp          Shows size of intermediate residuals  'no'%                           and the convergence history%   OPTIONS.NSigma        Take as target for the second and     'no'%                           following eigenvalues, the best  %                           approximate eigenvalues from the %                           test subspace.  %   OPTIONS.Pairs         Search for conjugated eigenpairs      'no'%   OPTIONS.LSolver       Linear solver                         'GMRES'%   OPTIONS.LS_Tol        Residual reduction linear solver      1,0.7,0.7^2,..%   OPTIONS.LS_MaxIt      Maximum number it.  linear solver     5%   OPTIONS.LS_ell        ell for BiCGstab(ell)                 4%   OPTIONS.Precond       Preconditioner  (see below)           identity.%   OPTIONS.Type_Precond  Way of using preconditioner           'left'%  %  For instance%  %    options=struct('Tol',1.0e-8,'LSolver','BiCGstab','LS_ell',4,'Precond',M);%  %  changes the convergence tolerance to 1.0e-8, takes BiCGstab as linear %  solver, and takes M as preconditioner (for ways of defining M, see below).%%%  PRECONDITIONING. The action M-inverse of the preconditioner M (an %  approximation of A-lamda*B) on an N-vector V can be defined in the %  OPTIONS%  %     OPTIONS.Precond%     OPTIONS.L_Precond     same as OPTIONS.Precond%     OPTIONS.U_Precond%     OPTIONS.P_Precond%%  If no preconditioner has been specified (or is []), then M\V=V (M is%  the identity).%  If Precond is an N by N matrix, say, K, then%        M\V = K\V.%  If Precond is an N by 2*N matrix, say, K, then%        M\V = U\L\V, where K=[L,U], and L and U are N by N matrices.%  If Precond is a string, say, 'Mi', then%        if Mi(V,'L') and Mi(V,'U') return N-vectors %               M\V = Mi(Mi(V,'L'),'U')%        otherwise %               M\V = Mi(V) or M\V=Mi(V,'preconditioner').%  Note that Precond and A can be the same string.%  If L_Precond and U_Precond are strings, say, 'Li' and 'Ui', %  respectively, then%        M\V=Ui(Li(V)).%  If (P_precond,) L_Precond, and U_precond are N by N matrices, say, %  (P,) L, and U, respectively, then%        M\V=U\L\(P*V)      (P*M=L*U)%%     OPTIONS.Type_Precond%  The preconditioner can be used as explicit left preconditioner%  ('left', default), as explicit right preconditioner ('right') or %  implicitly ('impl').%  %%  JDQZ without input arguments returns the options and its defaults.%%   Gerard Sleijpen.%   Copyright (c) 2002%%% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.4b.% The toolbox can be obtained from http://www.cs.unimaas.nl/l.vandermaaten% You are free to use, change, or redistribute this code in any way you% want for non-commercial purposes. However, it is appreciated if you % maintain the name of the original author.%% (C) Laurens van der Maaten% Maastricht University, 2007global Qschur Zschur Sschur Tschur ...       Operator_MVs Precond_Solves ...       MinvZ QastMinvZif nargin==0   possibilities, return,end%%% Read/set parameters[n,nselect,Sigma,kappa,SCHUR,...   jmin,jmax,tol0,maxit,V,AV,BV,TS,DISP,PAIRS,JDV0,FIX_tol,track,NSIGMA,...   lsolver,LSpar] = ReadOptions(varargin{1:nargin});Qschur = zeros(n,0);    Zschur=zeros(n,0);; MinvZ  = zeros(n,0);    QastMinvZ=zeros(0,0); Sschur = []; Tschur=[]; history = []; %%% Return if eigenvalueproblem is trivialif n<2  if n==1, Qschur=1; Zschur=1; [Sschur,Tschur]=MV(1); end  if nargout == 0, Lambda=Sschur/Tschur, else  [varargout{1:nargout}]=output(history,SCHUR,1,Sschur/Tschur); end, return, end%---------- SET PARAMETERS & STRINGS FOR OUTPUT -------------------------if     TS==0, testspace='sigma(1)''*Av+sigma(2)''*Bv';elseif TS==1, testspace='sigma(2)*Av-sigma(1)*Bv';elseif TS==2, testspace='v'; elseif TS==3, testspace='Bv';elseif TS==4, testspace='Av';endString=['\r#it=%i #MV=%3i, dim(V)=%2i, |r_%2i|=%6.1e  '];%------------------- JDQZ -----------------------------------------------% fprintf('Scaling with kappa=%6.4g.',kappa)k=0; nt=0; j=size(V,2); nSigma=size(Sigma,1);it=0; extra=0; Zero=[]; target=[]; tol=tol0/sqrt(nselect);INITIATE=1;  JDV=0; rKNOWN=0; EXPAND=0; USE_OLD=0; DETECTED=0;time=clock;if TS ~=2while (k<nselect & it<maxit)   %%% Initialize target, test space and interaction matrices   if INITIATE, % set new target      nt=min(nt+1,nSigma); sigma = Sigma(nt,:); nlit=0; lit=0;        if j<2        [V,AV,BV]=Arnoldi(V,AV,BV,sigma,jmin,nselect,tol);        rKNOWN=0; EXPAND=0; USE_OLD=0; DETECTED=0; target=[];        j=min(jmin,n-k);      end      if DETECTED & NSIGMA         [Ur,Ul,St,Tt] = SortQZ(WAV,WBV,sigma,kappa);         y=Ur(:,1); q=V*y; Av=AV*y; Bv=BV*y;          [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa);         sigma=ScaleEig(theta);         USE_OLD=NSIGMA; rKNOWN=1; lit=10;      end         NEWSHIFT= 1;       if DETECTED & TS<2, NEWSHIFT= ~min(target==sigma); end      target=sigma; ttarget=sigma;      if ischar(ttarget), ttrack=0; else, ttrack=track; end      if NEWSHIFT          v=V; Av=AV; Bv=BV; W=eval(testspace);         %%% V=RepGS(Qschur,V); [AV,BV]=MV(V); %%% more stability??         %%% W=RepGS(Zschur,eval(testspace));  %%% dangerous if sigma~lambda         if USE_OLD, W(:,1)=V(:,1); end,          W=RepGS(Zschur,W); WAV=W'*AV;  WBV=W'*BV;      end      INITIATE=0; DETECTED=0; JDV=0;   end % if INITIATE   %%% Solve the preconditioned correction equation   if rKNOWN,      if JDV, z=W; q=V; extra=extra+1;          if DISP,  fprintf('  %2i-d proj.\n',k+j-1), end       end      if FIX_tol*nr>1 & ~ischar(target), theta=target; else, FIX_tol=0; end      t=SolvePCE(theta,q,z,r,lsolver,LSpar,lit);       nlit=nlit+1; lit=lit+1; it=it+1;      EXPAND=1; rKNOWN=0; JDV=0;   end % if rKNOWN       %%% Expand the subspaces and the interaction matrices   if EXPAND      [v,zeta]=RepGS([Qschur,V],t);      V=[V,v];       [Av,Bv]=MV(v); AV=[AV,Av]; BV=[BV,Bv];       w=eval(testspace); w=RepGS([Zschur,W],w);      WAV=[WAV,W'*Av;w'*AV]; WBV=[WBV,W'*Bv;w'*BV]; W=[W,w];      j=j+1; EXPAND=0;      %%% Check for stagnation      if abs(zeta(size(zeta,1),1))/norm(zeta)<0.06, JDV=JDV0; end   end % if EXPAND    %%% Solve projected eigenproblem   if USE_OLD      [Ur,Ul,St,Tt]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin,y);    else      [Ur,Ul,St,Tt]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin);    end   %%% Compute approximate eigenpair and residual   y=Ur(:,1); q=V*y; Av=AV*y; Bv=BV*y;    [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa);    %%%=== an alternative, but less stable way of computing z =====   % beta=Tt(1,1); alpha=St(1,1); theta=[alpha,beta];   % r=RepGS(Zschur,beta*Av-alpha*Bv,0); nr=norm(r); z=W*Ul(:,1);   rKNOWN=1; if nr<ttrack, ttarget=ScaleEig(theta); end         if DISP,                                  %%% display history            fprintf(String,it,Operator_MVs,j,nlit,nr),          end          history=[history;nr,it,Operator_MVs];    %%% save history   %%% check convergence    if (nr<tol)      %%% EXPAND Schur form      Qschur=[Qschur,q]; Zschur=[Zschur,z];      Sschur=[[Sschur;zeros(1,k)],Zschur'*Av];       Tschur=[[Tschur;zeros(1,k)],Zschur'*Bv];  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;      %%% Expand preconditioned Schur matrix MinvZ=M\Zschur      UpdateMinvZ;      J=[2:j]; j=j-1; Ur=Ur(:,J); Ul=Ul(:,J);       V=V*Ur; AV=AV*Ur; BV=BV*Ur; W=W*Ul;       WAV=St(J,J); WBV=Tt(J,J);        rKNOWN=0; DETECTED=1;  USE_OLD=0;      %%% check for conjugate pair      if PAIRS & (abs(imag(theta(1)/theta(2)))>tol)          t=ImagVector(q); % t=conj(q); t=t-q*(q'*t);         if norm(t)>tol, t=RepGS([Qschur,V],t,0);             if norm(t)>200*tol               target=ScaleEig(conj(theta));               EXPAND=1; DETECTED=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); Ul=Ul(:,J);       V=V*Ur; AV=AV*Ur; BV=BV*Ur; W=W*Ul;       WAV=St(J,J); WBV=Tt(J,J);    end % if j==jmaxend % while kend % if TS~=2if TS==2Q0=Qschur; ZastQ=[];% WAV=V'*AV; WBV=V'*BV;while (k<nselect & it<maxit)   %%% Initialize target, test space and interaction matrices   if INITIATE & ( nSigma>k | NSIGMA), % set new target      nt=min(nt+1,nSigma); sigma = Sigma(nt,:); nlit=0; lit=0;                      if j<2        [V,AV,BV]=Arnoldi(V,AV,BV,sigma,jmin,nselect,tol);        rKNOWN=0; EXPAND=0; USE_OLD=0; DETECTED=0; target=[];        j=min(jmin,n-k);;       end      if DETECTED & NSIGMA         [Ur,Ul,St,Tt]=SortQZ(WAV,WBV,sigma,kappa,1);          q=RepGS(Zschur,V*Ur(:,1)); [Av,Bv]=MV(q);          [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa);          sigma=ScaleEig(theta);          USE_OLD=NSIGMA; rKNOWN=1; lit=10;      end         target=sigma; ttarget=sigma;      if ischar(ttarget), ttrack=0; else, ttrack=track; end      if ~DETECTED         %%% additional stabilisation. May not be needed         %%% V=RepGS(Zschur,V); [AV,BV]=MV(V);          %%% end add. stab.         WAV=V'*AV; WBV=V'*BV;      end      DETECTED=0; INITIATE=0; JDV=0;    end % if INITIATE   %%% Solve the preconditioned correction equation   if rKNOWN,      if JDV, z=V; q=V; extra=extra+1;          if DISP,  fprintf('  %2i-d proj.\n',k+j-1), end       end      if FIX_tol*nr>1 & ~ischar(target), theta=target; else, FIX_tol=0; end      t=SolvePCE(theta,q,z,r,lsolver,LSpar,lit);       nlit=nlit+1; lit=lit+1; it=it+1;      EXPAND=1; rKNOWN=0; JDV=0;   end % if rKNOWN   %%% expand the subspaces and the interaction matrices   if EXPAND      [v,zeta]=RepGS([Zschur,V],t); [Av,Bv]=MV(v);       WAV=[WAV,V'*Av;v'*AV,v'*Av]; WBV=[WBV,V'*Bv;v'*BV,v'*Bv];      V=[V,v]; AV=[AV,Av]; BV=[BV,Bv];       j=j+1; EXPAND=0;      %%% Check for stagnation      if abs(zeta(size(zeta,1),1))/norm(zeta)<0.06, JDV=JDV0; end    end % if EXPAND    %%% compute approximate eigenpair   if USE_OLD      [Ur,Ul]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin,Ur(:,1));    else      [Ur,Ul]=SortQZ(WAV,WBV,ttarget,kappa,(j>=jmax)*jmin);   end        %%% Compute approximate eigenpair and residual   q=V*Ur(:,1); Av=AV*Ur(:,1); Bv=BV*Ur(:,1);     [r,z,nr,theta]=Comp_rz(RepGS(Zschur,[Av,Bv],0),kappa);    rKNOWN=1; if nr<ttrack, ttarget=ScaleEig(theta); end          if DISP,                                 %%% display history             fprintf(String,it,Operator_MVs, j,nlit,nr),           end            history=[history;nr,it,Operator_MVs];   %%% save history      %%% check convergence    if (nr<tol)      %%% expand Schur form      [q,a]=RepGS(Q0,q); a1=a(k+1,1); a=a(1:k,1);      %%% ZastQ=Z'*Q0        Q0=[Q0,q]; %%% the final Qschur      ZastQ=[ZastQ,Zschur'*q;z'*Q0]; Zschur=[Zschur,z]; Qschur=[Qschur,z];       Sschur=[[Sschur;Zero],a1\(Zschur'*Av-[Sschur*a;0])];