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)';  endreturn%%%======== END sort QZ decomposition interaction matrices ==============%%%======================================================================%%%======== INITIALIZATION ==============================================%%%======================================================================function MyClearglobal 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 QastMinvZreturn%%%======================================================================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 defaultsglobal Operator_Form Operator_A Operator_B Operator_Params ...       Precond_Form Precond_L Precond_U Precond_P Precond_Params ...       CHORDALDISTANCEOperator_A = varargin{1};n=CheckMatrix(Operator_A,1);% defaults              %%%% search for 'xx' in fieldnamesnselect0= 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'%% initiationnselect=[]; 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];   endend%------- 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  endend%------- 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;   endelse   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;   endendif 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;   endelse   Precond_L = L;endif ~isempty(Precond_L)  Precond_U=findfield(options,fopts,'u_',[]);  Precond_P=findfield(options,fopts,'p_',[]);endPrecond_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; endendTP=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); endend 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)]; endendV=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; endPAIRS  = 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;   endelse   TS=max(0,min(4,round(TS0(1,1))));end%------- set targets ----------------if isempty(Sigma)   if TS==1, Sigma=[0,1]; else, Sigma = 'LM'; endelseif ischar(Sigma)   switch Sigma     case {'LM','LR','SR','BE','SM'}     if ~ok, TS=3; end   endelse  [k,l]=size(Sigma);  if l==1, Sigma=[Sigma,ones(k,1)]; l=2; end  Sigma=ScaleEig(Sigma);endif 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, endend%------- 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, endendif     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;endls_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');endreturn%------------------------------------------------------------------------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);  endreturn%------------------------------------------------------------------------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)  endreturn%------------------------------------------------------------------------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