summaryrefslogtreecommitdiffstats
path: root/SD-VBS/common/toolbox/MultiNcut/tim_eigs.m
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/common/toolbox/MultiNcut/tim_eigs.m')
-rwxr-xr-xSD-VBS/common/toolbox/MultiNcut/tim_eigs.m1084
1 files changed, 1084 insertions, 0 deletions
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 @@
1function varargout = tim_eigs(varargin)
2
3nombre_A_times_X = 0; %tim
4nombreIterations = 0; %tim
5
6%seule difference avec eigs :
7% arguments_Afun = varargin{7-Amatrix-Bnotthere:end};
8%(pour l'instant : n'accepte que 2 arguments dans le cas de Afun : Afun(W,X))
9%permet d'aller plus vite en minimisant les acces a varargin
10%(Timothee)
11
12%EIGS Find a few eigenvalues and eigenvectors of a matrix using ARPACK.
13% D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues.
14% A must be square and should be large and sparse.
15%
16% [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude
17% eigenvalues and a matrix V whose columns are the corresponding eigenvectors.
18%
19% [V,D,FLAG] = EIGS(A) also returns a convergence flag. If FLAG is 0
20% then all the eigenvalues converged; otherwise not all converged.
21%
22% EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D. B must
23% be symmetric (or Hermitian) positive definite and the same size as A.
24% EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D.
25%
26% EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues.
27%
28% EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues based on SIGMA:
29% 'LM' or 'SM' - Largest or Smallest Magnitude
30% For real symmetric problems, SIGMA may also be:
31% 'LA' or 'SA' - Largest or Smallest Algebraic
32% 'BE' - Both Ends, one more from high end if K is odd
33% For nonsymmetric and complex problems, SIGMA may also be:
34% 'LR' or 'SR' - Largest or Smallest Real part
35% 'LI' or 'SI' - Largest or Smallest Imaginary part
36% If SIGMA is a real or complex scalar including 0, EIGS finds the eigenvalues
37% closest to SIGMA. For scalar SIGMA, and also when SIGMA = 'SM' which uses
38% the same algorithm as SIGMA = 0, B need only be symmetric (or Hermitian)
39% positive semi-definite since it is not Cholesky factored as in the other cases.
40%
41% EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options:
42% OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{0} | 1]
43% OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [0 | {1}]
44% OPTS.tol: convergence: Ritz estimate residual <= tol*NORM(A) [scalar | {eps}]
45% OPTS.maxit: maximum number of iterations [integer | {300}]
46% OPTS.p: number of Lanczos vectors: K+1<p<=N [integer | {2K}]
47% OPTS.v0: starting vector [N-by-1 vector | {randomly generated by ARPACK}]
48% OPTS.disp: diagnostic information display level [0 | {1} | 2]
49% OPTS.cholB: B is actually its Cholesky factor CHOL(B) [{0} | 1]
50% OPTS.permB: sparse B is actually CHOL(B(permB,permB)) [permB | {1:N}]
51%
52% EIGS(AFUN,N) accepts the function AFUN instead of the matrix A.
53% Y = AFUN(X) should return
54% A*X if SIGMA is not specified, or is a string other than 'SM'
55% A\X if SIGMA is 0 or 'SM'
56% (A-SIGMA*I)\X if SIGMA is a nonzero scalar (standard eigenvalue problem)
57% (A-SIGMA*B)\X if SIGMA is a nonzero scalar (generalized eigenvalue problem)
58% N is the size of A. The matrix A, A-SIGMA*I or A-SIGMA*B represented by AFUN is
59% assumed to be real and nonsymmetric unless specified otherwise by OPTS.isreal
60% and OPTS.issym. In all these EIGS syntaxes, EIGS(A,...) may be replaced by
61% EIGS(AFUN,N,...).
62%
63% EIGS(AFUN,N,K,SIGMA,OPTS,P1,P2,...) and EIGS(AFUN,N,B,K,SIGMA,OPTS,P1,P2,...)
64% provide for additional arguments which are passed to AFUN(X,P1,P2,...).
65%
66% Examples:
67% A = delsq(numgrid('C',15)); d1 = eigs(A,5,'SM');
68% Equivalently, if dnRk is the following one-line function:
69% function y = dnRk(x,R,k)
70% y = (delsq(numgrid(R,k))) \ x;
71% then pass dnRk's additional arguments, 'C' and 15, to EIGS:
72% n = size(A,1); opts.issym = 1; d2 = eigs(@dnRk,n,5,'SM',opts,'C',15);
73%
74% See also EIG, SVDS, ARPACKC.
75
76% Copyright 1984-2002 The MathWorks, Inc.
77% $Revision: 1.45 $ $Date: 2002/05/14 18:50:58 $
78
79cputms = zeros(5,1);
80t0 = cputime; % start timing pre-processing
81
82if (nargout > 3)
83 error('Too many output arguments.')
84end
85
86% Process inputs and do error-checking
87if isa(varargin{1},'double')
88 A = varargin{1};
89 Amatrix = 1;
90else
91 A = fcnchk(varargin{1});
92 Amatrix = 0;
93end
94
95isrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma)
96if Amatrix
97 isrealprob = isreal(A);
98end
99
100issymA = 0;
101if Amatrix
102 issymA = isequal(A,A');
103end
104
105if Amatrix
106 [m,n] = size(A);
107 if (m ~= n)
108 error('A must be a square matrix or a function which computes A*x.')
109 end
110else
111 n = varargin{2};
112 nstr = 'Size of problem, ''n'', must be a positive integer.';
113 if ~isequal(size(n),[1,1]) | ~isreal(n)
114 error(nstr)
115 end
116 if (round(n) ~= n)
117 warning('MATLAB:eigs:NonIntegerSize',['%s\n ' ...
118 'Rounding input size.'],nstr)
119 n = round(n);
120 end
121 if issparse(n)
122 n = full(n);
123 end
124end
125
126Bnotthere = 0;
127Bstr = sprintf(['Generalized matrix B must be the same size as A and' ...
128 ' either a symmetric positive (semi-)definite matrix or' ...
129 ' its Cholesky factor.']);
130if (nargin < (3-Amatrix-Bnotthere))
131 B = [];
132 Bnotthere = 1;
133else
134 Bk = varargin{3-Amatrix-Bnotthere};
135 if isempty(Bk) % allow eigs(A,[],k,sigma,opts);
136 B = Bk;
137 else
138 if isequal(size(Bk),[1,1]) & (n ~= 1)
139 B = [];
140 k = Bk;
141 Bnotthere = 1;
142 else % eigs(9,8,...) assumes A=9, B=8, ... NOT A=9, k=8, ...
143 B = Bk;
144 if ~isa(B,'double') | ~isequal(size(B),[n,n])
145 error(Bstr)
146 end
147 isrealprob = isrealprob & isreal(B);
148 end
149 end
150end
151
152if Amatrix & ((nargin - ~Bnotthere)>4)
153 error('Too many inputs.')
154end
155
156if (nargin < (4-Amatrix-Bnotthere))
157 k = min(n,6);
158else
159 k = varargin{4-Amatrix-Bnotthere};
160end
161
162kstr = ['Number of eigenvalues requested, k, must be a' ...
163 ' positive integer <= n.'];
164if ~isa(k,'double') | ~isequal(size(k),[1,1]) | ~isreal(k) | (k>n)
165 error(kstr)
166end
167if issparse(k)
168 k = full(k);
169end
170if (round(k) ~= k)
171 warning('MATLAB:eigs:NonIntegerEigQty',['%s\n ' ...
172 'Rounding number of eigenvalues.'],kstr)
173 k = round(k);
174end
175
176whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string.']);
177if (nargin < (5-Amatrix-Bnotthere))
178 % default: eigs('LM') => ARPACK which='LM', sigma=0
179 eigs_sigma = 'LM';
180 whch = 'LM';
181 sigma = 0;
182else
183 eigs_sigma = varargin{5-Amatrix-Bnotthere};
184 if isstr(eigs_sigma)
185 % eigs(string) => ARPACK which=string, sigma=0
186 if ~isequal(size(eigs_sigma),[1,2])
187 whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
188 'LM','SM','LA','SA','BE')];
189 whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ...
190 'LM','SM','LR','SR','LI','SI')];
191 error(whchstr)
192 end
193 eigs_sigma = upper(eigs_sigma);
194 if isequal(eigs_sigma,'SM')
195 % eigs('SM') => ARPACK which='LM', sigma=0
196 whch = 'LM';
197 else
198 % eigs(string), where string~='SM' => ARPACK which=string, sigma=0
199 whch = eigs_sigma;
200 end
201 sigma = 0;
202 else
203 % eigs(scalar) => ARPACK which='LM', sigma=scalar
204 if ~isa(eigs_sigma,'double') | ~isequal(size(eigs_sigma),[1,1])
205 error('Eigenvalue shift sigma must be a scalar.')
206 end
207 sigma = eigs_sigma;
208 if issparse(sigma)
209 sigma = full(sigma);
210 end
211 isrealprob = isrealprob & isreal(sigma);
212 whch = 'LM';
213 end
214end
215
216tol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS)
217maxit = [];
218p = [];
219info = int32(0); % use a random starting vector
220display = 1;
221cholB = 0;
222
223if (nargin >= (6-Amatrix-Bnotthere))
224 opts = varargin{6-Amatrix-Bnotthere};
225 if ~isa(opts,'struct')
226 error('Options argument must be a structure.')
227 end
228
229 if isfield(opts,'issym') & ~Amatrix
230 issymA = opts.issym;
231 if (issymA ~= 0) & (issymA ~= 1)
232 error('opts.issym must be 0 or 1.')
233 end
234 end
235
236 if isfield(opts,'isreal') & ~Amatrix
237 if (opts.isreal ~= 0) & (opts.isreal ~= 1)
238 error('opts.isreal must be 0 or 1.')
239 end
240 isrealprob = isrealprob & opts.isreal;
241 end
242
243 if ~isempty(B) & (isfield(opts,'cholB') | isfield(opts,'permB'))
244 if isfield(opts,'cholB')
245 cholB = opts.cholB;
246 if (cholB ~= 0) & (cholB ~= 1)
247 error('opts.cholB must be 0 or 1.')
248 end
249 if isfield(opts,'permB')
250 if issparse(B) & cholB
251 permB = opts.permB;
252 if ~isequal(sort(permB),(1:n)) & ...
253 ~isequal(sort(permB),(1:n)')
254 error('opts.permB must be a permutation of 1:n.')
255 end
256 else
257 warning('MATLAB:eigs:IgnoredOptionPermB', ...
258 ['Ignoring opts.permB since B is not its sparse' ...
259 ' Cholesky factor.'])
260 end
261 else
262 permB = 1:n;
263 end
264 end
265 end
266
267 if isfield(opts,'tol')
268 if ~isequal(size(opts.tol),[1,1]) | ~isreal(opts.tol) | (opts.tol<=0)
269 error(['Convergence tolerance opts.tol must be a strictly' ...
270 ' positive real scalar.'])
271 else
272 tol = full(opts.tol);
273 end
274 end
275
276 if isfield(opts,'p')
277 p = opts.p;
278 pstr = ['Number of basis vectors opts.p must be a positive' ...
279 ' integer <= n.'];
280 if ~isequal(size(p),[1,1]) | ~isreal(p) | (p<=0) | (p>n)
281 error(pstr)
282 end
283 if issparse(p)
284 p = full(p);
285 end
286 if (round(p) ~= p)
287 warning('MATLAB:eigs:NonIntegerVecQty',['%s\n ' ...
288 'Rounding number of basis vectors.'],pstr)
289 p = round(p);
290 end
291 end
292
293 if isfield(opts,'maxit')
294 maxit = opts.maxit;
295 str = ['Maximum number of iterations opts.maxit must be' ...
296 ' a positive integer.'];
297 if ~isequal(size(maxit),[1,1]) | ~isreal(maxit) | (maxit<=0)
298 error(str)
299 end
300 if issparse(maxit)
301 maxit = full(maxit);
302 end
303 if (round(maxit) ~= maxit)
304 warning('MATLAB:eigs:NonIntegerIterationQty',['%s\n ' ...
305 'Rounding number of iterations.'],str)
306 maxit = round(maxit);
307 end
308 end
309
310 if isfield(opts,'v0')
311 if ~isequal(size(opts.v0),[n,1])
312 error('Start vector opts.v0 must be n-by-1.')
313 end
314 if isrealprob
315 if ~isreal(opts.v0)
316 error(['Start vector opts.v0 must be real for real problems.'])
317 end
318 resid = full(opts.v0);
319 else
320 resid(1:2:(2*n-1),1) = full(real(opts.v0));
321 resid(2:2:2*n,1) = full(imag(opts.v0));
322 end
323 info = int32(1); % use resid as starting vector
324 end
325
326 if isfield(opts,'disp')
327 display = opts.disp;
328 dispstr = 'Diagnostic level opts.disp must be an integer.';
329 if (~isequal(size(display),[1,1])) | (~isreal(display)) | (display<0)
330 error(dispstr)
331 end
332 if (round(display) ~= display)
333 warning('MATLAB:eigs:NonIntegerDiagnosticLevel', ...
334 '%s\n Rounding diagnostic level.',dispstr)
335 display = round(display);
336 end
337 end
338
339 if isfield(opts,'cheb')
340 warning('MATLAB:eigs:ObsoleteOptionCheb', ...
341 ['Ignoring polynomial acceleration opts.cheb' ...
342 ' (no longer an option).']);
343 end
344 if isfield(opts,'stagtol')
345 warning('MATLAB:eigs:ObsoleteOptionStagtol', ...
346 ['Ignoring stagnation tolerance opts.stagtol' ...
347 ' (no longer an option).']);
348 end
349
350end
351
352% Now we know issymA, isrealprob, cholB, and permB
353
354if isempty(p)
355 if isrealprob & ~issymA
356 p = min(max(2*k+1,20),n);
357 else
358 p = min(max(2*k,20),n);
359 end
360end
361if isempty(maxit)
362 maxit = max(300,ceil(2*n/max(p,1)));
363end
364if (info == int32(0))
365 if isrealprob
366 resid = zeros(n,1);
367 else
368 resid = zeros(2*n,1);
369 end
370end
371
372if ~isempty(B) % B must be symmetric (Hermitian) positive (semi-)definite
373 if cholB
374 if ~isequal(triu(B),B)
375 error(Bstr)
376 end
377 else
378 if ~isequal(B,B')
379 error(Bstr)
380 end
381 end
382end
383
384useeig = 0;
385if isrealprob & issymA
386 knstr = sprintf(['For real symmetric problems, must have number' ...
387 ' of eigenvalues k < n.\n']);
388else
389 knstr = sprintf(['For nonsymmetric and complex problems, must have' ...
390 ' number of eigenvalues k < n-1.\n']);
391end
392if isempty(B)
393 knstr = [knstr sprintf(['Using eig(full(A)) instead.'])];
394else
395 knstr = [knstr sprintf(['Using eig(full(A),full(B)) instead.'])];
396end
397if (k == 0)
398 useeig = 1;
399end
400if isrealprob & issymA
401 if (k > n-1)
402 if (n >= 6)
403 warning('MATLAB:eigs:TooManyRequestedEigsForRealSym', ...
404 '%s',knstr)
405 end
406 useeig = 1;
407 end
408else
409 if (k > n-2)
410 if (n >= 7)
411 warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ...
412 '%s',knstr)
413 end
414 useeig = 1;
415 end
416end
417
418if isrealprob & issymA
419 if ~isreal(sigma)
420 error(['For real symmetric problems, eigenvalue shift sigma must' ...
421 ' be real.'])
422 end
423else
424 if ~isrealprob & issymA & ~isreal(sigma)
425 warning('MATLAB:eigs:ComplexShiftForRealProblem', ...
426 ['Complex eigenvalue shift sigma on a Hermitian problem' ...
427 ' (all real eigenvalues).'])
428 end
429end
430
431if isrealprob & issymA
432 if strcmp(whch,'LR')
433 whch = 'LA';
434 warning('MATLAB:eigs:SigmaChangedToLA', ...
435 ['For real symmetric problems, sigma value ''LR''' ...
436 ' (Largest Real) is now ''LA'' (Largest Algebraic).'])
437 end
438 if strcmp(whch,'SR')
439 whch = 'SA';
440 warning('MATLAB:eigs:SigmaChangedToSA', ...
441 ['For real symmetric problems, sigma value ''SR''' ...
442 ' (Smallest Real) is now ''SA'' (Smallest Algebraic).'])
443 end
444 if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'})
445 whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ...
446 'LM','SM','LA','SA','BE')];
447 error(whchstr)
448 end
449else
450 if strcmp(whch,'BE')
451 warning('MATLAB:eigs:SigmaChangedToLM', ...
452 ['Sigma value ''BE'' is now only available for real' ...
453 ' symmetric problems. Computing ''LM'' eigenvalues instead.'])
454 whch = 'LM';
455 end
456 if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'})
457 whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ...
458 'LM','SM','LR','SR','LI','SI')];
459 error(whchstr)
460 end
461end
462
463% Now have enough information to do early return on cases eigs does not handle
464if useeig
465 if (nargout <= 1)
466 varargout{1} = eigs2(A,n,B,k,whch,sigma,cholB, ...
467 varargin{7-Amatrix-Bnotthere:end});
468 else
469 [varargout{1},varargout{2}] = eigs2(A,n,B,k,whch,sigma,cholB, ...
470 varargin{7-Amatrix-Bnotthere:end});
471 end
472 if (nargout == 3)
473 varargout{3} = 0; % flag indicates "convergence"
474 end
475 return
476end
477
478if isrealprob & ~issymA
479 sigmar = real(sigma);
480 sigmai = imag(sigma);
481end
482
483if isrealprob & issymA
484 if (p <= k)
485 error(['For real symmetric problems, must have number of' ...
486 ' basis vectors opts.p > k.'])
487 end
488else
489 if (p <= k+1)
490 error(['For nonsymmetric and complex problems, must have number of' ...
491 ' basis vectors opts.p > k+1.'])
492 end
493end
494
495if isequal(whch,'LM') & ~isequal(eigs_sigma,'LM')
496 % A*x = lambda*M*x, M symmetric (positive) semi-definite
497 % => OP = inv(A - sigma*M)*M and B = M
498 % => shift-and-invert mode
499 mode = 3;
500elseif isempty(B)
501 % A*x = lambda*x
502 % => OP = A and B = I
503 mode = 1;
504else % B is not empty
505 % Do not use mode=2.
506 % Use mode = 1 with OP = R'\(A*(R\x)) and B = I
507 % where R is B's upper triangular Cholesky factor: B = R'*R.
508 % Finally, V = R\V returns the actual generalized eigenvectors of A and B.
509 mode = 1;
510end
511
512if cholB
513 pB = 0;
514 RB = B;
515 RBT = B';
516end
517
518if (mode == 3) & Amatrix % need lu(A-sigma*B)
519 if issparse(A) & (isempty(B) | issparse(B))
520 if isempty(B)
521 AsB = A - sigma * speye(n);
522 else
523 if cholB
524 AsB = A - sigma * RBT * RB;
525 else
526 AsB = A - sigma * B;
527 end
528 end
529 [L,U,P,Q] = lu(AsB);
530 [perm,dummy] = find(Q);
531 else
532 if isempty(B)
533 AsB = A - sigma * eye(n);
534 else
535 if cholB
536 AsB = A - sigma * RBT * RB;
537 else
538 AsB = A - sigma * B;
539 end
540 end
541 [L,U,P] = lu(AsB);
542 end
543 dU = diag(U);
544 rcondestU = full(min(abs(dU)) / max(abs(dU)));
545 if (rcondestU < eps)
546 if isempty(B)
547 ds = sprintf(['(A-sigma*I) has small reciprocal condition' ...
548 ' estimate: %f\n'],rcondestU);
549 else
550 ds = sprintf(['(A-sigma*B) has small reciprocal condition' ...
551 ' estimate: %f\n'],rcondestU);
552 end
553 ds = [ds sprintf(['indicating that sigma is near an exact' ...
554 ' eigenvalue. The\nalgorithm may not converge unless' ...
555 ' you try a new value for sigma.\n'])];
556 warning('MATLAB:eigs:SigmaNearExactEig',ds)
557 end
558end
559
560if (mode == 1) & ~isempty(B) & ~cholB % need chol(B)
561 if issparse(B)
562 permB = symamd(B);
563 [RB,pB] = chol(B(permB,permB));
564 else
565 [RB,pB] = chol(B);
566 end
567 if (pB == 0)
568 RBT = RB';
569 else
570 error(Bstr)
571 end
572end
573
574% Allocate outputs and ARPACK work variables
575if isrealprob
576 if issymA % real and symmetric
577 prefix = 'ds';
578 v = zeros(n,p);
579 ldv = int32(size(v,1));
580 ipntr = int32(zeros(15,1));
581 workd = zeros(n,3);
582 lworkl = p*(p+8);
583 workl = zeros(lworkl,1);
584 lworkl = int32(lworkl);
585 d = zeros(k,1);
586 else % real but not symmetric
587 prefix = 'dn';
588 v = zeros(n,p);
589 ldv = int32(size(v,1));
590 ipntr = int32(zeros(15,1));
591 workd = zeros(n,3);
592 lworkl = 3*p*(p+2);
593 workl = zeros(lworkl,1);
594 lworkl = int32(lworkl);
595 workev = zeros(3*p,1);
596 d = zeros(k+1,1);
597 di = zeros(k+1,1);
598 end
599else % complex
600 prefix = 'zn';
601 zv = zeros(2*n*p,1);
602 ldv = int32(n);
603 ipntr = int32(zeros(15,1));
604 workd = complex(zeros(n,3));
605 zworkd = zeros(2*prod(size(workd)),1);
606 lworkl = 3*p^2+5*p;
607 workl = zeros(2*lworkl,1);
608 lworkl = int32(lworkl);
609 workev = zeros(2*2*p,1);
610 zd = zeros(2*(k+1),1);
611 rwork = zeros(p,1);
612end
613
614ido = int32(0); % reverse communication parameter
615if isempty(B) | (~isempty(B) & (mode == 1))
616 bmat = 'I'; % standard eigenvalue problem
617else
618 bmat = 'G'; % generalized eigenvalue problem
619end
620nev = int32(k); % number of eigenvalues requested
621ncv = int32(p); % number of Lanczos vectors
622iparam = int32(zeros(11,1));
623iparam([1 3 7]) = int32([1 maxit mode]);
624select = int32(zeros(p,1));
625
626cputms(1) = cputime - t0; % end timing pre-processing
627
628iter = 0;
629ariter = 0;
630
631
632%tim
633
634
635indexArgumentsAfun = 7-Amatrix-Bnotthere:length(varargin);
636nbArgumentsAfun = length(indexArgumentsAfun);
637if nbArgumentsAfun >=1
638 arguments_Afun = varargin{7-Amatrix-Bnotthere};
639end
640if nbArgumentsAfun >=2
641 arguments_Afun2 = varargin{7-Amatrix-Bnotthere+1};
642end
643if nbArgumentsAfun >=3
644 arguments_Afun3 = varargin{7-Amatrix-Bnotthere+2};
645end
646%fin tim
647
648
649
650while (ido ~= 99)
651
652 t0 = cputime; % start timing ARPACK calls **aupd
653
654 if isrealprob
655 arpackc( [prefix 'aupd'], ido, ...
656 bmat, int32(n), whch, nev, tol, resid, ncv, ...
657 v, ldv, iparam, ipntr, workd, workl, lworkl, info);
658 else
659 zworkd(1:2:end-1) = real(workd);
660 zworkd(2:2:end) = imag(workd);
661 arpackc( 'znaupd', ido, ...
662 bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
663 ldv, iparam, ipntr, zworkd, workl, lworkl, ...
664 rwork, info );
665 workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]);
666 end
667
668 if (info < 0)
669 es = sprintf('Error with ARPACK routine %saupd: info = %d',...
670 prefix,full(info));
671 error(es)
672 end
673
674 cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd
675 t0 = cputime; % start timing MATLAB OP(X)
676
677 % Compute which columns of workd ipntr references
678
679
680
681
682
683 %[row,col1] = ind2sub([n,3],double(ipntr(1)));
684 %tim
685 row = mod(double(ipntr(1))-1,n) + 1 ;
686 col1 = floor((double(ipntr(1))-1)/n) + 1;
687
688
689 if (row ~= 1)
690 str = sprintf(['ipntr(1)=%d does not refer to the start of a' ...
691 ' column of the %d-by-3 array workd.'],full(ipntr(1)),n);
692 error(str)
693 end
694
695
696
697 %[row,col2] = ind2sub([n,3],double(ipntr(2)));
698 %tim
699 row = mod(double(ipntr(2))-1,n) + 1 ;
700 col2 = floor((double(ipntr(2))-1)/n) + 1;
701
702
703
704 if (row ~= 1)
705 str = sprintf(['ipntr(2)=%d does not refer to the start of a' ...
706 ' column of the %d-by-3 array workd.'],full(ipntr(2)),n);
707 error(str)
708 end
709 if ~isempty(B) & (mode == 3) & (ido == 1)
710 [row,col3] = ind2sub([n,3],double(ipntr(3)));
711 if (row ~= 1)
712 str = sprintf(['ipntr(3)=%d does not refer to the start of a' ...
713 ' column of the %d-by-3 array workd.'],full(ipntr(3)),n);
714 error(str)
715 end
716 end
717
718 switch (ido)
719 case {-1,1}
720 if Amatrix
721 if (mode == 1)
722 if isempty(B)
723 % OP = A*x
724 workd(:,col2) = A * workd(:,col1);
725 else
726 % OP = R'\(A*(R\x))
727 if issparse(B)
728 workd(permB,col2) = RB \ workd(:,col1);
729 workd(:,col2) = A * workd(:,col2);
730 workd(:,col2) = RBT \ workd(permB,col2);
731 else
732 workd(:,col2) = RBT \ (A * (RB \ workd(:,col1)));
733 end
734 end
735 elseif (mode == 3)
736 if isempty(B)
737 if issparse(A)
738 workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
739 else
740 workd(:,col2) = U \ (L \ (P * workd(:,col1)));
741 end
742 else % B is not empty
743 if (ido == -1)
744 if cholB
745 workd(:,col2) = RBT * (RB * workd(:,col1));
746 else
747 workd(:,col2) = B * workd(:,col1);
748 end
749 if issparse(A) & issparse(B)
750 workd(perm,col2) = U \ (L \ (P * workd(:,col1)));
751 else
752 workd(:,col2) = U \ (L \ (P * workd(:,col1)));
753 end
754 elseif (ido == 1)
755 if issparse(A) & issparse(B)
756 workd(perm,col2) = U \ (L \ (P * workd(:,col3)));
757 else
758 workd(:,col2) = U \ (L \ (P * workd(:,col3)));
759 end
760 end
761 end
762 else % mode is not 1 or 3
763 error(['Unknown mode returned from ' prefix 'aupd.'])
764 end
765 else % A is not a matrix
766 if (mode == 1)
767 if isempty(B)
768 % OP = A*x
769 %workd(:,col2) = feval(A,workd(:,col1),varargin{7-Amatrix-Bnotthere:end});
770
771
772
773
774
775 nombre_A_times_X = nombre_A_times_X + 1;
776
777
778
779 pause(0); %voir
780
781 if nbArgumentsAfun == 1
782 workd(:,col2) = feval(A,workd(:,col1),arguments_Afun);
783 %workd(:,col2) = max(workd(:,col2),0);
784 elseif nbArgumentsAfun == 2
785 workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2);
786 elseif nbArgumentsAfun == 3
787 workd(:,col2) = feval(A,workd(:,col1),arguments_Afun,arguments_Afun2,arguments_Afun3);
788 else
789 workd(:,col2) = feval(A,workd(:,col1),varargin{indexArgumentsAfun});
790 end
791 %workd(:,col2) = tim_w_times_x_c(workd(:,col1),arguments_Afun); %slower
792
793 else
794 % OP = R'\(A*(R\x))
795 if issparse(B)
796 workd(permB,col2) = RB \ workd(:,col1);
797 workd(:,col2) = feval(A,workd(:,col2),arguments_Afun);
798 workd(:,col2) = RBT \ workd(permB,col2);
799
800 else
801 workd(:,col2) = RBT \ feval(A,(RB\workd(:,col1)),arguments_Afun);
802 end
803 end
804 elseif (mode == 3)
805 if isempty(B)
806 workd(:,col2) = feval(A,workd(:,col1), arguments_Afun);
807 else
808 if (ido == -1)
809 if cholB
810 workd(:,col2) = RBT * (RB * workd(:,col1));
811 else
812 workd(:,col2) = B * workd(:,col1);
813 end
814 workd(:,col2) = feval(A,workd(:,col2), arguments_Afun);
815 elseif (ido == 1)
816 workd(:,col2) = feval(A,workd(:,col3), arguments_Afun);
817 end
818 end
819 else % mode is not 1 or 3
820 error(['Unknown mode returned from ' prefix 'aupd.'])
821 end
822 end % if Amatrix
823 case 2
824 if (mode == 3)
825 if cholB
826 workd(:,col2) = RBT * (RB * workd(:,col1));
827 else
828 workd(:,col2) = B * workd(:,col1);
829 end
830 else
831 error(['Unknown mode returned from ' prefix 'aupd.'])
832 end
833 case 3
834 % setting iparam(1) = ishift = 1 ensures this never happens
835 warning('MATLAB:eigs:WorklShiftsUnsupported', ...
836 ['EIGS does not yet support computing the shifts in workl.' ...
837 ' Setting reverse communication parameter to 99 and returning.'])
838 ido = int32(99);
839 case 99
840 otherwise
841 error(['Unknown value of reverse communication parameter' ...
842 ' returned from ' prefix 'aupd.'])
843 end
844
845 cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X)
846
847 %tim
848 if nombreIterations ~= double(ipntr(15))
849 nombreIterations = double(ipntr(15));
850 end
851
852 if display >= 1 && display <=2
853 iter = double(ipntr(15));
854 if (iter > ariter) & (ido ~= 99)
855 ariter = iter;
856 ds = sprintf(['Iteration %d: a few Ritz values of the' ...
857 ' %d-by-%d matrix:'],iter,p,p);
858 disp(ds)
859 if isrealprob
860 if issymA
861 dispvec = [workl(double(ipntr(6))+(0:p-1))];
862 if isequal(whch,'BE')
863 % roughly k Large eigenvalues and k Small eigenvalues
864 disp(dispvec(max(end-2*k+1,1):end))
865 else
866 % k eigenvalues
867 disp(dispvec(max(end-k+1,1):end))
868 end
869 else
870 dispvec = [complex(workl(double(ipntr(6))+(0:p-1)), ...
871 workl(double(ipntr(7))+(0:p-1)))];
872 % k+1 eigenvalues (keep complex conjugate pairs together)
873 disp(dispvec(max(end-k,1):end))
874 end
875 else
876 dispvec = [complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ...
877 workl(2*double(ipntr(6))+(0:2:2*(p-1))))];
878 disp(dispvec(max(end-k+1,1):end))
879 end
880 end
881 end
882
883end % while (ido ~= 99)
884
885t0 = cputime; % start timing post-processing
886
887flag = 0;
888if (info < 0)
889 es = sprintf('Error with ARPACK routine %saupd: info = %d',prefix,full(info));
890 error(es)
891else
892 if (nargout >= 2)
893 rvec = int32(1); % compute eigenvectors
894 else
895 rvec = int32(0); % do not compute eigenvectors
896 end
897
898 if isrealprob
899 if issymA
900 arpackc( 'dseupd', rvec, 'A', select, ...
901 d, v, ldv, sigma, ...
902 bmat, int32(n), whch, nev, tol, resid, ncv, ...
903 v, ldv, iparam, ipntr, workd, workl, lworkl, info );
904 if isequal(whch,'LM') | isequal(whch,'LA')
905 d = flipud(d);
906 if (rvec == 1)
907 v(:,1:k) = v(:,k:-1:1);
908 end
909 end
910 if ((isequal(whch,'SM') | isequal(whch,'SA')) & (rvec == 0))
911 d = flipud(d);
912 end
913 else
914 arpackc( 'dneupd', rvec, 'A', select, ...
915 d, di, v, ldv, sigmar, sigmai, workev, ...
916 bmat, int32(n), whch, nev, tol, resid, ncv, ...
917 v, ldv, iparam, ipntr, workd, workl, lworkl, info );
918 d = complex(d,di);
919 if rvec
920 d(k+1) = [];
921 else
922 zind = find(d == 0);
923 if isempty(zind)
924 d = d(k+1:-1:2);
925 else
926 d(max(zind)) = [];
927 d = flipud(d);
928 end
929 end
930 end
931 else
932 zsigma = [real(sigma); imag(sigma)];
933 arpackc( 'zneupd', rvec, 'A', select, ...
934 zd, zv, ldv, zsigma, workev, ...
935 bmat, int32(n), whch, nev, tol, resid, ncv, zv, ...
936 ldv, iparam, ipntr, zworkd, workl, lworkl, ...
937 rwork, info );
938 if issymA
939 d = zd(1:2:end-1);
940 else
941 d = complex(zd(1:2:end-1),zd(2:2:end));
942 end
943 v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]);
944 end
945
946 if (info ~= 0)
947 es = ['Error with ARPACK routine ' prefix 'eupd: '];
948 switch double(info)
949 case 2
950 ss = sum(select);
951 if (ss < k)
952 es = [es ...
953 ' The logical variable select was only set with ' int2str(ss) ...
954 ' 1''s instead of nconv=' int2str(double(iparam(5))) ...
955 ' (k=' int2str(k) ').' ...
956 ' Please contact the ARPACK authors at arpack@caam.rice.edu.'];
957 else
958 es = [es ...
959 'The LAPACK reordering routine ' prefix(1) ...
960 'trsen did not return all ' int2str(k) ' eigenvalues.']
961 end
962 case 1
963 es = [es ...
964 'The Schur form could not be reordered by the LAPACK routine ' ...
965 prefix(1) 'trsen.' ...
966 ' Please contact the ARPACK authors at arpack@caam.rice.edu.'];
967 case -14
968 es = [es prefix ...
969 'aupd did not find any eigenvalues to sufficient accuracy.'];
970 otherwise
971 es = [es sprintf('info = %d',full(info))];
972 end
973 error(es)
974 else
975 nconv = double(iparam(5));
976 if (nconv == 0)
977 if (nargout < 3)
978 warning('MATLAB:eigs:NoEigsConverged', ...
979 'None of the %d requested eigenvalues converged.',k)
980 else
981 flag = 1;
982 end
983 elseif (nconv < k)
984 if (nargout < 3)
985 warning('MATLAB:eigs:NotAllEigsConverged', ...
986 'Only %d of the %d requested eigenvalues converged.',nconv,k)
987 else
988 flag = 1;
989 end
990 end
991 end % if (info ~= 0)
992end % if (info < 0)
993
994if (issymA) | (~isrealprob)
995 if (nargout <= 1)
996 if isrealprob
997 varargout{1} = d;
998 else
999 varargout{1} = d(k:-1:1,1);
1000 end
1001 else
1002 varargout{1} = v(:,1:k);
1003 varargout{2} = diag(d(1:k,1));
1004 if (nargout >= 3)
1005 varargout{3} = flag;
1006 end
1007 end
1008else
1009 if (nargout <= 1)
1010 varargout{1} = d;
1011 else
1012 cplxd = find(di ~= 0);
1013 % complex conjugate pairs of eigenvalues occur together
1014 cplxd = cplxd(1:2:end);
1015 v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ...
1016 complex(v(:,cplxd),-v(:,cplxd+1))];
1017 varargout{1} = v(:,1:k);
1018 varargout{2} = diag(d);
1019 if (nargout >= 3)
1020 varargout{3} = flag;
1021 end
1022 end
1023end
1024
1025if (nargout >= 2) & (mode == 1) & ~isempty(B)
1026 if issparse(B)
1027 varargout{1}(permB,:) = RB \ varargout{1};
1028 else
1029 varargout{1} = RB \ varargout{1};
1030 end
1031end
1032
1033cputms(4) = cputime-t0; % end timing post-processing
1034
1035cputms(5) = sum(cputms(1:4)); % total time
1036
1037if (display >= 2) %tim
1038 if (mode == 1)
1039 innerstr = sprintf(['Compute A*X:' ...
1040 ' %f\n'],cputms(3));
1041 elseif (mode == 2)
1042 innerstr = sprintf(['Compute A*X and solve B*X=Y for X:' ...
1043 ' %f\n'],cputms(3));
1044 elseif (mode == 3)
1045 if isempty(B)
1046 innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ...
1047 ' %f\n'],cputms(3));
1048 else
1049 innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ...
1050 ' %f\n'],cputms(3));
1051 end
1052 end
1053 if ((mode == 3) & (Amatrix))
1054 if isempty(B)
1055 prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ...
1056 ' %f\n'],cputms(1));
1057 else
1058 prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ...
1059 ' %f\n'],cputms(1));
1060 end
1061 elseif ((mode == 2) & (~cholB))
1062 prepstr = sprintf(['Pre-processing, including chol(B):' ...
1063 ' %f\n'],cputms(1));
1064 else
1065 prepstr = sprintf(['Pre-processing:' ...
1066 ' %f\n'],cputms(1));
1067 end
1068 sstr = sprintf(['***********CPU Timing Results in seconds***********']);
1069 ds = sprintf(['\n' sstr '\n' ...
1070 prepstr ...
1071 'ARPACK''s %saupd: %f\n' ...
1072 innerstr ...
1073 'Post-processing with ARPACK''s %seupd: %f\n' ...
1074 '***************************************************\n' ...
1075 'Total: %f\n' ...
1076 sstr '\n'], ...
1077 prefix,cputms(2),prefix,cputms(4),cputms(5));
1078 disp(ds)
1079 if nombre_A_times_X > 0 %tim
1080 disp(sprintf('Number of iterations : %f\n',nombreIterations));
1081 disp(sprintf('Number of times A*X was computed : %f\n',nombre_A_times_X));
1082 disp(sprintf('Average time for A*X : %f\n',cputms(3)/nombre_A_times_X));
1083 end
1084end
generated by cgit v1.2.2 (git 2.25.0) at 2025-04-30 19:40:19 -0400