jdqz.m

数据降维工具箱

    q(1,1)=q(1,1)+1; q=q/norm(q); q=[zeros(i-1,1);q]; 
    z(1,1)=z(1,1)+1; z=z/norm(z); z=[zeros(i-1,1);z];
    S=S-(S*q)*(2*q)'; S=S-(2*z)*(z'*S);
    T=T-(T*q)*(2*q)'; T=T-(2*z)*(z'*T);
    Q=Q-(Q*q)*(2*q)'; Z=Z-(Z*z)*(2*z)';
  end
return

%%%======== END sort QZ decomposition interaction matrices ==============

%%%======================================================================
%%%======== INITIALIZATION ==============================================
%%%======================================================================
function MyClear

global Operator_Form Operator_A Operator_B Operator_Params ...
       Precond_L Precond_U Precond_P Precond_Params ...
       Precond_Form Precond_Type ...
       Operator_MVs Precond_Solves ...
       CHORDALDISTANCE ...
       Qschur Zschur Sschur Tschur ...
       MinvZ QastMinvZ

return
%%%======================================================================
function [n,nselect,Sigma,kappa,SCHUR,...
          jmin,jmax,tol,maxit,V,AV,BV,TS,DISP,PAIRS,JDV,FIX,track,NSIGMA,...
          lsolver,par] = ReadOptions(varargin)
% Read options and set defaults

global Operator_Form Operator_A Operator_B Operator_Params ...
       Precond_Form Precond_L Precond_U Precond_P Precond_Params ...
       CHORDALDISTANCE


Operator_A = varargin{1};

n=CheckMatrix(Operator_A,1);

% defaults              %%%% search for 'xx' in fieldnames
nselect0= 5; 
maxit   = 200;          %%%% 'ma'
SCHUR   = 0;            %%%% 'sch'
tol     = 1e-8;         %%%% 'to'
DISP    = 0;            %%%% 'di'
p0      = 5; %%% jmin=nselect+p0 %%%% 'jmi'
p1      = 5; %%% jmax=jmin+p1    %%%% 'jma'
TS      = 1;            %%%% 'te'
PAIRS   = 0;            %%%% 'pai'
JDV     = 0;            %%%% 'av'
track   = 1e-4;         %%%% 'tr'
FIX     = 1000;         %%%% 'fix'
NSIGMA  = 0;            %%%% 'ns'
CHORD   = 1;            %%%% 'ch'
lsolver = 'gmres';      %%%% 'lso'
ls_maxit= 200;          %%%% 'ls_m'
ls_tol  = [0.7,0.49];   %%%% 'ls_t' 
ell     = 4;            %%%% 'ls_e'
TP      = 0;            %%%% 'ty'
L       = [];           %%%% 'l_'
U       = [];           %%%% 'u_'
P       = [];           %%%% 'p_'
kappa   = 1;            %%%% 'sca'
V0      = 'ones(n,1)+rand(n,1)';  %%%% 'v0'

%% initiation
nselect=[]; Sigma=[]; options=[]; Operator_B=[];
jmin=-1; jmax=-1; V=[]; AV=[]; BV=[]; par=[];

%------------------------------------------------
%------- Find quantities ------------------------
%------------------------------------------------
jj=[];
for j = 2:nargin
   if isstruct(varargin{j})
      options = varargin{j};
   elseif ischar(varargin{j})
      s=varargin{j}; 
      if exist(s)==2 & isempty(Operator_B)
         Operator_B=s;
      elseif  length(s) == 2 & isempty(Sigma)
        s=upper(s);
        switch s
          case {'LM','SM','LR','SR','BE'}, Sigma=s;
          otherwise
            jj=[jj,j];
        end
      else
        jj=[jj,j];
      end
   elseif min([n,n]==size(varargin{j})) & isempty(Operator_B)
      Operator_B=varargin{j}; 
   elseif length(varargin{j}) == 1
      s = varargin{j};
      if isempty(nselect) & isreal(s) & (s == fix(s)) & (s > 0)
         nselect = min(n,s);
      elseif isempty(Sigma)
         Sigma = s; 
      else
        jj=[jj,j];
      end
   elseif min(size(varargin{j}))==1  & isempty(Sigma)
      Sigma = varargin{j}; if size(Sigma,1)==1, Sigma=Sigma'; end 
   elseif min(size(varargin{j}))==2  & isempty(Sigma)
      Sigma = varargin{j}; if size(Sigma,2)>2 , Sigma=Sigma'; end 
   else
      jj=[jj,j];
   end
end

%------- find parameters for operators -----------
Operator_Params=[]; Operator_Params{2}='';
k=length(jj);
if k>0
  Operator_Params(3:k+2)=varargin(jj);
  if ~ischar(Operator_A)
    msg=sprintf(', %i',jj);
    msg=sprintf('Input argument, number%s, not recognized.',msg);
    button=questdlg(msg,'Input arguments','Ignore','Stop','Ignore');
    if strcmp(button,'Stop'), n=-1; return, end
  end
end
%------- operator B -----------------------------
if isempty(Operator_B)
   if ischar(Operator_A)
      Operator_Form=2; % or Operator_Form=3, or Operator_Form=5;
   else
      Operator_Form=8;
   end
else
   if ischar(Operator_B)
      if ischar(Operator_A), Operator_Form=1; else, Operator_Form=6; end
   elseif ischar(Operator_A)
      Operator_Form=4;
   else
      Operator_Form=7;
   end
end
if n<2, return, end

%------- number of eigs to be computed ----------
if isempty(nselect), nselect=min(n,nselect0); end

%------------------------------------------------
%------- Analyse Options ------------------------
%------------------------------------------------
fopts = [];
if ~isempty(options), fopts = fieldnames(options); end

%------- preconditioner -------------------------
Precond_L=findfield(options,fopts,'pr',[]);
[L,ok]=findfield(options,fopts,'l_',Precond_L);
if ok & ~isempty(Precond_L),
   msg =sprintf('A preconditioner is defined in');
   msg =[msg,sprintf('\n''Precond'', but also in ''L_precond''.')];
   msg=[msg,sprintf('\nWhat is the correct one?')];
   button=questdlg(msg,'Preconditioner','L_Precond','Precond','L_Precond');
   if strcmp(button,'L_Precond'), 
     Precond_L = L;
   end
else
   Precond_L = L;
end

if ~isempty(Precond_L)
  Precond_U=findfield(options,fopts,'u_',[]);
  Precond_P=findfield(options,fopts,'p_',[]);
end

Precond_Params=[]; Precond_Params{2}='';
Params=findfield(options,fopts,'par',[]);
[l,k]=size(Params);
if k>0, 
  if iscell(Params), Precond_Params(3:k+2)=Params; 
  else, Precond_Params{3}=Params; end
end
TP=findfield(options,fopts,'ty',TP);
n=SetPrecond(n,TP); if n<2, return, end

%------- max, min dimension search subspace ------
jmin=min(n,findfield(options,fopts,'jmi',jmin));
jmax=min(n,findfield(options,fopts,'jma',jmax));
if jmax < 0
   if jmin<0, jmin=min(n,nselect+p0); end
   jmax=min(n,jmin+p1); 
else
   if jmin<0, jmin=max(1,jmax-p1); end
end 
maxit=findfield(options,fopts,'ma',maxit);

%------- initial search subspace ----------------
V=findfield(options,fopts,'v',[]);
[m,d]=size(V); 
if m~=n 
  if m>n, V = V(1:n,:); end
  if m<n, V = [V;0.001*rand(n-m-1,d)]; end
end
V=orth(V); [m,d]=size(V);
if d==0, nr=0; while nr==0, V=eval(V0); nr=norm(V); V=V/nr; end, end


%------- Check definition B, Compute AV, BV -----
[AV,BV,n]=CheckDimMV(V); if n<2, return, end

%------- Other options --------------------------
tol=findfield(options,fopts,'to',tol);
kappa  = findfield(options,fopts,'sca',kappa);
kappa  = abs(kappa(1,1)); if kappa==0, kappa=1; end
PAIRS  = findfield(options,fopts,'pai',PAIRS,[0,1]);
SCHUR  = findfield(options,fopts,'sch',SCHUR,[0,1,0.5]);
DISP   = findfield(options,fopts,'di',DISP,[0,1]);
JDV    = findfield(options,fopts,'av',JDV,[0,1]);
track  = max(abs(findfield(options,fopts,'tr',track,[0,track,inf])),0);
NSIGMA = findfield(options,fopts,'ns',NSIGMA,[0,1]);
FIX    = max(abs(findfield(options,fopts,'fix',0,[0,FIX,inf])),0);
CHORDALDISTANCE = findfield(options,fopts,'ch',CHORD,[0,1]);

[TS0,ok] = findfield(options,fopts,'te',TS);

if ok & ischar(TS0)
   if     strncmpi(TS0,'st',2),  TS=0; %% 'standard'
   elseif strncmpi(TS0,'ha',2),  TS=1; %% 'harmonic'
   elseif strncmpi(TS0,'se',2),  TS=2; %% 'searchspace'
   elseif strncmpi(TS0,'bv',2),  TS=3;
   elseif strncmpi(TS0,'av',2),  TS=4;
   end
else
   TS=max(0,min(4,round(TS0(1,1))));
end

%------- set targets ----------------
if isempty(Sigma)
   if TS==1, Sigma=[0,1]; else, Sigma = 'LM'; end
elseif ischar(Sigma)
   switch Sigma
     case {'LM','LR','SR','BE','SM'}
     if ~ok, TS=3; end
   end
else
  [k,l]=size(Sigma);
  if l==1, Sigma=[Sigma,ones(k,1)]; l=2; end
  Sigma=ScaleEig(Sigma);
end

if ischar(Sigma) & TS<2 
   msg1=sprintf('   The choice sigma = ''%s'' does not match the\n',Sigma);
   msg2=sprintf('   selected test subspace. Specify a numerical\n');
   msg3=sprintf('   value for sigma (e.g. sigma = '); msg4='';
   switch Sigma
      case {'LM','LR'}
         msg4=sprintf('[1,0]');     
      case {'SM','SR','BE'}
         msg4=sprintf(' [0,1]');  
   end
   msg5=sprintf('),\n   or continue with ''TestSpace''=''B*V''.');
   msg=[msg1,msg2,msg3,msg4,msg5];
   button=questdlg(msg,'Targets and test subspaces','Continue','Stop','Continue');
   if strcmp(button,'Continue'), TS=3; else, n=-1; return, end
end

%------- linear solver --------------------------
lsolver  = findfield(options,fopts,'lso',lsolver);

method=strvcat('exact','olsen','iluexact','gmres','cgstab','bicgstab');
j=strmatch(lower(lsolver),method);
if isempty(j), 
  msg=['The linear solver ''',lsolver,''' is not included.'];
  msg=[msg,sprintf('\nIs GMRES ok?')];
  button=questdlg(msg,'Linear solver','Yes','No','Yes');
  if strcmp(button,'Yes'), j=4; ls_maxit=5; else, n=-1; return, end
end
if     j==1, lsolver='exact'; Precond_Form = 0;
elseif j==2 | j==3, lsolver='olsen';
elseif j==4, lsolver='gmres';  ls_maxit=5;
else,        lsolver='cgstab'; ls_tol=1.0e-10;
end

ls_maxit= findfield(options,fopts,'ls_m',ls_maxit);
ls_tol  = findfield(options,fopts,'ls_t',ls_tol);
ell     = findfield(options,fopts,'ls_e',ell);

par=[ls_tol,ls_maxit,ell];


%----- Display the parameters that are used ----------------
if DISP 

  fprintf('\n'),fprintf('PROBLEM\n')
  switch Operator_Form
     case {1,4,6,7} 
  fprintf('%13s: %s\n','A',StrOp(Operator_A));
  fprintf('%13s: %s\n','B',StrOp(Operator_B));
     case 2
  fprintf('%13s: ''%s''  ([Av,Bv] = %s(v))\n','A,B',Operator_A,Operator_A);
     case 3
  fprintf('%13s: ''%s''\n','A,B',Operator_A);
  fprintf('%15s(Av = %s(v,''A''), Bv = %s(v,''B''))\n',...
                                '',Operator_A,Operator_A);
     case {5,8}
  fprintf('%13s: %s\n','A',StrOp(Operator_A));
  fprintf('%13s: %s\n','B','Identity   (B*v = v)');
  end

  fprintf('%13s: %i\n','dimension',n);
  fprintf('%13s: %i\n\n','nselect',nselect);

  if length(jj)>0 & (ischar(Operator_A) | ischar(Operator_B))
    msgj=sprintf(', %i',jj); msgo=''; 
    if ischar(Operator_A), msgo=sprintf(' ''%s''',Operator_A); end
    if ischar(Operator_B), msgo=sprintf('%s ''%s''.',msgo,Operator_B); end
    fprintf(' The JDQZ input arguments, number%s, are\n',msgj)
    fprintf(' taken as input parameters 3:%i for%s.\n\n',length(jj)+2,msgo);
  end

  fprintf('TARGET\n')
  if ischar(Sigma)
     fprintf('%13s: ''%s''\n','sigma',Sigma)
  else 
     fprintf('%13s: %s\n','sigma',mydisp(Sigma(1,:)))
     for j=2:size(Sigma,1), 
        fprintf('%13s: %s\n','',mydisp(Sigma(j,:)))
     end
  end

  fprintf('\nOPTIONS\n')
  fprintf('%13s: %g\n','Schur',SCHUR)
  fprintf('%13s: %g\n','Tol',tol)
  fprintf('%13s: %i\n','Disp',DISP)
  fprintf('%13s: %i\n','jmin',jmin)
  fprintf('%13s: %i\n','jmax',jmax)
  fprintf('%13s: %i\n','MaxIt',maxit)
  fprintf('%13s: %s\n','v0',StrOp(V))
  fprintf('%13s: %i\n','Pairs',PAIRS)
  fprintf('%13s: %i\n','AvoidStag',JDV)
  fprintf('%13s: %i\n','NSigma',NSIGMA)
  fprintf('%13s: %g\n','Track',track)
  fprintf('%13s: %g\n','FixShift',FIX)
  fprintf('%13s: %i\n','Chord',CHORDALDISTANCE)
  fprintf('%13s: ''%s''\n','LSolver',lsolver)
  str=sprintf('%g ',ls_tol);
  fprintf('%13s: [ %s]\n','LS_Tol',str)
  fprintf('%13s: %i\n','LS_MaxIt',ls_maxit)
  if strcmp(lsolver,'cgstab')
    fprintf('%13s: %i\n','LS_ell',ell)
  end

  DisplayPreconditioner(n); fprintf('\n')
   
  switch TS
     case 0, str='Standard, W = alpha''*A*V + beta''*B*V';
     case 1, str='Harmonic, W = beta*A*V - alpha*B*V';
     case 2, str='SearchSpace, W = V';
     case 3, str='Petrov, W = B*V';
     case 4, str='Petrov, W = A*V';
   end
   fprintf('%13s: %s\n','TestSpace',str); fprintf('\n\n');

end

return
%------------------------------------------------------------------------
function msg=StrOp(Op)

  if ischar(Op)
    msg=sprintf('''%s''',Op);
  elseif issparse(Op), [n,k]=size(Op);
    msg=sprintf('[%ix%i sparse]',n,k);
  else, [n,k]=size(Op);
    msg=sprintf('[%ix%i double]',n,k);
  end

return
%------------------------------------------------------------------------
function DisplayPreconditioner(n)
  global Precond_Form Precond_Type ...
         Precond_L Precond_U Precond_P

  FP=Precond_Form;
  switch Precond_Form
    case 0,    
      fprintf('%13s: %s\n','Precond','No preconditioner'); 
    case {1,2,5}
      fprintf('%13s: %s','Precond',StrOp(Precond_L));
      if FP==2, fprintf('  (M\\v = %s(v,''preconditioner''))',Precond_L); end 
      fprintf('\n')
    case {3,4,6,7} 
      fprintf('%13s: %s\n','L precond',StrOp(Precond_L));
      fprintf('%13s: %s\n','U precond',StrOp(Precond_U));
      if FP==7, fprintf('%13s: %s\n','P precond',StrOp(Precond_P)); end
      if FP==4,fprintf('%15s(M\\v = %s(%s(v,''L''),''U''))\n',...
          '',Precond_L,Precond_L); end 
  end

  if FP
  switch Precond_Type
    case 0, str='explicit left';
    case 1, str='explicit right';
    case 2, str='implicit';
  end
  fprintf('%15sTo be used as %s preconditioner.\n','',str)
  end

return
%------------------------------------------------------------------------
function possibilities

 fprintf('\n')
 fprintf('PROBLEM\n')
 fprintf('            A: [ square matrix | string ]\n');
 fprintf('            B: [ square matrix {identity} | string ]\n');
 fprintf('      nselect: [ positive integer {5} ]\n\n');
 fprintf('TARGET\n')
 fprintf('        sigma: [ vector of scalars | \n');
 fprintf('                 pair of vectors o