jdqz.m

数据降维工具箱

global Precond_Form Precond_L Precond_U Precond_P Precond_Params Precond_Solves

    switch Precond_Form
      case 0,     
      case 1,     Precond_Params{1}=y; 
                  y=feval(Precond_L,Precond_Params{:}); 
      case 2,     Precond_Params{1}=y; Precond_Params{2}='preconditioner';
                  y=feval(Precond_L,Precond_Params{:});
      case 3,     Precond_Params{1}=y; 
                  Precond_Params{1}=feval(Precond_L,Precond_Params{:}); 
                  y=feval(Precond_U,Precond_Params{:});
      case 4,     Precond_Params{1}=y; Precond_Params{2}='L'; 
                  Precond_Params{1}=feval(Precond_L,Precond_Params{:}); 
                  Precond_Params{2}='U';
                  y=feval(Precond_L,Precond_Params{:});
      case 5,     y=Precond_L\y;
      case 6,     y=Precond_U\(Precond_L\y);
      case 7,     y=Precond_U\(Precond_L\(Precond_P*y));
    end

    if Precond_Form
      Precond_Solves = Precond_Solves +size(y,2);
    end
%% y(1:5,1), pause
return
%------------------------------------------------------------------------
function [v,u]=PreMV(theta,Q,Z,M,v)
% v=Atilde*v

global Precond_Type

  if Precond_Type
    u=SkewProj(Q,Z,M,SolvePrecond(v)); 
    [v,w]=MV(u); v=theta(2)*v-theta(1)*w;
  else
    [v,u]=MV(v); u=theta(2)*v-theta(1)*u;
    v=SkewProj(Q,Z,M,SolvePrecond(u));
  end
  
return
%------------------------------------------------------------------------
function  r=SkewProj(Q,Z,M,r);

   if ~isempty(Q), 
      r=r-Z*(M\(Q'*r));
   end 

return
%------------------------------------------------------------------------
function ppar=spar(par,nit)

k=size(par,2)-2;
ppar=par(1,k:k+2);

if k>1
   if nit>k
      ppar(1,1)=par(1,k)*((par(1,k)/par(1,k-1))^(nit-k));
   else
      ppar(1,1)=par(1,max(nit,1));
   end
end

ppar(1,1)=max(ppar(1,1),1.0e-8);

return
%------------------------------------------------------------------------
function u=ImagVector(u)
% computes "essential" imaginary part of a vector
  maxu=max(u); maxu=maxu/abs(maxu); u=imag(u/maxu);
return
%------------------------------------------------------------------------
function Sigma=ScaleEig(Sigma)
% 
%
  n=sign(Sigma(:,2)); n=n+(n==0);
  d=sqrt((Sigma.*conj(Sigma))*[1;1]).*n;
  Sigma=diag(d)\Sigma;

return
%%%======== COMPUTE r AND z =============================================
function [r,z,nrm,theta]=Comp_rz(E,kappa)
%
% [r,z,nrm,theta]=Comp_rz(E)
%   computes the direction r of the minimal residual,
%   the left projection vector z,
%   the approximate eigenvalue theta
%
% [r,z,nrm,theta]=Comp_rz(E,kappa)
%   kappa is a scaling factor.
%

% coded by Gerard Sleijpen, version Januari 7, 1998

  if nargin == 1
    kappa=norm(E(:,1))/norm(E(:,2)); kappa=2^(round(log2(kappa)));
  end

  if kappa ~=1, E(:,1)=E(:,1)/kappa; end

  [Q,sigma,u]=svd(E,0);   %%% E*u=Q*sigma, sigma(1,1)>sigma(2,2)
  r=Q(:,2); z=Q(:,1); nrm=sigma(2,2); 
  % nrm=nrm/sigma(1,1); nrmz=sigma(1,1)

  u(1,:)=u(1,:)/kappa; theta=[-u(2,2),u(1,2)];
  
return
%%%======== END computation r and z =====================================
%%%======================================================================
%%%======== Orthogonalisation ===========================================
%%%======================================================================
function [V,R]=RepGS(Z,V,gamma)
%
% Orthonormalisation using repeated Gram-Schmidt
%   with the Daniel-Gragg-Kaufman-Stewart (DGKS) criterion
%
% Q=RepGS(V)
%   The n by k matrix V is orthonormalized, that is,
%   Q is an n by k orthonormal matrix and 
%   the columns of Q span the same space as the columns of V
%   (in fact the first j columns of Q span the same space
%   as the first j columns of V for all j <= k).
%
% Q=RepGS(Z,V)
%  Assuming Z is n by l orthonormal, V is orthonormalized against Z:
%  [Z,Q]=RepGS([Z,V])
%
% Q=RepGS(Z,V,gamma)
%  With gamma=0, V is only orthogonalized against Z
%  Default gamma=1 (the same as Q=RepGS(Z,V))
%
% [Q,R]=RepGS(Z,V,gamma)
%  if gamma == 1, V=[Z,Q]*R; else, V=Z*R+Q; end
 
% coded by Gerard Sleijpen, March, 2002

% if nargin == 1, V=Z; Z=zeros(size(V,1),0); end   
if nargin <3, gamma=1; end

[n,dv]=size(V); [m,dz]=size(Z);

if gamma, l0=min(dv+dz,n); else, l0=dz; end
R=zeros(l0,dv); 

if dv==0, return, end
if dz==0 & gamma==0, return, end

% if m~=n
%   if m<n, Z=[Z;zeros(n-m,dz)]; end
%   if m>n, V=[V;zeros(m-n,dv)]; n=m; end
% end

if (dz==0 & gamma)
   j=1; l=1; J=1;
   q=V(:,1); nr=norm(q); R(1,1)=nr;
   while nr==0, q=rand(n,1); nr=norm(q); end, V(:,1)=q/nr;
   if dv==1, return, end    
else
   j=0; l=0; J=[];
end

while j<dv,
   j=j+1; q=V(:,j); nr_o=norm(q); nr=eps*nr_o;
   if dz>0, yz=Z'*q;     q=q-Z*yz;      end
   if l>0,  y=V(:,J)'*q; q=q-V(:,J)*y;  end
   nr_n=norm(q);
  
   while (nr_n<0.5*nr_o & nr_n > nr)
      if dz>0, sz=Z'*q;     q=q-Z*sz;     yz=yz+sz; end
      if l>0,  s=V(:,J)'*q; q=q-V(:,J)*s; y=y+s;    end
      nr_o=nr_n; nr_n=norm(q);                
   end
   if dz>0, R(1:dz,j)=yz; end
   if l>0,  R(dz+J,j)=y;  end

   if ~gamma
     V(:,j)=q;
   elseif l+dz<n, l=l+1; 
     if nr_n <= nr % expand with a random vector
       % if nr_n==0
       V(:,l)=RepGS([Z,V(:,J)],rand(n,1));
       % else % which can be numerical noice
       %   V(:,l)=q/nr_n;
       % end
     else
       V(:,l)=q/nr_n; R(dz+l,j)=nr_n; 
     end
     J=[1:l];
   end 

end % while j

if gamma & l<dv, V=V(:,J); end

return
%%%======== END  Orthogonalisation ======================================

%%%======================================================================
%%%======== Sorts Schur form ============================================
%%%======================================================================
function [Q,Z,S,T]=SortQZ(A,B,tau,kappa,gamma,u)
%
% [Q,Z,S,T]=SortQZ(A,B,tau)
%   A and B are k by k matrices, tau is a complex pair [alpha,beta].
%   SortQZ computes the qz-decomposition of (A,B) with prescribed
%   ordering: A*Q=Z*S, B*Q=Z*T; 
%             Q and Z are unitary k by k matrices,
%             S and T are upper triangular k by k matrices.
%   The ordering is as follows:
%   (diag(S),diag(T)) are the eigenpairs of (A,B) ordered
%   with increasing "chordal distance" w.r.t. tau.
%
%   If tau is a scalar then [tau,1] is used.
%   Default value for tau is tau=0, i.e., tau=[0,1].
%
% [Q,Z,S,T]=SortQZ(A,B,tau,kappa)
%   kappa scales A first: A/kappa. Default kappa=1.
%
% [Q,Z,S,T]=SortQZ(A,B,tau,kappa,gamma)
%   Sorts the first MAX(gamma,1) elements. Default: gamma=k
%
% [Q,Z,S,T]=SortQZ(A,B,tau,kappa,gamma,u)
%   Now, with ordering such that angle u and Q(:,1) is less than 45o,
%   and, except for the first pair, (diag(S),diag(T)) are 
%   with increasing "chordal distance" w.r.t. tau.
%   If such an ordering does not exist, ordering is as without u.
%

% coded by Gerard Sleijpen, version April, 2002

  k=size(A,1); 
  if k==1, Q=1; Z=1; S=A;T=B; return, end
  kk=k-1;
  if nargin < 3, tau=[0,1]; end 
  if nargin < 4, kappa=1;   end 
  if nargin > 4, kk=max(1,min(gamma,k-1)); end

  %%    kappa=max(norm(A,inf)/max(norm(B,inf),1.e-12),1);
  %%    kappa=2^(round(log2(kappa)));

%%%------ compute the qz factorization -------
  [S,T,Z,Q]=qz(A,B); Z=Z';
                 % kkappa=max(norm(A,inf),norm(B,inf));
                 % Str='norm(A*Q-Z*S)';texttest(Str,eval(Str))
                 % Str='norm(B*Q-Z*T)';texttest(Str,eval(Str))

%%%------ scale the eigenvalues --------------
  t=diag(T); n=sign(t); n=n+(n==0); D=diag(n);  
  Q=Q/D; S=S/D; T=T/D; 

%%%------ sort the eigenvalues ---------------
  I=SortEig(diag(S),real(diag(T)),tau,kappa);
               
%%%------ swap the qz form -------------------
  [Q,Z,S,T]=SwapQZ(Q,Z,S,T,I(1:kk)); 
                % Str='norm(A*Q-Z*S)';texttest(Str,eval(Str))
                % Str='norm(B*Q-Z*T)';texttest(Str,eval(Str))



  if nargin < 6 | size(u,2) ~= 1
     return
  else
%%% repeat SwapQZ if angle is too small
     kk=min(size(u,1),k); J=1:kk; u=u(J,1)'*Q(J,:); 
     if abs(u(1,1))>0.7, return, end
     for j=2:kk
        J=1:j;
        if norm(u(1,J))>0.7
          J0=[j,1:j-1];  
          [Qq,Zz,Ss,Tt]=SwapQZ(eye(j),eye(j),S(J,J),T(J,J),J0);
          if abs(u(1,J)*Qq(:,1))>0.71
             Q(:,J)=Q(:,J)*Qq; Z(:,J)=Z(:,J)*Zz; 
             S(J,J)=Ss; S(J,j+1:k)=Zz'*S(J,j+1:k);
             T(J,J)=Tt; T(J,j+1:k)=Zz'*T(J,j+1:k);
                 % Str='norm(A*Q-Z*S)';texttest(Str,eval(Str))
                 % Str='norm(B*Q-Z*T)';texttest(Str,eval(Str))
                 fprintf('  Took %2i:%6.4g\n',j,S(1,1)./T(1,1))
             return
          end
        end
     end
     disp(['  Selection problem:  took ',num2str(1)])
     return
  end

%%%======================================================================
function I=SortEig(s,t,sigma,kappa);
%
% I=SortEig(S,T,SIGMA) sorts the indices of [S,T] as prescribed by SIGMA
%   S and T are K-vectors.
%
%   If SIGMA=[ALPHA,BETA] is a complex pair then 
%     if CHORDALDISTANCE
%       sort [S,T] with increasing chordal distance w.r.t. SIGMA.
%     else 
%       sort S./T with increasing distance w.r.t. SIGMA(1)/SIGMA(2)
%
%  The chordal distance D between a pair A and a pair B is 
%  defined as follows.
%    Scale A by a scalar F such that NORM(F*A)=1.
%    Scale B by a scalar G such that NORM(G*B)=1.
%    Then D(A,B)=SQRT(1-ABS((F*A)*(G*B)')).
%
% I=SortEig(S,T,SIGMA,KAPPA). Kappa is a caling that effects the
% chordal distance: [S,KAPPA*T] w.r.t. [SIGMA(1),KAPPA*SIGMA(2)].

% coded by Gerard Sleijpen, version April 2002

global CHORDALDISTANCE

if ischar(sigma)
  warning off, s=s./t; warning on
  switch sigma
    case 'LM'
      [s,I]=sort(-abs(s));
    case 'SM'
      [s,I]=sort(abs(s));
    case 'LR';
      [s,I]=sort(-real(s));
    case 'SR';
      [s,I]=sort(real(s));
    case 'BE';
      [s,I]=sort(real(s)); I=twistdim(I,1);
  end
elseif CHORDALDISTANCE
  if kappa~=1, t=kappa*t; sigma(2)=kappa*sigma(2); end
  n=sqrt(s.*conj(s)+t.*t);
  [s,I]=sort(-abs([s,t]*sigma')./n);
else
  warning off, s=s./t; warning on
  if sigma(2)==0; [s,I]=sort(-abs(s));
  else, [s,I]=sort(abs(s-sigma(1)/sigma(2)));
  end 
end

return
%------------------------------------------------------------------------
function t=twistdim(t,k)

  d=size(t,k); J=1:d; J0=zeros(1,2*d);
  J0(1,2*J)=J; J0(1,2*J-1)=flipdim(J,2); I=J0(1,J);
  if k==1, t=t(I,:); else, t=t(:,I); end

return
%%%======================================================================
function [Q,Z,S,T]=SwapQZ(Q,Z,S,T,I)
% [Q,Z,S,T]=SwapQZ(QQ,ZZ,SS,TT,P)
%    QQ and ZZ are K by K unitary,  SS and TT are K by K uper triangular.
%    P is the first part of a permutation of (1:K)'.
%
%    Then Q and Z are K by K unitary, S and T are K by K upper triangular,
%    such that, for A = ZZ*SS*QQ' and B = ZZ*T*QQ', we have 
%    A*Q = Z*S, B*Q = Z*T  and LAMBDA(1:LENGTH(P))=LLAMBDA(P) where 
%    LAMBDA=DIAG(S)./DIAGg(T) and LLAMBDA=DIAG(SS)./DIAG(TT).
%
%    Computation uses Givens rotations. 
%

% coded by Gerard Sleijpen, version October 12, 1998
  

  kk=min(length(I),size(S,1)-1);
  j=1; while (j<=kk & j==I(j)), j=j+1; end
  while j<=kk
    i=I(j);
    for k = i-1:-1:j, 
      %%% i>j, move ith eigenvalue to position j 
      J = [k,k+1]; 
      q = T(k+1,k+1)*S(k,J) - S(k+1,k+1)*T(k,J);
      if q(1) ~= 0 
        q = q/norm(q);
        G = [[q(2);-q(1)],q'];
        Q(:,J) = Q(:,J)*G; 
        S(:,J) = S(:,J)*G; T(:,J) = T(:,J)*G;
      end 
      if abs(S(k+1,k))<abs(T(k+1,k)), q=T(J,k); else q=S(J,k); end
      if q(2) ~= 0
        q=q/norm(q);
        G = [q';q(2),-q(1)];
        Z(:,J) = Z(:,J)*G'; 
        S(J,:) = G*S(J,:); T(J,:) = G*T(J,:);
      end 
      T(k+1,k) = 0;
      S(k+1,k) = 0; 
    end
    I=I+(I<i); 
    j=j+1; while (j<=kk & j==I(j)), j=j+1; end
  end

return
%------------------------------------------------------------------------
function [Q,Z,S,T]=SwapQZ0(Q,Z,S,T,I)
%
% [Q,Z,S,T]=sortqz0(A,B,s,t,k)
%    A and B are k by k matrices, t and s are k vectors such that
%    (t(i),s(i)) eigenpair (A,B), i.e. t(i)*A-s(i)*B singular.
%    Computes the Schur form with a ordering prescribed by (t,s):
%       A*Q=Z*S, B*Q=Z*T such that diag(S)./diag(T)=s./t.
%    Computation uses Householder reflections. 
%

% coded by Gerard Sleijpen, version October 12, 1997
  
  k=size(S,1); s=diag(S); t=diag(T); s=s(I,1); t=t(I,1);

  for i=1:k-1
    %%% compute q s.t. C*q=(t(i,1)*S-s(i,1)*T)*q=0
    C=t(i,1)*S(i:k,i:k)-s(i,1)*T(i:k,i:k);
    [q,r,p]=qr(C);          %% C*P=Q*R
    %% check whether last but one diag. elt r nonzero
    j=k-i; while abs(r(j,j))<eps*norm(r); j=j-1; end; j=j+1;
    r(j,j)=1; e=zeros(j,1); e(j,1)=1;
    q=p*([r(1:j,1:j)\e;zeros(k-i+1-j,1)]); q=q/norm(q);%% C*q
    %%% end computation q
    z=conj(s(i,1))*S(i:k,i:k)*q+conj(t(i,1))*T(i:k,i:k)*q; z=z/norm(z);
    a=q(1,1); if a ~=0, a=abs(a)/a; q=a*q; end
    a=z(1,1); if a ~=0, a=abs(a)/a; z=a*z; end