jdqr.m

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

return

%=======================================================================
%======= SET PARAMETERS ================================================
%=======================================================================
function [n,nselect,sigma,SCHUR,...
         jmin,jmax,tol,maxit,V,INTERIOR,SHOW,PAIRS,JDV,OLD,...
         lsolver,par] = ReadOptions(varargin)
% Read options and set defaults

global A_operator L_precond U_precond

A_operator = varargin{1};  

%%% determine dimension
if ischar(A_operator)
  n=-1;
  if exist(A_operator) ~=2
    msg=sprintf('  Can not find the M-file ''%s.m''  ',A_operator);
    errordlg(msg,'MATRIX'),n=-2;
  end
  if n==-1, eval('n=feval(A_operator,[],''dimension'');','n=-1;'), end
else
  [n,n] = size(A_operator);
  if any(size(A_operator) ~= n)
    msg=sprintf('  The operator must be a square matrix or a string.  ');
    errordlg(msg,'MATRIX'),n=-3;
  end
end

%%% defaults
SCHUR   = 0;
jmin    = -1;
jmax    = -1;
p0      = 5; % jmin=nselect+p0
p1      = 5; % jmax=jmin+p1
tol     = 1e-8; 
maxit   = 200;
V       = zeros(0,0);
INTERIOR= 0;
SHOW    = 0;
PAIRS   = 0;
JDV     = 0;
OLD     = 1e-4;
lsolver = 'gmres';
ls_maxit= 200; 
ls_tol  = [1,0.7];  
ell     = 4;
par     = [ls_tol,ls_maxit,ell];
 

options=[]; sigma=[]; varg=[]; L_precond = []; U_precond = [];
for j = 2:nargin
  if isstruct(varargin{j})
    options = varargin{j};
  elseif ischar(varargin{j}) 
    if length(varargin{j}) == 2 & isempty(sigma)
      sigma = varargin{j};
    elseif isempty(L_precond)
      L_precond=varargin{j};
    elseif isempty(U_precond)
      U_precond=varargin{j};
    end
  elseif length(varargin{j}) == 1
    varg = [varg,varargin{j}];
  elseif min(size(varargin{j}))==1 
    sigma = varargin{j}; if size(sigma,1)==1, sigma=conj(sigma'); end 
  elseif isempty(L_precond)
    L_precond=varargin{j};
  elseif isempty(U_precond)
    U_precond=varargin{j};
  end
end

if ischar(sigma)
  sigma0=sigma; sigma=upper(sigma);
  switch sigma
    case {'LM','LR','SR','BE','SM'}
    otherwise
      if exist(sigma0)==2 & isempty(L_precond)
        ok=1; eval('v=feval(sigma,zeros(n,1));','ok=0')
        if ok, L_precond=sigma0; sigma=[]; end
      end
  end
end

[s,I]=sort(varg); I=flipdim(I,2);  
J=[]; j=0; 
while j<length(varg)
  j=j+1; jj=I(j); s=varg(jj);
  if isreal(s) & (s == fix(s)) & (s > 0)
    if n==-1
      n=s; eval('v=feval(A_operator,zeros(n,0));','n=-1;')
      if n>-1, J=[J,jj]; end
    end
  else
    if isempty(sigma), sigma=s; 
    elseif ischar(sigma) & isempty(L_precond)
      ok=1; eval('v=feval(sigma0,zeros(n,1));','ok=0')
      if ok, L_precond=sigma0; sigma=s; end
    end, J=[J,jj];  
  end 
end
varg(J)=[];

if n==-1,
  msg1=sprintf('  Cannot find the dimension of ''%s''.  \n',A_operator);
  msg2=sprintf('  Put the dimension n in the parameter list:  \n  like');
  msg3=sprintf('\t\n\n\t jdqr(''%s'',n,..),  \n\n',A_operator);
  msg4=sprintf('  or let\n\n\t n = %s(',A_operator);
  msg5=sprintf('[],''dimension'')\n\n  give n.');
  msg=[msg1,msg2,msg3,msg4,msg5];
  errordlg(msg,'MATRIX')
end

nselect=[]; 
if n<2, return, end

if length(varg) == 1
   nselect=min(n,varg);
elseif length(varg)>1
   if isempty(sigma), sigma=varg(end); varg(end)=[]; end
   nselect=min(n,min(varg));
end

fopts = []; if ~isempty(options), fopts=fields(options); end

if isempty(L_precond) 
  if strmatch('Precond',fopts) 
     L_precond = options.Precond;
  elseif strmatch('L_Precond',fopts)
     L_precond = options.L_Precond;
  end
end

if isempty(U_precond) & strmatch('U_Precond',fopts)
  U_precond = options.U_Precond;
end

if isempty(L_precond), ls_tol  = [0.7,0.49]; end
if ~isempty(L_precond) & ischar(L_precond)
  if exist(L_precond) ~=2
    msg=sprintf('  Can not find the M-file ''%s.m''  ',L_precond); n=-1;
  elseif ~isempty(U_precond) & ~ischar(U_precond) & n>0
    msg=sprintf('  L and U should both be strings or matrices'); n=-1;
  elseif strcmp(L_precond,U_precond)
    eval('v=feval(L_precond,zeros(n,1),''L'');','n=-1;')
    eval('v=feval(L_precond,zeros(n,1),''U'');','n=-1;')
    if n<0
       msg='L and U use the same M-file';
       msg1=sprintf(' %s.m   \n',L_precond);
       msg2='Therefore L and U are called';
       msg3=sprintf(' as\n\n\tw=%s(v,''L'')',L_precond); 
       msg4=sprintf(' \n\tw=%s(v,''U'')\n\n',L_precond); 
       msg5=sprintf('Check the dimensions and/or\n');
       msg6=sprintf('put this "switch" in %s.m.',L_precond);
       msg=[msg,msg1,msg2,msg3,msg4,msg5,msg6];
    else
       U_precond=0;
    end 
  elseif ischar(A_operator) & strcmp(A_operator,L_precond) 
    U_precond='preconditioner';
    eval('v=feval(L_precond,zeros(n,1),U_precond);','n=-1;')
    if n<0
       msg='Preconditioner and matrix use the same M-file';
       msg1=sprintf(' %s.   \n',L_precond);
       msg2='Therefore the preconditioner is called';
       msg3=sprintf(' as\n\n\tw=%s(v,''preconditioner'')\n\n',L_precond); 
       msg4='Put this "switch" in the M-file.';
       msg=[msg,msg1,msg2,msg3,msg4];
    end 
  else
    eval('v=feval(L_precond,zeros(n,1));','n=-1')
    if n<0
       msg=sprintf('''%s'' should produce %i-vectors',L_precond,n); 
    end
  end
end

Ud=1;
if ~isempty(L_precond) & ~ischar(L_precond) & n>0
  if ~isempty(U_precond) & ischar(U_precond)
    msg=sprintf('  L and U should both be strings or matrices'); n=-1;
  elseif ~isempty(U_precond)
    if ~min([n,n]==size(L_precond) &  [n,n]==size(U_precond))
      msg=sprintf('Both L and U should be %iX%i.',n,n); n=-1; 
    end
  elseif min([n,n]==size(L_precond))
    U_precond=speye(n); Ud=0;
  elseif min([n,2*n]==size(L_precond)) 
    U_precond=L_precond(:,n+1:2*n); L_precond=L_precond(:,1:n);
  else 
    msg=sprintf('The preconditioning matrix\n');
    msg2=sprintf('should be %iX%i or %ix%i ([L,U]).\n',n,n,n,2*n); 
    msg=[msg,msg2]; n=-1;
  end
end
if n<0, errordlg(msg,'PRECONDITIONER'), return, end

ls_tol0=ls_tol;
if strmatch('Tol',fopts),   tol = options.Tol;       end

if isempty(nselect), nselect=min(n,5); end

if strmatch('jmin',fopts), jmin=min(n,options.jmin); end 
if strmatch('jmax',fopts)
  jmax=min(n,options.jmax);
  if jmin<0, jmin=max(1,jmax-p1); end
else
  if jmin<0, jmin=min(n,nselect+p0); end
  jmax=min(n,jmin+p1); 
end 

if strmatch('MaxIt',fopts), maxit = abs(options.MaxIt);   end

if strmatch('v0',fopts);
  V = options.v0; 
  [m,d]=size(V); 
  if m~=n | d==0 
    if m>n, V = V(1:n,:); end
    nrV=norm(V); if nrV>0, V=V/nrV; end
    d=max(d,1);
    V = [V; ones(n-m,d) +0.1*rand(n-m,d)]; 
  end
else
  V = ones(n,1) +0.1*rand(n,1); d=1;
end

if strmatch('TestSpace',fopts), INTERIOR = boolean(options.TestSpace,...
         [INTERIOR,0,1,1.1,1.2],strvcat('standard','harmonic'));
end

if isempty(sigma)
  if INTERIOR, sigma=0; else, sigma = 'LM'; end
elseif ischar(sigma)
  switch sigma
    case {'LM','LR','SR','BE'}
    case {'SM'}
      sigma=0;
    otherwise
      if INTERIOR, sigma=0; else, sigma='LM'; end
   end
end

if strmatch('Schur',fopts),    SCHUR = boolean(options.Schur,SCHUR);    end
if strmatch('Disp',fopts),     SHOW  = boolean(options.Disp,[SHOW,2]);  end
if strmatch('Pairs',fopts),    PAIRS = boolean(options.Pairs,PAIRS);    end
if strmatch('AvoidStag',fopts),JDV   = boolean(options.AvoidStag,JDV);  end
if strmatch('Track',fopts)
  OLD = boolean(options.Track,[OLD,0,OLD,inf],strvcat('no','yes'));     end 
OLD=max(abs(OLD),10*tol);

if strmatch('LSolver',fopts), lsolver = lower(options.LSolver);         end

switch lsolver
  case {'exact'}
    L_precond=[]; 
    if ischar(A_operator)
  msg=sprintf('The operator must be a matrix for ''exact''.');
  msg=[msg,sprintf('\nDo you want to solve the correction equation')];
  msg=[msg,sprintf('\naccurately with an iterative solver (BiCGstab)?')];
      button=questdlg(msg,'Solving exactly','Yes','No','Yes');
      if strcmp(button,'No'), n=-1; return, end
    end
  case {'iluexact'}
    if ischar(L_precond)
  msg=sprintf('The preconditioner must be matrices for ''iluexact''.');
      errordlg(msg,'Solving with ''iluexact''')
      n=-1; return
    end
  case {'olsen'}
  case {'cg','minres','symmlq'}
    ls_tol=1.0e-12;
  case 'gmres'
    ls_maxit=5;
  case 'bicgstab'
    ls_tol=1.0e-10;
  otherwise
    error(['Unknown method ''' lsolver '''.']);
end

if strmatch('LS_MaxIt',fopts), ls_maxit=abs(options.LS_MaxIt); end
if strmatch('LS_Tol',fopts),   ls_tol=abs(options.LS_Tol);     end
if strmatch('LS_ell',fopts),   ell=round(abs(options.LS_ell)); end

par=[ls_tol,ls_maxit,ell];

if SHOW

  fprintf('\n'),fprintf('PROBLEM\n')
  if ischar(A_operator)
    fprintf('            A: ''%s''\n',A_operator);
  elseif issparse(A_operator)
    fprintf('            A: [%ix%i sparse]\n',n,n);
  else
    fprintf('            A: [%ix%i double]\n',n,n);
  end
    fprintf('    dimension: %i\n',n);
    fprintf('      nselect: %i\n\n',nselect);

  fprintf('TARGET\n')
  if ischar(sigma)
    fprintf('        sigma: ''%s''',sigma)
  else 
    Str=ShowLambda(sigma);
    fprintf('        sigma: %s',Str)
  end
    fprintf('\n\n')

  fprintf('OPTIONS\n');
    fprintf('        Schur: %i\n',SCHUR);
    fprintf('          Tol: %g\n',tol);
    fprintf('         Disp: %i\n',SHOW);
    fprintf('         jmin: %i\n',jmin);
    fprintf('         jmax: %i\n',jmax);
    fprintf('        MaxIt: %i\n',maxit);
    fprintf('           v0: [%ix%i double]\n',size(V));
    fprintf('    TestSpace: %g\n',INTERIOR);
    fprintf('        Pairs: %i\n',PAIRS);
    fprintf('    AvoidStag: %i\n',JDV);
    fprintf('        Track: %g\n',OLD);

    fprintf('      LSolver: ''%s''\n',lsolver);


  switch lsolver
    case {'exact','iluexact','olsen'}
    case {'cg','minres','symmlq','gmres','bicgstab'}
    if length(ls_tol)>1
      fprintf('       LS_Tol: ['); fprintf(' %g',ls_tol); fprintf(' ]\n');
    else
      fprintf('       LS_Tol: %g\n',ls_tol); 
    end
      fprintf('     LS_MaxIt: %i\n',ls_maxit);
  end
  if strcmp(lsolver,'bicgstab')
    fprintf('       LS_ell: %i\n',ell);
  end

  if isempty(L_precond)
    fprintf('      Precond: []\n');
  else
    StrL='Precond'; StrU='U precond'; Us='double'; Ls=Us; ok=0;
    if issparse(L_precond), Ls='sparse'; end
    if ~isempty(U_precond) & Ud, StrL='L precond'; ok=1;
      if issparse(U_precond), Us='sparse'; end
    end
    if ischar(L_precond)
      fprintf('%13s: ''%s''\n',StrL,L_precond);
      if ok 
        if U_precond~=0 & ~strcmp(U_precond,'preconditioner')
          fprintf('%13s: ''%s''\n',StrU,U_precond);
        else
          fprintf('%13s: ''%s''\n',StrU,L_precond);
        end
      end
    else
      fprintf('%13s: [%ix%i %s]\n',StrL,n,n,Ls);
      if ok & Ud, fprintf('%13s: [%ix%i %s]\n',StrU,n,n,Us); end
    end 
  end  
  fprintf('\n')

  string1='%13s: ''%s'''; string2='\n\t       %s';
  switch INTERIOR
  case 0
      fprintf(string1,'TestSpace','Standard, W = V ')
      fprintf(string2,'V W: V orthogonal')
      fprintf(string2,'W=A*V')
   case 1
      fprintf(string1,'TestSpace','Harmonic, W = A*V - sigma*V')
      fprintf(string2,'V W: V and W orthogonal')
      fprintf(string2,'AV-Q*E=W*R where AV=A*V-sigma*V and E=Q''*AV')
   case 1.1
      fprintf(string1,'TestSpace','Harmonic, W = A*V - sigma*V')
      fprintf(string2,'V W AV: V and W orthogonal')
      fprintf(string2,'AV=A*V-sigma*V, AV-Q*Q''*AV=W*R')
   case 1.2
      fprintf(string1,'TestSpace','Harmonic, W = A*V - sigma*V')
      fprintf(string2,'V W: W orthogonal')
      fprintf(string2,'W=AV-Q*E where AV=A*V-sigma*V and E=Q''*AV')
  otherwise
      fprintf(string1,'TestSpace','Experimental')
  end % switch INTERIOR
  fprintf('\n\n')
end % if SHOW

if ischar(sigma) & INTERIOR >=1
   msg1=sprintf('\n   The choice sigma = ''%s'' does not match',sigma);
   msg2=sprintf('\n   the search for INTERIOR eigenvalues.');
   msg3=sprintf('\n   Specify a numerical value for sigma');
   msg4=sprintf(',\n   for instance, a value that is ');
   switch sigma
      case {'LM'}
         % sigma='SM';
         msg5=sprintf('absolute large.\n');
         msg4=[msg4,msg5];
      case {'LR'}
         % sigma='SP'; % smallest positive real
         msg5=sprintf('positive large.\n');
         msg4=[msg4,msg5];
      case {'SR'}
         % sigma='LN'; % largest negative real
         msg5=sprintf('negative and absolute large.\n');
         msg4=[msg4,msg5];
      case {'BE'}
         % sigma='AS'; % alternating smallest pos., largest neg
         msg4=sprintf('.\n');
   end
   msg=[msg1,msg2,msg3,msg4];
   msg5=sprintf('   Do you want to continue with sigma=0?');
   msg=[msg,msg5];
   button=questdlg(msg,'Finding Interior Eigenvalues','Yes','No','Yes');
   if strcmp(button,'Yes'), sigma=0, else, n=-1; end
end

return

%-------------------------------------------------------------------
function x = boolean(x,gamma,string)
%Y = BOOLEAN(X,GAMMA,STRING)
%  GAMMA(1) is the default. 
%  If GAMMA is not specified, GAMMA = 0.
%  STRING is a matrix of accepted strings. 
%  If STRING is not specified STRING = ['no ';'yes']
%  STRING(I,:) and GAMMA(I) are accepted expressions for X 
%  If X=GAMMA(I) then Y=X. If X=STRING(I,:), then Y=GAMMA(I+1).
%  For other values of X, Y=GAMMA(1);

if nargin < 2, gamma=0; end
if nargin < 3, string=strvcat('no','yes'); gamma=[gamma,0,1]; end

if ischar(x)
  i=strmatch(lower(x),string,'exact'); 
  if isempty(i),i=1; else, i=i+