jdqr.m

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

  end

return

%===========================================================================
function CheckSortSchur(sigma)
global Qschur Rschur

k=size(Rschur,1); if k==0, return, end

I=SortEig(diag(Rschur),sigma);

if ~min((1:k)'==I)
   [U,Rschur]=SwapSchur(eye(k),Rschur,I);
   Qschur=Qschur*U;
end

return

%%%=========== COMPUTE SORTED JORDAN FORM ==================================
function [X,Jordan]=Jordan(S)
% [X,J]=JORDAN(S)
%   For S  k by k upper triangular matrix with ordered diagonal elements,
%   JORDAN computes the Jordan decomposition.
%   X is a k by k matrix of vectors spanning invariant spaces.
%   If J(i,i)=J(i+1,i+1) then X(:,i)'*X(:,i+1)=0.
%   J is a k by k matrix, J is Jordan such that S*X=X*J.
%   diag(J)=diag(S) are the eigenvalues

% coded by Gerard Sleijpen, Januari 14, 1998

k=size(S,1); X=zeros(k);  
if k==0, Jordan=[]; return, end

%%% accepted separation between eigenvalues:
delta=2*sqrt(eps)*norm(S,inf); delta=max(delta,10*eps);

T=eye(k); s=diag(S); Jordan=diag(s);
for i=1:k
  I=[1:i]; e=zeros(i,1); e(i,1)=1;
  C=S(I,I)-s(i,1)*T(I,I); C(i,i)=1;
  j=i-1; q=[]; jj=0; 
  while j>0
    if abs(C(j,j))<delta, jj=jj+1; j=j-1; else, j=0; end
  end
  q=X(I,i-jj:i-1);
  C=[C,T(I,I)*q;q',zeros(jj)]; 
  q=C\[e;zeros(jj,1)]; nrm=norm(q(I,1));
  Jordan(i-jj:i-1,i)=-q(i+1:i+jj,1)/nrm;
  X(I,i)=q(I,1)/nrm;
end

return
%========== OUTPUT =========================================================
function varargout=output(history,X,Lambda)

global Qschur Rschur

if nargout == 1, varargout{1}=diag(Rschur);             return, end
if nargout > 2,  varargout{nargout}=history;                    end
if nargout > 3,  varargout{3} = Qschur;  varargout{4} = Rschur; end
if nargout < 4 & size(X,2)<2
  varargout{1}=Qschur; varargout{2}=Rschur; return
else
  varargout{1}=X; varargout{2}=Lambda;
end

return

%===========================================================================
%===== UPDATE PRECONDITIONED SCHUR VECTORS =================================
%===========================================================================
function   solved=UpdateMinv(u,solved)

global Qschur PinvQ Pinv_u Pu L_precond

   if ~isempty(L_precond)
      if ~solved, Pinv_u=SolvePrecond(u); end
      Pu=[[Pu;u'*PinvQ],Qschur'*Pinv_u]; 
      PinvQ=[PinvQ,Pinv_u]; 
      solved=0;
   end

return
%===========================================================================
%===== SOLVE CORRECTION EQUATION ===========================================
%===========================================================================
function t=Solve_pce(theta,u,r,lsolver,par,nit)

global Qschur PinvQ Pinv_u Pu L_precond

  switch lsolver
    case 'exact'
      t = exact(theta,[Qschur,u],r);

    case 'iluexact'
      t = iluexact(theta,[Qschur,u],r);

    case {'gmres','bicgstab','olsen'}
      if isempty(L_precond) %%% no preconditioning

         t = feval(lsolver,theta,...
            [Qschur,u],[Qschur,u],1,r,spar(par,nit));

      else  %%% solve left preconditioned system
               
         %%% compute vectors and matrices for skew projection
         Pinv_u=SolvePrecond(u); 
         mu=[Pu,Qschur'*Pinv_u;u'*PinvQ,u'*Pinv_u];

         %%% precondion and project r
         r=SolvePrecond(r);
         r=SkewProj([Qschur,u],[PinvQ,Pinv_u],mu,r);

         %%% solve preconditioned system
         t = feval(lsolver,theta,...
            [Qschur,u],[PinvQ,Pinv_u],mu,r,spar(par,nit));
      end

    case {'cg','minres','symmlq'}
      if isempty(L_precond) %%% no preconditioning

         t = feval(lsolver,theta,...
            [Qschur,u],[Qschur,u],1,r,spar(par,nit));

      else  %%% solve two-sided expl. precond. system

         %%% compute vectors and matrices for skew projection
         Pinv_u=SolvePrecond(u,'L'); 
         mu=[Pu,Qschur'*Pinv_u;u'*PinvQ,u'*Pinv_u];

         %%% precondion and project r
         r=SolvePrecond(r,'L');
         r=SkewProj([Qschur,u],[PinvQ,Pinv_u],mu,r);

         %%% solve preconditioned system
         t = feval(lsolver,theta,...
            [Qschur,u],[PinvQ,Pinv_u],mu,r,spar(par,nit));

         %%% "unprecondition" solution
         t=SkewProj([PinvQ,Pinv_u],[Qschur,u],mu,t);
         t=SolvePrecond(t,'U');
      end
      
  end
    
return
%=======================================================================
%======= LINEAR SOLVERS ================================================
%=======================================================================
function x = exact(theta,Q,r)

global A_operator

   [n,k]=size(Q); [n,l]=size(r);

   if ischar(A_operator)
      [x,xtol]=bicgstab(theta,Q,Q,1,r,[5.0e-14/norm(r),200,4]);
      return
   end

   x = [A_operator-theta*speye(n,n),Q;Q',zeros(k,k)]\[r;zeros(k,l)];
   x = x(1:n,1:l);  

return
%----------------------------------------------------------------------
function x = iluexact(theta,Q,r)

global L_precond U_precond

   [n,k]=size(Q); [n,l]=size(r);

   y = L_precond\[r,Q]; 
   x = [U_precond,y(:,l+1:k+1);Q',zeros(k,k)]\[y(:,1:l);zeros(k,l)]; 
   x = x(1:n,1:l); 
 
return
%----------------------------------------------------------------------
function r = olsen(theta,Q,Z,M,r,par)
return
%
%======= Iterative methods =============================================
%
function x = bicgstab(theta,Q,Z,M,r,par)
% BiCGstab(ell)
% [x,rnrm] = bicgstab(theta,Q,Z,M,r,par)
% Computes iteratively an approximation to the solution 
% of the linear system Q'*x = 0 and Atilde*x=r 
% where Atilde=(I-Z*M^(-1)*Q')*(U\L\(A-theta)).
%
% This function is specialized for use in JDQZ.
% integer nmv: number of matrix multiplications
% rnrm: relative residual norm
%
%  par=[tol,mxmv,ell] where 
%    integer m: max number of iteration steps
%    real tol: residual reduction
%
% nmv:  number of MV with Atilde
% rnrm: obtained residual reduction
%
% -- References: ETNA

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

% -- Initialization --
%

tol=par(1); max_it=par(2); l=par(3); n=size(r,1);
rnrm=1; nmv=0;
 
if max_it==0 | tol>=1, x=r; return, end

rnrm=norm(r); snrm=rnrm; tol=tol*rnrm;

sigma=1; omega=1; 
x=zeros(n,1); u=zeros(n,1); tr=r;
%   hist=rnrm;
% -- Iteration loop
while (rnrm > tol) & (nmv <= max_it)

  sigma=-omega*sigma;
  for j = 1:l,
    rho=tr'*r(:,j);  bet=rho/sigma;
    u=r-bet*u;
    %%%%%% u(:,j+1)=Atilde*u(:,j)
    u(:,j+1)=mvp(theta,Q,Z,M,u(:,j)); 
    sigma=tr'*u(:,j+1);  alp=rho/sigma;
    x=x+alp*u(:,1);
    r=r-alp*u(:,2:j+1);
    %%%%%% r(:,j+1)=Atilde*r(:,j)
    r(:,j+1)=mvp(theta,Q,Z,M,r(:,j));
  end

  gamma=r(:,2:l+1)\r(:,1); omega=gamma(l,1);
  x=x+r*[gamma;0]; u=u*[1;-gamma]; r=r*[1;-gamma];

  rnrm = norm(r); nmv = nmv+2*l;
    %  hist=[hist,rnrm];
end
    %      figure(3),
    %      plot([0:length(hist)-1]*2*l,log10(hist/snrm),'-*'),
    %      drawnow, 

rnrm = rnrm/snrm;

return
%----------------------------------------------------------------------
function v = gmres(theta,Q,Z,M,v,par)
% GMRES
% [x,rnrm] = gmres(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: Saad

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

% -- Initialization

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

if max_it==0 | tol>=1, return, end
 
H = zeros(max_it +1,max_it); Rot=[ones(1,max_it);zeros(1,max_it)];

rnrm = norm(v); v = v/rnrm;  V = [v];
tol = tol * rnrm; snrm = rnrm;
y = [ rnrm ; zeros(max_it,1) ];
j=0;          %  hist=rnrm;
while (nmv < max_it) & (rnrm > tol),
  j=j+1; nmv=nmv+1;
  v=mvp(theta,Q,Z,M,v); 
  % [v, H(1:j+1,j)] = RepGS(V,v);                      
  [v,h] = RepGS(V,v); H(1:size(h,1),j)=h; 
  V = [V, v]; 
  for i = 1:j-1,
    a = Rot(:,i);
    H(i:i+1,j) = [a'; -a(2) a(1)]*H(i:i+1,j);
  end
  J=[j, j+1];
  a=H(J,j);
  if a(2) ~= 0
    cs = norm(a); 
    a = a/cs; Rot(:,j) = a;
    H(J,j) = [cs; 0];
    y(J) = [a'; -a(2) a(1)]*y(J);
  end 
  rnrm = abs(y(j+1)); 
             %  hist=[hist,rnrm];
end
             %    figure(3)
             %    plot([0:length(hist)-1],log10(hist/snrm),'-*')
             %    drawnow, pause

J=[1:j];
v = V(:,J)*(H(J,J)\y(J));
rnrm = rnrm/snrm;

return
%======================================================================
%========== BASIC OPERATIONS ==========================================
%======================================================================
function v=MV(v)

global A_operator nm_operations

if ischar(A_operator)
  v = feval(A_operator,v);
else
  v = A_operator*v;
end

nm_operations = nm_operations+1;

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

   v = MV(v) - theta*v;
   v = SolvePrecond(v);
   v = SkewProj(Q,Z,M,v);
  
return
%----------------------------------------------------------------------
function u=SolvePrecond(u,flag);

global L_precond U_precond

if isempty(L_precond), return, end

if nargin<2

   if ischar(L_precond)
      if ischar(U_precond)
         u=feval(L_precond,u,U_precond);  
      elseif isempty(U_precond)
         u=feval(L_precond,u); 
      else
         u=feval(L_precond,u,'L'); u=feval(L_precond,u,'U'); 
      end
   else
      u=U_precond\(L_precond\u); 
   end

else

   switch flag
     case 'U'
       if ischar(L_precond), u=feval(L_precond,u,'U'); else, u=U_precond\u; end

     case 'L'
       if ischar(L_precond), u=feval(L_precond,u,'L'); else, u=L_precond\u; end

   end
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)
% Changes par=[tol(:),max_it,ell]  to
% ppap=[TOL,max_it,ell] where 
% if lenght(tol)==1
%   TOL=tol
% else
%   red=tol(end)/told(end-1); tole=tol(end);
%   tol=[tol,red*tole,red^2*tole,red^3*tole,...]
%   TOL=tol(nit);
% end

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
%
%======= Iterative methods for symmetric systems =======================
%
function x = cg(theta,Q,Z,M,r,par)
% CG
% [x,rnrm] = cg(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: Hestenes and Stiefel

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

           b=r;

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