jdqr.m

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

     j=jmin; J=[1:j]; 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 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 EXPAND
     %%% Expand the subspaces of the interaction matrix  
     v=RepGS([Qschur,V],v); 
     if size(v,2)>0
       w=MV(v); if tau ~=0, w=w-tau*v;end
       AV=[AV,w]; R=[R,W'*w];
       w=RepGS([Qschur,W],w);
       R=[R;w'*AV]; M=[M,W'*v;w'*V,w'*v]; 
       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.2

%%% The JD loop (Harmonic Ritz values)
%%%    W orthonormal, V and W orthogonal to Qschur, 
%%%    W'*W=eye(j), Qschur'*V=0, Qschur'*W=0
%%%    W=(A*V-tau*V)-Qschur*E, E=Qschur'*(A*V-tau*V), 
%%%    M=W'*V
V=V/R; M=M/R; temptarget='LM';            E=E/R;
while (k<nselect) & (nit<maxit) 

   %%% Compute approximate eigenpair and residual
   [UR,S]=SortSchur(M,temptarget,j==jmax,jmin);
   y=UR(:,1); u=V*y; nrm=norm(u); y=y/nrm; u=u/nrm;
   theta=S(1,1)'/(nrm*nrm); w=W*y; r=w-theta*u; nr=norm(r); r_KNOWN=1;
   if nr<t_tol, temptarget=[S(1,1);inf]; 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=1./diag(S)+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;
     y=E*y; 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
       V=V/R; j=size(V,2);  M=M/R;         E=E/R;
     else
       J=[2:j]; j=j-1; UR=UR(:,J); M=S(J,J);
       V=V*UR; W=W*UR; [r,a]=RepGS(u,r,0); 
       s=u'*V; V=V-u*s; W=W-r*s; M=M-s'*(r'*V)-(W'*u)*s;     
                                 E=[E*UR-y*s;(tau-theta-a)*s];
         
       if (nr*norm(s))^2>eps, [W,R]=qr(W,0); V=V/R; M=(R'\M)/R; E=E/R; 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=[1/(conj(theta)-tau);inf];
     else, nlit=0; temptarget='LM';
       if nt<n_tar
         nt=nt+1; tau0=tau; tau=sigma(nt,1); 
         [W,R]=qr(W+(tau0-tau)*V,0); V=V/R; M=W'*V; E=E/R; 
       end
     end

   end

   %%% 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;               E=E*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 EXPAND
     %%% Expand the subspaces of the interaction matrix  
     v=RepGS(Qschur,v,0);
     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); 
       nrw=y(j+1,1); y=y(1:j,:);   
       v=v-V*y; v=v/nrw;                     e=e-E*y; e=e/nrw;
       M=[M,W'*v;w'*V,w'*v]; 
       V=[V,v]; W=[W,w]; j=j+1;              E=[E,e];

       if 1/cond(M)<10*tol
         [V,W,R,E,M]=SetInitialSpaces(V,nselect,tau,jmin,tol,W,E);    
         k=size(Rschur,1); if k>=nselect, break, end
         V=V/R; M=M/R; j=size(V,2);temptarget='LM'; E=E/R;
       end 

       EXPAND=0; tol=tol0;
     else
       tol=2*tol;
     end
   end 

end % while (nit<maxit)


end % case


time_needed=etime(clock,time);

Refine([Qschur,V],1);% 2-SCHUR);
CheckSortSchur(sigma);

Lambda=[]; X=zeros(n,0); 
if ~SCHUR & k>0, [z,Lambda]=Jordan(Rschur); X=Qschur*z; end

%-------------- display results ----------------------------
if SHOW == 2, MovieTheta, figure(FIG), end
if SHOW & size(history,1)>0

   switch INTERIOR
      case 0
        testspace='V, V orthonormal';     
      case 1
        testspace='A*V-sigma*V, V and W orthonormal';    
      case 1.1
        testspace='A*V-sigma*V, V and W orthonormal, AV';   
      case 1.2
        testspace='A*V-sigma*V, W orthogonal';     
      otherwise 
        testspace='Experimental'; 
   end

   StringT=sprintf('The test subspace W is computed as  W = %s.',testspace);
   StringX=sprintf('JDQZ with jmin=%g, jmax=%g, residual tolerance %g.',...
           jmin,jmax,tol); 
   StringY=sprintf('Correction equation solved with %s.',lsolver);   
   date=fix(clock);
   String=sprintf('\n%2i-%2i-%2i, %2i:%2i:%2i',date(3:-1:1),date(4:6));

   StringL='log_{10} || r_{#it} ||_2';
   for pl=1:SHOW
     subplot(SHOW,1,pl), t=history(:,pl+1);
     plot(t,log10(history(:,1)),'*-',t,log10(tol)+0*t,':')
     legend(StringL), title(StringT)
     StringL='log_{10} || r_{#MV} ||_2'; StringT=StringX; 
   end 
   if SHOW==2, xlabel([StringY,String])
   else,  xlabel([StringX,String]), ylabel(StringY), end
   drawnow
end

if SHOW
   str1=num2str(abs(k-nselect)); str='s';
   if k>nselect, 
     if k==nselect+1, str1='one'; str=''; end
     fprintf('\n\nDetected %s additional eigenpair%s.',str1,str)
   end
   if k<nselect, 
     if k==0, str1='any'; str=''; elseif k==nselect-1, str1='one'; str=''; end
     fprintf('\n\nFailed detection of %s eigenpair%s.',str1,str)
   end
   if k>0, ShowLambda(diag(Rschur)); else, fprintf('\n'); end

   Str='time_needed';                              DispResult(Str,eval(Str))
   if (k>0)
      if ~SCHUR
        Str='norm(MV(X)-X*Lambda)';                DispResult(Str,eval(Str))
      end
      Str='norm(MV(Qschur)-Qschur*Rschur)';        DispResult(Str,eval(Str))
      I=eye(k); Str='norm(Qschur''*Qschur-I)';     DispResult(Str,eval(Str))  
   end
   fprintf('\n\n')

end

if nargout == 0, if ~SHOW, eigenvalues=diag(Rschur), end, return, end
[varargout{1:nargout}]=output(history,X,Lambda);

return

%===========================================================================
%======= PREPROCESSING =====================================================
%===========================================================================

%======= INITIALIZE SUBSPACE ===============================================
function [V,W,R,E,M]=SetInitialSpaces(V,nselect,tau,jmin,tol,W,E);
%[V,W,R,E,M]=SetInitialSpaces(VV,nselect,tau,jmin,tol);
%  Output: V(:,1:SIZE(VV,2))=ORTH(VV),
%          V'*V=W'*W=EYE(JMIN), M=W'*V;
%          such that A*V-tau*V=W*R+Qschur*E, 
%          with R upper triangular, and E=Qschur'*(A*V-tau*V).
%
%[V,W,R,E,M]=SetInitialSpaces(VV,nselect,tau,jmin,tol,AV,EE);
%  Input such that
%  A*VV-tau*VV=AV+Qschur*EE, EE=Qschur'*(A*VV-tau*VV);
%
%  Output: V(:,1:SIZE(VV,2))=ORTH(VV),
%          V'*V=W'*W=EYE(JMIN), M=W'*V;
%          such that A*V-tau*V=W*R+Qschur*E, 
%          with R upper triangular, and E=Qschur'*(A*V-tau*V).

global Qschur Rschur

[n,j]=size(V); k=size(Qschur,2);

if j>1, 
  [V,R]=qr(V,0);
  if nargin <6,   
    W=MV(V); R=eye(j); if tau~=0, W=W-tau*V; end
    if k>0, E=Qschur'*W; W=W-Qschur*E; else, E=zeros(0,j); end
  end
  [V,W,R,E,M]=CheckForNullSpace(V,nselect,tau,tol,W,E,R);
  l=size(Qschur,2); j=size(V,2);
  if l>=nselect, if size(V,2)==0; R=1; M=1; return, end, end
  if l>k, UpdateMinv(Qschur(:,k+1:l),0); end, k=l;
end
if j==0, nr=0;
  while nr==0
    V = ones(n,1)+0.1*rand(n,1); V=RepGS(Qschur,V); nr=norm(V);
  end, j=1;
end
if j==1
  [V,H,E]=Arnoldi(V,tau,jmin,nselect,tol); 
  l=size(Qschur,2); j=max(size(H,2),1);
  if l>=nselect, W=V; R=eye(j); M=R; return, end
  if l>k, UpdateMinv(Qschur(:,k+1:l),0); end
  [Q,R]=qr(full(H),0);
  W=V*Q; V(:,j+1)=[]; M=Q(1:j,:)';
  %% W=V*Q; V=V(:,1:j)/R; E=E/R; R=eye(j); M=Q(1:j,:)'/R;
  %% W=V*H; V(:,j+1)=[];R=R'*R;   M=H(1:j,:)';
end

return
%%%======== ARNOLDI (for initializing spaces) ===============================
function [V,H,E]=Arnoldi(v,tau,jmin,nselect,tol)
%
%[V,AV,H,nMV,tau]=ARNOLDI(A,V0,TAU,JMIN,NSELECT,TOL)
%    ARNOLDI computes the Arnoldi factorization of dimenison JMIN+1:
%    (A-tau)*V(:,1:JMIN)=V*H where V is n by JMIN+1 orthonormal with
%    first column a multiple of V0, and H is JMIN+1 by JMIN Hessenberg.
%
%    If an eigenvalue if H(1:j,1:j) is an eigenvalue of A
%    within the required tolerance TOL then the Schurform
%    A*Qschur=Qschur*Rschur is expanded and the Arnoldi factorization
%    (A-tau)*V(:,1:j)=V(:,1:j+1)*H(1:j+1,1:j) is deflated.
%    Returns if size(Qschur,2) = NSELECT or size(V,2) = JMIN+1

%    (A-tau)*V(:,1:JMIN)=V*H+Qschur*E, Qschur'*V=0

%   Coded November 5, 1998, G. Sleijpen

global Qschur Rschur

k=size(Qschur,2); [n,j]=size(v);

if ischar(tau), tau=0; end

H=zeros(1,0); V=zeros(n,0); E=[];
j=0; nr=norm(v);

while j<jmin & k<nselect & j+k<n
   if nr>=tol
      v=v/nr; V=[V,v]; j=j+1;
      Av=MV(v);
   end
   if j==0 
      H=zeros(1,0); j=1;
      nr=0; while nr==0, v=RepGS(Qschur,rand(n,1)); nr=norm(v); end
      v=v/nr; V=v; Av=MV(v);
   end
   if tau~=0; Av=Av-tau*v; end, [v,e] = RepGS(Qschur,Av,0);
   if k==0, E=zeros(0,j); else, E = [E,e(1:k,1)]; end
   [v,y] = RepGS(V,v,0);        H = [H,y(1:j,1)];
   nr = norm(v);                H = [H;zeros(1,j-1),nr];
   [Q,U,H1] = DeflateHess(full(H),tol); 
   j=size(U,2); l=size(Q,2);
   if l>0 %--- expand Schur form ------
      Qschur=[Qschur,V*Q]; 
      Rschur=[Rschur,E*Q; zeros(l,k),H1(1:l,1:l)+tau*eye(l)]; k=k+l;
      E=[E*U;H1(1:l,l+1:l+j)];
      if j>0, V=V*U; H=H1(l+1:l+j+1,l+1:l+j);
      else, V=zeros(n,0); H=zeros(1,0); end
   end
end % while

if nr>=tol
   v=v/nr; V=[V,v]; 
end

return
%----------------------------------------------------------------------
function [Q,U,H]=DeflateHess(H,tol)
% H_in*[Q,U]=[Q,U]*H_out such that H_out(K,K) upper triangular
% where K=1:SIZE(Q,2) and ABS(Q(end,2)*H_in(j+1,j))<TOL, j=SIZE(H,2),

[j1,j]=size(H); 
if j1==j, [Q,H]=schur(H); U=zeros(j,0); return, end

nr=H(j+1,j);
U=eye(j); i=1; J=i:j;
for l=1:j
   [X,Lambda]=eig(H(J,J));
   I=find(abs(X(size(X,1),:)*nr)<tol);
   if isempty(I), break, end
   q=X(:,I(1)); q=q/norm(q); 
   q(1,1)=q(1,1)+sign(q(1,1)); q=q/norm(q);
   H(:,J)=H(:,J)-(2*H(:,J)*q)*q'; 
   H(J,:)=H(J,:)-2*q*(q'*H(J,:));
   U(:,J)=U(:,J)-(2*U(:,J)*q)*q'; 
   i=i+1; J=i:j;
end

[Q,HH]=RestoreHess(H(i:j+1,J)); 
H(:,J)=H(:,J)*Q; H(J,:)=Q'*H(J,:);
U(:,J)=U(:,J)*Q; 

Q=U(:,1:i-1); U=U(:,i:j);

return
%----------------------------------------------------------------------
function [Q,M]=RestoreHess(M)

[j1,j2]=size(M); Q=eye(j2);

for j=j1:-1:2
   J=1:j-1;
   q=M(j,J)'; q=q/norm(q);
   q(j-1,1)=q(j-1,1)+sign(q(j-1,1));
   q=q/norm(q);
   M(:,J)=M(:,J)-2*(M(:,J)*q)*q';
   M(J,:)=M(J,:)-2*q*q'*M(J,:);
   Q(:,J)=Q(:,J)-2*Q(:,J)*q*q';
end

return
%%%=========== END ARNOLDI ============================================  
function [V,W,R,E,M]=CheckForNullSpace(V,nselect,tau,tol,W,E,Rv);
% V,W orthonormal, A*V-tau*V=W*R+Qschur'*E

global Qschur Rschur

  k=size(Rschur,1); j=size(V,2); 
 
  [W,R]=qr(W,0); E=E/Rv; R=R/Rv; M=W'*V; 
  %%% not accurate enough M=Rw'\(M/Rv);

  if k>=nselect, return, end

  CHECK=1; l=k;

  [S,T,Z,Q]=qz(R,M); Z=Z';
  while CHECK 
    I=SortEigPairVar(S,T,2); [Q,Z,S,T]=SwapQZ(Q,Z,S,T,I(1)); 
    s=abs(S(1,1)); t=min(abs(T(1,1)),1); CHECK=(s*sqrt(1-t*t)<tol);
    if CHECK 
      V=V*Q; W=W*Z; E=E*Q; 
      u=V(:,1); [r,a]=RepGS(u,(W(:,1)-T(1,1)'*u)*S(1,1),0);
      Qschur=[Qschur,u]; t=(T(1,1)'-a/S(1,1))*S(1,:);
      Rschur=[Rschur,E(:,1);zeros(1,k),tau+t(1,1)]; k=k+1;
      J=[2:j]; j=j-1; 
      V=V(:,J); W=W(:,J); E=[E(:,J);t(1,J)];  s=S(1,J)/S(1,1);
      R=S(J,J); M=T(J,J); Q=eye(j); Z=eye(j); 

      s=s/R; nrs=norm(r)*norm(s);
      if nrs>=tol,
        W=W+r*s; M=M+s'*(r'*V); 
        if nrs^2>eps, [W,R0]=qr(W,0); R=R0*R; M=R0'\M; end
      end

      S=R; T=M;

      CHECK=(k<nselect & j>0);
    end
  end

return


%===========================================================================
%======= POSTPROCESSING ====================================================
%===========================================================================
function Refine(V,gamma);

if gamma==0, return, end

global Qschur Rschur

  J=1:size(Rschur,1); 

  if gamma==1, 
    [V,R]=qr(V(:,J),0); W=MV(V); M=V'*W;
    [U,Rschur]=schur(M);
    [U,Rschur]=rsf2csf(U,Rschur); Qschur=V*U;
    return 
  elseif gamma==2
    [V,R]=qr(V,0); W=MV(V); M=V'*W;
    [U,S]=schur(M); [U,S]=rsf2csf(U,S); 
    R=R*U; F=R'*R-S'*S;
    [X,Lambda]=Jordan(S); 
    % Xinv=inv(X); D=sqrt(diag(Xinv*Xinv')); X=X*diag(D);
    [d,I]=sort(abs(diag(X'*F*X)));
    [U,S]=SwapSchur(U,S,I(J));
    Qschur=V*U(:,J); Rschur=S(J,J);