对“视觉机器学习20讲配套仿真代码”的研究心得---流形学习

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%功能:演示流形学习算法在计算机视觉中的应用
%基于流形学习实现目标分类;
%环境:Win7,Matlab2012b
%Modi: NUDT-VAP
%时间:2014-02-04
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% SWISS ROLL DATASET


  N=2000;
  K=12;
  d=2; 


clf; colordef none; colormap jet; set(gcf,'Position',[200,400,620,200]);


% PLOT TRUE MANIFOLD
  tt0 = (3*pi/2)*(1+2*[0:0.02:1]); hh = [0:0.125:1]*30;
  xx = (tt0.*cos(tt0))'*ones(size(hh));
  yy = ones(size(tt0))'*hh;
  zz = (tt0.*sin(tt0))'*ones(size(hh));
  cc = tt0'*ones(size(hh));


  subplot(1,3,1); cla;
  surf(xx,yy,zz,cc);
  view([12 20]); grid off; axis off; hold on;
  lnx=-5*[3,3,3;3,-4,3]; lny=[0,0,0;32,0,0]; lnz=-5*[3,3,3;3,3,-3];
  lnh=line(lnx,lny,lnz);
  set(lnh,'Color',[1,1,1],'LineWidth',2,'LineStyle','-','Clipping','off');
  axis([-15,20,0,32,-15,15]);


% GENERATE SAMPLED DATA
  tt = (3*pi/2)*(1+2*rand(1,N));  height = 21*rand(1,N);
  X = [tt.*cos(tt); height; tt.*sin(tt)];


% SCATTERPLOT OF SAMPLED DATA
  subplot(1,3,2); cla;
  scatter3(X(1,:),X(2,:),X(3,:),12,tt,'+');
  view([12 20]); grid off; axis off; hold on;
  lnh=line(lnx,lny,lnz);
  set(lnh,'Color',[1,1,1],'LineWidth',2,'LineStyle','-','Clipping','off');
  axis([-15,20,0,32,-15,15]); drawnow;


% RUN LLE ALGORITHM
Y=lle(X,K,d);


% SCATTERPLOT OF EMBEDDING
  subplot(1,3,3); cla;
  scatter(Y(1,:),Y(2,:),12,tt,'+');
  grid off;
  set(gca,'XTick',[]); set(gca,'YTick',[]); 



&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

function  [v,d,flag] = eigs(varargin)
%EIGS   Find a few eigenvalues and eigenvectors.
%   EIGS solves the eigenvalue problem A*v = lambda*v or the generalized
%   eigenvalue problem A*v = lambda*B*v.  Only a few selected eigenvalues,
%   or eigenvalues and eigenvectors, are computed.
%
%   [V,D,FLAG] = EIGS(A)
%   [V,D,FLAG] = EIGS('Afun',N)
%
%   The first input argument is either a square matrix (which can be
%   full or sparse, symmetric or nonsymmetric, real or complex), or a
%   string containing the name of an M-file which applies a linear
%   operator to the columns of a given matrix.  In the latter case,
%   the second input argument must be N, the order of the problem.
%   For example, EIGS('fft',...) is much faster than EIGS(F,...)
%   where F is the explicit FFT matrix.
%
%   The remaining input arguments are optional and can be given in
%   practically any order:
%
%   [V,D,FLAG] = EIGS(A,B,K,SIGMA,OPTIONS)
%   [V,D,FLAG] = EIGS('Afun',N,B,K,SIGMA,OPTIONS)
%
%   where
%
%       B         A symmetric positive definite matrix the same size as A.
%       K         An integer, the number of eigenvalues desired.
%       SIGMA     A scalar shift or a two letter string.
%       OPTIONS   A structure containing additional parameters.
%
%   With one output argument, D is a vector containing K eigenvalues.
%   With two output arguments, D is a K-by-K diagonal matrix and V is
%   a matrix with K columns so that A*V = V*D or A*V = B*V*D.
%   With three output arguments, FLAG indicates whether or not the
%   eigenvalues converged to the desired tolerance.  FLAG = 0 indicated
%   convergence, FLAG = 1 not.  FLAG = 2 indicates that EIGS stagnated
%   i.e. two consecutive iterates were the same.
%
%   If B is not specified, B = eye(size(A)) is used.  Note that if B is
%   specified, then its Cholesky factorization is computed.
%
%   If K is not specified, K = MIN(N,6) eigenvalues are computed.
%
%   If SIGMA is not specified, the K-th eigenvalues largest in magnitude
%   are computed.  If SIGMA is zero, the K-th eigenvalues smallest in
%   magnitude are computed.  If SIGMA is a real or complex scalar, the
%   "shift", the K-th eigenvalues nearest SIGMA are computed using
%   shift-invert.  Note that this requires computation of the LU
%   factorization of A-SIGMA*B.  If SIGMA is one of the following strings,
%   it specifies the desired eigenvalues.
%
%       SIGMA            Specified eigenvalues
%
%       'LM'             Largest Magnitude  (the default)
%       'SM'             Smallest Magnitude (same as sigma = 0)
%       'LR'             Largest Real part
%       'SR'             Smallest Real part
%       'BE'             Both Ends.  Computes k/2 eigenvalues
%                        from each end of the spectrum (one more
%                        from the high end if k is odd.)
%
%   The OPTIONS structure specifies certain parameters in the algorithm.
%
%     Field name       Parameter                             Default
%
%     OPTIONS.tol      Convergence tolerance:                1e-10 (symmetric)
%                      norm(A*V-V*D,1) <= tol * norm(A,1)    1e-6 (nonsymm.)
%     OPTIONS.stagtol  Stagnation tolerance: quit when       1e-6
%                      consecutive iterates are the same.
%     OPTIONS.p        Dimension of the Arnoldi basis.       2*k
%     OPTIONS.maxit    Maximum number of iterations.         300
%     OPTIONS.disp     Number of eigenvalues displayed       20
%                      at each iteration.  Set to 0 for
%                      no intermediate output.
%     OPTIONS.issym    Positive if Afun is symmetric.        0
%     OPTIONS.cheb     Positive if A is a string,            0
%                      sigma is 'LR','SR', or a shift, and
%                      polynomial acceleration should be
%                      applied.
%     OPTIONS.v0       Starting vector for the Arnoldi       rand(n,1)-.5
%                      factorization
%
%   See also EIG, SVDS.


%   Richard J. Radke and Dan Sorensen.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.19 $  $Date: 1998/08/14 18:50:27 $


global ChebyStruct   % Variables used by accpoly.m 


%  ======   Check input values and set defaults where needed


A = varargin{1};
if isstr(A)
   n = varargin{2};
else
   [n,n] = size(A);
   if any(size(A) ~= n)
      error('Matrix A must be a square matrix or a string.')
   end
end


if n == 0
   v = [];
   d = [];
   flag = 0;
   return
end


% No need to iterate - also elimiates confusion with
% remaining input arguments (they are all n-by-n).


if n == 1
   B = 1;
   for j = (2+isstr(A)):nargin
      if isequal(size(varargin{j}),[n,n]) & ~isstruct(varargin{j})
         B = varargin{j};
         if ~isreal(B)
            error('B must be real and symmetric positive definite.')
      end
         break
      end
   end
   if isstr(A)
      if nargout < 2
         v = feval(A,1) / B;
      else
         v = 1;
         d = feval(A,1) / B;
         flag = 0;
      end
      return
   else
      if nargout < 2
         v = A / B;
      else
         v = 1;
         d = A / B;
         flag = 0;
      end
      return
   end
end


B = [];
b2 = 1;
k = [];
sigma = [];
options = [];
for j = (2+isstr(A)):nargin
   if isequal(size(varargin{j}),[n,n])
      B = varargin{j};
      if ~isreal(B)
         error('B must be real and symmetric positive definite.')
      end
   elseif isstruct(varargin{j})
      options = varargin{j};
   elseif isstr(varargin{j}) & length(varargin{j}) == 2
      sigma = varargin{j};
   elseif max(size(varargin{j})) == 1
      s = varargin{j};
      if isempty(k) & isreal(s) & (s == fix(s)) & (s > 0)
         k = min(n,s);
      else
         sigma = s;
      end
   else
      error(['Input argument number ' int2str(j) ' not recognized.'])
   end
end


% Set defaults


if isempty(k)
   k = min(n,6);
end
if isempty(sigma)
   sigma = 'LM';
end
if isempty(options)
   fopts = [];
else
   fopts = fieldnames(options);
end


% A is the matrix of all zeros (not detectable if A is a string)


if ~isstr(A)
   if nnz(A) == 0
      if nargout < 2
         v = zeros(k,1);
      else
         v = eye(n,k);
         d = zeros(k,k);
         flag = 0;
      end
      return
   end
end


%  Trick to get faster convergence.  The tail end (closest to cut off
%  of the sort) will not converge as fast as the leading end of the
%  "wanted" Ritz values.


ksave = k;
k = min(n-1,k+3);


% Set issym, specifying a symmetric matrix or operator.


if ~isstr(A)
   issym = isequal(A,A');
elseif strmatch('issym',fopts)
   issym = options.issym;
else
   issym = 0;
end


% Set tol, the convergence tolerance.


if strmatch('tol',fopts);
   tol = options.tol;
else
   if issym
      tol = 1e-10;         % Default tol for symmetric problems.
   else
      tol = 1e-6;          % Default tol for nonsymmetric problems.
   end      
end


% Set stagtol, the stagnation tolerance.


if strmatch('stagtol',fopts);
   stagtol = options.stagtol;
else
   stagtol = 1e-6;
end


% Set p, the dimension of the Arnoldi basis.


if strmatch('p',fopts);
   p = options.p;
else
   p = min(max(2*k,20),n);
end


% Set maxit, the maximum number of iterations.


if strmatch('maxit',fopts);
   maxit = options.maxit;
else
   maxit = max(300,2*n/p);
end


% Set display option.


if strmatch('disp',fopts);
   dispn = options.disp;
else
   dispn = 20;
end


if strmatch('v0',fopts);
   v0 = options.v0;
else
   v0 = rand(n,1)-.5;
end


% Check if Chebyshev iteration is requested.


if strmatch('cheb',fopts);
   cheb = options.cheb;
else
   cheb = 0;
end


if cheb & issym
   if isstr(A) & strcmp(sigma,'LR')
      sigma = 'LO';
      ChebyStruct.sig = 'LR';
      ChebyStruct.filtpoly = [];
   elseif isstr(A) & strcmp(sigma,'SR')
      sigma = 'SO';
      ChebyStruct.sig = 'SR';
      ChebyStruct.filtpoly = [];
   elseif isstr(A) & ~isstr(sigma)
      ChebyStruct.sig = sigma;
      sigma = 'LO';
      ChebyStruct.filtpoly = [];
   elseif isstr(A)
      warning('Cheb option only available for symmetric operators and when sigma is ''LR'', ''SR'', or numeric.  Proceeding without cheb.')
      cheb = 0;
   end
elseif cheb
   warning('Cheb option only available for symmetric operators and when sigma is ''LR'', ''SR'', or numeric.  Proceeding without cheb.')
   cheb = 0;
end


if isstr(sigma)
   sigma = upper(sigma);
end


if strcmp(sigma,'SM') & ~isstr(A)
   sigma = 0;
end


if k > n | k < 0
   error('k must be an integer between 0 and n.')
elseif isstr(A) & ~isstr(sigma) & ~cheb
   error('Operator A and numeric sigma not compatible (unless options.cheb > 0).')
elseif fix(k) ~= k
   error('k must be an integer between 0 and n.')
elseif p > n
   error('p must satisfy p <= n.')
elseif issym & ~isreal(sigma)
   error('Use real sigma for a symmetric matrix.')
end


if ~isstr(A) & ~isempty(B);
   pmmd = symmmd(spones(A) | spones(B));
   A = A(pmmd,pmmd);
   B = B(pmmd,pmmd);
elseif ~isstr(A)
   pmmd = symmmd(A);
   A = A(pmmd,pmmd);
end   


B2 = B;


if ~isequal(B,B')
   error('B must be real and symmetric positive definite.')
end
if ~isempty(B)
   [BL,b2] = chol(B);
   if b2
      error('B must be real and symmetric positive definite.')
   end
end


info = 0;
v = v0;
p = p - k;
psave = p;
op = A;
ChebyStruct.op = op;


if ~isempty(B)
   beta = 1.0/sqrt(v(:,1)'*B*v(:,1));
else
   beta = 1.0/norm(v(:,1));
end


v(:,1) = v(:,1)*beta;


%  Compute w = Av;


if ~isempty(B)
   if ~isstr(A) & ~isstr(sigma),
      [L,U] = lu(A - sigma*B);
      condU = condest(U);
      dsigma = n * full(max(max(abs(A)))) * eps;
      if sigma < 0
         sgnsig = -1;
      else
         sgnsig = 1;
      end
      sigitr = 1;
      while condU > 1/eps & ((dsigma <= 1 & sigitr <= 10) | ~isfinite(condU))
         disps1 = sprintf(['sigma = %10e is near an exact solution ' ...
               'to the generalized eigenvalue problem of A and B,\n' ...
               'so we cannot use the LU factorization of (A-sigma*B): ' ...
               'condest(U) = %10e.\n'], ...
            sigma,condU);
         if abs(sigma) < 1
            sigma = sigma + sgnsig * dsigma;
            disps2 = sprintf('We are trying sigma + %10e = %10e instead.\n', ...
               sgnsig*dsigma,sigma);
         else
            sigma = sigma * (1 + dsigma);
            disps2 = sprintf('We are trying sigma * (1 + %10e) = %10e instead.\n', ...
               dsigma,sigma);
         end
         if nargout < 3 & dispn ~= 0             
            disp([disps1 disps2])
         end   
         [L,U] = lu(A - sigma*B);
         condU = condest(U);
         dsigma = 10 * dsigma;
         sigitr = sigitr + 1;
      end
      if sigitr > 1 & nargout < 3 & dispn ~= 0
         disps = sprintf(['LU factorization of (A-sigma*I) for ' ...
               'sigma = %10e: condest(U) = %10e\n'],sigma,condU);
         disp(disps)             
      end
      v(:,1) = U \ (L \ (B*v(:,1)));
      if sum(~isfinite(v(:,1)))
         d = Inf;
         if nargout <= 1
            v = d;
         end
         flag = 1;
         return
      end
      beta = 1.0/sqrt(v(:,1)'*B*v(:,1));
      v(:,1) = v(:,1)*beta;
      v(:,2) = U \ (L \ (B*v(:,1)));
      if sum(~isfinite(v(:,2)))
         d = Inf;
         if nargout <= 1
            v = d;
         end
         flag = 1;
         return
      end
   elseif isstr(A) & b2
      v(:,2) = B \ feval(A,v(:,1));
   elseif isstr(A) & ~b2
      v(:,2) = BL \ feval(A,(BL'\v(:,1)));
   elseif ~b2
      v(:,2) = BL \ (A*(BL'\v(:,1)));
   else
      v(:,2) = B \ (A*v(:,1));
   end
else
   if ~isstr(A) & ~isstr(sigma),
      [L,U] = lu(A - sigma*speye(n));
      condU = condest(U);
      dsigma = n * full(max(max(abs(A)))) * eps;
      if sigma < 0
         sgnsig = -1;
      else
         sgnsig = 1;
      end
      sigitr = 1;
      while condU > 1/eps & ((dsigma <= 1 & sigitr <= 10) | ~isfinite(condU))
         disps1 = sprintf(['sigma = %10e is near an exact eigenvalue of A,\n' ...
               'so we cannot use the LU factorization of (A-sigma*I): ' ...
               ' condest(U) = %10e.\n'],sigma,condU);
         if abs(sigma) < 1
            sigma = sigma + sgnsig * dsigma;
            disps2 = sprintf('We are trying sigma + %10e = %10e instead.\n', ...
               sgnsig*dsigma,sigma);
         else
            sigma = sigma * (1 + dsigma);
            disps2 = sprintf('We are trying sigma * (1 + %10e) = %10e instead.\n', ...
               dsigma,sigma);
         end
         if nargout < 3 & dispn ~= 0
            disp([disps1 disps2])
         end
         [L,U] = lu(A - sigma*speye(n));
         condU = condest(U);
         dsigma = 10 * dsigma;
         sigitr = sigitr + 1;
      end
      if sigitr > 1 & nargout < 3 & dispn ~= 0
         disps = sprintf(['LU factorization of (A-sigma*I) for ' ...
               'sigma = %10e: condest(U) = %10e\n'],sigma,condU);
         disp(disps)             
      end         
      v(:,1) = U \ (L \ v(:,1));
      if sum(~isfinite(v(:,1)))
         d = Inf;
         if nargout <= 1
            v = d;
         end
         flag = 1;
         return
      end
      beta = 1.0/norm(v(:,1));
      v(:,1) = v(:,1)*beta;
      v(:,2) = U \ (L \ v(:,1));
      if sum(~isfinite(v(:,2)))
         d = Inf;
         if nargout <= 1
            v = d;
         else
            v = v(:,1);
         end
         flag = 1;
         return
      end
   elseif isstr(A)
      v(:,2) = feval(A,v(:,1));
   else
      v(:,2) = A*v(:,1);
   end
end


if ~isempty(B)
   beta = 1.0/sqrt(v(:,1)'*B*v(:,1));
else
   beta = 1.0/norm(v(:,1));
end


v(:,1) = v(:,1)*beta;


if ~isempty(B)
   alpha = v(:,1)'*B*v(:,2);
else
   alpha = v(:,1)'*v(:,2);
end


h(1,1) = alpha;


%  Compute the residual in the second column of V


v(:,2) = v(:,2) - alpha*v(:,1);


%  Perform one step of iterative refinement to correct any
%  orthogonality problems


if ~isempty(B)
   alpha = v(:,1)'*B*v(:,2);
else
   alpha = v(:,1)'*v(:,2);
end


v(:,2) = v(:,2) - alpha*v(:,1);


h(1,1) = h(1,1) + alpha;


%  Compute k steps of the Arnoldi sequence


kstart = 1;
ritz = 1.0;
kp1 = k + 1;
kend = k + p;
k1 = 1;


if ~isempty(B)
   if ~isstr(A) & ~isstr(sigma)
      [v,h,info] =  arnold2(k1,k,L,U,B,v,h,tol);
      if info == -1
         d = h(1,1);
         if nargout <= 1
            v = d;
         else
            v = v(:,1);
         end
         flag = 1;
         return
      end
   elseif ~b2
      [v,h] =  arnold(k1,k,A,BL,v,h,1);
   else
      [v,h] =  arnold(k1,k,A,B,v,h,0);
   end
else
   if ~isstr(A) & ~isstr(sigma)
      [v,h,info] =  arnold2nob(k1,k,L,U,v,h,tol);
      if info == -1
         d = h(1,1);
         if nargout <= 1
            v = d;
         else
            v = v(:,1);
         end
         flag = 1;
         return
      end
   else
      [v,h] =  arnoldnob(k1,k,A,v,h);
   end
end


%  Now update the Arnoldi sequence in place


iter = 0;
knew = k;
ritzes = zeros(ksave,1);
ritzests = ones(ksave,1);
stopcrit = 1;
beta = 1;
betanew = 1;
residest = 1;
stagnation = 0;


%  MAIN LOOP


while (((stopcrit > tol) & (iter < maxit) & ~stagnation) | iter < 2)
   
   iter = iter + 1;
   if dispn
      iter
   end
   
   ChebyStruct.iter = iter;
   
   %     Compute p additional steps of the Arnoldi sequence
   
   kold = k;
   k = knew;
   
   if ~isempty(B)
      if ~isstr(A) & ~isstr(sigma)
         [v,h,info] = arnold2(k,kend,L,U,B,v,h,tol);
         if info == -1
            iter = iter - 1;
            flag = 1;
            break
         end
      elseif ~b2
         [v,h,info] = arnold(k,kend,A,BL,v,h,1);
      else
         [v,h,info] =  arnold(k,kend,A,B,v,h,0);
      end
   else
      if ~isstr(A) & ~isstr(sigma)
         [v,h,info] = arnold2nob(k,kend,L,U,v,h,tol);
         if info == -1
            iter = iter - 1;
            flag = 1;
            break
         end
      else
         [v,h,info] = arnoldnob(k,kend,A,v,h);
  end
   end
   
   %     If we're doing Chebyshev iteration, when the Ritz estimates
   %     on the extreme values of the spectrum are good enough, then
   %     apply the interpolant polynomial instead.
   
   if cheb == 1 & (ritzests(1) < .1) & (ritzests(2) < .1)
      
      if dispn, disp('Starting polynomial acceleration'), end
      
      A = 'accpoly';
      eigh = eig(h);       
      sig = sigma;
      iter = 1;
      ChebyStruct.iter = iter;
      
      %  Restart the basis with the largest/smallest
      %  ritz vector, as appropriate.
      
      v = v(:,2);
      if ~isempty(B)
         beta = 1.0/sqrt(v(:,1)'*B*v(:,1));
      else
         beta = 1.0/norm(v(:,1));
      end
      
      v(:,1) = v(:,1)*beta;
      
      if ~isempty(B)
         if isstr(A) & b2
            v(:,2) = B \ feval(A,v(:,1));
         elseif isstr(A) & ~b2
            v(:,2) = BL \ feval(A,(BL'\v(:,1)));
         elseif ~b2
            v(:,2) = BL \ (A*(BL'\v(:,1)));
         else
            v(:,2) = B \ (A*v(:,1));
         end
         
         beta = 1.0/sqrt(v(:,1)'*B*v(:,1));
         v(:,1) = v(:,1)*beta;
         alpha = v(:,1)'*B*v(:,2);
         h = alpha;
         v(:,2) = v(:,2) - alpha*v(:,1);
         alpha = v(:,1)'*B*v(:,2);
         v(:,2) = v(:,2) - alpha*v(:,1);
         h = h + alpha;
      else
         if isstr(A)
            v(:,2) = feval(A,v(:,1));
         else
            v(:,2) = A*v(:,1);
         end
         
         beta = 1.0/norm(v(:,1));
         v(:,1) = v(:,1)*beta;
         alpha = v(:,1)'*v(:,2);
         h = alpha;
         v(:,2) = v(:,2) - alpha*v(:,1);
         alpha = v(:,1)'*v(:,2);
         v(:,2) = v(:,2) - alpha*v(:,1);
         h = h + alpha;
      end
      
      if ~b2
         [v,h,info] =  arnold(1,kend,A,BL,v,h,1);
      elseif ~isempty(B)
         [v,h,info] =  arnold(1,kend,A,B,v,h,0);
      else
         [v,h,info] =  arnoldnob(1,kend,A,v,h);
      end
      
      k1 = 1;
      knew = k;
      ritzes = zeros(ksave,1);
      cheb = -1;
      residest = 1;
      sigma = 'LR';
      
   end
   
   k = kold;
   
   %     Compute p shifts based on sigma
   
   %     If A is symmetric, keep H tridiagonal (to avoid spurious
   %     imaginary parts).
   
   if issym
      for i=1:kend-2
         h(i,i+2:kend) = zeros(1,kend-i-1);
         hav = mean([h(i,i+1),h(i+1,i)]);
         h(i,i+1) = hav;
         h(i+1,i) = hav;
      end
      hav = mean([h(kend,kend-1),h(kend-1,kend)]);
      h(kend,kend-1) = hav;
      h(kend-1,kend) = hav;
   end   
   
   [w,q1] = shftit(h,kstart,kend,sigma);
   
   %     Update the command window with current eigenvalue estimates
   
   ritzesold = ritzes;
   
   if ~isstr(sigma)
      ritzes = sigma + ...
         1./w((kend-kstart+1):-1:(kend-ksave+1));
      ritzes = [(sigma+1./eig(h(1:kstart-1,1:kstart-1)));...
            ritzes];
   else
      ritzes = w((kend-kstart+1):-1:(kend-ksave+1));
      ritzes = [eig(h(1:kstart-1,1:kstart-1));ritzes];
   end
   
   if dispn
      eigs = ritzes(1:min(dispn,length(ritzes)))
   end
   
   if iter == 1 & cheb > 0
      ChebyStruct.lbd = min(ritzes);
      ChebyStruct.ubd = max(ritzes);
   elseif cheb > 0
      ChebyStruct.lbd = min(ChebyStruct.lbd, min(ritzes));
      ChebyStruct.ubd = max(ChebyStruct.ubd, max(ritzes));
   end
   
   [m1,m2]=size(q1);
   ritz = norm(q1(m1,p+2:m1));
   
   if ~isempty(B)
      betanew = sqrt(v(:,kend+1)'*B*v(:,kend+1));
   else
      betanew = norm(v(:,kend+1));
   end
   
   ritznew = betanew*q1(m1,1:m1);
   jj = m1;
   kconv = 0;
   while(jj > 0),
      if (abs(ritznew(jj)) <= tol),
         jj = jj - 1;
         kconv = kconv+1;
      else
         jj = -1;
      end
   end
   kkconv = kconv;
   
   %     The while loop counts the number of converged ritz values.
   
   ritzests = [w abs(ritznew)'];
   
   ritzests = ritzests(size(ritzests,1) : -1 : ...
      size(ritzests,1)-ksave+1,2);
   
   stopcritold = stopcrit;
   stopcrit = max(ritzests);
   residest = norm(ritzests);    
   
   %  At first, we use the Ritz estimates to estimate convergence.
   %  However, as the algorithm converges, the Ritz estimates become
   %  poor estimates of the actual error in each Ritz pair.  So when the
   %  Ritz estimates become too small, we actually form the matrix of
   %  errors || AV - VD || where V are the estimates for the eigenvectors
   %  and eigenvalues.  This is expensive computationally but ensures
   %  that the user gets the desired eigenpairs to the requested
   %  tolerance.
   
   if max(ritzests) <= tol*max(norm(h),1)
      
      if ~b2 & isstr(sigma)
         vee = BL \ v(:,1:kend);
      else
         vee = v(:,1:kend);
      end
      
      if ~isempty(B2)
         if isstr(A)
            vee = vee * q1(1:kend,kend:-1:kend-ksave+1);
            for i = 1 : ksave
               nvi = norm(v(:,i));
               if isfinite(nvi) & nvi ~= 0
                  vee(:,i) = vee(:,i) / nvi;
               end
            end
            Avee = feval(A,vee);
            errmat = Avee - B2 * vee * diag(ritzes);
            residest = norm(errmat,1) / norm(Avee,1);
         else
            vee = vee * q1(1:kend, kend:-1:kend-ksave+1);
            for i = 1 : ksave
               nvi = norm(vee(:,i));
               if isfinite(nvi) & nvi ~= 0
                  vee(:,i) = vee(:,i) / nvi;
               end
            end
            errmat = A * vee - B2 * vee * diag(ritzes);
            residest = norm(errmat,1) / norm(A,1);
         end
      else
         if isstr(A)
            vee = vee * q1(1:kend,kend:-1:kend-ksave+1);
            for i = 1 : ksave
               nvi = norm(v(:,i));
               if isfinite(nvi) & nvi ~= 0
                  vee(:,i) = vee(:,i) / nvi;
               end
            end
            Avee = feval(A,vee);
            errmat = Avee - vee * diag(ritzes);
            residest = norm(errmat,1) / norm(Avee,1);
         else
            vee = vee * q1(1:kend, kend:-1:kend-ksave+1);
            for i = 1 : ksave
               nvi = norm(vee(:,i));
               if isfinite(nvi) & nvi ~= 0
                  vee(:,i) = vee(:,i) / nvi;
               end
            end
            errmat = A * vee - vee * diag(ritzes);
            residest = norm(errmat,1) / norm(A,1);
         end
      end
      
      for ii = 1:length(ritzes);
         ritzests(ii) = norm(errmat(:,ii));
      end
      
      stopcrit = residest;
      
   end
   
   if (abs(stopcritold-stopcrit) < stagtol*abs(stopcrit))
      stagind = (ritzesold~=0) & (ritzes~=0);
      stagn = abs(ritzesold - ritzes);
      stagn(stagind) = stagn(stagind) ./ ritzes(stagind);
      stagnation = max(stagn) < stagtol;
end


   if (((stopcrit > tol) & ~stagnation) | iter < 2)
      
      %     Apply the p implicit shifts if convergence has not yet 
      %     happened.  Otherwise don't apply them and get out of the
      %     loop on next loop test. We need to keep same test here as
      %     in the main loop test to   avoid applying shifts and then
      %     quitting, which would lead to a wrong size factorization
      %     on return.
      
      %        If some ritz values have converged then
      %        adjust k and p to move the "boundary"
      %        of the filter cutoff.
      
      if kconv > 0
         kk = ksave + 3 + kconv;
         p = max(ceil(psave/3),kend-kk);
         k = kend - p;
      end
      
      if any(any(imag(v))) | any(any(imag(h)))
         [v,h,knew] = apshft1(v,h,w,k,p);
      else
         [v,h,knew] = apshft2(v,h,real(w),imag(w),k,p);
      end
      
      if ~isempty(B)
         betanew = sqrt(v(:,kp1)'*B*v(:,kp1));
      else
         betanew = norm(v(:,kp1));
      end   
      
   end
   
   if dispn
      stopcrit
      disp('==========================')
   end
   
end         %  End of Arnoldi iteration  (MAIN LOOP)


if stopcrit <= tol
   flag = 0;
else
   if nargout < 3 & dispn ~= 0
      if iter >= maxit
         disp('Exiting: Maximum number of iterations exceeded.')
      elseif stagnation
         disp('Exiting: Two consecutive iterates were the same.')
      else
         disp('Exiting: Eigenvalues did not converge.')        
      end
   else
      if stagnation
         flag = 2;
      else
         flag = 1;
      end
   end
end


%  Compute the eigenvalues and eigenvectors of h
kend = min(kend,size(h,1));
[w,q1] = shftit(h,1,kend,sigma);


k = ksave;
p = psave;


%  Transform the converged eigenvalues back to 
%  the original problem and return them in the diagonal
%  k by k matrix d.


%  Set v equal to the wanted eigenvectors


v = v(:,1:kend) * q1(1:kend, kend:-1:kend-k+1);


if ~b2 & isstr(sigma)
   v = BL \ v;
end


for i = 1:k
   nvi = norm(v(:,i));
   if isfinite(nvi) & nvi ~= 0
      v(:,i) = v(:,i) / nvi;
   end   
end


%  In Chebyshev iteration, we recover the eigenvalues by Rayleigh
%  quotients.  Otherwise, the eigenvalues are recovered from the
%  set of shifts w.


if cheb == -1
   d = [];
   if ~isempty(B2)
      for i=1:ksave
         d = [d; ( (v(:,i)' * feval(op,v(:,i)) ) ./ ...
               (v(:,i)' * B2 * v(:,i)))];
      end
   else
      for i=1:ksave
         d = [d; ( (v(:,i)' * feval(op,v(:,i)) ) ./ ...
               norm(v(:,i)))];
      end
   end
   d = diag(d);
elseif ~isstr(sigma)
   t = sigma + 1./w;
   d = diag(t(kend:-1:kend-k+1));
else
   d = diag(w(kend:-1:kend-k+1));
end


if ~isstr(A)
   v(pmmd,:) = v;
end


if nargout <= 1
   v = diag(d);
end


% ====================================================


function [v,h,k] =  apshft1(v,h,w,k,p)


%  APSHFT1  Apply shifts to update an Arnoldi factorization
%
%  APSHFT1 was designed to be called by EIGS when V or H is complex.
%
%  [V,H,k] = APSHFT1(V,H,w,k,p) implicitly applies the
%     p real shifts held in w to update the existing Arnoldi
%     factorization AV - VH = re'  .
%                               k+p
%  The routine results in
%
%     A(VQ) - (VQ)(Q'HQ) = re'  Q
%                            k+p
%
%  where the orthogonal matrix Q is the product of the Givens
%  rotations resulting from p bulge chase sweeps.
%
%  The updated residual is placed in V(:,k+1) and the updated
%  Arnoldi factorization V <- VQ, H <- Q'HQ is returned.


%  Dan Sorensen and Richard J. Radke, 11/95.


k1 = 1;
kend = k+p;
kp1 = k + 1;
q = zeros(1,kend);
q(kend) = 1.0;


dh = diag(h,-1);
ix = find(dh(1:end-1)==0);     % Find the column indices of 0
ix = [0 ; ix ; kend-1];       % subdiagonals in H.
nx = size(ix,1);           


for jj = 1:p               % Loop over shifts
   for ii = 1:nx-1         % Loop over blocks in H
      k1 = ix(ii)+1; k2 = ix(ii+1);
      c = h(k1,k1) - w(jj);
      s = h(k1+1,k1);
      [G,R] = qr([c;s]);
      for i = k1:k2        % Loop over rows in the block
         if i > k1 
            [G,R] = qr(h(i:i+1,i-1));
            h(i:i+1,i-1) = R(:,1);
         end
         
         %           apply rotation from left to rows of H
         
         h(i:i+1,i:kend) = G'* h(i:i+1,i:kend);
         
         
         %           apply rotation from right to columns of H
         
         ip2 = i+2;
         if ip2 > kend
            ip2 = kend;
         end
         
         h(1:ip2,i:i+1) = h(1:ip2,i:i+1)*G;
         
         %           apply rotation from right to columns of V
         
         v(:,i:i+1) = v(:,i:i+1)*G;
         
         %           accumulate e'  Q so residual can be updated
         %                       k+p
         
         q(i:i+1) = q(i:i+1)*G;
      end
   end
end


%  Update the residual and store in the k+1 -st column of v


v(:,kend+1) = v(:,kend+1)*q(k);
v(:,kp1) = v(:,kend+1) + v(:,kp1)*h(kp1,k);


% ==========================================================


function [v,h,k] =  apshft2(v,h,wr,wi,k,p)


%  APSHFT2  Apply shifts to update an Arnoldi factorization
%
%  APSHFT2 was designed to be called by EIGS when V and H are real.
%
%  [V,H,k] = APSHFT2(V,H,wr,wi,k,p) implicitly applies
%     the p complex shifts given by w = wr + i*wi, to update the
%     existing Arnoldi factorization AV - VH = re'  .
%                                                k+p
%
%  The routine results in
%
%     A(VQ) - (VQ)(Q'HQ) = re'  Q
%                            k+p
%
%  where the orthogonal matrix Q is the product of the Givens
%  rotations resulting from p bulge chase sweeps.
%
%  The updated residual is placed in V(:,k+1) and the updated
%  Arnoldi factorization V <- VQ, H <- Q'HQ is returned.


%  Dan Sorensen and Richard J. Radke, 11/95.


k1 = 1;
kend = k+p;
kp1 = k + 1;
q = zeros(kend,1);
q(kend) = 1.0;
mark = 0;
num = 0;


dh = diag(h,-1);
ix1 = find(dh(1:end-1)==0);    % Find the column indices of 0
ix = [0 ; ix1 ; kend-1];         % subdiagonals in H.
jx = [0 ; ix1 ; kend-2];
nx = size(ix,1);           


for jj = 1:p
   
   %     compute and apply a bulge chase sweep initiated by the
   %     implicit shift held in w(jj)
   
   if (abs(wi(jj)) == 0.0) 
      
      %        apply a real shift using 2 by 2 Givens rotations
      
      for ii = 1:nx-1      % Loop over blocks in H
         k1 = ix(ii)+1;
         k2 = ix(ii+1);
         c = h(k1,k1) - wr(jj);
         s = h(k1+1,k1) ;
         t = norm([c s]);
         if t == 0.0
            c = 1.0;
         else
            c = c/t;
            s = s/t;
         end
         for i = k1:k2
            if i > k1
               t = norm(h(i:i+1,i-1));
               if t == 0.0
                  c = 1.0;
                  s = 0.0;
               else
                  c = h(i,i-1)/t;
                  s = h(i+1,i-1)/t;
                  h(i,i-1) = t;
                  h(i+1,i-1) = 0.0;
               end
            end
            
            %              apply rotation from left to rows of H
            
            G = [c s ; -s c];
            h(i:i+1,i:kend) = G* h(i:i+1,i:kend);
            
            %              apply rotation from right to columns of H
            
            ip2 = i+2;
            if ip2 > kend
               ip2 = kend;
            end
            
            h(1:ip2,i:i+1) =  h(1:ip2,i:i+1)*G';
            
            %              apply rotation from right to columns of V
            
            v(:,i:i+1) =  v(:,i:i+1)*G';
            
            %              accumulate e'  Q so residual can be updated
            %                          k+p
            
            q(i:i+1) =  G*q(i:i+1);
         end
      end
      num = num + 1;
   else
      
      %        Apply a double complex shift using 3 by 3 Householder 
      %        transformations
      
      if (jj == p | mark == 1) ,
         mark = 0;       % skip application of conjugate shift
      else
         num = num + 2;    % mark that a complex conjugate
         mark = 1;         % pair has been applied
         
         for ii = 1:nx-1   % Loop over blocks in H
            k1 = jx(ii)+1;
            k2 = k1+1;
            k3 = jx(ii+1);
            c = h(k1,k1)*h(k1,k1) + h(k1,k2)*h(k2,k1) ...
               - 2.0*wr(jj)*h(k1,k1);
            c = c + wr(jj)^2 + wi(jj)^2;
            s = h(k2,k1)*(h(k1,k1) + h(k2,k2) - 2.0*wr(jj));
            g = h(k2+1,k2)*h(k2,k1);
            t = norm([c s g]);
            sig = -sign(c);
            c = c -t*sig;
            for i = k1:k3
               if i > k1
                  t = norm(h(i:i+2,i-1));
                  sig = -sign(h(i,i-1));
                  c = h(i,i-1) - t*sig;
                  s = h(i+1,i-1);
                  g = h(i+2,i-1);
                  h(i,i-1) = t;
                  h(i+1,i-1) = 0.0;
                  h(i+2,i-1) = 0.0;
               end
               t = norm([c s g]);
               if t ~= 0.0
                  c = c/t;
                  s = s/t;
                  g = g/t;
               end
               z = [c s g]';
               
               %                 apply transformation from left to rows of H
               
               t =  sig*2.0*(z'*h(i:i+2,i:kend));
               h(i:i+2,i:kend) = sig*h(i:i+2,i:kend) - z*t;
               
               %                 apply transformation from right to columns of H
               
               ip3 = i+3;
               if ip3 > kend
                  ip3 = kend;
               end
               t =  sig*2.0*h(1:ip3,i:i+2)*z;
               h(1:ip3,i:i+2) = sig*h(1:ip3,i:i+2) - t*z';
               
               %                 apply transformation from right to columns of V
               
               t =  sig*2.0*v(:,i:i+2)*z;
               v(:,i:i+2) = sig*v(:,i:i+2) - t*z';
               
               %                 accumulate e'  Q so residual can be updated
               %                             k+p
               
               t =  sig*2.0*z'*q(i:i+2);
               q(i:i+2) = sig*q(i:i+2) - z*t;
            end      
         end
         
         %           clean up step with Givens rotation
         
         i = kend-1;
         if i > k1
            t = norm([h(i,i-1) h(i+1,i-1)]);
            if t ~= 0.0
               c = h(i,i-1)/t;
               s = h(i+1,i-1)/t;
            else
               c = 1.0;
               s = 0.0;
            end
            h(i,i-1) = t;
            h(i+1,i-1) = 0.0;
         end
         
         %           apply rotation from left to rows of H
         
         G = [c s ; -s c];
         h(i:i+1,i:kend) = G* h(i:i+1,i:kend);
         
         %           apply rotation from right to columns of H
         
         ip2 = i+2;
         if ip2 > kend
            ip2 = kend;
         end
         h(1:ip2,i:i+1) =  h(1:ip2,i:i+1)*G';
         
         %           apply rotation from right to columns of V
         
         v(:,i:i+1) =  v(:,i:i+1)*G';
         
         %           accumulate e'  Q so residual can be updated
         %                       k+p
         
         q(i:i+1) =  G*q(i:i+1);
      end
   end
end


%  update residual and store in the k+1 -st column of v


k = kend - num;
v(:,kend+1) = v(:,kend+1)*q(k);
if k < size(h,1)
   v(:,k+1) = v(:,kend+1) + v(:,k+1)*h(k+1,k);
end


% ===========================================================


function [v,h,defloc] =  arnold(k1,k2,A,B,v,h,bfactor)


%  ARNOLD  Computes or extends an Arnoldi factorization.
%
%  ARNOLD was designed to be called by EIGS in cases where the
%  shift parameter sigma is non-numeric.  This routine uses power
%  methods and not inverse iteration. 
%
%  Let V and H be the Arnoldi factors from a k1-step Arnoldi
%  process such that
%
%     AV - BVH = re' ,           where r is stored in
%            k1            the k1+1-st column of V.

%  [V,H,defloc] = ARNOLD(k1,k2,A,B,V,H,bfactor)
%
%  extends the existing factorization by k2 - k1 steps, resulting in
%
%          AV - V H = r e' .          where r is stored in
%            +   + +   + k2         the k2+1-st column of V.
%
%  On input and output, the columns of V should be B-orthogonal,
%  i.e. V'BV = I.
%
%  If B can be Cholesky factored as B = LL', where L is lower
%  triangular, then the factor L should be passed to ARNOLD in
%  the B position and the parameter bfactor should be set to true.
%
%  See also EIGS.


%  Dan Sorensen and Richard J. Radke, 11/95


if nargin < 7
   bfactor = 0;
end


if bfactor
   BL = B;
   B = [];
end


defloc = 0;


for j = k1+1:k2,
   jm1 = j-1; 
   jp1 = j+1;
   if bfactor
      vtBL = v(:,j)' * BL;
      beta = norm(vtBL);
    else
      beta = sqrt(v(:,j)'*B*v(:,j));
   end
   h(j,jm1) = beta;
   if (beta <= 10*eps*norm(h(1:jm1,jm1)))     % If beta is "small"
      v(:,j) = rand(size(v,1),1);             % we deflate by
      if bfactor
         s = v(:,1:jm1)'*BL*(v(:,j)'*BL)';
      else
         s = v(:,1:jm1)'*B*v(:,j);               % setting the
      end
      v(:,j) = v(:,j) - v(:,1:jm1)*s;         % corresponding
      if bfactor
         s = v(:,1:jm1)'*BL*(v(:,j)'*BL)';
      else
         s = v(:,1:jm1)'*B*v(:,j);               % subdiagonal of H
      end
      v(:,j) = v(:,j) - v(:,1:jm1)*s;         % equal to 0 and
      if bfactor
         vtBL = v(:,j)'*BL;
         beta = norm(vtBL);
      else
         beta = sqrt(v(:,j)'*B*v(:,j));          % starting the
      end
      beta = 1.0/beta;                        % basis for a new 
      h(j,jm1) = 0.0;                         % invariant subspace.
      defloc = j;
   else
      beta = 1.0/beta;
   end
   v(:,j) = v(:,j)*beta;
   
   %     Compute w = Av and store w in j+1 -st col of V
   
   if bfactor
      v2 = BL \ v(:,j);
   else
      v2 = v(:,j);
   end
   
   if isstr(A)
      w = feval(A,v2);
   else
      w = A*v2;
   end
   
   if bfactor
      v(:,jp1) = BL' \ w;
   else
      v(:,jp1) = B \ w;
   end
   
   %     Compute the next (j-th) column of H
   
   if bfactor
      h(1:j,j) = v(:,1:j)'*BL*(v(:,jp1)'*BL)';
   else
      h(1:j,j) = v(:,1:j)'*B*v(:,jp1);
   end
   
   %     Compute the residual in the j+1 -st col of V
   
   v(:,jp1) = v(:,jp1) - v(:,1:j)*h(1:j,j);
   
   %     Perform one step of iterative refinement to correct the orthogonality
   
   if bfactor
      s = v(:,1:j)'*BL*(v(:,jp1)'*BL)';
     else
      s = v(:,1:j)'*B*v(:,jp1);
   end
   v(:,jp1) = v(:,jp1) - v(:,1:j)*s;
   h(1:j,j) = h(1:j,j) + s;
end


% ==============================================================


function [v,h,defloc] =  arnoldnob(k1,k2,A,v,h)


%  ARNOLDNOB  Computes or extends an Arnoldi factorization.
%
%  ARNOLDNOB is exactly the same as ARNOLD when B = I.
%


defloc = 0;


for j = k1+1:k2,
   jm1 = j-1; 
   jp1 = j+1;
   beta = norm(v(:,j)); 
   h(j,jm1) = beta;
   if (beta <= 10*eps*norm(h(1:jm1,jm1)))    % If beta is "small"
      v(:,j) = rand(size(v,1),1);            % we deflate by
      s = v(:,1:jm1)'*v(:,j);                % setting the
      v(:,j) = v(:,j) - v(:,1:jm1)*s;        % corresponding
      s = v(:,1:jm1)'*v(:,j);                % subdiagonal of H
      v(:,j) = v(:,j) - v(:,1:jm1)*s;        % equal to 0 and
      beta = norm(v(:,j));                   % starting the
      beta = 1.0/beta;                       % basis for a new 
      h(j,jm1) = 0.0;                        % invariant subspace.
      defloc = j;
   else
      beta = 1.0/beta;
   end
   v(:,j) = v(:,j)*beta;
   
   %     Compute w = Av and store w in j+1 -st col of V
   
   v2 = v(:,j);
   
   if isstr(A)
      v(:,jp1) = feval(A,v2);
   else
      v(:,jp1) = A*v2;
   end
   
   %     Compute the next (j-th) column of H
   
   h(1:j,j) = v(:,1:j)'*v(:,jp1);
   
   %     Compute the residual in the j+1 -st col of V
   
   v(:,jp1) = v(:,jp1) - v(:,1:j)*h(1:j,j);
   
   %     Perform one step of iterative refinement to correct the orthogonality
   
   s = v(:,1:j)'*v(:,jp1);
   v(:,jp1) = v(:,jp1) - v(:,1:j)*s;
   h(1:j,j) = h(1:j,j) + s;
end


% ===============================================================


function [v,h,defloc] =  arnold2(k1,k2,L,U,B,v,h,tol);


%  ARNOLD2  Computes or extends an Arnoldi factorization.
%
%  ARNOLD2 was designed to be called by EIGS in cases where the
%  shift parameter sigma is numeric.  This routine uses inverse
%  iteration and the LU factors of (A - sigma*B).
%
%  Let V and H be the Arnoldi factors from a k1-step Arnoldi
%  process such that
%
%    inv(L*U) B V - V H = re'  ,     where r is stored in
%                 k1       the k1+1-st column of V.

%  [V,H,defloc] = ARNOLD2(k1,k2,L,U,B,V,H,tol)
%  extends the existing factorization by k2 - k1 steps, resulting in
%
%         inv(L*U) B V - V H = r e' .     where r is stored in
%                     +   + +   + k2      the k2+1-st column of V.
%
%  On input and output, the columns of V should be B-orthogonal,
%  i.e. V'BV = I.
%
%  The parameter tol provides a tolerance for deciding when to
%  deflate the problem.  If the problem can be deflated, a pointer
%  to the appropriate location in H of the start of the active basis
%  is returned in the output argument defloc.
%
%  If the L|U factors of A-sigma*B are ill-conditioned and return
%  unusable U\(L\v), then v and h remain unchanged and defloc is
%  assigned -1.
%
%  Dan Sorensen and Richard J. Radke, 11/95


defloc = 0;
vold = v;
hold = h;


for j = k1+1:k2,
   jm1 = j-1; 
   jp1 = j+1;
   beta = sqrt(v(:,j)'*B*v(:,j));
   h(j,jm1) = beta;
   if (beta <= tol)
      v(j,j) = 1.0;
      s = v(:,1:jm1)'*B*v(:,j);
      v(:,j) = v(:,j) - v(:,1:jm1)*s;
      s = v(:,1:jm1)'*B*v(:,j); 
      v(:,j) = v(:,j) - v(:,1:jm1)*s;
      beta = sqrt(v(:,j)'*B*v(:,j));
      beta = 1.0/beta;
      defloc = j;
   else
      beta = 1.0/beta;
   end
   v(:,j) = v(:,j)*beta;
   
   %     Compute w = Av and store w in j+1 -st col of V
   
   v(:,jp1) = U \ (L \ (B*v(:,j)));
   
   if sum(~isfinite(v(:,jp1)))
      v = vold;
      h = hold;
      defloc = -1;
      return
   end
   
   %     Compute the next (j-th) column of H
   
   h(1:j,j) = v(:,1:j)'*B*v(:,jp1);
   
   %     Compute the residual in the j+1 -st col of V
   
   v(:,jp1) = v(:,jp1) - v(:,1:j)*h(1:j,j);
   
   %     One step of iterative refinement to correct the orthogonality
   
   s = v(:,1:j)'*B*v(:,jp1);
   v(:,jp1) = v(:,jp1) - v(:,1:j)*s;
   h(1:j,j) = h(1:j,j) + s;
end


% ================================================================


function [v,h,defloc] =  arnold2nob(k1,k2,L,U,v,h,tol);


%  ARNOLD2NOB  Computes or extends an Arnoldi factorization.
%
%  ARNOLD2NOB is the same as ARNOLD2 when B = I.
%


defloc = 0;
vold = v;
hold = h;


for j = k1+1:k2,
   jm1 = j-1; 
   jp1 = j+1;
   beta = norm(v(:,j));
   h(j,jm1) = beta;
   if (beta <= tol)
      v(j,j) = 1.0;
      s = v(:,1:jm1)'*v(:,j);
      v(:,j) = v(:,j) - v(:,1:jm1)*s;
      s = v(:,1:jm1)'*v(:,j); 
      v(:,j) = v(:,j) - v(:,1:jm1)*s;
      beta = sqrt(v(:,j)'*v(:,j));
      beta = 1.0/beta;
      defloc = j;
   else
      beta = 1.0/beta;
   end
   v(:,j) = v(:,j)*beta;
   
   %     Compute w = Av and store w in j+1 -st col of V
   
   v(:,jp1) = U \ (L \ v(:,j));
   
   if sum(~isfinite(v(:,jp1)))
      v = vold;
      h = hold;
      defloc = -1;
      return
   end
   
   %     Compute the next (j-th) column of H
   
   h(1:j,j) = v(:,1:j)'*v(:,jp1);
   
   %     Compute the residual in the j+1 -st col of V
   
   v(:,jp1) = v(:,jp1) - v(:,1:j)*h(1:j,j);
   
   %     One step of iterative refinement to correct the orthogonality
   
   s = v(:,1:j)'*v(:,jp1);
   v(:,jp1) = v(:,jp1) - v(:,1:j)*s;
   h(1:j,j) = h(1:j,j) + s;
end


% ==============================================================


function w = accpoly(v)


%   ACCPOLY applies an accelerant polynomial in the operator
%   specified by the global variable Chebystruct.op 
%       to the input vector v.
%
%   ACCPOLY was designed to be called by EIGS when the input argument
%       A to EIGS is an m-file specifying a symmetric matrix-vector
%       product.
%
%   Approximate lower and upper bounds on the eigenvalues of op
%   should be present in the global variables lbd and ubd.
%   The global variable sig is 'SO' or 'LO' depending on whether
%   the smallest or largest eigenvalues of op are desired.
%   sig can also be a numeric shift.
%
%   The accelerant polynomial in ChebyStruct.op is either a scaled 
%       Chebyshev polynomial over [lbd,ubd] or an equiripple single-peak
%       polynomial obtained by the filter design code contained in bp_fer.


%   Richard J. Radke, 3/96.


global ChebyStruct


op = ChebyStruct.op;
lbd = ChebyStruct.lbd;
ubd = ChebyStruct.ubd;
sig = ChebyStruct.sig;
filtpoly = ChebyStruct.filtpoly;
iter = ChebyStruct.iter;


p = .25;  %  1-p is the fraction of the spectrum that gets mapped
%  into the equiripple region (sig = string)
m = 10;


if ~isstr(sig)
   slope = 2/(ubd-lbd);
   intercept = 1 - (ubd*slope);
   if sig > ubd
      sig = ubd-.01*(ubd-lbd);
   elseif sig < lbd
      sig = lbd+.01*(ubd-lbd);
   end
elseif strcmp(sig,'SO') | strcmp(sig,'SR')
   slope = 2/(ubd-lbd-p*(ubd-lbd));slope = 2/(ubd-p*(ubd-lbd)-lbd);
   intercept = -1 - (lbd*slope);
   intercept = 1 - (ubd*slope);
elseif strcmp(sig,'LO') | strcmp(sig,'LR')
   slope = 2/(ubd-p*(ubd-lbd)-lbd);
   intercept = -1 - (lbd*slope);
end


if isstr(sig)
   
   w0 = v;
   w1 = slope*feval(op,v)+intercept*v;
   
   for jj = 2:m
      w = 2*(slope*feval(op,w1)+intercept*w1)-w0;
      w0 = w1;
      w1 = w;
   end
   
else
   
   if iter == 1
      
      wp = acos(slope*sig+intercept);
      m = 10;
      N = 2*m+1;
      L = 4;
      del = 0.01;
      h = bp_fer(N,L,wp,del,-del,2^7);
      ChebyStruct.filtpoly = h2x(h);
      filtpoly = ChebyStruct.filtpoly;
      iter = 2;
      
   end
   
   w = filtpoly(1)*v;
   
   for i=2:(m+1)
      w = (slope*feval(op,w)+intercept*w) + filtpoly(i)*v;
   end
   
end


% ========================================================


function [w,qq] = shftit(h, kstart, kend, sigma)


%  SHFTIT  Calculate shifts to update an Arnoldi factorization.
%
%  SHFTIT was designed to be called by EIGS.
%
%  Syntax: [w,q] = SHFTIT(H, kstart, kend, sigma) where
%
%  H is an upper Hessenberg matrix (from the Arnoldi factorization
%     AV = VH + fe_k'),
%
%  kstart points to the start of the active block in H,
%
%  kend points to the end of the active block in H,
%
%  sigma is a numeric shift or one of 'LR','SR','LM','SM','BE'.
%     (If sigma is a number, sigma = 'LM' is used, since
%     the problem has already been shifted by sigma.  The
%     real eigenvalues are recovered later.)
%
%  SHFTIT calculates [q,w] = eig(The active block of H), where
%  the eigenvalues and eigenvectors are reordered according to sigma,
%  with the eigenvalues to use as shifts put first.


%  Dan Sorensen and Richard J. Radke, 7/95 


[q,ww] = eig(h(kstart:kend,kstart:kend));


w = diag(ww);


k = kend-kstart+1;


%  select filter mechanism by activating appropriate choice below


if ~isstr(sigma)
   
   [s,ir] = sort(abs(w));
   
elseif strcmp(sigma,'BE')
   
   %     sort on real part, alternating from high to low end of
   %     the spectrum
   
   par = rem(k,2);
   
   [s,ir] = sort(-real(w));
   
   ix(1:2:k-1+par) = 1:ceil(k/2);
   ix(2:2:k-par) = k:-1:(ceil(k/2)+1);
   
   ir = flipud(ir(ix));
   
elseif strcmp(sigma,'LO')
   
   %     keep the k-1 largest and 1 smallest e-vals
   %     at the end of the sort
   
   [s,ir] = sort(real(w));
   ir = ir([2:k,1]);
   
elseif strcmp(sigma,'SO')
   
   %     keep the k-1 smallest and 1 largest e-vals
   %     at the end of the sort
   
   [s,ir] = sort(-real(w));
   ir = ir([2:k,1]);
   
elseif strcmp(sigma,'LR')
   
   %     sort for largest real part (shifts are smallest real part)
   
   [s,ir] = sort(real(w));
   
elseif strcmp(sigma,'SR')
   
   %     sort for smallest real part (shifts are largest real part)
   
   [s,ir] = sort(-real(w));
   
elseif strcmp(sigma,'SM')
   
   %     sort for smallest absolute val (shifts are largest abs val)
   
   [s,ir] = sort(-abs(w));
   
else
   
   %     sort for largest absolute val (shifts are smallest abs val)
   
   [s,ir] = sort(abs(w));
   
end


%  apply sort to w


w = w(ir);
qq = q(:,ir);


% =====================================================


function p = add_poly(p1,p2)
%
%   function p = add_poly(p1,p2);
%   Adds 2 polynomials.
%
l1 = length(p1); l2 = length(p2);
if l1 == 0, p = rlz(p2); break; end
if l2 == 0, p = rlz(p1); break; end
p1 = rlz(p1(:)); p2 = rlz(p2(:));
l1 = length(p1); l2 = length(p2);
if l1 > l2
   p2 = [zeros(l1-l2,1); p2];
else
   p1 = [zeros(l2-l1,1); p1]; 
end
p = rlz((p1 + p2).');


% =======================================================


function [h,h2,rs] = bp_fer(N,L,wp,us,ls,g)
% A program for the design of linear phase bandpass FIR filters with a
% Flat monotonically decreasing Pass band (on either side of wp)
% and  an EquiRipple Stop band.
%       With this program, the user can specify the stop band ripple size
%       but not the passband frequency of flatness.
%
% N  : length of total filter
% L  : degree of flatness
% wp : passband frequency of flatness
% us : upper bound in stop band
% ls : lower boudn in stop band
% g  : grid size
%
% EXAMPLE
%    N  =  25;
%    L  =  8;
%    wp  = .4*pi;
%    us =  0.01;
%    ls = -us;
%    g  =  2^8;
%    [h,h2,rs] = bp_fer(N,L,wp,us,ls,g);
% or
%    N = 31; L = 16; 
%    bp_fer(N,L,wp,us,ls,g);




if (rem(N,2) == 0) | (rem(L,4) ~= 0)
   disp('N must be odd and L must be divisible by 4')
   return
else
   SN = 1e-9;                                      % SN : SMALL NUMBER
   q  = (N-L+1)/2;
   w  = [0:g]'*pi/g;
   ws1 = wp*0.8;
   a = 5; b = 1;
   ws2 = (a*wp+b*pi)/(a+b);
   
   d = ws1/(pi-ws2);            % q1 : number of ref. set freq. in 1st stopband
   q1 = round(q/(1+1/d));       % q2 : number of ref. set freq. in 2nd stopband
   if q1 == 0
      q1 = 1;
   elseif q1 == q
      q1 = q - 1;
   end
   
   q2 = q - q1;
   
   if q1 == 1;
      rs1 = ws1;
   else
      rs1 = [0:q1-1]'*(ws1/(q1-1));
   end
   if q2 == 1
      rs2 = ws2;
   else
      rs2 = [0:q2-1]'*(pi-ws2)/(q2-1)+ws2;
   end
   rs = [rs1; rs2];
   
   Y1 = [ls*(1-(-1).^(q1:-1:1))/2 + us*((-1).^(q1:-1:1)+1)/2]';
   Y2 = [ls*(1-(-1).^(1:q2))/2 + us*((-1).^(1:q2)+1)/2]';
   Y = [Y1; Y2];
   
   n  = 0:q-1;
   Z   = zeros(2*(g-q)+1,1);
   % A1  = (-1)^(L/2) * (sin(w/2-wp/2).*sin(w/2+wp/2)).^(L/2);
   % A1r = (-1)^(L/2) * (sin(rs/2-wp/2).*sin(rs/2+wp/2)).^(L/2);
   A1  = (-1)^(L/2) * ((cos(wp)-cos(w))/2).^(L/2);
   A1r = (-1)^(L/2) * ((cos(wp)-cos(rs))/2).^(L/2);
   it = 0;
   while 1 & (it < 35)
      if length(rs) ~= length(n)
         rs ,n
         error('Filter design code error')
      end
      a2 = cos(rs*n)\[(Y-1)./A1r];
      A2 = real(fft([a2(1);a2(2:q)/2;Z;a2(q:-1:2)/2])); A2 = A2(1:g+1);
      A  = 1 + A1.*A2;
      %   figure(1), plot(w/pi,A),
      %        hold on, plot(rs/pi,Y,'o'), hold off, % pause(.1)
      %   figure(2), plot(w/pi,20*log10(abs(A))),
      %        hold on, plot(rs/pi,20*log10(abs(Y)),'o'), hold off, % pause(.1)
      %        pause(.5)
      ri = sort([localmax(A); localmax(-A)]);
      lri = length(ri);
      % ri(1:length(ri)-q) = [];
      if lri ~= q+1
         if abs(A(ri(lri))-A(ri(lri-1))) < abs(A(ri(1))-A(ri(2)))
            ri(lri) = [];
         else
            ri(1) = [];
         end
      end
      rs = (ri-1)*pi/g;
      [temp, k] = min(abs(rs - wp)); rs(k) = [];
      q1 = sum(rs < wp);
      q2 = sum(rs > wp);
      
      Y1 = [ls*(1-(-1).^(q1:-1:1))/2 + us*((-1).^(q1:-1:1)+1)/2]';
      Y2 = [ls*(1-(-1).^(1:q2))/2 + us*((-1).^(1:q2)+1)/2]';
      Y = [Y1; Y2];
      
      % rs = refine2(a2,L/2,rs);
      % A1r = (-1)^(L/2) * (sin(rs/2-wp/2).*sin(rs/2+wp/2)).^(L/2);
      A1r = (-1)^(L/2) * ((cos(wp)-cos(rs))/2).^(L/2);
      Ar = 1 + (cos(rs*n)*a2) .* A1r;
      Err = max([max(Ar)-us, ls-min(Ar)]);
      %   fprintf(1,'    Err = %18.15f\n',Err);
      if Err < SN, break, end
      it = it + 1;
   end
   
   h2 = [a2(q:-1:2)/2; a2(1); a2(2:q)/2];
   h = h2;
   for k = 1:L/2
      h = conv(h,[1 -2*cos(wp) 1]')/4;
   end
   h((N+1)/2) = h((N+1)/2) + 1;
   
end


% ========================================================


function C = chebpoly(n)


% C = cheb_poly(n)
% Chebychev polynomial
%


if n == 0
   C = 1;
elseif n == 1
   C = [1 0];
else
   A = 1;
   B = [1 0];
   for k = 2:n
      C = 2*[B 0] - [0 0 A];
      A = B;
      B = C;
   end
end


% =========================================================


function b = cos2x(a)
%
% converts the cos polynomial, 
%       a(f) = a(1) + a(2)*cos(w) + ... + a(m+1)*cos(m*w)
% over [0,1/2],
% to the polynomial
%       b(x) = b(m+1) + b(m)*x + ... + b(1)*x^m
% over [-1,1]
%
% the transformation is : x = cos(w)
%
% (x == 1) and (w == 0) are mapped together
% (x == -1) and (w == pi) are mapped together
%


m = length(a)-1;
c = a(1);
for k = 2:m+1
   c = add_poly(c,a(k)*chebpoly(k-1));
end
b = zeros(size(a));
b(:) = c;


% =============================================


function a = h2cos(h)
% a = h2cos(h);


N = length(h);
if sum(abs(h-h(N:-1:1))) > 0.00001
   disp('for symmetric h only')
   a = [];
   return
end
if rem(N,2)==1
   % N even
   a = 2*h((N+1)/2:N);
   a(1) = a(1)/2;
else
   disp('for odd length only')
   return
end


% ==================================================


function x = h2x(h)
%
% x = h2x(h)
%
x = cos2x(h2cos(h));


% ====================================================


function k = localmax(x)
% k = localmax(x)
% finds location of local maxima
%
s = size(x); x = [x(:)].'; N = length(x);
b1 = x(1:N-1)<=x(2:N); b2 = x(1:N-1)>x(2:N);
k = find(b1(1:N-2)&b2(2:N-1))+1;
if x(1)>x(2), k = [k, 1]; end
if x(N)>x(N-1), k = [k, N]; end
k = sort(k); if s(2) == 1, k = k'; end


% =============================================


function p = rlz(p)


if isempty(p)
   break
end


while p(1) == 0 
   p(1) = [];
   if isempty(p)
      break
   end
end


% =================================================

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&7

function  varargout = eigs(varargin)
%EIGS  Find a few eigenvalues and eigenvectors of a matrix using ARPACK.
%   D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues.
%   A must be square and should be large and sparse.
%
%   [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude
%   eigenvalues and a matrix V whose columns are the corresponding eigenvectors.
%
%   [V,D,FLAG] = EIGS(A) also returns a convergence flag.  If FLAG is 0
%   then all the eigenvalues converged; otherwise not all converged.
%
%   EIGS(AFUN,N) accepts the function AFUN instead of the matrix A.
%   Y = AFUN(X) should return Y = A*X.  N is the size of A.  The matrix A
%   represented by AFUN is assumed to be real and nonsymmetric.  In all these
%   calling sequences, EIGS(A,...) may be replaced by EIGS(AFUN,N,...).
%
%   EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D.  B must
%   be symmetric (or Hermitian) positive definite and the same size as A.
%   EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D.
%
%   EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues.
%
%   EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues based on SIGMA:
%      'LM' or 'SM' - Largest or Smallest Magnitude
%   For real symmetric problems, SIGMA may also be:
%      'LA' or 'SA' - Largest or Smallest Algebraic
%      'BE' - Both Ends, one more from high end if K is odd
%   For nonsymmetric and complex problems, SIGMA may also be:
%      'LR' or 'SR' - Largest or Smallest Real part
%      'LI' or 'SI' - Largest or Smallest Imaginary part
%   If SIGMA is a real or complex scalar, EIGS finds the eigenvalues closest to
%   SIGMA.  In this case, B need only be symmetric (or Hermitian) positive
%   semi-definite and the function Y = AFUN(X) must return Y = (A-SIGMA*B)\X.
%
%   EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options:
%   OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{0} | 1]
%   OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [0 | {1}]
%   OPTS.tol: convergence: abs(d_comp-d_true) < tol*abs(d_comp) [scalar | {eps}]
%   OPTS.maxit: maximum number of iterations [integer | {300}]
%   OPTS.p: number of Lanczos vectors: K+1<p<=N [integer | {2K}]
%   OPTS.v0: starting vector [N-by-1 vector | {randomly generated by ARPACK}]
%   OPTS.disp: diagnostic information display level [0 | {1} | 2]
%   OPTS.cholB: B is actually its Cholesky factor CHOL(B) [{0} | 1]
%   OPTS.permB: sparse B is actually CHOL(B(permB,permB)) [permB | {1:N}]
%
%   EIGS(AFUN,N,K,SIGMA,OPTS,P1,...) and EIGS(AFUN,N,B,K,SIGMA,OPTS,P1,...)
%   provide for additional arguments which are passed to AFUN(X,P1,...).
%
%   Examples:
%      A = delsq(numgrid('C',15));  d1 = eigs(A,5,'SM');
%   Equivalently, if dnRk is the following one-line function:
%      function y = dnRk(x,R,k)
%      y = (delsq(numgrid(R,k))) * x;
%   then pass dnRk's additional arguments, 'C' and 15, to eigs:
%      n = size(A,1);  opts.issym = 1;  d2 = eigs(@dnRk,n,5,'SM',opts,'C',15);
%
%   See also EIG, SVDS, ARPACKC.


%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 1.38 $  $Date: 2000/08/02 21:50:54 $


cputms = zeros(5,1);
t0 = cputime; % start timing pre-processing


if (nargout > 3)
    error('Too many output arguments')
end


% Process inputs and do error-checking
if isa(varargin{1},'double')
    A = varargin{1};
    Amatrix = 1;
else
A = fcnchk(varargin{1});
    Amatrix = 0;
end


isrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma)
if Amatrix
    isrealprob = isreal(A);
end


issymA = 0;
if Amatrix
    issymA = isequal(A,A');
end


if Amatrix
    [m,n] = size(A);
    if (m ~= n)
        error('A must be a square matrix or a function which computes A*x')
    end
else
    n = varargin{2};
    nstr = sprintf('Size of problem n must be a positive integer');
    if ~isequal(size(n),[1,1]) | ~isreal(n)
        error(nstr)
    end
    if (round(n) ~= n)
        warning(nstr)
        n = round(n);
    end
    if issparse(n)
        n = full(n);
    end
end


Bnotthere = 0;
Bstr = sprintf(['Generalized matrix B must be the same size as A and' ...
        ' either a symmetric positive (semi-)definite matrix or' ...
        ' its Cholesky factor']);
if (nargin < (3-Amatrix-Bnotthere))
    B = [];
    Bnotthere = 1;
else
    Bk = varargin{3-Amatrix-Bnotthere};
    if isempty(Bk) % allow eigs(A,[],k,sigma,opts);
        B = Bk;
    else
        if isequal(size(Bk),[1,1]) & (n ~= 1)
            B = [];
            k = Bk;
            Bnotthere = 1;
        else % eigs(9,8,...) assumes A=9, B=8, ... NOT A=9, k=8, ...
            B = Bk;
            if ~isa(B,'double') | ~isequal(size(B),[n,n])
                error(Bstr)
            end
            isrealprob = isrealprob & isreal(B);
        end
    end
end


if Amatrix & ((nargin - ~Bnotthere)>4)
    error('Too many inputs')
end


if (nargin < (4-Amatrix-Bnotthere))
    k = min(n,6);
else
    k = varargin{4-Amatrix-Bnotthere};
end


kstr = sprintf(['Number of eigenvalues requested k must be a' ...
        ' positive integer <= n']);
if ~isa(k,'double') | ~isequal(size(k),[1,1]) | ~isreal(k) | (k>n)
    error(kstr)
end
if issparse(k)
    k = full(k);
end
if (round(k) ~= k)
    warning(kstr)
    k = round(k);
end


whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string']);
if (nargin < (5-Amatrix-Bnotthere))
    whch = 'LM';
    sigma = 0;
else
    whch = varargin{5-Amatrix-Bnotthere};
    if isstr(whch)
        if ~isequal(size(whch),[1,2])
            error(whchstr)
        end
        whch = upper(whch);
        sigma = 0;
    else
        sigma = whch;
        if ~isa(sigma,'double') | ~isequal(size(sigma),[1,1])
            error('Eigenvalue shift sigma must be a scalar')
        end
        if issparse(sigma)
            sigma = full(sigma);
        end
        isrealprob = isrealprob & isreal(sigma);
        if (sigma == 0)
            whch = 'SM';
        else
            whch = 'LM';
        end
    end
end


tol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS)
maxit = [];
p = [];
info = int32(0); % use a random starting vector
display = 1;
cholB = 0;


if (nargin >= (6-Amatrix-Bnotthere))
    opts = varargin{6-Amatrix-Bnotthere};
    if ~isa(opts,'struct')
        error('Options argument must be a structure')
    end
    
    if isfield(opts,'issym') & ~Amatrix
        issymA = opts.issym;
        if (issymA ~= 0) & (issymA ~= 1)
            error('opts.issym must be 0 or 1')
        end
    end
    
    if isfield(opts,'isreal') & ~Amatrix
        if (opts.isreal ~= 0) & (opts.isreal ~= 1)
            error('opts.isreal must be 0 or 1')
        end
        isrealprob = isrealprob & opts.isreal;
    end
    
    if ~isempty(B) & (isfield(opts,'cholB') | isfield(opts,'permB'))
        if isfield(opts,'cholB')
            cholB = opts.cholB;
            if (cholB ~= 0) & (cholB ~= 1)
                error('opts.cholB must be 0 or 1')
            end
            if isfield(opts,'permB')
                if issparse(B) & cholB
                    permB = opts.permB;
                    if ~isequal(sort(permB),(1:n)) & ...
                            ~isequal(sort(permB),(1:n)')
                        error('opts.permB must be a permutation of 1:n')
                    end
                else
                    warning(['Ignoring opts.permB since B is not its sparse' ...
                            ' Cholesky factor'])
                end
            else
                permB = 1:n;
            end
        end
    end
    
    if isfield(opts,'tol')
        if ~isequal(size(opts.tol),[1,1]) | ~isreal(opts.tol) | (opts.tol<=0)
            error(['Convergence tolerance opts.tol must be a strictly' ...
                    ' positive real scalar'])
        else
            tol = full(opts.tol);
        end
    end
    
    if isfield(opts,'p')
        p = opts.p;
        pstr = sprintf(['Number of basis vectors opts.p must be a positive' ...
                ' integer <= n']);
        if ~isequal(size(p),[1,1]) | ~isreal(p) | (p<=0) | (p>n)
            error(pstr)
        end
        if issparse(p)
            p = full(p);
        end
        if (round(p) ~= p)
            warning(pstr)
            p = round(p);
        end
    end
    
    if isfield(opts,'maxit')
        maxit = opts.maxit;
        str = sprintf(['Maximum number of iterations opts.maxit must be' ...
                ' a positive integer']);
        if ~isequal(size(maxit),[1,1]) | ~isreal(maxit) | (maxit<=0)
            error(str)
        end
        if issparse(maxit)
            maxit = full(maxit);
        end
        if (round(maxit) ~= maxit)
            warning(str)
            maxit = round(maxit);
        end
    end
    
    if isfield(opts,'v0')
        if ~isequal(size(opts.v0),[n,1])
            error('Start vector opts.v0 must be n-by-1')
        end
        if isrealprob
            if ~isreal(opts.v0)
                error(['Start vector opts.v0 must be real for real problems'])
            end
            resid = full(opts.v0);
        else
            resid(1:2:(2*n-1),1) = full(real(opts.v0));
            resid(2:2:2*n,1) = full(imag(opts.v0));
        end
        info = int32(1); % use resid as starting vector
    end
    
    if isfield(opts,'disp')
        display = opts.disp;
        dispsyt = sprintf('Diagnostic level opts.disp must be an integer');
        if (~isequal(size(display),[1,1])) | (~isreal(display)) | (display<0)
            error(dispstr)
        end
        if (round(display) ~= display)
            warning(dispstr)
            display = round(display);
        end
    end
    
    if isfield(opts,'cheb')
        warning(['Ignoring polynomial acceleration opts.cheb' ...
                ' (no longer an option)']);
    end
    if isfield(opts,'stagtol')
        warning(['Ignoring stagnation tolerance opts.stagtol' ...
                ' (no longer an option)']);
    end
    
end


% Now we know issymA, isrealprob, cholB, and permB


if isempty(p)
    if isrealprob & ~issymA
        p = min(max(2*k+1,20),n);
    else
        p = min(max(2*k,20),n);
    end
end
if isempty(maxit)
    maxit = max(300,ceil(2*n/max(p,1)));
end
if (info == int32(0))
    if isrealprob
        resid = zeros(n,1);
    else
        resid = zeros(2*n,1);
    end
end


if ~isempty(B) % B must be symmetric (Hermitian) positive (semi-)definite
    if cholB
        if ~isequal(triu(B),B)
            error(Bstr)
        end
    else
        if ~isequal(B,B')
            error(Bstr)
        end
    end
end


useeig = 0;
if isrealprob & issymA
    knstr = sprintf(['For real symmetric problems, must have number' ...
            ' of eigenvalues k < n.\n']);
else
    knstr = sprintf(['For nonsymmetric and complex problems, must have' ...
            ' number of eigenvalues k < n-1.\n']);
end
if isempty(B)
    knstr = [knstr sprintf(['Using eig(full(A)) instead.'])];
else
    knstr = [knstr sprintf(['Using eig(full(A),full(B)) instead.'])];
end
if (k == 0)
    useeig = 1;
end
if isrealprob & issymA
    if (k > n-1)
        if (n >= 6)
            warning(knstr)
        end
        useeig = 1;
    end
else
    if (k > n-2)
        if (n >= 7)
            warning(knstr)
        end
        useeig = 1;
    end
end


if isrealprob & issymA
    if ~isreal(sigma)
        error(['For real symmetric problems, eigenvalue shift sigma must' ...
                ' be real'])
    end
else
    if ~isrealprob & issymA & ~isreal(sigma)
        warning(['Complex eigenvalue shift sigma on a Hermitian problem' ...
                ' (all real eigenvalues)'])
    end
end


if isrealprob & issymA
    if strcmp(whch,'LR')
        whch = 'LA';
        warning(['For real symmetric problems, sigma value ''LR''' ...
                ' (Largest Real) is now ''LA'' (Largest Algebraic)'])
    end
    if strcmp(whch,'SR')
        whch = 'SA';
        warning(['For real symmetric problems, sigma value ''SR''' ...
                ' (Smallest Real) is now ''SA'' (Smallest Algebraic)'])
    end
    if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'})
        error(whchstr)
    end
else
    if strcmp(whch,'BE')
        warning(['Sigma value ''BE'' is now only available for real' ...
                ' symmetric problems.  Computing ''LM'' eigenvalues instead.'])
        whch = 'LM';
    end
    if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'})
        error(whchstr)
    end
end


% Now have enough information to do early return on cases eigs does not handle
if useeig
    if (nargout <= 1)
        varargout{1} = eigs2(A,n,B,k,whch,sigma,cholB, ...
            varargin{7-Amatrix-Bnotthere:end});
    else
        [varargout{1},varargout{2}] = eigs2(A,n,B,k,whch,sigma,cholB, ...
            varargin{7-Amatrix-Bnotthere:end});
    end
    if (nargout == 3)
        varargout{3} = 0; % flag indicates "convergence"
    end
    return
end


if isrealprob & ~issymA
    sigmar = real(sigma);
    sigmai = imag(sigma);
end


if isrealprob & issymA
    if (p <= k)
        error(['For real symmetric problems, must have number of' ...
                ' basis vectors opts.p > k'])
    end
else
    if (p <= k+1)
        error(['For nonsymmetric and complex problems, must have number of' ...
                ' basis vectors opts.p > k+1'])
    end
end


% Determine mode
if isequal(whch,'LM') & ~isequal(sigma,0)
    mode = 3;
elseif isempty(B)
    mode = 1;
else
    mode = 2;
end


if cholB
    pB = 0;
    RB = B;
    RBT = B';
end


if (mode == 3) & Amatrix % need lu(A-sigma*B)
    if issparse(A) & (isempty(B) | issparse(B))
        if isempty(B)
            AsB = A - sigma * speye(n);
        else
            if cholB
                AsB = A - sigma * RBT * RB;
            else
                AsB = A - sigma * B;
            end
        end
        perm = colmmd(AsB);
        [L,U,P] = lu(AsB(:,perm));
    else
        if isempty(B)
            AsB = A - sigma * eye(n);
        else
            if cholB
                AsB = A - sigma * RBT * RB;
            else
                AsB = A - sigma * B;
            end
        end
        [L,U,P] = lu(AsB);
    end
    dU = diag(U);
    rcondestU = double(min(abs(dU)) / max(abs(dU)));
    if (rcondestU < eps)
        if isempty(B)
            ds = sprintf(['(A-sigma*I) has small reciprocal condition' ...
                    ' estimate: %f\n'],rcondestU);
        else
            ds = sprintf(['(A-sigma*B) has small reciprocal condition' ...
                    ' estimate: %f\n'],rcondestU);
        end
        ds = [ds sprintf(['indicating that sigma is near an exact' ...
                    ' eigenvalue. The\nalgorithm may not converge unless' ...
                    ' you try a new value for sigma.\n'])];
        disp(ds)
        pause(2)
    end
end


if (mode == 2) & ~cholB % need chol(B)
    if issparse(B)
        permB = symmmd(B);
        [RB,pB] = chol(B(permB,permB));
    else
        [RB,pB] = chol(B);
    end
    if (pB == 0)
        RBT = RB';
    else
        error(Bstr)
    end
end


% Allocate outputs and ARPACK work variables
if isrealprob
    if issymA % real and symmetric
        prefix = 'ds';
        v = zeros(n,p);
        ldv = int32(size(v,1));
        ipntr = int32(zeros(15,1));
        workd = zeros(n,3);
        lworkl = p*(p+8);
        workl = zeros(lworkl,1);
        lworkl = int32(lworkl);
        d = zeros(k,1);
    else % real but not symmetric
        prefix = 'dn';
        v = zeros(n,p);
        ldv = int32(size(v,1));
        ipntr = int32(zeros(15,1));
        workd = zeros(n,3);
        lworkl = 3*p*(p+2);
        workl = zeros(lworkl,1);
        lworkl = int32(lworkl);
        workev = zeros(3*p,1);
        d = zeros(k+1,1);
        di = zeros(k+1,1);
    end
else % complex
    prefix = 'zn';
    zv = zeros(2*n*p,1);
    ldv = int32(n);
    ipntr = int32(zeros(15,1));
    workd = complex(zeros(n,3));
    zworkd = zeros(2*prod(size(workd)),1);
    lworkl = 3*p^2+5*p;
    workl = zeros(2*lworkl,1);
    lworkl = int32(lworkl);
    workev = zeros(2*2*p,1);
    zd = zeros(2*(k+1),1);
    rwork = zeros(p,1);
end


ido = int32(0); % reverse communication parameter
if isempty(B)
    bmat = 'I'; % standard eigenvalue problem
else
    bmat = 'G'; % generalized eigenvalue problem
end
nev = int32(k); % number of eigenvalues requested
ncv = int32(p); % number of Lanczos vectors
iparam = int32(zeros(11,1));
iparam([1 3 7]) = int32([1 maxit mode]);
select = int32(zeros(p,1));


cputms(1) = cputime - t0; % end timing pre-processing


iter = 0;
ariter = 0;


while (ido ~= 99)
        
    t0 = cputime; % start timing ARPACK calls **aupd
        
    if isrealprob
        arpackc( [prefix 'aupd'], ido, ...
            bmat, int32(n), whch, nev, tol, resid, ncv, ...
            v, ldv, iparam, ipntr, workd, workl, lworkl, info);
    else
        zworkd(1:2:end-1) = real(workd);
        zworkd(2:2:end) = imag(workd);
        arpackc( 'znaupd', ido, ...
            bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
            ldv, iparam, ipntr, zworkd, workl, lworkl, ...
            rwork, info );
        workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]);
    end


    if (info < 0)
        es = sprintf('Error with ARPACK routine %saupd: info = %d',...
           prefix,double(info));
        error(es)
    end
     
    cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd
    t0 = cputime; % start timing MATLAB OP(X)
    
    % Compute which columns of workd ipntr references
    [row,col1] = ind2sub([n,3],double(ipntr(1)));
    if (row ~= 1)
        str = sprintf(['ipntr(1)=%d does not refer to the start of a' ...
                ' column of the %d-by-3 array workd'],double(ipntr(1)),n);
        error(str)
    end
    [row,col2] = ind2sub([n,3],double(ipntr(2)));
    if (row ~= 1)
        str = sprintf(['ipntr(2)=%d does not refer to the start of a' ...
                ' column of the %d-by-3 array workd'],double(ipntr(2)),n);
        error(str)
    end
    if ~isempty(B) & (mode == 3) & (ido == 1)
        [row,col3] = ind2sub([n,3],double(ipntr(3)));
        if (row ~= 1)
            str = sprintf(['ipntr(3)=%d does not refer to the start of a' ...
                    ' column of the %d-by-3 array workd'],double(ipntr(3)),n);
            error(str)
        end
    end
    
    if ((ido == -1) | (ido == 1))
        if Amatrix
            if (mode == 1)
                workd(:,col2) = A * workd(:,col1);
            elseif (mode == 2)
                workd(:,col1) = A * workd(:,col1);
                if issparse(B)
                    workd(permB,col2) = RB \ (RBT \workd(permB,col1));
                else
                    workd(:,col2) = RB \ (RBT \workd(:,col1));
                end
            elseif (mode == 3)
                if isempty(B)
                    if issparse(A)
                        workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
                    else
                        workd(:,col2) = U \ (L \ (P * workd(:,col1)));
                    end
                else % B is not empty
                    if (ido == -1)
                        if cholB
                            workd(:,col2) = RBT * (RB * workd(:,col1));
                        else
                            workd(:,col2) = B * workd(:,col1);
                        end
                        if issparse(A) & issparse(B)
                            workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
                        else
                            workd(:,col2) = U \ (L \ (P * workd(:,col1)));
                        end
                    elseif (ido == 1)
                        if issparse(A) & issparse(B)
                            workd(perm,col2) = U \ (L \ (P * workd(:,col3)));
                        else
                            workd(:,col2) = U \ (L \ (P * workd(:,col3)));
                        end
                    end
                end
            else % mode is not 1,2 or 3
                error(['Unknown mode returned from ' prefix 'aupd'])
            end
        else % A is not a matrix
            if (mode == 1)
                workd(:,col2) = feval(A,workd(:,col1), ...
                    varargin{7-Amatrix-Bnotthere:end});
            elseif (mode == 2)
                workd(:,col1) = feval(A,workd(:,col1), ...
                    varargin{7-Amatrix-Bnotthere:end});
                if issparse(B)
                    workd(permB,col2) = RB \ (RBT \workd(permB,col1));
                else
                    workd(:,col2) = RB \ (RBT \workd(:,col1));
                end
            elseif (mode == 3)
                if isempty(B)
                    workd(:,col2) = feval(A,workd(:,col1), ...
                        varargin{7-Amatrix-Bnotthere:end});
                else
                    if (ido == -1)
                        if cholB
                            workd(:,col2) = RBT * (RB * workd(:,col1));
                        else
                            workd(:,col2) = B * workd(:,col1);
                        end
                        workd(:,col2) = feval(A,workd(:,col2), ...
                            varargin{7-Amatrix-Bnotthere:end});
                    elseif (ido == 1)
                        workd(:,col2) = feval(A,workd(:,col3), ...
                            varargin{7-Amatrix-Bnotthere:end});
                    end
                end
            else % mode is not 1,2 or 3
                error(['Unknown mode returned from ' prefix 'aupd'])
            end
        end % if Amatrix
    elseif (ido == 2)
        if (mode == 2) | (mode == 3)
            if cholB
                workd(:,col2) = RBT * (RB * workd(:,col1));
            else
                workd(:,col2) = B * workd(:,col1);
            end
        else
            error(['Unknown mode returned from ' prefix 'aupd'])
        end
    elseif (ido == 3)
        % setting iparam(1) = ishift = 1 ensures this never happens
        warning(['eigs does not yet support computing the shifts in workl.' ...
                ' Setting reverse communication parameter to 99 and returning'])
        ido = int32(99);
    elseif (ido ~= 99)
        error(['Unknown value of reverse communication parameter' ...
                ' returned from ' prefix 'aupd'])      
    end
    
    cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X)


    if display
        iter = double(ipntr(15));
        if (iter > ariter) & (ido ~= 99)
            ariter = iter;
            ds = sprintf(['Iteration %d: a few Ritz values of the' ...
                    ' %d-by-%d matrix:'],iter,p,p);
            disp(ds)
            if isrealprob
                if issymA
                    dispvec = [workl(double(ipntr(6))+(0:p-1))];
                    if isequal(whch,'BE')
                        % roughly k Large eigenvalues and k Small eigenvalues
                        disp(dispvec(end-2*k+1:end))
                    else
                        % k eigenvalues
                        disp(dispvec(end-k+1:end))
                    end
                else
                    dispvec = [complex(workl(double(ipntr(6))+(0:p-1)), ...
                            workl(double(ipntr(7))+(0:p-1)))];
                    % k+1 eigenvalues (keep complex conjugate pairs together)
                    disp(dispvec(end-k:end))
                end
            else
                dispvec = [complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ...
                        workl(2*double(ipntr(6))+(0:2:2*(p-1))))];
                disp(dispvec(end-k+1:end))
            end
        end
    end
    
end % while (ido ~= 99)


t0 = cputime; % start timing post-processing


flag = 0;
if (info < 0)
    es = sprintf('Error with ARPACK routine %saupd: info = %d',prefix,double(info));
    error(es)
else
    if (nargout >= 2)
        rvec = int32(1); % compute eigenvectors
    else
        rvec = int32(0); % do not compute eigenvectors
    end
    
    if isrealprob
        if issymA
            arpackc( 'dseupd', rvec, 'A', select, ...
                d, v, ldv, sigma, ...
                bmat, int32(n), whch, nev, tol, resid, ncv, ...
                v, ldv, iparam, ipntr, workd, workl, lworkl, info );
            if isequal(whch,'LM') | isequal(whch,'LA')
                d = flipud(d);
                if (rvec == 1)
                    v(:,1:k) = v(:,k:-1:1);
                end
            end
            if ((isequal(whch,'SM') | isequal(whch,'SA')) & (rvec == 0))
                d = flipud(d);
            end
        else
            arpackc( 'dneupd', rvec, 'A', select, ...
                d, di, v, ldv, sigmar, sigmai, workev, ...
                bmat, int32(n), whch, nev, tol, resid, ncv, ...
                v, ldv, iparam, ipntr, workd, workl, lworkl, info );
            d = complex(d,di);
            if rvec
                d(k+1) = [];
            else
                zind = find(d == 0);
                if isempty(zind)
                    d = d(k+1:-1:2);
                else
                    d(max(zind)) = [];
                    d = flipud(d);
                end
            end
        end
    else
        zsigma = [real(sigma); imag(sigma)];
        arpackc( 'zneupd', rvec, 'A', select, ...
            zd, zv, ldv, zsigma, workev, ...
            bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
            ldv, iparam, ipntr, zworkd, workl, lworkl, ...
            rwork, info );
        if issymA
            d = zd(1:2:end-1);
        else
            d = complex(zd(1:2:end-1),zd(2:2:end));
        end
        v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]);
    end
    
    if (info ~= 0)
        es = ['Error with ARPACK routine ' prefix 'eupd: '];
        switch double(info)
            case 2
            ss = sum(select);
            if (ss < k)
            es = [es ...
                    '  The logical variable select was only set with ' int2str(ss) ...
                    ' 1''s instead of nconv=' int2str(double(iparam(5))) ...
                    ' (k=' int2str(k) ').' ...
                    ' Please contact the ARPACK authors at arpack@caam.rice.edu'];
            else
            es = [es ...
                    'The LAPACK reordering routine ' prefix(1) ...
                    'trsen did not return all ' int2str(k) ' eigenvalues.']
            end
            case 1
            es = [es ...
                    'The Schur form could not be reordered by the LAPACK routine ' ...
                    prefix(1) 'trsen.' ...
                    ' Please contact the ARPACK authors at arpack@caam.rice.edu'];
            case -14
            es = [es prefix ...
                    'aupd did not find any eigenvalues to sufficient accuracy'];
            otherwise
            es = [es sprintf('info = %d',double(info))];
        end
        error(es)
    else
        nconv = double(iparam(5));
        if (nconv == 0)
            if (nargout < 3)
                ws = sprintf(['None of the %d requested eigenvalues' ...
                        ' converged'],k);
                warning(ws)
            else
                flag = 1;
            end
        elseif (nconv < k)
            if (nargout < 3)
                ws = sprintf(['Only %d of the %d requested eigenvalues' ...
                        ' converged'],nconv,k);
                warning(ws)
            else
                flag = 1;
            end
        end
    end % if (info ~= 0)
end % if (info < 0)


if (issymA) | (~isrealprob)
    if (nargout <= 1)
        if isrealprob
            varargout{1} = d;
        else
            varargout{1} = d(k:-1:1,1);
        end
    else
        varargout{1} = v(:,1:k);
        varargout{2} = diag(d(1:k,1));
        if (nargout >= 3)
            varargout{3} = flag;
        end
    end
else
    if (nargout <= 1)
        varargout{1} = d;
    else
        cplxd = find(di ~= 0);
        % complex conjugate pairs of eigenvalues occur together
        cplxd = cplxd(1:2:end);
        v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ...
                complex(v(:,cplxd),-v(:,cplxd+1))];
        varargout{1} = v(:,1:k);
        varargout{2} = diag(d);
        if (nargout >= 3)
            varargout{3} = flag;
        end
    end
end


cputms(4) = cputime-t0; % end timing post-processing


cputms(5) = sum(cputms(1:4)); % total time


if (display == 2)
    if (mode == 1)
        innerstr = sprintf(['Compute A*X:' ...
                '                               %f\n'],cputms(3));
    elseif (mode == 2)
        innerstr = sprintf(['Compute A*X and solve B*X=Y for X:' ...
                '         %f\n'],cputms(3));
    elseif (mode == 3)
        if isempty(B)
            innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ...
                    '               %f\n'],cputms(3));
        else
            innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ...
                    '             %f\n'],cputms(3));
        end
    end
    if ((mode == 3) & (Amatrix))
        if isempty(B)
            prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ...
                    '   %f\n'],cputms(1));
        else
            prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ...
                    '   %f\n'],cputms(1));
        end
    elseif ((mode == 2) & (~cholB))
        prepstr = sprintf(['Pre-processing, including chol(B):' ...
                '         %f\n'],cputms(1));
    else
        prepstr = sprintf(['Pre-processing:' ...
                '                            %f\n'],cputms(1));
    end
    sstr = sprintf(['***********CPU Timing Results in seconds***********']);
    ds = sprintf(['\n' sstr '\n' ...
            prepstr ...
            'ARPACK''s %saupd:                           %f\n' ...
            innerstr ...
            'Post-processing with ARPACK''s %seupd:      %f\n' ...
            '***************************************************\n' ...
            'Total:                                     %f\n' ...
            sstr '\n'], ...
        prefix,cputms(2),prefix,cputms(4),cputms(5));
    disp(ds)
end

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

% LLE ALGORITHM (using K nearest neighbors)
%
% [Y] = lle(X,K,dmax)
%
% X = data as D x N matrix (D = dimensionality, N = #points)
% K = number of neighbors
% dmax = max embedding dimensionality
% Y = embedding as dmax x N matrix


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [Y] = lle(X,K,d)


[D,N] = size(X);
fprintf(1,'LLE running on %d points in %d dimensions\n',N,D);




% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS 
fprintf(1,'-->Finding %d nearest neighbours.\n',K);


X2 = sum(X.^2,1);
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;


[sorted,index] = sort(distance);
neighborhood = index(2:(1+K),:);






% STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS
fprintf(1,'-->Solving for reconstruction weights.\n');


if(K>D) 
  fprintf(1,'   [note: K>D; regularization will be used]\n'); 
  tol=1e-3; % regularlizer in case constrained fits are ill conditioned
else
  tol=0;
end


W = zeros(K,N);
for ii=1:N
   z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % shift ith pt to origin
   C = z'*z;                                        % local covariance
   C = C + eye(K,K)*tol*trace(C);                   % regularlization (K>D)
   W(:,ii) = C\ones(K,1);                           % solve Cw=1
   W(:,ii) = W(:,ii)/sum(W(:,ii));                  % enforce sum(w)=1
end;




% STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF COST MATRIX M=(I-W)'(I-W)
fprintf(1,'-->Computing embedding.\n');


% M=eye(N,N); % use a sparse matrix with storage for 4KN nonzero elements
M = sparse(1:N,1:N,ones(1,N),N,N,4*K*N); 
for ii=1:N
   w = W(:,ii);
   jj = neighborhood(:,ii);
   M(ii,jj) = M(ii,jj) - w';
   M(jj,ii) = M(jj,ii) - w;
   M(jj,jj) = M(jj,jj) + w*w';
end;


% CALCULATION OF EMBEDDING
options.disp = 0; options.isreal = 1; options.issym = 1; 
[Y,eigenvals] = eigs(M,d+1,0,options);
Y = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0




fprintf(1,'Done.\n');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




% other possible regularizers for K>D
%   C = C + tol*diag(diag(C));                       % regularlization
%   C = C + eye(K,K)*tol*trace(C)*K;                 % regularlization

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页