jdqz.m

数据降维工具箱

      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 jmax

end % while k
Qschur=Q0;
end

time_needed=etime(clock,time);

if JDV0 & extra>0 & DISP
  fprintf('\n\n# j-dim. proj.: %2i\n\n',extra)
end

I=CheckSortSchur(Sigma,kappa); Target(1:length(I),:)=Target(I,:);

XKNOWN=0;

if nargout == 0 
   if ~DISP
   eigenvalues=diag(Sschur)./diag(Tschur)
   % Result(eigenvalues)
   return, end
else
  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)
end
if k<nselect & DISP
   fprintf('\nFailed to detect %i eigenpairs.\n',nselect-k)
end

if (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')
end
if k==0
  disp('no eigenvalue could be detected with the required precision')
end

return
%%%======== 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 Qschur

if 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; 
end

V=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 % while

return
%%%======== 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 Tschur

k=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));
end

if ~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;
end

return
%========================================================================
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, 2002

k=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
  return
end

tol=k*(norm(S,1)+norm(T,1))*eps;
[X,Jor,I]=PseudoJordan(S,T,tol);


if SCHUR == 0
for 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;
  end
end
end

Jor=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
  end
end
J=[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;
end

J=[];
for i=1:l2
  for k=l:-1:1
    j=jj(k)+i; if j<=jj(k+1), J=[J,j]; end
  end
end

X=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 Tschur

if nargout == 1, varargout{1}=diag(Sschur)./diag(Tschur); return, end
if nargout > 2,  varargout{nargout}=history;        end
if 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; end
end

%-------------- 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; end
end

return
%%%======================================================================
%%%======== UPDATE PRECONDITIONED SCHUR VECTORS =========================
%%%======================================================================
function   UpdateMinvZ

global 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'}