jdqr.m

一个很好的Matlab编制的数据降维处理软件

function varargout=jdqr(varargin)
%JDQR computes a partial Schur decomposition of a square matrix or operator. 
%  Lambda = JDQR(A) returns the absolute largest eigenvalues in a K vector
%  Lambda. Here K=min(5,N) (unless K has been specified), where N=size(A,1).
%  JDQR(A) (without output argument) displays the K eigenvalues.
%
%  [X,Lambda] = JDQR(A) returns the eigenvectors in the N by K matrix X and 
%  the eigenvalues in the K by K diagonal matrix Lambda. Lambda contains the 
%  Jordan structure if there are multiple eigenvalues.
%
%  [X,Lambda,HISTORY] = JDQR(A) returns also the convergence history
%  (that is, norms of the subsequential residuals).
%
%  [X,Lambda,Q,S] = JDQR(A) returns also a partial Schur decomposition:
%  S is an K by K upper triangular matrix and Q is an N by K orthonormal matrix 
%  such that A*Q = Q*S. The diagonal elements of S are eigenvalues of A.
%
%  [X,Lambda,Q,S,HISTORY] = JDQR(A) returns also the convergence history.
%
%  [X,Lambda,HISTORY] = JDQR('Afun')
%  [X,Lambda,HISTORY] = JDQR('Afun',N)
%  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 a given column vector. In the latter case, the M-file must 
%  return the the order N of the problem with N = Afun([],'dimension') or 
%  N must be specified in the list of input arguments.
%  For example, EIGS('fft',...) is much faster than EIGS(F,...)
%  where F is the explicit FFT matrix.
%
%  The remaining input arguments are optional and can be given in
%  practically any order:
%  ... = JDQR(A,K,SIGMA,OPTIONS) 
%  ... = JDQR('Afun',K,SIGMA,OPTIONS)
%  where
%
%      K         An integer, the number of eigenvalues desired.
%      SIGMA     A scalar shift or a two letter string.
%      OPTIONS   A structure containing additional parameters.
%
%  With one output argument, S is a vector containing K eigenvalues.
%  With two output arguments, S is a K-by-K upper triangular matrix 
%  and Q is a matrix with K columns so that A*Q = Q*S and Q'*Q=I.
%  With three output arguments, HISTORY contains the convergence history.
%
%  If K is not specified, then K = MIN(N,5) eigenvalues are computed.
%
%  If SIGMA is not specified, then the K-th eigenvalues largest in magnitude
%  are computed.  If SIGMA is zero, then the K-th eigenvalues smallest in
%  magnitude are computed.  If SIGMA is a real or complex scalar then the 
%  K-th eigenvalues nearest SIGMA are computed.  If SIGMA is one of the
%  following strings, then it specifies the desired eigenvalues.
%
%    SIGMA             Location wanted eigenvalues
%
%    'LM'              Largest Magnitude  (the default)
%    '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.)
%
%
%  The OPTIONS structure specifies certain parameters in the algorithm.
%
%   Field name         Parameter                            Default
%
%   OPTIONS.Tol        Convergence tolerance:               1e-8 
%                      norm(A*Q-Q*S,1) <= tol * norm(A,1)   
%   OPTIONS.jmin       minimum dimension search subspace    k+5
%   OPTIONS.jmax       maximum dimension search subspace    jmin+5
%   OPTIONS.MaxIt      Maximum number of iterations.        100
%   OPTIONS.v0         Starting space                       ones+0.1*rand
%   OPTIONS.Schur      Gives schur decomposition            'no'
%                      also in case of 2 or 3 output 
%                      arguments (X=Q, Lambda=R).
%
%   OPTIONS.TestSpace  For using harmonic Ritz values       'Standard'
%                      If 'TestSpace'='Harmonic' then 
%                      sigma=0 is the default value for SIGMA
%
%   OPTIONS.Disp       Shows size of intermediate residuals  1
%                      and displays the appr. eigenvalues.
%
%   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                       LU=[[],[]].
%
% For instance
%
%   OPTIONS=STRUCT('Tol',1.0e-10,'LSolver','BiCGstab','LS_ell',2,'Precond',M);
%
% changes the convergence tolerance to 1.0e-10, takes BiCGstab(2) as linear 
% solver and the preconditioner defined in M.m if M is the string 'M', 
% or M = L*U if M is an n by 2*n matrix: M = [L,U].
%
% The preconditoner can be specified in the OPTIONS structure,
% but also in the argument list:
%  ... = JDQR(A,K,SIGMA,M,OPTIONS) 
%  ... = JDQR(A,K,SIGMA,L,U,OPTIONS) 
%  ... = JDQR('Afun',K,SIGMA,'M',OPTIONS)
%  ... = JDQR('Afun',K,SIGMA,'L','U',OPTIONS)
% as an N by N matrix M (then M is the preconditioner), or an N by 2*N 
% matrix M (then  L*U is the preconditioner, where  M = [L,U]),
% or as N by N matrices L and U (then  L*U is the preconditioner),
% or as one or two strings containing the name of  M-files ('M', or 
% 'L' and 'U') which apply a linear operator to a given column vector.
%
%  JDQR (without input arguments) lists the options and the defaults.


%   Gerard Sleijpen.
%   Copyright (c) 98
%
%

% 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, 2007

global Qschur Rschur PinvQ Pinv_u Pu nm_operations

if nargin==0, possibilities, return, end

%%% Read/set parameters
[n,nselect,sigma,SCHUR,...
   jmin,jmax,tol,maxit,V,INTERIOR,SHOW,PAIRS,JDV0,t_tol,...
   lsolver,LSpar] = ReadOptions(varargin{1:nargin});
LSpar0=LSpar; JDV=0; tol0=tol; LOCK0=~ischar(sigma); 
if nargout>3, SCHUR=0; end
tau=0; if INTERIOR>=1 & LOCK0, tau=sigma(1); end
n_tar=size(sigma,1); nt=1; FIG=gcf;

%%% Initiate global variables
Qschur = zeros(n,0); Rschur = []; 
PinvQ  = zeros(n,0); Pinv_u = zeros(n,1); Pu = [];
nm_operations = 0; history = [];

%%% Return if eigenvalueproblem is trivial
if n<2
  if n==1, Qschur=1; Rschur=MV(1); end
  if nargout == 0, eigenvalue=Rschur, else
  [varargout{1:nargout}]=output(history,Qschur,Rschur); end, 
return, end

String = ['\r#it=%i #MV=%i dim(V)=%i |r_%i|=%6.1e  '];
StrinP = '--- Checking for conjugate pair ---\n';
time = clock;

%%% Initialize V, W:
%%%   V,W orthonormal, A*V=W*R+Qschur*E, R upper triangular
[V,W,R,E,M]=SetInitialSpaces(V,nselect,tau,jmin,tol); 
j=size(V,2); k=size(Rschur,1);
nit=0; nlit=0; SOLVED=0;

switch INTERIOR

case 0

%%% The JD loop (Standard)
%%%    V orthogonal, V orthogonal to Qschur
%%%    V*V=eye(j), Qschur'*V=0, 
%%%    W=A*V, M=V'*W
%%%
W=W*R; if tau ~=0; W=W+tau*V; end, M=M'*R; temptarget=sigma(nt,:);  
while (k<nselect) & (nit < maxit) 

   %%% Compute approximate eigenpair and residual
   [UR,S]=SortSchur(M,temptarget,j==jmax,jmin); 
   y=UR(:,1); theta=S(1,1); u=V*y; w=W*y; 
   r=w-theta*u; [r,s]=RepGS(Qschur,r,0); nr=norm(r); r_KNOWN=1;
   if LOCK0 & nr<t_tol, temptarget=[theta;sigma(nt,:)]; end

          % defekt=abs(norm(RepGS(Qschur,MV(u)-theta*u,0))-nr); 
          % DispResult('defekt',defekt,3)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   history=[history;nr,nit,nm_operations];                         %%%
   if SHOW, fprintf(String,nit,nm_operations,j,nlit,nr)            %%%
     if SHOW == 2, LOCK =  LOCK0 & nr<t_tol;                       %%%
       if MovieTheta(n,nit,diag(S),jmin,sigma(nt,:),LOCK,j==jmax)  %%%
   break, end, end, end                                            %%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


   %%% Check for convergence
   if nr<tol

      %%% Expand the partial Schur form
      Qschur=[Qschur,u]; 
      %% Rschur=[[Rschur;zeros(1,k)],Qschur'*MV(u)]; k=k+1; 
      Rschur=[Rschur,s;zeros(1,k),theta];  k=k+1;

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      if SHOW, ShowLambda(theta,k), end %%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      if k>=nselect, break, end, r_KNOWN=0; 
      

      %%% Expand preconditioned Schur matrix PinvQ
      SOLVED=UpdateMinv(u,SOLVED);

      if j==1, 
         [V,W,R,E,M]=SetInitialSpaces(zeros(n,0),nselect,tau,jmin,tol); 
         k=size(Rschur,1); if k>=nselect, break, end
         W=W*R; if tau ~=0; W=W+tau*V; end; M=M'*R; j=size(V,2);
      else,
         J=[2:j]; j=j-1; UR=UR(:,J); 
         M=S(J,J); V=V*UR; W=W*UR; 
      end

     if PAIRS & abs(imag(theta))>tol, v=imag(u/sign(max(u)));
       if norm(v)>tol, v=RepGS(Qschur,v,0); EXPAND=(norm(v)>sqrt(tol)); end
     end

     if EXPAND, temptarget=conj(theta); if SHOW, fprintf(StrinP), end
     else, nlit=0; nt=min(nt+1,n_tar); temptarget=sigma(nt,:); end

   end % nr<tol

   %%% Check for shrinking the search subspace
   if j>=jmax
      j=jmin; J=[1:j]; UR=UR(:,J);
      M=S(J,J); V=V*UR; W=W*UR;
   end % if j>=jmax
 
   if r_KNOWN
      %%% Solve correction equation
      v=Solve_pce(theta,u,r,lsolver,LSpar,nlit); SOLVED=1;
      nlit=nlit+1; nit=nit+1; r_KNOWN=0; EXPAND=1;
   end % if r_KNOWN

   if EXPAND
      %%% Expand the subspaces of the interaction matrix  
      v=RepGS([Qschur,V],v); 
      if size(v,2)>0
        w=MV(v);
        M=[M,V'*w;v'*W,v'*w]; 
        V=[V,v]; W=[W,w]; j=j+1; EXPAND=0; tol=tol0;
      else
        tol=2*tol;
      end
   end % if EXPAND

end % while (nit<maxit)


case 1

%%% The JD loop (Harmonic Ritz values)
%%%    Both V and W orthonormal and orthogonal w.r.t. Qschur
%%%    V*V=eye(j), Qschur'*V=0, W'*W=eye(j), Qschur'*W=0
%%%    (A*V-tau*V)=W*R+Qschur*E, E=Qschur'*(A*V-tau*V), M=W'*V
%%%
temptarget=0; FIXT=1; lsolver0=lsolver;
while (k<nselect) & (nit<maxit) 

   %%% Compute approximate eigenpair and residual
   [UR,UL,S,T]=SortQZ(R,M,temptarget,j>=jmax,jmin);
   y=UR(:,1); theta=T(1,1)'*S(1,1); 
   u=V*y; w=W*(R*y); r=w-theta*u; nr=norm(r); r_KNOWN=1; 
   if nr<t_tol, temptarget=[theta;0]; end, theta=theta+tau;

           % defekt=abs(norm(RepGS(Qschur,MV(u)-theta*u,0))-nr); 
           % DispResult('defect',defekt,3)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   history=[history;nr,nit,nm_operations];                           %%%
   if SHOW, fprintf(String,nit,nm_operations,j,nlit,nr)              %%%
     if SHOW == 2, Lambda=diag(S)./diag(T)+tau; Lambda(1)=theta;     %%%
       if MovieTheta(n,nit,Lambda,jmin,sigma(nt,:),nr<t_tol,j==jmax) %%%
   break, end, end, end                                              %%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   %%% Check for convergence
   if nr<tol

     %%% Expand the partial Schur form
     Qschur=[Qschur,u]; 
     %% Rschur=[[Rschur;zeros(1,k)],Qschur'*MV(u)]; k=k+1;
     Rschur=[Rschur,E*y;zeros(1,k),theta];   k=k+1;  

     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     if SHOW, ShowLambda(theta,k), end %%
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     if k>=nselect, break, end, r_KNOWN=0; JDV=0;

     %%% Expand preconditioned Schur matrix PinvQ
     SOLVED=UpdateMinv(u,SOLVED);

     if j==1,
       [V,W,R,E,M]=...
         SetInitialSpaces(zeros(n,0),nselect,tau,jmin,tol); 
       k=size(Rschur,1); if k>=nselect, break, end, j=size(V,2);
     else
       J=[2:j]; j=j-1; UR=UR(:,J); UL=UL(:,J);
       R=S(J,J); M=T(J,J); V=V*UR; W=W*UL;    
       [r,a]=RepGS(u,r,0);      E=[E*UR;(T(1,1)'-a/S(1,1))*S(1,J)];
                                     
       s=(S(1,J)/S(1,1))/R; W=W+r*s; M=M+s'*(r'*V);
       if (nr*norm(s))^2>eps, [W,R0]=qr(W,0); R=R0*R; M=R0'\M; end

     end

     if PAIRS & abs(imag(theta))>tol, v=imag(u/sign(max(u)));
       if norm(v)>tol, v=RepGS(Qschur,v,0); EXPAND=(norm(v)>sqrt(tol)); end
     end

     if EXPAND, if SHOW, fprintf(StrinP), end
       temptarget=[conj(theta)-tau;0];
     else, nlit=0; temptarget=0;
       if nt<n_tar
         nt=nt+1; tau0=tau; tau=sigma(nt,1); tau0=tau0-tau;
         [W,R]=qr(W*R+tau0*V,0); M=W'*V;
       end
     end

   end

   %%% Check for shrinking the search subspace
   if j>=jmax
      j=jmin; J=[1:j]; UR=UR(:,J); UL=UL(:,J);
      R=S(J,J); M=T(J,J); V=V*UR; W=W*UL;           E=E*UR;
   end % if j>=jmax

   if r_KNOWN
      %%% Solve correction equation
      if JDV, disp('Stagnation'),
        LSpar(end-1)=(LSpar(end-1)+15)*2;
        % lsolver='bicgstab'; LSpar=[1.e-2,300,4];
      else
        LSpar=LSpar0; JDV=0; lsolver=lsolver0;
      end
      if nr>0.001 & FIXT, theta=tau; else, FIXT=0; end
      v=Solve_pce(theta,u,r,lsolver,LSpar,nlit);
      nlit=nlit+1; nit=nit+1; r_KNOWN=0; EXPAND=1; SOLVED=1; JDV=0;
   end
 
   if EXPAND
      %%% Expand the subspaces of the interaction matrix  
      [v,zeta]=RepGS([Qschur,V],v);
      if JDV0 & abs(zeta(end,1))/norm(zeta)<0.06, JDV=JDV+1; end
      if size(v,2)>0 
        w=MV(v); if tau ~=0, w=w-tau*v; end
        [w,e]=RepGS(Qschur,w,0); [w,y]=RepGS(W,w); 
        R=[[R;zeros(1,j)],y]; M=[M,W'*v;w'*V,w'*v];  E=[E,e];
        V=[V,v]; W=[W,w]; j=j+1; EXPAND=0; tol=tol0;
      else
        tol=2*tol;
      end
   end
     
end % while (nit<maxit)

case 1.1

%%% The JD loop (Harmonic Ritz values)
%%%    V W AV.
%%%    Both V and W orthonormal and orthogonal w.r.t. Qschur, AV=A*V-tau*V
%%%    V*V=eye(j),  W'*W=eye(j), Qschur'*V=0, Qschur'*W=0, 
%%%    (I-Qschur*Qschur')*AV=W*R, M=W'*V; R=W'*AV;
%%%
AV=W*R; temptarget=0;
while (k<nselect) & (nit<maxit) 

   %%% Compute approximate eigenpair and residual
   [UR,UL,S,T]=SortQZ(R,M,temptarget,j>=jmax,jmin);
   y=UR(:,1); u=V*y; w=AV*y; theta=u'*w;  
   r=w-theta*u; [r,y]=RepGS(Qschur,r,0); nr=norm(r); r_KNOWN=1;
   if nr<t_tol, temptarget=[theta;0]; end, theta=theta+tau;

           % defekt=abs(norm(RepGS(Qschur,MV(u)-theta*u,0))-nr); 
           % DispResult('defect',defekt,3)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   history=[history;nr,nit,nm_operations];                           %%%
   if SHOW, fprintf(String,nit,nm_operations,j,nlit,nr)              %%%
     if SHOW == 2,Lambda=diag(S)./diag(T)+tau; Lambda(1)=theta;      %%%
       if MovieTheta(n,nit,Lambda,jmin,sigma(nt,:),nr<t_tol,j==jmax) %%%
   break, end, end, end                                              %%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
   %%% Check for convergence
   if nr<tol

     %%% Expand the partial Schur form
     Qschur=[Qschur,u]; 
     %% Rschur=[[Rschur;zeros(1,k)],Qschur'*MV(u)]; k=k+1;
     Rschur=[Rschur,y;zeros(1,k),theta];  k=k+1; 

     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     if SHOW, ShowLambda(theta,k), end %%
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     if k>=nselect, break, end, r_KNOWN=0;

     %%% Expand preconditioned Schur matrix PinvQ
     SOLVED=UpdateMinv(u,SOLVED);

     if j==1
       [V,W,R,E,M]=SetInitialSpaces(zeros(n,0),nselect,tau,jmin,tol); 
       k=size(Rschur,1); if k>=nselect, break, end
       AV=W*R; j=size(V,2); 
     else
       J=[2:j]; j=j-1; UR=UR(:,J); UL=UL(:,J);
       AV=AV*UR; R=S(J,J); M=T(J,J); V=V*UR; W=W*UL;
     end
      
     if PAIRS & abs(imag(theta))>tol, v=imag(u/sign(max(u)));
       if norm(v)>tol, v=RepGS(Qschur,v,0); EXPAND=(norm(v)>sqrt(tol)); end
     end

     if EXPAND,if SHOW, fprintf(StrinP), end
       temptarget=[conj(theta)-tau;0];
     else, nlit=0; temptarget=0;
       if nt<n_tar
         nt=nt+1; tau0=tau; tau=sigma(nt,1); tau0=tau0-tau;
         AV=AV+tau0*V; [W,R]=qr(W*R+tau0*V,0); M=W'*V;
       end
     end

   end

   %%% Check for shrinking the search subspace
   if j>=jmax