jdqr.m

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

x= zeros(n,1); u=zeros(n,1);

rho = norm(r);  snrm = rho; rho = rho*rho; 
tol = tol*tol*rho; 
sigma=1;

while  ( rho > tol & nmv < max_it )

  beta=rho/sigma;   
  u=r-beta*u;
  y=smvp(theta,Q,Z,M,u); nmv=nmv+1;  
  sigma=y'*r; alpha=rho/sigma;  
  x=x+alpha*u; 
  r=r-alpha*y; sigma=-rho; rho=r'*r;  

end % while

rnrm=sqrt(rho)/snrm;

return
%----------------------------------------------------------------------
function x = minres(theta,Q,Z,M,r,par)
% MINRES
% [x,rnrm] = minres(theta,Q,Z,M,b,par)
% Computes iteratively an approximation to the solution 
% of the linear system Q'*x = 0 and Atilde*x=b 
% where Atilde=(I-Z*M^(-1)*Q')*(U\L\(A-theta)).
%
% par=[tol,m] where
%  integer m: degree of the minimal residual polynomial
%  real tol: residual reduction
%
% nmv:  number of MV with Atilde
% rnrm: obtained residual reduction
%
% -- References: Paige and Saunders

% Gerard Sleijpen ([e-mail])
% Copyright (c) 1998, Gerard Sleijpen

% -- Initialization


tol=par(1); max_it=par(2); n = size(r,1);
rnrm = 1; nmv=0; 

if max_it ==0 | tol>=1, x=r; return, end

x=zeros(n,1); rho = norm(r); v = r/rho; snrm=rho;
beta = 0; v_old = zeros(n,1); 
beta_t = 0; c = -1; s = 0;
w = zeros(n,1); www = v;

tol=tol*rho;

while  ( nmv < max_it  &  abs(rho) > tol )

   wv =smvp(theta,Q,Z,M,v)-beta*v_old;  nmv=nmv+1;  
   alpha = v'*wv; wv = wv-alpha*v;
   beta = norm(wv); v_old = v; v = wv/beta;

   l1 = s*alpha - c*beta_t; l2 = s*beta;

   alpha_t = -s*beta_t - c*alpha;  beta_t = c*beta;
   l0 = sqrt(alpha_t*alpha_t+beta*beta); 
   c = alpha_t/l0; s = beta/l0;

   ww = www - l1*w; www = v - l2*w; w = ww/l0;

   x =  x + (rho*c)*w; rho =  s*rho; 

end % while

rnrm=abs(rho)/snrm;

return

%----------------------------------------------------------------------
function x = symmlq(theta,Q,Z,M,r,par)
% SYMMLQ
% [x,rnrm] = symmlq(theta,Q,Z,M,b,par)
% Computes iteratively an approximation to the solution 
% of the linear system Q'*x = 0 and Atilde*x=b 
% where Atilde=(I-Z*M^(-1)*Q')*(U\L\(A-theta)).
%
% par=[tol,m] where
%  integer m: degree of the minimal residual polynomial
%  real tol: residual reduction
%
% nmv:  number of MV with Atilde
% rnrm: obtained residual reduction
%
% -- References: Paige and Saunders

% Gerard Sleijpen ([e-mail])
% Copyright (c) 1998, Gerard Sleijpen

% -- Initialization


tol=par(1); max_it=par(2); n = size(r,1);
rnrm = 1; nmv=0; 

if max_it ==0 | tol>=1, x=r; return, end

x=zeros(n,1); rho = norm(r); v = r/rho; snrm=rho;
beta = 0;  beta_t = 0; c = -1;  s = 0;
v_old = zeros(n,1); w = v; gtt = rho; g = 0;
      
tol=tol*rho;
   
while  ( nmv < max_it  &  rho > tol )

   wv = smvp(theta,Q,Z,M,v) - beta*v_old; nmv=nmv+1;
   alpha = v'*wv; wv = wv - alpha*v;
   beta = norm(wv); v_old = v;  v = wv/beta;

   l1 = s*alpha - c*beta_t; l2 = s*beta; 
      
   alpha_t = -s*beta_t - c*alpha; beta_t = c*beta;
   l0 = sqrt(alpha_t*alpha_t+beta*beta); 
   c = alpha_t/l0; s = beta/l0;

   gt = gtt - l1*g; gtt = -l2*g; g = gt/l0;

   rho = sqrt(gt*gt+gtt*gtt); 

   x = x + (g*c)*w + (g*s)*v;
   w =  s*w - c*v;  
     
end % while

rnrm=rho/snrm; 

return

%----------------------------------------------------------------------
function v=smvp(theta,Q,Z,M,v)
% v=Atilde*v

   v = SkewProj(Z,Q,M,v);
   v = SolvePrecond(v,'U');
   v = MV(v) - theta*v; 
   v = SolvePrecond(v,'L');
   v = SkewProj(Q,Z,M,v);
  
return
%=======================================================================
%========== Orthogonalisation ==========================================
%=======================================================================
function [v,y]=RepGS(V,v,gamma)
% [v,y]=REP_GS(V,w)
% If V orthonormal then [V,v] orthonormal and w=[V,v]*y;
% If size(V,2)=size(V,1) then w=V*y;
%
% The orthonormalisation uses repeated Gram-Schmidt
% with the Daniel-Gragg-Kaufman-Stewart (DGKS) criterion.
%
% [v,y]=REP_GS(V,w,GAMMA)
% GAMMA=1 (default) same as [v,y]=REP_GS(V,w)
% GAMMA=0, V'*v=zeros(size(V,2)) and  w = V*y+v (v is not normalized).

 
% coded by Gerard Sleijpen, August 28, 1998

if nargin < 3, gamma=1; end

[n,d]=size(V);

if size(v,2)==0, y=zeros(d,0); return, end

nr_o=norm(v); nr=eps*nr_o; y=zeros(d,1);
if d==0
  if gamma, v=v/nr_o; y=nr_o; else, y=zeros(0,1); end, return
end

y=V'*v; v=v-V*y; nr_n=norm(v); ort=0;

while (nr_n<0.5*nr_o & nr_n > nr)
  s=V'*v; v=v-V*s; y=y+s; 
  nr_o=nr_n; nr_n=norm(v);     ort=ort+1; 
end

if nr_n <= nr, if ort>2, disp(' dependence! '), end
  if gamma  % and size allows, expand with a random vector
    if d<n, v=RepGS(V,rand(n,1)); y=[y;0]; else, v=zeros(n,0); end
  else, v=0*v; end
elseif gamma, v=v/nr_n; y=[y;nr_n]; end

return
%=======================================================================
%============== Sorts Schur form =======================================
%=======================================================================
function [Q,S]=SortSchur(A,sigma,gamma,kk)
%[Q,S]=SortSchur(A,sigma)
%  A*Q=Q*S with diag(S) in order prescribed by sigma.
%  If sigma is a scalar then with increasing distance from sigma.
%  If sigma is string then according to string
%  ('LM' with decreasing modulus, etc)
%
%[Q,S]=SortSchur(A,sigma,gamma,kk)
%  if gamma==0, sorts only for the leading element
%  else, sorts for the kk leading elements

  l=size(A,1);
  if l<2, Q=1;S=A; return, 
  elseif nargin==2, kk=l-1; 
  elseif gamma, kk=min(kk,l-1); 
  else, kk=1; sigma=sigma(1,:); end

%%%------ compute schur form -------------
  [Q,S]=schur(A); %% A*Q=Q*S, Q'*Q=eye(size(A));
%%% transform real schur form to complex schur form
  if norm(tril(S,-1),1)>0, [Q,S]=rsf2csf(Q,S); end

%%%------ find order eigenvalues ---------------
  I = SortEig(diag(S),sigma); 

%%%------ reorder schur form ----------------
  [Q,S] = SwapSchur(Q,S,I(1:kk)); 

return
%----------------------------------------------------------------------
function I=SortEig(t,sigma);
%I=SortEig(T,SIGMA) sorts the indices of T.
%
% T is a vector of scalars, 
% SIGMA is a string or a vector of scalars.
% I is a permutation of (1:LENGTH(T))' such that:
%   if SIGMA is a vector of scalars then
%   for K=1,2,...,LENGTH(T) with KK = MIN(K,SIZE(SIGMA,1))
%      ABS( T(I(K))-SIGMA(KK) ) <= ABS( T(I(J))-SIGMA(KK) ) 
%      SIGMA(kk)=INF: ABS( T(I(K)) ) >= ABS( T(I(J)) ) 
%         for all J >= K

if ischar(sigma)
  switch sigma
    case 'LM'
      [s,I]=sort(-abs(t));
    case 'SM'
      [s,I]=sort(abs(t));
    case 'LR';
      [s,I]=sort(-real(t));
    case 'SR';
      [s,I]=sort(real(t));
    case 'BE';
      [s,I]=sort(real(t)); I=twistdim(I,1);
  end
else

  [s,I]=sort(abs(t-sigma(1,1))); 
  ll=min(size(sigma,1),size(t,1)-1);
  for j=2:ll
    if sigma(j,1)==inf
      [s,J]=sort(abs(t(I(j:end)))); J=flipdim(J,1);
    else
      [s,J]=sort(abs(t(I(j:end))-sigma(j,1)));
    end
    I=[I(1:j-1);I(J+j-1)];
  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,S]=SwapSchur(Q,S,I)
% [Q,S]=SwapSchur(QQ,SS,P)
%    QQ and SS are square matrices of size K by K
%    P is the first part of a permutation of (1:K)'.
%
%    If    M = QQ*SS*QQ'  and  QQ'*QQ = EYE(K), SS upper triangular
%    then  M*Q = Q*S      with   Q'*Q = EYE(K),  S upper triangular
%    and   D(1:LENGTH(P))=DD(P) where D=diag(S), DD=diag(SS)
%
%    Computations uses Givens rotations.

  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
      q = [S(k,k)-S(k+1,k+1),S(k,k+1)];
      if q(1) ~= 0
        q = q/norm(q);
        G = [[q(2);-q(1)],q'];
        J = [k,k+1];
        Q(:,J) = Q(:,J)*G;
        S(:,J) = S(:,J)*G;
        S(J,:) = G'*S(J,:);
      end
      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]=SortQZ(A,B,sigma,gamma,kk)
%
% [Q,Z,S,T]=SORTQZ(A,B,SIGMA)
%   A and B are K by K matrices, SIGMA is a complex scalar or string.
%   SORTQZ computes the qz-decomposition of (A,B) with prescribed
%   ordering: A*Q=Z*S, B*Q=Z*T; 
%             Q and Z are K by K unitary,
%             S and T are K by K upper triangular.
%   The ordering is as follows:
%   (DAIG(S),DAIG(T)) are the eigenpairs of (A,B) ordered
%   as prescribed by SIGMA.
%

% coded by Gerard Sleijpen, version Januari 12, 1998


  l=size(A,1); 
  if l<2; Q=1; Z=1; S=A; T=B; return
  elseif nargin==3, kk=l-1; 
  elseif gamma, kk=min(kk,l-1); 
  else, kk=1; sigma=sigma(1,:); end

%%%------ compute qz form ----------------
  [S,T,Z,Q]=qz(A,B); Z=Z'; S=triu(S);
  
%%%------ sort eigenvalues ---------------

  I=SortEigPair(diag(S),diag(T),sigma); 
 
%%%------ sort qz form -------------------
  [Q,Z,S,T]=SwapQZ(Q,Z,S,T,I(1:kk)); 

return
%----------------------------------------------------------------------
function I=SortEigPair(s,t,sigma)
% I=SortEigPair(S,T,SIGMA)
%   S is a complex K-vectors, T a positive real K-vector
%   SIGMA is a string or a vector of pairs of complex scalars.
%   SortEigPair gives the index set I that sorts the pairs (S,T).
%
%   If SIGMA is a pair of scalars then the sorting is 
%   with increasing "chordal distance" w.r.t. SIGMA.
%
%  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 and F*A(2)>=0,
%  scale B by a scalar G such that NORM(G*B)=1 and G*B(2)>=0,
%  then D(A,B)=ABS((F*A)*RROT(G*B)) where RROT(alpha,beta)=(beta,-alpha)


% coded by Gerard Sleijpen, version Januari 14, 1998

n=sign(t); n=n+(n==0); t=abs(t./n); s=s./n; 

if ischar(sigma)
  switch sigma
    case {'LM','SM'}
    case {'LR','SR','BE'}
      s=real(s);
  end
  [s,I]=sort((-t./sqrt(s.*conj(s)+t.*t)));
  switch sigma
    case {'LM','LR'}
      I=flipdim(I,1);
    case {'SM','SR'}
    case 'BE'
      I=twistdim(I,1);  
  end
else

  n=sqrt(sigma.*conj(sigma)+1); ll=size(sigma,1); 
  tau=[ones(ll,1)./n,-sigma./n]; tau=tau.';

  n=sqrt(s.*conj(s)+t.*t); s=[s./n,t./n];

  [t,I]=sort(abs(s*tau(:,1))); 
  ll = min(ll,size(I,1)-1); 
  for j=2:ll
    [t,J]=sort(abs(s(I(j:end),:)*tau(:,j))); 
    I=[I(1:j-1);I(J+j-1)];
  end 

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