THESE/codes/graphe/Ncut_9/eigs_new.m

   1 function  varargout = eigs(varargin)
   2 %EIGS  Find a few eigenvalues and eigenvectors of a matrix using ARPACK
   3 %   D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues.
   4 %   A must be square and should be large and sparse.
   5 %
   6 %   [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude
   7 %   eigenvalues and a matrix V whose columns are the corresponding
   8 %   eigenvectors.
   9 %
  10 %   [V,D,FLAG] = EIGS(A) also returns a convergence flag. If FLAG is 0 then
  11 %   all the eigenvalues converged; otherwise not all converged.
  12 %
  13 %   EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D. B
  14 %   must be symmetric (or Hermitian) positive definite and the same size as
  15 %   A. EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D.
  16 %
  17 %   EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues.
  18 %
  19 %   EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues. If SIGMA is:
  20 %      'LM' or 'SM' - Largest or Smallest Magnitude
  21 %   For real symmetric problems, SIGMA may also be:
  22 %      'LA' or 'SA' - Largest or Smallest Algebraic
  23 %      'BE' - Both Ends, one more from high end if K is odd
  24 %   For nonsymmetric and complex problems, SIGMA may also be:
  25 %      'LR' or 'SR' - Largest or Smallest Real part
  26 %      'LI' or 'SI' - Largest or Smallest Imaginary part
  27 %   If SIGMA is a real or complex scalar including 0, EIGS finds the
  28 %   eigenvalues closest to SIGMA. For scalar SIGMA, and when SIGMA = 'SM',
  29 %   B need only be symmetric (or Hermitian) positive semi-definite since it
  30 %   is not Cholesky factored as in the other cases.
  31 %
  32 %   EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options:
  33 %   OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{false} |
  34 %   true]
  35 %   OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [false | {true}]
  36 %   OPTS.tol: convergence: Ritz estimate residual <= tol*NORM(A) [scalar | {eps}]
  37 %   OPTS.maxit: maximum number of iterations [integer | {300}]
  38 %   OPTS.p: number of Lanczos vectors: K+1<p<=N [integer | {2K}]
  39 %   OPTS.v0: starting vector [N-by-1 vector | {randomly generated}]
  40 %   OPTS.disp: diagnostic information display level [0 | {1} | 2]
  41 %   OPTS.cholB: B is actually its Cholesky factor CHOL(B) [{false} | true]
  42 %   OPTS.permB: sparse B is actually CHOL(B(permB,permB)) [permB | {1:N}]
  43 %   Use CHOL(B) instead of B when SIGMA is a string other than 'SM'.
  44 %
  45 %   EIGS(AFUN,N) accepts the function AFUN instead of the matrix A. AFUN is
  46 %   a function handle and Y = AFUN(X) should return
  47 %      A*X            if SIGMA is unspecified, or a string other than 'SM'
  48 %      A\X            if SIGMA is 0 or 'SM'
  49 %      (A-SIGMA*I)\X  if SIGMA is a nonzero scalar (standard problem)
  50 %      (A-SIGMA*B)\X  if SIGMA is a nonzero scalar (generalized problem)
  51 %   N is the size of A. The matrix A, A-SIGMA*I or A-SIGMA*B represented by
  52 %   AFUN is assumed to be real and nonsymmetric unless specified otherwise
  53 %   by OPTS.isreal and OPTS.issym. In all these EIGS syntaxes, EIGS(A,...)
  54 %   may be replaced by EIGS(AFUN,N,...).
  55 %
  56 %   Example:
  57 %      A = delsq(numgrid('C',15));  d1 = eigs(A,5,'SM');
  58 %
  59 %   Equivalently, if dnRk is the following one-line function:
  60 %      %----------------------------%
  61 %      function y = dnRk(x,R,k)
  62 %      y = (delsq(numgrid(R,k))) \ x;
  63 %      %----------------------------%
  64 %
  65 %      n = size(A,1);  opts.issym = 1;
  66 %      d2 = eigs(@(x)dnRk(x,'C',15),n,5,'SM',opts);
  67 %
  68 %   See also EIG, SVDS, ARPACKC, FUNCTION_HANDLE.
  69
  70 %   Copyright 1984-2008 The MathWorks, Inc.
  71 %   $Revision: 1.45.4.11 $  $Date: 2008/12/01 07:19:19 $
  72
  73 %   EIGS provides the reverse communication interface to ARPACK library
  74 %   routines. EIGS attempts to provide an interface for as many different
  75 %   algorithms as possible. The reverse communication interfaces are
  76 %   documented in the ARPACK Users' Guide, ISBN 0-89871-407-9.
  77
  78 cputms = zeros(5,1);
  79 t0 = cputime; % start timing pre-processing
  80
  81 % Process inputs and do error-checking
  82 if (nargout > 3)
  83     error('MATLAB:eigs:TooManyOutputs', 'Too many output arguments.')
  84 end
  85
  86 % Platform dependent integer type
  87 if strfind(computer, '64')
  88     intconvert = @(arraytoconvert) int64(arraytoconvert);
  89     inttype = 'int64';
  90 else
  91     intconvert = @(arraytoconvert) int32(arraytoconvert);
  92     inttype = 'int32';
  93 end
  94
  95 % Error check inputs and derive some information from them
  96 [A,Amatrix,isrealprob,issymA,n,B,classAB,k,eigs_sigma,whch, ...
  97     sigma,tol,maxit,p,info,eigs_display,cholB,permB,resid,useeig, ...
  98     afunNargs,style,mode,Afactors,Bfactors] = ...
  99     checkInputs(varargin{:});
 100
 101 % Now have enough information to do early return on cases EIGS does not
 102 % handle. For these cases, use the full EIG code.
 103 if useeig
 104     fullEig(nargout);
 105     return
 106 end
 107
 108 if ~isempty(Afactors)
 109     L = Afactors.L;     Afactors.L = [];
 110     U = Afactors.U;     Afactors.U = [];
 111     pp = Afactors.pp;   Afactors.pp = [];
 112     qq = Afactors.qq;   Afactors.qq = [];
 113     dgAsB = Afactors.dgAsB;    Afactors.dgAsB = [];
 114     clear Afactors;
 115 end
 116
 117 if ~isempty(Bfactors)
 118     BisHpd = Bfactors.BisHpd;
 119     if BisHpd
 120         RB = Bfactors.RB;       Bfactors.RB = [];
 121         RBT = Bfactors.RBT;     Bfactors.RBT = [];
 122         permB = Bfactors.permB; Bfactors.permB = [];
 123     else
 124         LB = Bfactors.LB;       Bfactors.LB = [];
 125         UB = Bfactors.UB;       Bfactors.UB = [];
 126         ppB = Bfactors.ppB;     Bfactors.ppB = [];
 127         qqB = Bfactors.qqB;     Bfactors.qqB = [];
 128         dgB = Bfactors.dgB;     Bfactors.dgB = [];
 129     end
 130     clear Bfactors;
 131 end
 132
 133 if isempty(B)
 134     if mode ~= 3
 135         % OP = A
 136         applyOP = @(v)Amtimes(v);
 137     else
 138         % OP = (A-\sigma*I)^{-1}
 139         applyOP = @(v)AminusSigmaBsolve(v);
 140     end
 141 else
 142     if mode ~= 3 && BisHpd == true
 143         % OP = L^{-1}AL^{-T} (B = LL^T)
 144         applyOP = @(v)RBTsolve(Amtimes(RBsolve(v)));
 145     elseif mode ~= 3 && BisHpd == false
 146         % OP = U^{-1}L^{-1}A (B = LU)
 147         applyOP = @(v)Bsolve(Amtimes(v));
 148     else
 149         % OP = (A-\sigma*B)^{-1}B
 150         applyOP = @(v)AminusSigmaBsolve(Bmtimes(v));
 151     end
 152 end
 153 if strcmp(style,'G')
 154     if ~isempty(B) && BisHpd == true && mode == 3
 155         applyM = @(v) Bmtimes(v);
 156     else
 157         applyM = @(v) v;
 158     end
 159 end
 160
 161 % Allocate outputs and ARPACK work variables
 162 if isrealprob
 163     if issymA % real and symmetric
 164         if strcmp(classAB,'single')
 165             aupdfun = 'ssaupd';
 166             eupdfun = 'sseupd';
 167         else
 168             aupdfun = 'dsaupd';
 169             eupdfun = 'dseupd';
 170         end
 171       lworkl = intconvert(p*(p+8));
 172         d = zeros(k,1,classAB);
 173     else % real but not symmetric
 174         if strcmp(classAB,'single')
 175             aupdfun = 'snaupd';
 176             eupdfun = 'sneupd';
 177         else
 178             aupdfun = 'dnaupd';
 179             eupdfun = 'dneupd';
 180         end
 181       lworkl = intconvert(3*p*(p+2));
 182         workev = zeros(3*p,1,classAB);
 183         d = zeros(k+1,1,classAB);
 184         di = zeros(k+1,1,classAB);
 185     end
 186     v = zeros(n,p,classAB);
 187     workd = zeros(n,3,classAB);
 188     workl = zeros(lworkl,1,classAB);
 189 else % complex
 190     if strcmp(classAB,'single')
 191         aupdfun = 'cnaupd';
 192         eupdfun = 'cneupd';
 193     else
 194         aupdfun = 'znaupd';
 195         eupdfun = 'zneupd';
 196     end
 197     zv = zeros(2*n*p,1,classAB);
 198     workd = complex(zeros(n,3,classAB));
 199     zworkd = zeros(2*numel(workd),1,classAB);
 200    lworkl = intconvert(2*(3*p^2+5*p));
 201    workl = zeros(lworkl,1,classAB);
 202     workev = zeros(2*2*p,1,classAB);
 203     zd = zeros(2*(k+1),1,classAB);
 204     rwork = zeros(p,1,classAB);
 205 end
 206 ldv = intconvert(n);
 207 ipntr = zeros(15,1,inttype);
 208 ido = intconvert(0); % reverse communication parameter, initial value
 209 if strcmp(style,'S')
 210     bmat = 'I'; % standard eigenvalue problem
 211 else
 212     bmat = 'G'; % generalized eigenvalue problem
 213 end
 214 nev = intconvert(k); % number of eigenvalues requested
 215 ncv = intconvert(p); % number of Lanczos vectors
 216 iparam = zeros(11,1,inttype);
 217 % iparam(1) = ishift = 1 ensures we are never asked to handle ido=3
 218 iparam([1 3 7]) = [1 maxit mode];
 219 select = zeros(p,1,inttype);
 220
 221 % To Do: Remove this error when ARPACKC supports singles
 222 if strcmp(classAB,'single')
 223     error('MATLAB:eigs:single', ...
 224         'EIGS does not support single precision inputs.')
 225 end
 226
 227 % The ARPACK routines return to EIGS many times per each iteration but we
 228 % only want to display the Ritz values once per iteration (if opts.disp>0).
 229 % Keep track of whether we've displayed this iteration yet in eigs_iter.
 230 eigs_iter = 0;
 231
 232 cputms(1) = cputime - t0; % end timing pre-processing
 233
 234 % Iterate until ARPACK's reverse communication parameter ido says to stop
 235 while (ido ~= 99)
 236
 237     t0 = cputime; % start timing ARPACK calls **aupd
 238
 239     if isrealprob
 240         arpackc( aupdfun, ido, ...
 241          bmat, intconvert(n), whch, nev, tol, resid, ncv, ...
 242             v, ldv, iparam, ipntr, workd, workl, lworkl, info );
 243     else
 244         % The FORTRAN ARPACK routine expects the complex input zworkd to have
 245         % real and imaginary parts interleaved, but the OP about to be
 246         % applied to workd expects it in MATLAB's complex representation with
 247         % separate real and imaginary parts. Thus we need both.
 248         zworkd(1:2:end-1) = real(workd);
 249         zworkd(2:2:end) = imag(workd);
 250         arpackc( aupdfun, ido, ...
 251          bmat, intconvert(n), whch, nev, tol, resid, ncv, ...
 252             zv, ldv, iparam, ipntr, zworkd, workl, lworkl, rwork, info );
 253         workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]);
 254     end
 255
 256     if (info < 0)
 257         error('MATLAB:eigs:ARPACKroutineError', ...
 258             'Error with ARPACK routine %s: info = %d', ...
 259             aupdfun,full(double(info)))
 260     end
 261
 262     cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd
 263     t0 = cputime; % start timing MATLAB OP(X)
 264
 265     % Compute which columns of workd ipntr references
 266     cols = checkIpntr;
 267
 268     % The ARPACK reverse communication parameter ido tells EIGS what to do
 269
 270     switch ido
 271         case -1
 272             workd(:,cols(2)) = applyOP(workd(:,cols(1)));
 273             if strcmp(style,'G') && mode ~= 3
 274                 workd(:,cols(1)) = workd(:,cols(2));
 275             end
 276         case 1
 277             if strcmp(style,'G') && mode == 3 && (isempty(B) || BisHpd)
 278                 % use with B-inner product for mode 3; see applyM
 279                 workd(:,cols(2)) = AminusSigmaBsolve(workd(:,cols(3)));
 280             else
 281                 % same work as case -1
 282                 workd(:,cols(2)) = applyOP(workd(:,cols(1)));
 283                 if strcmp(style,'G') && mode ~= 3
 284                     % use with std inner product; see applyM
 285                     workd(:,cols(1)) = workd(:,cols(2));
 286                 end
 287             end
 288         case 2
 289             workd(:,cols(2)) = applyM(workd(:,cols(1)));
 290         case 99
 291             % ARPACK has converged
 292         otherwise
 293             error('MATLAB:eigs:UnknownIdo','Unknown ido.');
 294     end
 295
 296     cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X)
 297
 298     if eigs_display
 299         displayRitzValues;
 300     end
 301
 302 end % while (ido ~= 99)
 303
 304 t0 = cputime; % start timing post-processing
 305
 306 if (info < 0)
 307     error('MATLAB:eigs:ARPACKroutineError', ...
 308         'Error with ARPACK routine %s: info = %d',aupdfun,full(info));
 309 end % if (info < 0)
 310
 311 if (nargout >= 2)
 312    rvec = intconvert(true); % compute eigenvectors
 313 else
 314    rvec = intconvert(false); % do not compute eigenvectors
 315 end
 316
 317 if isrealprob
 318     if issymA
 319         arpackc( eupdfun, rvec, 'A', select, ...
 320             d, v, ldv, sigma, ...
 321          bmat, intconvert(n), whch, nev, tol, resid, ncv, ...
 322             v, ldv, iparam, ipntr, workd, workl, lworkl, info );
 323         if strcmp(whch,'LM') || strcmp(whch,'LA')
 324             d = flipud(d);
 325             if (rvec == 1)
 326                 v(:,1:k) = v(:,k:-1:1);
 327             end
 328         end
 329         if ((strcmp(whch,'SM') || strcmp(whch,'SA')) && (rvec == 0))
 330             d = flipud(d);
 331         end
 332     else
 333         % If sigma is complex, isrealprob=true and we use [c,z]neupd.
 334         % So use sigmar=sigma and sigmai=0 here in dneupd.
 335         arpackc( eupdfun, rvec, 'A', select, ...
 336             d, di, v, ldv, sigma, 0, workev, ...
 337          bmat, intconvert(n), whch, nev, tol, resid, ncv, ...
 338             v, ldv, iparam, ipntr, workd, workl, lworkl, info );
 339         d = complex(d,di);
 340         if rvec
 341             d(k+1) = [];
 342         else
 343             zind = find(d == 0);
 344             if isempty(zind)
 345                 d = d(k+1:-1:2);
 346             else
 347                 d(max(zind)) = [];
 348                 d = flipud(d);
 349             end
 350         end
 351     end
 352 else
 353     zsigma = [real(sigma); imag(sigma)];
 354     arpackc( eupdfun, rvec, 'A', select, ...
 355         zd, zv, ldv, zsigma, workev, ...
 356       bmat, intconvert(n), whch, nev, tol, resid, ncv, zv, ...
 357         ldv, iparam, ipntr, zworkd, workl, lworkl, ...
 358         rwork, info );
 359     if issymA
 360         d = zd(1:2:end-1);
 361     else
 362         d = complex(zd(1:2:end-1),zd(2:2:end));
 363     end
 364     v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]);
 365 end
 366
 367 flag = processEUPDinfo(nargin<3);
 368
 369 if (issymA) || (~isrealprob)
 370     if (nargout <= 1)
 371         if isrealprob
 372             varargout{1} = d;
 373         else
 374             varargout{1} = d(k:-1:1,1);
 375         end
 376     else
 377         varargout{1} = v(:,1:k);
 378         varargout{2} = diag(d(1:k,1));
 379         if (nargout >= 3)
 380             varargout{3} = flag;
 381         end
 382     end
 383 else
 384     if (nargout <= 1)
 385         varargout{1} = d;
 386     else
 387         cplxd = find(di ~= 0);
 388         % complex conjugate pairs of eigenvalues occur together
 389         cplxd = cplxd(1:2:end);
 390         v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ...
 391             complex(v(:,cplxd),-v(:,cplxd+1))];
 392         varargout{1} = v(:,1:k);
 393         varargout{2} = diag(d);
 394         if (nargout >= 3)
 395             varargout{3} = flag;
 396         end
 397     end
 398 end
 399
 400 if (nargout >= 2) && (mode ~= 3) && ~isempty(B) && BisHpd == true
 401     varargout{1} = RBsolve(varargout{1});
 402 end
 403
 404 if ~isempty(B)
 405     if nargout >= 2
 406         varargout{2} = varargout{2}/scaleB;
 407     else
 408         varargout{1} = varargout{1}/scaleB;
 409     end
 410 end
 411
 412 if BisHpd == true && (nargout >= 2)
 413     vnorms = zeros(k,1);
 414     for ii = 1 : k
 415         vnorms(ii) = scaleB*(varargout{1}(:,ii)'*(Bmtimes(varargout{1}(:,ii))));
 416         varargout{1}(:,ii) = varargout{1}(:,ii)/sqrt(vnorms(ii));
 417     end
 418 end
 419
 420 cputms(4) = cputime-t0; % end timing post-processing
 421
 422 cputms(5) = sum(cputms(1:4)); % total time
 423
 424 if (eigs_display == 2)
 425     printTimings;
 426 end
 427
 428 %-------------------------------------------------------------------------%
 429 % Nested functions
 430 %-------------------------------------------------------------------------%
 431
 432 % checkInputs error checks the inputs to EIGS and also derives some
 433 %   variables from them:
 434 % A may be a matrix or a function applying OP.
 435 % Amatrix is true if A is a matrix, false if A is a function.
 436 % isrealprob is true if all of A, B and sigma are real, false otherwise.
 437 % issymA is true if A is symmetric, false otherwise.
 438 % n is the size of (square) A and B.
 439 % B is [] for the standard problem. Otherwise it may be one of B, CHOL(B)
 440 %   or CHOL(B(permB,permB)).
 441 % classAB is single if either A or B is single, otherwise double.
 442 % k is the number of eigenvalues to be computed.
 443 % eigs_sigma is the value for sigma passed in by the user, 'LM' if it was
 444 %   unspecified. eigs_sigma may be either a string or a scalar value.
 445 % whch is the ARPACK string corresponding to eigs_sigma and mode.
 446 % sigma is the ARPACK scalar corresponding to eigs_sigma and mode.
 447 % tol is the convergence tolerance.
 448 % maxit is the maximum number of iterations.
 449 % p is the number of Lanczos vectors.
 450 % info is the start value, initialized to 1 or 0 to indicate whether to use
 451 % resid as the start vector or not.
 452 % eigs_display is true if Ritz values should be displayed, false otherwise.
 453 % cholB is true if CHOL(B) was passed in instead of B, false otherwise.
 454 % permB may be [], otherwise it is the permutation in CHOL(B(permB,permB)).
 455 % resid is the start vector if specified and info=1, otherwise all zero.
 456 % useeig is true if we need to use EIG instead of ARPACK, otherwise false.
 457 % afunNargs is the range of EIGS' varargin that are to be passed as
 458 %   trailing parameters to the function as in afun(X,P1,P2,...).
 459     function [A,Amatrix,isrealprob,issymA,n,B,classAB,k, ...
 460             eigs_sigma,whch,sigma,tol,maxit,p,info,eigs_display,cholB,...
 461             permB,resid,useeig,afunNargs,style,mode,Afactors,Bfactors] = checkInputs(varargin)
 462         % Process inputs and do error-checking
 463
 464         % Process the input A or the inputs AFUN and N
 465         % Start to derive some qualities (real, symmetric) about the problem
 466         if isfloat(varargin{1})
 467             A = varargin{1};
 468             Amatrix = true;
 469         else
 470             % By checking the function A with fcnchk, we can now use direct
 471             % function evaluation on the result, without resorting to feval
 472             A = fcnchk(varargin{1});
 473             Amatrix = false;
 474         end
 475         isrealprob = true;
 476         issymA = false;
 477         if Amatrix
 478             isrealprob = isreal(A);
 479             issymA = ishermitian(A);
 480             [m,n] = size(A);
 481             if (m ~= n)
 482                 error('MATLAB:eigs:NonSquareMatrixOrFunction',...
 483                     'A must be a square matrix or a function.')
 484             end
 485         else
 486             n = varargin{2};
 487             nstr = 'Size of problem, ''n'', must be a positive integer.';
 488             if ~isscalar(n) || ~isreal(n)
 489                 error('MATLAB:eigs:NonPosIntSize', nstr)
 490             end
 491             if issparse(n)
 492                 n = full(n);
 493             end
 494             if (round(n) ~= n)
 495                 warning('MATLAB:eigs:NonPosIntSize',['%s\n         ' ...
 496                     'Rounding input size.'],nstr)
 497                 n = round(n);
 498             end
 499         end
 500
 501         % Process the input B and derive the class of the problem.
 502         % Is B present in the eigs call or not?
 503         Bpresent = true;
 504         Bstr = 'Generalized matrix B must be the same size as A.';
 505         if (nargin < (3-Amatrix))
 506             B = [];
 507             Bpresent = false;
 508         else
 509             % Is the next input B or K?
 510             B = varargin{3-Amatrix};
 511             if ~isempty(B) % allow eigs(A,[],k,sigma,opts);
 512                 if isscalar(B)
 513                     if n ~= 1
 514                         % this input is really K and B is not specified
 515                         B = [];
 516                         Bpresent = false;
 517                     else
 518                         % This input could be B or K.
 519                         % If A is scalar, then the only valid value for k is 1.
 520                         % So if this input is scalar, let it be B, namely
 521                         % eigs(4,2,...) assumes A=4, B=2, NOT A=4, k=2
 522                         if ~isnumeric(B)
 523                             error('MATLAB:eigs:BsizeMismatchA', Bstr);
 524                         end
 525                         % Unless, of course, the scalar is 1, in which case
 526                         % assume the that it is meant to be K.
 527                         if (B == 1) && ((Amatrix && nargin <= 3) || ...
 528                                 (~Amatrix && nargin <= 4))
 529                             B = [];
 530                             Bpresent = false;
 531                         elseif ~isfloat(B)
 532                             error('MATLAB:eigs:BsizeMismatchA', Bstr);
 533                         end
 534                     end
 535                 else
 536                     % B is a not a scalar.
 537                     if ~isfloat(B) || ~isequal(size(B),[n,n])
 538                         error('MATLAB:eigs:BsizeMismatchA', Bstr);
 539                     end
 540                     isrealprob = isrealprob && isreal(B);
 541                 end
 542             end
 543         end
 544         % ARPACK can only handle homogeneous inputs
 545         if Amatrix
 546             classAB = superiorfloat(A,B);
 547             A = cast(A,classAB);
 548             B = cast(B,classAB);
 549         else
 550             if ~isempty(B)
 551                 classAB = class(B);
 552             else
 553                 classAB = 'double';
 554             end
 555         end
 556
 557         % argOffset tells us where to get the eigs inputs K, SIGMA and OPTS.
 558         % If A is really the function afun, then it also helps us find the
 559         % trailing parameters in eigs(afun,n,[B],k,sigma,opts,P1,P2,...)
 560         % Values of argOffset:
 561         %  0: Amatrix is false and Bpresent is true:
 562         %     eigs(afun,n,B,k,sigma,opts,P1,P2,...)
 563         %  1: Amatrix and Bpresent are both true, or both false
 564         %     eigs(A,B,k,sigma,opts)
 565         %     eigs(afun,n,k,sigma,opts,P1,P2,...)
 566         %  2: Amatrix is true and Bpresent is false:
 567         %     eigs(A,k,sigma,opts)
 568         argOffset = Amatrix + ~Bpresent;
 569
 570         if Amatrix && ((nargin - Bpresent)>4)
 571             error('MATLAB:eigs:TooManyInputs', 'Too many inputs.')
 572         end
 573
 574         % Process the input K.
 575         if (nargin < (4-argOffset))
 576             k = min(n,6);
 577         else
 578             k = varargin{4-argOffset};
 579             kstr = ['Number of eigenvalues requested, k, must be a' ...
 580                 ' positive integer <= n.'];
 581             if ~isnumeric(k) || ~isscalar(k) || ~isreal(k) || (k>n)
 582                 error('MATLAB:eigs:NonIntegerEigQty', kstr)
 583             end
 584             if issparse(k)
 585                 k = full(k);
 586             end
 587             if (round(k) ~= k)
 588                 warning('MATLAB:eigs:NonIntegerEigQty',['%s\n         ' ...
 589                     'Rounding number of eigenvalues.'],kstr)
 590                 k = round(k);
 591             end
 592         end
 593
 594         % Process the input SIGMA and derive ARPACK values whch and sigma.
 595         % eigs_sigma is the value documented in the help as "SIGMA" that is
 596         % passed in to EIGS. eigs_sigma may be either a scalar, including 0,
 597         % or a string, including 'SM'.
 598         % In ARPACK, eigs_sigma corresponds to two variables:
 599         % 1.  which, called "whch" to avoid conflict with MATLAB's function
 600         % 2.  sigma
 601         % whch is always a string. sigma is always a scalar.
 602         % Valid combinations are shown below. Note eigs_sigma = 0/'SM' has
 603         % the same sigma/whch values as eigs_sigma='LM' (default) so these
 604         % must be distinguished by the mode.
 605         % eigs_sigma = 'SM' or 0 => sigma = 0, whch = 'LM' (mode=3)
 606         % eigs_sigma is a string not 'SM' => sigma = 0, whch = eigs_sigma (mode=1)
 607         % eigs_sigma is any scalar => sigma = eigs_sigma, whch = 'LM'
 608         % (mode=1)
 609         whchstr = 'Eigenvalue range sigma must be a valid 2-element string.';
 610         if (nargin < (5-argOffset))
 611             % default: eigs 'LM' => ARPACK which='LM', sigma=0
 612             eigs_sigma = 'LM';
 613             whch = 'LM';
 614             sigma = 0;
 615         else
 616             eigs_sigma = varargin{5-argOffset};
 617             if ischar(eigs_sigma)
 618                 % eigs(string) => ARPACK which=string, sigma=0
 619                 if ~isequal(size(eigs_sigma),[1,2])
 620                     error('MATLAB:eigs:EigenvalueRangeNotValid', ...
 621                         [whchstr '\nFor real symmetric A, the' ...
 622                         ' choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.' ...
 623                         '\nFor non-symmetric or complex' ...
 624                         ' A, the choices are ''%s'', ''%s'', ''%s'', ''%s'',' ...
 625                         ' ''%s'' or ''%s''.\n'], ...
 626                         'LM','SM','LA','SA','BE','LM','SM','LR','SR','LI','SI')
 627                 end
 628                 eigs_sigma = upper(eigs_sigma);
 629                 if strcmp(eigs_sigma,'SM')
 630                     % eigs('SM') => ARPACK which='LM', sigma=0
 631                     whch = 'LM';
 632                 else
 633                     % eigs(string), where string~='SM' => ARPACK which=string, sigma=0
 634                     whch = eigs_sigma;
 635                 end
 636                 sigma = zeros(classAB);
 637             else
 638                 % eigs(scalar) => ARPACK which='LM', sigma=scalar
 639                 if ~isfloat(eigs_sigma) || ~isscalar(eigs_sigma)
 640                     error('MATLAB:eigs:EigenvalueShiftNonScalar',...
 641                         'Eigenvalue shift sigma must be a scalar.')
 642                 end
 643                 sigma = eigs_sigma;
 644                 if issparse(sigma)
 645                     sigma = full(sigma);
 646                 end
 647                 sigma = cast(sigma,classAB);
 648                 isrealprob = isrealprob && isreal(sigma);
 649                 whch = 'LM';
 650             end
 651         end
 652
 653         % Process the input OPTS and derive some ARPACK values.
 654         % ARPACK's minimum tolerance is eps/2 ([S/D]LAMCH's EPS)
 655         tol = eps(classAB);
 656         maxit = [];
 657         p = [];
 658         style = [];
 659         % Always use resid as the start vector, whether it is OPTS.v0 or
 660         % randomly generated within eigs.  We default resid to empty here.
 661         % If the user does not initialize it, we provide a random residual
 662         % below.
 663       info = intconvert(1);
 664         resid = [];
 665         eigs_display = 1;
 666         cholB = false; % do we have B or its Cholesky factor?
 667         permB = []; % if cholB, is it chol(B), or chol(B(permB,permB))?
 668         if (nargin >= (6-argOffset))
 669             opts = varargin{6-argOffset};
 670             if ~isa(opts,'struct')
 671                 error('MATLAB:eigs:OptionsNotStructure',...
 672                     'Options argument must be a structure.')
 673             end
 674             if isfield(opts,'issym') && ~Amatrix
 675                 issymA = opts.issym;
 676                 if (issymA ~= false) && (issymA ~= true)
 677                     error('MATLAB:eigs:InvalidOptsIssym', ...
 678                         'opts.issym must be true or false.')
 679                 end
 680             end
 681             if isfield(opts,'isreal') && ~Amatrix
 682                 if (opts.isreal ~= false) && (opts.isreal ~= true)
 683                     error('MATLAB:eigs:InvalidOptsIsreal', ...
 684                         'opts.isreal must be true or false.')
 685                 end
 686                 isrealprob = isrealprob && opts.isreal;
 687             end
 688             if ~isempty(B) && (isfield(opts,'cholB') || isfield(opts,'permB'))
 689                 if isfield(opts,'cholB')
 690                     cholB = opts.cholB;
 691                     if (cholB ~= false) && (cholB ~= true)
 692                         error('MATLAB:eigs:InvalidOptsCholB', ...
 693                             'opts.cholB must be true or false.')
 694                     end
 695                     if isfield(opts,'permB')
 696                         if issparse(B) && cholB
 697                             permB = opts.permB;
 698                             if ~isvector(permB) || ~isequal(sort(permB(:)),(1:n)')
 699                                 error('MATLAB:eigs:InvalidOptsPermB',...
 700                                     'opts.permB must be a permutation of 1:n.')
 701                             end
 702                         else
 703                             warning('MATLAB:eigs:IgnoredOptionPermB', ...
 704                                 ['Ignoring opts.permB since B is not its sparse' ...
 705                                 ' Cholesky factor.'])
 706                         end
 707                     end
 708                 end
 709             end
 710             if isfield(opts,'tol')
 711                 if ~isfloat(tol) || ~isscalar(opts.tol) || ~isreal(opts.tol) || (opts.tol<=0)
 712                     error('MATLAB:eigs:InvalidOptsTol',...
 713                         ['Convergence tolerance opts.tol must be a strictly' ...
 714                         ' positive real scalar.'])
 715                 end
 716                 tol = cast(full(opts.tol),classAB);
 717             end
 718             if isfield(opts,'p')
 719                 p = opts.p;
 720                 pstr = ['Number of basis vectors opts.p must be a positive' ...
 721                     ' integer <= n.'];
 722                 if ~isnumeric(p) || ~isscalar(p) || ~isreal(p) || (p<=0) || (p>n)
 723                     error('MATLAB:eigs:InvalidOptsP', pstr)
 724                 end
 725                 if issparse(p)
 726                     p = full(p);
 727                 end
 728                 if (round(p) ~= p)
 729                     warning('MATLAB:eigs:NonIntegerVecQty',['%s\n         ' ...
 730                         'Rounding number of basis vectors.'],pstr)
 731                     p = round(p);
 732                 end
 733             end
 734             if isfield(opts,'maxit')
 735                 maxit = opts.maxit;
 736                 str = ['Maximum number of iterations opts.maxit must be' ...
 737                     ' a positive integer.'];
 738                 if ~isnumeric(maxit) || ~isscalar(maxit) || ~isreal(maxit) || (maxit<=0)
 739                     error('MATLAB:eigs:OptsMaxitNotPosInt', str)
 740                 end
 741                 if issparse(maxit)
 742                     maxit = full(maxit);
 743                 end
 744                 if (round(maxit) ~= maxit)
 745                     warning('MATLAB:eigs:NonIntegerIterationQty',['%s\n         ' ...
 746                         'Rounding number of iterations.'],str)
 747                     maxit = round(maxit);
 748                 end
 749             end
 750             if isfield(opts,'v0')
 751                 if ~isfloat(opts.v0) || ~isequal(size(opts.v0),[n,1])
 752                     error('MATLAB:eigs:WrongSizeOptsV0',...
 753                         'Start vector opts.v0 must be n-by-1.')
 754                 end
 755                 if isrealprob
 756                     if ~isreal(opts.v0)
 757                         error('MATLAB:eigs:NotRealOptsV0',...
 758                             'Start vector opts.v0 must be real for real problems.')
 759                     end
 760                     resid(1:n,1) = full(opts.v0);
 761                 else
 762                     resid(2:2:2*n,1) = full(imag(opts.v0));
 763                     resid(1:2:(2*n-1),1) = full(real(opts.v0));
 764                 end
 765             end
 766             if isfield(opts,'disp')
 767                 eigs_display = opts.disp;
 768                 dispstr = 'Diagnostic level opts.disp must be an integer.';
 769                 if ~isnumeric(eigs_display) || ~isscalar(eigs_display) || ...
 770                         ~isreal(eigs_display) || (eigs_display<0)
 771                     error('MATLAB:eigs:NonIntegerDiagnosticLevel', dispstr)
 772                 end
 773                 if (round(eigs_display) ~= eigs_display)
 774                     warning('MATLAB:eigs:NonIntegerDiagnosticLevel', ...
 775                         '%s\n         Rounding diagnostic level.',dispstr)
 776                     eigs_display = round(eigs_display);
 777                 end
 778             end
 779             if isfield(opts,'cheb')
 780                 error('MATLAB:eigs:ObsoleteOptionCheb', ...
 781                     'Polynomial acceleration opts.cheb is an obsolete option.');
 782             end
 783             if isfield(opts,'stagtol')
 784                 error('MATLAB:eigs:ObsoleteOptionStagtol', ...
 785                     'Stagnation tolerance opts.stagtol is an obsolete option.');
 786             end
 787             if isfield(opts,'style')
 788                 style = opts.style;
 789                 str = 'style must be ''S'' or ''G''.';
 790                 if ~ischar(style)
 791                     error('MATLAB:eigs:InvalidStyle',str);
 792                 elseif ~(strcmp(style,'S') || strcmp(style,'G'))
 793                     error('MATLAB:eigs:InvalidStyle',str);
 794                 end
 795             end
 796         end
 797         if (isempty(resid))
 798             if isrealprob
 799                 resid = cast(rand(n,1),classAB);
 800             else
 801                 resid = cast(rand(2*n,1),classAB);
 802             end
 803         end
 804
 805         afunNargs = zeros(1,0);
 806         if ~Amatrix
 807             % The trailing parameters for afun start at varargin{7-argOffset}
 808             % in eigs(afun,n,[B],k,sigma,opts,P1,P2,...). If there are no
 809             % trailing parameters in eigs, then afunNargs is a 1-by-0 empty
 810             % and no trailing parameters are passed to afun(x)
 811             afunNargs = 7-argOffset:nargin;
 812         end
 813
 814         % Now that OPTS has been processed, do final error checking and
 815         % assign ARPACK variables
 816
 817         % Extra check on input B
 818         if ~isempty(B)
 819             % B must be symmetric (Hermitian) positive (semi-)definite
 820             if cholB
 821                 if ~isequal(triu(B),B)
 822                     error('MATLAB:eigs:BNotChol', ...
 823                         'opts.CHOLB specified, but B is not upper triangular.');
 824                 end
 825             end
 826         end
 827
 828         if isempty(style)
 829             if strcmp(eigs_sigma,'SM') || isscalar(eigs_sigma) || ~isempty(B)
 830                 style = 'G';
 831             else
 832                 style = 'S';
 833             end
 834         end
 835
 836         if ~isempty(B)
 837             scaleB = norm(B,'fro')/sqrt(n);
 838             scaleB = 2^(floor(log2(scaleB+1)));
 839             B = B/scaleB;
 840             if cholB
 841                 scaleB = scaleB^2;
 842             end
 843             if isscalar(eigs_sigma)
 844                 sigma = scaleB*eigs_sigma;
 845             end
 846         end
 847
 848
 849         if strcmp(eigs_sigma,'SM') || ~ischar(eigs_sigma)
 850             % eigs(A,B,k,scalarSigma) or eigs(A,B,k,'SM'), B may be []
 851             % Note: sigma must be real for [s,d]saupd and [s,d]naupd
 852             % If sigma is complex, even if A and B are both real, we use
 853             % [c,z]naupd.
 854             % This means that mode=3 in [s,d]naupd, which has
 855             % OP = real(inv(A - sigma*M)*M) and B = M
 856             % reduces to the same OP as [s,d]saupd and [c,z]naupd.
 857             % A*x = lambda*M*x, M symmetric (positive) semi-definite
 858             % => OP = inv(A - sigma*M)*M and B = M
 859             % => shift-and-invert mode
 860             mode = 3;
 861         elseif strcmp(style,'S')
 862             % eigs(A,k,stringSigma) or eigs(A,[],k,stringSigma),
 863             % stringSigma~='SM'
 864             % A*x = lambda*x
 865             % => OP = A and B = I
 866             mode = 1;
 867         else
 868             % eigs(A,B,k,stringSigma), stringSigma~='SM'
 869             % A*x = lambda*B*x
 870             % => OP = inv(B)*A and use standard inner product.
 871             mode = 2;
 872         end
 873
 874         BisHpd = false;
 875         if cholB || (~isempty(B) && ishermitian(B))
 876             % The reordering permutation permB is [] unless B is sparse
 877             [RB,RBT,permB,BisHpd] = CHOLfactorB;
 878             if mode == 3 && ~cholB
 879                 RB = [];    RBT = [];   permB = [];
 880             end
 881         end
 882
 883         qqB = [];
 884         if BisHpd == false && (mode == 1 || mode == 2)
 885             [LB,UB,ppB,qqB,dgB] = LUfactorB;
 886         end
 887
 888         Bfactors = [];
 889         if ~isempty(B)
 890             if BisHpd == true
 891                 Bfactors = struct('RB',RB,'RBT',RBT,'permB',permB,'BisHpd',BisHpd);
 892             elseif (mode == 1 || mode == 2)
 893                 Bfactors = struct('LB',LB,'UB',UB,'ppB',ppB,'qqB',qqB,'dgB',dgB,'BisHpd',BisHpd);
 894             end
 895         end
 896
 897         Afactors = [];
 898         qq = [];
 899         if (mode == 3) && Amatrix % need lu(A-sigma*B)
 900             % The reordering permutation permAsB is [] unless A-sigma*B is sparse
 901             [L,U,pp,qq,dgAsB] = LUfactorAminusSigmaB;
 902             Afactors = struct('L',L,'U',U,'pp',pp,'qq',qq,'dgAsB',dgAsB);
 903         end % if (mode == 3) && Amatrix
 904
 905         % under these conditions, OP must be unsymmetric
 906         % note that OP = inv(A-\sigma*B)*B IS symmetric if A is symmetric
 907         % and B-inner product is used!
 908         if ~isempty(B) && (BisHpd == false || (strcmp(style,'S') && mode == 3))
 909             issymA = false;
 910         end
 911         % Extra check on input K
 912         % We fall back on using the full EIG code if K is too large.
 913         useeig = false;
 914         if isrealprob && issymA
 915             knstr = sprintf(['For real symmetric problems, must have' ...
 916                 ' number of eigenvalues k < n.\n']);
 917         else
 918             knstr = sprintf(['For nonsymmetric and complex problems,' ...
 919                 ' must have number of eigenvalues k < n-1.\n']);
 920         end
 921         if isempty(B)
 922             knstr = [knstr 'Using eig(full(A)) instead.'];
 923         else
 924             knstr = [knstr 'Using eig(full(A),full(B)) instead.'];
 925         end
 926         if (k == 0)
 927             useeig = true;
 928         end
 929         if isrealprob && issymA
 930             if (k > n-1)
 931                 if (n >= 6)
 932                     warning('MATLAB:eigs:TooManyRequestedEigsForRealSym', ...
 933                         '%s',knstr)
 934                 end
 935                 useeig = true;
 936             end
 937         else
 938             if (k > n-2)
 939                 if (n >= 7)
 940                     warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ...
 941                         '%s',knstr)
 942                 end
 943                 useeig = true;
 944             end
 945         end
 946
 947         % Extra check on input SIGMA
 948         if isrealprob && issymA
 949             if ~isreal(sigma)
 950                 error('MATLAB:eigs:ComplexShiftForRealSymProblem',...
 951                     ['For real symmetric problems, eigenvalue shift sigma must' ...
 952                     ' be real.'])
 953             end
 954         else
 955             if ~isrealprob && issymA && ~isreal(sigma)
 956                 warning('MATLAB:eigs:ComplexShiftForHermitianProblem', ...
 957                     ['Complex eigenvalue shift sigma on a Hermitian problem' ...
 958                     ' (all real eigenvalues).'])
 959             end
 960         end
 961         if isrealprob && issymA
 962             if strcmp(whch,'LR')
 963                 whch = 'LA';
 964                 warning('MATLAB:eigs:SigmaChangedToLA', ...
 965                     ['For real symmetric problems, sigma value ''LR''' ...
 966                     ' (Largest Real) is now ''LA'' (Largest Algebraic).'])
 967             end
 968             if strcmp(whch,'SR')
 969                 whch = 'SA';
 970                 warning('MATLAB:eigs:SigmaChangedToSA', ...
 971                     ['For real symmetric problems, sigma value ''SR''' ...
 972                     ' (Smallest Real) is now ''SA'' (Smallest Algebraic).'])
 973             end
 974             if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'})
 975                 error('MATLAB:eigs:EigenvalueRangeNotValid', ...
 976                     [whchstr '\nFor real symmetric A, the' ...
 977                     ' choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
 978                     'LM','SM','LA','SA','BE');
 979             end
 980         else
 981             if strcmp(whch,'BE')
 982                 warning('MATLAB:eigs:SigmaChangedToLM', ...
 983                     ['Sigma value ''BE'' is now only available for real' ...
 984                     ' symmetric problems.  Computing ''LM'' eigenvalues instead.'])
 985                 whch = 'LM';
 986             end
 987             if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'})
 988                 error('MATLAB:eigs:EigenvalueRangeNotValid', ...
 989                     [whchstr '\nFor non-symmetric or complex' ...
 990                     ' A, the choices are ''%s'', ''%s'', ''%s'', ''%s'',' ...
 991                     ' ''%s'' or ''%s''.\n'],'LM','SM','LR','SR','LI','SI');
 992             end
 993         end
 994
 995         % The remainder of the error checking does not apply for the large
 996         % values of K that force us to use full EIG instead of ARPACK.
 997         if useeig
 998             return
 999         end
1000
1001         % Extra check on input OPTS.p
1002         if isempty(p)
1003             if isrealprob && ~issymA
1004                 p = min(max(2*k+1,20),n);
1005             else
1006                 p = min(max(2*k,20),n);
1007             end
1008         else
1009             if isrealprob && issymA
1010                 if (p <= k)
1011                     error('MATLAB:eigs:InvalidOptsPforRealSymProb',...
1012                         ['For real symmetric problems, must have number of' ...
1013                         ' basis vectors opts.p > k.'])
1014                 end
1015             else
1016                 if (p <= k+1)
1017                     error('MATLAB:eigs:InvalidOptsPforComplexOrNonSymProb',...
1018                         ['For nonsymmetric and complex problems, must have number of' ...
1019                         ' basis vectors opts.p > k+1.'])
1020                 end
1021             end
1022         end
1023
1024         % Extra check on input OPTS.maxit
1025         if isempty(maxit)
1026             maxit = max(300,ceil(2*n/max(p,1)));
1027         end
1028
1029
1030         % Nested functions in checkInputs
1031         %-------------------------------------------------------------------------%
1032         function [RB,RBT,ppB,BisHpd] = CHOLfactorB
1033             % permB may be [] (from checkInputs) if the problem is not sparse
1034             % or if it was not passed in as opts.permB
1035             ppB = permB;
1036             if cholB
1037                 % CHOL(B) was passed in as B
1038                 RB = B;
1039                 RBT = B';
1040                 BisHpd = true;
1041             else
1042                 % CHOL(B) was not passed into EIGS
1043                 if ~isempty(B)
1044                     % Algorithm requires CHOL(B) to be computed
1045                     if issparse(B)
1046                         ppB = symamd(B);
1047                         [RB,idxB] = chol(B(ppB,ppB));
1048                     else
1049                         [RB,idxB] = chol(B);
1050                     end
1051                     if mode == 3
1052                         RB = [];    ppB = [];
1053                     end
1054                     RBT = RB';
1055                     if (idxB == 0)
1056                         BisHpd = true;
1057                     elseif mode == 3 && isreal(B)
1058                         % LDL decomposition is only for 'SM' eigenvalues of the
1059                         % pair (A,B) where B is Hermitian positive
1060                         % semi-definite; in this case, as ARPACK users' guide
1061                         % suggests, one should still use B-(semi)inner product
1062                         [LB,DB,pvB] = ldl(B,'vector'); %#ok
1063                         BisHpd = checkTridiagForHSD(diag(DB), diag(DB,1));
1064                     else
1065                         BisHpd = false;
1066                     end
1067                 end
1068             end
1069             if ~isempty(B) && issparse(B)
1070                tmp = speye(length(B));
1071                ppB = tmp(ppB,1:length(B));
1072             end
1073         end % CHOLfactorB
1074         %-------------------------------------------------------------------------%
1075         function [LB,UB,ppB,qqB,dgB] = LUfactorB
1076             % LU factor B, including a reordering perm if it is sparse
1077             if issparse(B)
1078                 [LB,UB,ppB,qqB,dgB] = lu(B);
1079             else
1080                 [LB,UB,ppB] = lu(B,'vector');
1081                 qqB = []; dgB = [];
1082             end
1083             % Warn if lu(B) is ill-conditioned
1084             dUB = diag(UB);
1085             if any(dUB == 0) || any(diag(LB) == 0)
1086                 error('MATLAB:eigs:SingularB', ...
1087                     ['B is singular.\nUnable to compute the specified eigenvalues'...
1088                      ' because infinite eigenvalue(s) exist']);
1089             end
1090             rcondestUB = full(min(abs(dUB)) / max(abs(dUB)));
1091             if (rcondestUB < eps)
1092                 warning('MATLAB:eigs:IllConditionedB',...
1093                     ['B has small reciprocal condition' ...
1094                     ' estimate: %f\n' ...
1095                     '         indicating that B^{-1}A might not be applied accurately.\n' ...
1096                     '         and there may be infinite eigenvalue(s).\n'],...
1097                     rcondestUB);
1098             end
1099         end % LUfactorB
1100
1101         %-------------------------------------------------------------------------%
1102         function [L,U,pp,qq,dgAsB] = LUfactorAminusSigmaB
1103             % LU factor A-sigma*B, including a reordering perm if it is sparse
1104             if sigma == 0
1105                 AsB = A;
1106             elseif isempty(B)
1107                 if issparse(A)
1108                     AsB = A - sigma * speye(n);
1109                 else
1110                     AsB = A - sigma * eye(n);
1111                 end
1112             else
1113                 if cholB
1114                     if issparse(B)
1115                         AsB = A - sigma * checkInputBmtimes(speye(n));
1116                     else
1117                         AsB = A - sigma * checkInputBmtimes(eye(n));
1118                     end
1119                 else
1120                     AsB = A - sigma * B;
1121                 end
1122             end
1123             if issparse(AsB)
1124                 [L,U,pp,qq,dgAsB] = lu(AsB);
1125             else
1126                 [L,U,pp] = lu(AsB,'vector');
1127                 qq = []; dgAsB = [];
1128             end
1129             % Warn if lu(A-sigma*B) is ill-conditioned
1130             % => sigma is close to an exact eigenvalue of (A,B)
1131             if isempty(B)
1132                 ds = '(A-sigma*I)';
1133             else
1134                 ds = '(A-sigma*B)';
1135             end
1136             dU = diag(U);
1137             if any(dU == 0) || any(diag(L) == 0)
1138                 error('MATLAB:eigs:AminusBSingular', ...
1139                     [ds 'is singular. The shift is an eigenvalue.\n' ...
1140                     'Try to use some other shift please.\n']);
1141             end
1142             rcondestU = full(min(abs(dU)) / max(abs(dU)));
1143             if (rcondestU < eps)
1144                 warning('MATLAB:eigs:SigmaNearExactEig',...
1145                     [ds ' has small reciprocal condition' ...
1146                     ' estimate: %f\n' ...
1147                     '         indicating that sigma is near an exact' ...
1148                     ' eigenvalue.\n         The algorithm may not converge unless' ...
1149                     ' you try a new value for sigma.\n'], ...
1150                     rcondestU);
1151             end
1152         end % LUfactorAminusSigmaB
1153
1154 %-------------------------------------------------------------------------%
1155         function v = checkInputBmtimes(u)
1156             % Matrix-vector multiply v = B*u
1157             if cholB % use B's cholesky factor and its transpose
1158                 if ~isempty(permB)
1159                     v = permB'*(RBT * (RB * (permB*u)));
1160                 else
1161                     v = RBT * (RB * u);
1162                 end
1163             else
1164                 v = B * u;
1165             end
1166         end
1167
1168     end % checkInputs
1169
1170 %-------------------------------------------------------------------------%
1171     function fullEig(nOutputs)
1172         % Use EIG(FULL(A)) or EIG(FULL(A),FULL(B)) instead of ARPACK
1173         if ~isempty(B)
1174             B = Bmtimes(eye(n));
1175             B = B*scaleB;
1176         end
1177         if isfloat(A)
1178             if issparse(A);
1179                 A = full(A);
1180             end
1181         else
1182             % A is specified by a function.
1183             % Form the matrix A by applying the function.
1184             if ischar(eigs_sigma) && ~strcmp(eigs_sigma,'SM')
1185                 % A is a function multiplying A*x
1186                 AA = eye(n);
1187                 for i = 1:n
1188                     AA(:,i) = A(AA(:,i),varargin{afunNargs});
1189                 end
1190                 A = AA;
1191             else
1192                 if (isfloat(eigs_sigma) && eigs_sigma == 0) || strcmp(eigs_sigma,'SM')
1193                     % A is a function solving A\x
1194                     invA = eye(n);
1195                     for i = 1:n
1196                         invA(:,i) = A(invA(:,i),varargin{afunNargs});
1197                     end
1198                     A = eye(n) / invA;
1199                 else
1200                     % A is a function solving (A-sigma*B)\x
1201                     % B may be [], indicating the identity matrix
1202                     % U = (A-sigma*B)\sigma*B
1203                     % => (A-sigma*B)*U = sigma*B
1204                     % => A*U = sigma*B(U + eye(n))
1205                     % => A = sigma*B(U + eye(n)) / U
1206                     if isempty(B)
1207                         sB = eigs_sigma*eye(n);
1208                     else
1209                         sB = eigs_sigma*B;
1210                     end
1211                     U = zeros(n,n);
1212                     for i = 1:n
1213                         U(:,i) = A(sB(:,i),varargin{afunNargs});
1214                     end
1215                     A = sB*(U+eye(n)) / U;
1216                 end
1217             end
1218         end
1219
1220         if isempty(B)
1221             eigInputs = {A};
1222         else
1223             eigInputs = {A,B};
1224         end
1225         % Now with full floating point matrices A and B, use EIG:
1226         if (nOutputs <= 1)
1227             d = eig(eigInputs{:});
1228         else
1229             [V,D] = eig(eigInputs{:});
1230             d = diag(D);
1231         end
1232
1233         % Grab the eigenvalues we want, based on sigma
1234         firstKindices = 1:k;
1235         lastKindices = n:-1:n-k+1;
1236         if ischar(eigs_sigma)
1237             switch eigs_sigma
1238                 case 'LM'
1239                     [ignore,ind] = sort(abs(d));
1240                     range = lastKindices;
1241                 case 'SM'
1242                     [ignore,ind] = sort(abs(d));
1243                     range = firstKindices;
1244                 case 'LA'
1245                     [ignore,ind] = sort(d);
1246                     range = lastKindices;
1247                 case 'SA'
1248                     [ignore,ind] = sort(d);
1249                     range = firstKindices;
1250                 case 'LR'
1251                     [ignore,ind] = sort(abs(real(d)));
1252                     range = lastKindices;
1253                 case 'SR'
1254                     [ignore,ind] = sort(abs(real(d)));
1255                     range = firstKindices;
1256                 case 'LI'
1257                     [ignore,ind] = sort(abs(imag(d)));
1258                     range = lastKindices;
1259                 case 'SI'
1260                     [ignore,ind] = sort(abs(imag(d)));
1261                     range = firstKindices;
1262                 case 'BE'
1263                     [ignore,ind] = sort(abs(d));
1264                     range = [1:floor(k/2), n-ceil(k/2)+1:n];
1265                 otherwise
1266                     error('MATLAB:eigs:fullEigSigma','Unknown value of sigma');
1267             end
1268         else
1269             % sigma is a scalar
1270             [ignore,ind] = sort(abs(d-eigs_sigma));
1271             range = 1:k;
1272         end
1273
1274         if (nOutputs <= 1)
1275             varargout{1} = d(ind(range));
1276         else
1277             varargout{1} = V(:,ind(range));
1278             varargout{2} = D(ind(range),ind(range));
1279             if (nOutputs == 3)
1280                 % flag indicates "convergence"
1281                 varargout{3} = 0;
1282             end
1283         end
1284
1285     end % FULLEIG
1286
1287
1288 %-------------------------------------------------------------------------%
1289     function cols = checkIpntr
1290         % Check that ipntr returned from ARPACK refers to the start of a
1291         % column of workd.
1292         if (ido == 1) && (mode == 3) && strcmp(style,'G')
1293             inds = double(ipntr(1:3));
1294         else
1295             inds = double(ipntr(1:2));
1296         end
1297         rows = mod(inds-1,n)+1;
1298         cols = (inds-rows)/n+1;
1299         if ~all(rows==1)
1300             error('MATLAB:eigs:ipntrMismatchWorkdColumn', ...
1301                 ['One of ipntr(1:3) does not refer to the start' ...
1302                 ' of a column of the %d-by-3 array workd.'],n)
1303         end
1304     end % checkIpntr
1305
1306 %-------------------------------------------------------------------------%
1307     function v = Amtimes(u)
1308         % Matrix-vector multiply v = A*u
1309         if Amatrix
1310             v = A * u;
1311         else % A is a function
1312             v = A(u,varargin{afunNargs});
1313             if isrealprob && ~isreal(v)
1314                 error('MATLAB:eigs:complexFunction', ...
1315                     'AFUN is complex; set opts.isreal = false.');
1316             end
1317         end
1318     end
1319
1320 %-------------------------------------------------------------------------%
1321     function v = Bmtimes(u)
1322         % Matrix-vector multiply v = B*u
1323         if cholB % use B's cholesky factor and its transpose
1324             if ~isempty(permB)
1325                 v = permB'*(RBT * (RB * (permB*u)));
1326             else
1327                 v = RBT * (RB * u);
1328             end
1329         else
1330             v = B * u;
1331         end
1332     end
1333
1334 %-------------------------------------------------------------------------%
1335     function v = RBsolve(u)
1336         % Solve v = RB\u for v
1337         if issparse(B)
1338             if ~isempty(permB)
1339                 v = permB'*(RB \ u);
1340             else
1341                 v = RB \ u;
1342             end
1343         else
1344             RBopts.UT = true;
1345             v = linsolve(RB,u,RBopts);
1346         end
1347     end
1348
1349 %-------------------------------------------------------------------------%
1350     function v = RBTsolve(u)
1351         % Solve v = RB'\u for v
1352         if issparse(B)
1353             if ~isempty(permB)
1354                 v = RBT \ (permB*u);
1355             else
1356                 v = RBT \ u;
1357             end
1358         else
1359             RBTopts.LT = true;
1360             v = linsolve(RBT,u,RBTopts);
1361         end
1362     end
1363
1364 %-------------------------------------------------------------------------%
1365     function v = AminusSigmaBsolve(u)
1366         % Solve v = (A-sigma*B)\u for v
1367         if Amatrix
1368             if ~isempty(dgAsB)
1369                 % use LU reordering permAsB
1370                 v = qq*(U \ (L \ (pp*(dgAsB\u))));
1371             else
1372                 v = U \ (L \ u(pp));
1373             end
1374         else % A is a function
1375             v = A(u,varargin{afunNargs});
1376             if isrealprob && ~isreal(v)
1377                 error('MATLAB:eigs:complexFunction', ...
1378                     'AFUN is complex; set opts.isreal = false.');
1379             end
1380         end
1381     end % AminusSigmaBsolve
1382 %-------------------------------------------------------------------------%
1383     function v = Bsolve(u)
1384         % Solve v = (A-sigma*B)\u for v
1385         if ~isempty(dgB)
1386             % use LU reordering permAsB
1387             v = qqB*(UB \ (LB \ (ppB*(dgB\u))));
1388         else
1389             v = UB \ (LB \ u(ppB));
1390         end
1391     end % AminusSigmaBsolve
1392
1393 %-------------------------------------------------------------------------%
1394     function displayRitzValues
1395         % Display a few Ritz values at the current iteration
1396         iter = double(ipntr(15));
1397         if (iter > eigs_iter) && (ido ~= 99)
1398             eigs_iter = iter;
1399             ds = sprintf(['Iteration %d: a few Ritz values of the' ...
1400                 ' %d-by-%d matrix:'],iter,p,p);
1401             disp(ds)
1402             if isrealprob
1403                 if issymA
1404                     dispvec = workl(double(ipntr(6))+(0:p-1));
1405                     if strcmp(whch,'BE')
1406                         % roughly k Large eigenvalues and k Small eigenvalues
1407                         disp(dispvec(max(end-2*k+1,1):end))
1408                     else
1409                         % k eigenvalues
1410                         disp(dispvec(max(end-k+1,1):end))
1411                     end
1412                 else
1413                     dispvec = complex(workl(double(ipntr(6))+(0:p-1)), ...
1414                         workl(double(ipntr(7))+(0:p-1)));
1415                     % k+1 eigenvalues (keep complex conjugate pairs together)
1416                     disp(dispvec(max(end-k,1):end))
1417                 end
1418             else
1419                 dispvec = complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ...
1420                     workl(2*double(ipntr(6))+(0:2:2*(p-1))));
1421                 disp(dispvec(max(end-k+1,1):end))
1422             end
1423         end
1424     end
1425
1426 %-------------------------------------------------------------------------%
1427     function flag = processEUPDinfo(warnNonConvergence)
1428         % Process the info flag returned by the ARPACK routine **eupd
1429         flag = 0;
1430         if (info ~= 0)
1431             es = ['Error with ARPACK routine ' eupdfun ':\n'];
1432             switch double(info)
1433                 case 2
1434                     ss = sum(select);
1435                     if (ss < k)
1436                         error('MATLAB:eigs:ARPACKroutineError02ssLTk', ...
1437                             [es 'The logical variable select was only set' ...
1438                             ' with %d 1''s instead of nconv=%d (k=%d).\n' ...
1439                             'Please report this to the ARPACK authors at' ...
1440                             ' arpack@caam.rice.edu.'], ...
1441                             ss,double(iparam(5)),k)
1442                     else
1443                         error('MATLAB:eigs:ARPACKroutineError02', ...
1444                             [es 'The LAPACK reordering routine %strsen' ...
1445                             ' did not return all %d eigenvalues.'], ...
1446                             aupdfun(1),k);
1447                     end
1448                 case 1
1449                     error('MATLAB:eigs:ARPACKroutineError01', ...
1450                         [es 'The Schur form could not be reordered by the' ...
1451                         ' LAPACK routine %strsen.\nPlease report this to the' ...
1452                         ' ARPACK authors at arpack@caam.rice.edu.'], ...
1453                         aupdfun(1))
1454                 case -14
1455                     error('MATLAB:eigs:ARPACKroutineErrorMinus14', ...
1456                         [es aupdfun ...
1457                         ' did not find any eigenvalues to sufficient accuracy.']);
1458                 otherwise
1459                     error('MATLAB:eigs:ARPACKroutineError', ...
1460                         [es 'info = %d. Please consult the ARPACK Users''' ...
1461                         ' Guide for more information.'],full(info));
1462             end
1463         else
1464             nconv = double(iparam(5));
1465             if (nconv == 0)
1466                 if (warnNonConvergence)
1467                     warning('MATLAB:eigs:NoEigsConverged', ...
1468                         'None of the %d requested eigenvalues converged.',k)
1469                 else
1470                     flag = 1;
1471                 end
1472             elseif (nconv < k)
1473                 if (warnNonConvergence)
1474                     warning('MATLAB:eigs:NotAllEigsConverged', ...
1475                         'Only %d of the %d requested eigenvalues converged.', ...
1476                         nconv,k)
1477                 else
1478                     flag = 1;
1479                 end
1480             end
1481         end
1482     end % processEUPDinfo
1483
1484 %-------------------------------------------------------------------------%
1485     function printTimings
1486         % Print the time taken for each major stage of the EIGS algorithm
1487         if (mode == 1)
1488             innerstr = sprintf(['Compute A*X:' ...
1489                 '                               %f\n'],cputms(3));
1490         elseif (mode == 3)
1491             if isempty(B)
1492                 innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ...
1493                     '               %f\n'],cputms(3));
1494             else
1495                 innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ...
1496                     '             %f\n'],cputms(3));
1497             end
1498         end
1499         if ((mode == 3) && (Amatrix))
1500             if isempty(B)
1501                 prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ...
1502                     '   %f\n'],cputms(1));
1503             else
1504                 prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ...
1505                     '   %f\n'],cputms(1));
1506             end
1507         else
1508             prepstr = sprintf(['Pre-processing:' ...
1509                 '                            %f\n'],cputms(1));
1510         end
1511         sstr = sprintf('***********CPU Timing Results in seconds***********');
1512         ds = sprintf(['\n' sstr '\n' ...
1513             prepstr ...
1514             'ARPACK''s %s:                           %f\n' ...
1515             innerstr ...
1516             'Post-processing with ARPACK''s %s:      %f\n' ...
1517             '***************************************************\n' ...
1518             'Total:                                     %f\n' ...
1519             sstr '\n'], ...
1520             aupdfun,cputms(2),eupdfun,cputms(4),cputms(5));
1521         disp(ds)
1522     end % printTimings
1523
1524 %-------------------------------------------------------------------------%
1525 % End of nested functions
1526 %-------------------------------------------------------------------------%
1527
1528 end % EIGS
1529
1530 %-------------------------------------------------------------------------%
1531 % Subfunctions
1532 %-------------------------------------------------------------------------%
1533 function tf = ishermitian(A)
1534 %ISHERMITIAN
1535 tf = isequal(A,A');
1536 end % ishermititan
1537
1538 function BisHpd = checkTridiagForHSD(alpha, beta)
1539 % CHECKTRIDIAGFORHSD
1540 %   Uses Sturm sequence on alpha (diagonal) and beta (superdiagonal) to
1541 %   determine if the matrix diag(alpha,0) + diag(beta,1) + diag(beta,-1) is
1542 %   Positive Semi-definite.
1543 n = length(alpha);
1544 BisHpd = true;
1545 d = alpha(1);
1546 if d < 0
1547     BisHpd = false;
1548     return;
1549 end
1550 for k = 1:(n-1)
1551     if d == 0
1552         d = eps*(abs(beta(k))+eps);
1553     end
1554     d = alpha(k+1) - beta(k)*(beta(k)/d);
1555     if d < 0
1556         BisHpd = false;
1557         return;
1558     end
1559 end
1560 end % checkTridiagForHSD
1561 %-------------------------------------------------------------------------%
1562 % End of subfunctions
1563 %-------------------------------------------------------------------------%