From d17b33131c14864bd1eae275f49a3f148e21cf29 Mon Sep 17 00:00:00 2001 From: Leo Chan Date: Thu, 22 Oct 2020 01:53:21 -0400 Subject: Squashed commit of the sb-vbs branch. Includes the SD-VBS benchmarks modified to: - Use libextra to loop as realtime jobs - Preallocate memory before starting their main computation - Accept input via stdin instead of via argc Does not include the SD-VBS matlab code. Fixes libextra execution in LITMUS^RT. --- SD-VBS/common/toolbox/MultiNcut/tim_eigs.m | 1084 ++++++++++++++++++++++++++++ 1 file changed, 1084 insertions(+) create mode 100755 SD-VBS/common/toolbox/MultiNcut/tim_eigs.m (limited to 'SD-VBS/common/toolbox/MultiNcut/tim_eigs.m') diff --git a/SD-VBS/common/toolbox/MultiNcut/tim_eigs.m b/SD-VBS/common/toolbox/MultiNcut/tim_eigs.m new file mode 100755 index 0000000..4b80262 --- /dev/null +++ b/SD-VBS/common/toolbox/MultiNcut/tim_eigs.m @@ -0,0 +1,1084 @@ +function varargout = tim_eigs(varargin) + +nombre_A_times_X = 0; %tim +nombreIterations = 0; %tim + +%seule difference avec eigs : +% arguments_Afun = varargin{7-Amatrix-Bnotthere:end}; +%(pour l'instant : n'accepte que 2 arguments dans le cas de Afun : Afun(W,X)) +%permet d'aller plus vite en minimisant les acces a varargin +%(Timothee) + +%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(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 including 0, EIGS finds the eigenvalues +% closest to SIGMA. For scalar SIGMA, and also when SIGMA = 'SM' which uses +% the same algorithm as SIGMA = 0, B need only be symmetric (or Hermitian) +% positive semi-definite since it is not Cholesky factored as in the other cases. +% +% 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: Ritz estimate residual <= tol*NORM(A) [scalar | {eps}] +% OPTS.maxit: maximum number of iterations [integer | {300}] +% OPTS.p: number of Lanczos vectors: K+1 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 = '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('MATLAB:eigs:NonIntegerSize',['%s\n ' ... + 'Rounding input size.'],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 = ['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('MATLAB:eigs:NonIntegerEigQty',['%s\n ' ... + 'Rounding number of eigenvalues.'],kstr) + k = round(k); +end + +whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string.']); +if (nargin < (5-Amatrix-Bnotthere)) + % default: eigs('LM') => ARPACK which='LM', sigma=0 + eigs_sigma = 'LM'; + whch = 'LM'; + sigma = 0; +else + eigs_sigma = varargin{5-Amatrix-Bnotthere}; + if isstr(eigs_sigma) + % eigs(string) => ARPACK which=string, sigma=0 + if ~isequal(size(eigs_sigma),[1,2]) + whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ... + 'LM','SM','LA','SA','BE')]; + whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ... + 'LM','SM','LR','SR','LI','SI')]; + error(whchstr) + end + eigs_sigma = upper(eigs_sigma); + if isequal(eigs_sigma,'SM') + % eigs('SM') => ARPACK which='LM', sigma=0 + whch = 'LM'; + else + % eigs(string), where string~='SM' => ARPACK which=string, sigma=0 + whch = eigs_sigma; + end + sigma = 0; + else + % eigs(scalar) => ARPACK which='LM', sigma=scalar + if ~isa(eigs_sigma,'double') | ~isequal(size(eigs_sigma),[1,1]) + error('Eigenvalue shift sigma must be a scalar.') + end + sigma = eigs_sigma; + if issparse(sigma) + sigma = full(sigma); + end + isrealprob = isrealprob & isreal(sigma); + whch = 'LM'; + 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('MATLAB:eigs:IgnoredOptionPermB', ... + ['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 = ['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('MATLAB:eigs:NonIntegerVecQty',['%s\n ' ... + 'Rounding number of basis vectors.'],pstr) + p = round(p); + end + end + + if isfield(opts,'maxit') + maxit = opts.maxit; + str = ['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('MATLAB:eigs:NonIntegerIterationQty',['%s\n ' ... + 'Rounding number of iterations.'],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; + dispstr = '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('MATLAB:eigs:NonIntegerDiagnosticLevel', ... + '%s\n Rounding diagnostic level.',dispstr) + display = round(display); + end + end + + if isfield(opts,'cheb') + warning('MATLAB:eigs:ObsoleteOptionCheb', ... + ['Ignoring polynomial acceleration opts.cheb' ... + ' (no longer an option).']); + end + if isfield(opts,'stagtol') + warning('MATLAB:eigs:ObsoleteOptionStagtol', ... + ['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('MATLAB:eigs:TooManyRequestedEigsForRealSym', ... + '%s',knstr) + end + useeig = 1; + end +else + if (k > n-2) + if (n >= 7) + warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ... + '%s',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('MATLAB:eigs:ComplexShiftForRealProblem', ... + ['Complex eigenvalue shift sigma on a Hermitian problem' ... + ' (all real eigenvalues).']) + end +end + +if isrealprob & issymA + if strcmp(whch,'LR') + whch = 'LA'; + warning('MATLAB:eigs:SigmaChangedToLA', ... + ['For real symmetric problems, sigma value ''LR''' ... + ' (Largest Real) is now ''LA'' (Largest Algebraic).']) + end + if strcmp(whch,'SR') + whch = 'SA'; + warning('MATLAB:eigs:SigmaChangedToSA', ... + ['For real symmetric problems, sigma value ''SR''' ... + ' (Smallest Real) is now ''SA'' (Smallest Algebraic).']) + end + if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'}) + whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ... + 'LM','SM','LA','SA','BE')]; + error(whchstr) + end +else + if strcmp(whch,'BE') + warning('MATLAB:eigs:SigmaChangedToLM', ... + ['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'}) + whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ... + '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 + +if isequal(whch,'LM') & ~isequal(eigs_sigma,'LM') + % A*x = lambda*M*x, M symmetric (positive) semi-definite + % => OP = inv(A - sigma*M)*M and B = M + % => shift-and-invert mode + mode = 3; +elseif isempty(B) + % A*x = lambda*x + % => OP = A and B = I + mode = 1; +else % B is not empty + % Do not use mode=2. + % Use mode = 1 with OP = R'\(A*(R\x)) and B = I + % where R is B's upper triangular Cholesky factor: B = R'*R. + % Finally, V = R\V returns the actual generalized eigenvectors of A and B. + mode = 1; +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 + [L,U,P,Q] = lu(AsB); + [perm,dummy] = find(Q); + 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 = full(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'])]; + warning('MATLAB:eigs:SigmaNearExactEig',ds) + end +end + +if (mode == 1) & ~isempty(B) & ~cholB % need chol(B) + if issparse(B) + permB = symamd(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) | (~isempty(B) & (mode == 1)) + 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; + + +%tim + + +indexArgumentsAfun = 7-Amatrix-Bnotthere:length(varargin); +nbArgumentsAfun = length(indexArgumentsAfun); +if nbArgumentsAfun >=1 + arguments_Afun = varargin{7-Amatrix-Bnotthere}; +end +if nbArgumentsAfun >=2 + arguments_Afun2 = varargin{7-Amatrix-Bnotthere+1}; +end +if nbArgumentsAfun >=3 + arguments_Afun3 = varargin{7-Amatrix-Bnotthere+2}; +end +%fin tim + + + +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,full(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))); + %tim + row = mod(double(ipntr(1))-1,n) + 1 ; + col1 = floor((double(ipntr(1))-1)/n) + 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.'],full(ipntr(1)),n); + error(str) + end + + + + %[row,col2] = ind2sub([n,3],double(ipntr(2))); + %tim + row = mod(double(ipntr(2))-1,n) + 1 ; + col2 = floor((double(ipntr(2))-1)/n) + 1; + + + + if (row ~= 1) + str = sprintf(['ipntr(2)=%d does not refer to the start of a' ... + ' column of the %d-by-3 array workd.'],full(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.'],full(ipntr(3)),n); + error(str) + end + end + + switch (ido) + case {-1,1} + if Amatrix + if (mode == 1) + if isempty(B) + % OP = A*x + workd(:,col2) = A * workd(:,col1); + else + % OP = R'\(A*(R\x)) + if issparse(B) + workd(permB,col2) = RB \ workd(:,col1); + workd(:,col2) = A * workd(:,col2); + workd(:,col2) = RBT \ workd(permB,col2); + else + workd(:,col2) = RBT \ (A * (RB \ workd(:,col1))); + end + 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 or 3 + error(['Unknown mode returned from ' prefix 'aupd.']) + end + else % A is not a matrix + if (mode == 1) + if isempty(B) + % OP = A*x + %workd(:,col2) = feval(A,workd(:,col1),varargin{7-Amatrix-Bnotthere:end}); + + + + + + nombre_A_times_X = nombre_A_times_X + 1; + + + + pause(0); %voir + + if nbArgumentsAfun == 1 + workd(:,col2) = feval(A,workd(:,col1),arguments_Afun); + %workd(:,col2) = max(workd(:,col2),0); + elseif nbArgumentsAfun == 2 + workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2); + elseif nbArgumentsAfun == 3 + workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2,arguments_Afun3); + else + workd(:,col2) = feval(A,workd(:,col1),varargin{indexArgumentsAfun}); + end + %workd(:,col2) = tim_w_times_x_c(workd(:,col1),arguments_Afun); %slower + + else + % OP = R'\(A*(R\x)) + if issparse(B) + workd(permB,col2) = RB \ workd(:,col1); + workd(:,col2) = feval(A,workd(:,col2),arguments_Afun); + workd(:,col2) = RBT \ workd(permB,col2); + + else + workd(:,col2) = RBT \ feval(A,(RB\workd(:,col1)),arguments_Afun); + end + end + elseif (mode == 3) + if isempty(B) + workd(:,col2) = feval(A,workd(:,col1), arguments_Afun); + 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), arguments_Afun); + elseif (ido == 1) + workd(:,col2) = feval(A,workd(:,col3), arguments_Afun); + end + end + else % mode is not 1 or 3 + error(['Unknown mode returned from ' prefix 'aupd.']) + end + end % if Amatrix + case 2 + if (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 + case 3 + % setting iparam(1) = ishift = 1 ensures this never happens + warning('MATLAB:eigs:WorklShiftsUnsupported', ... + ['EIGS does not yet support computing the shifts in workl.' ... + ' Setting reverse communication parameter to 99 and returning.']) + ido = int32(99); + case 99 + otherwise + error(['Unknown value of reverse communication parameter' ... + ' returned from ' prefix 'aupd.']) + end + + cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X) + + %tim + if nombreIterations ~= double(ipntr(15)) + nombreIterations = double(ipntr(15)); + end + + if display >= 1 && display <=2 + 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(max(end-2*k+1,1):end)) + else + % k eigenvalues + disp(dispvec(max(end-k+1,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(max(end-k,1):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(max(end-k+1,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,full(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',full(info))]; + end + error(es) + else + nconv = double(iparam(5)); + if (nconv == 0) + if (nargout < 3) + warning('MATLAB:eigs:NoEigsConverged', ... + 'None of the %d requested eigenvalues converged.',k) + else + flag = 1; + end + elseif (nconv < k) + if (nargout < 3) + warning('MATLAB:eigs:NotAllEigsConverged', ... + 'Only %d of the %d requested eigenvalues converged.',nconv,k) + 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 + +if (nargout >= 2) & (mode == 1) & ~isempty(B) + if issparse(B) + varargout{1}(permB,:) = RB \ varargout{1}; + else + varargout{1} = RB \ varargout{1}; + end +end + +cputms(4) = cputime-t0; % end timing post-processing + +cputms(5) = sum(cputms(1:4)); % total time + +if (display >= 2) %tim + 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) + if nombre_A_times_X > 0 %tim + disp(sprintf('Number of iterations : %f\n',nombreIterations)); + disp(sprintf('Number of times A*X was computed : %f\n',nombre_A_times_X)); + disp(sprintf('Average time for A*X : %f\n',cputms(3)/nombre_A_times_X)); + end +end -- cgit v1.2.2