* Buf #12772 fixed- eigs() failed when trying to solve a sparse matrix eigen value...
[scilab.git] / scilab / modules / arnoldi / macros / eigs.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
3 //
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10 function [d, v] = eigs(varargin)
11     lhs = argn(1);
12     rhs = argn(2);
13     if(rhs == 0 | rhs > 6)
14         error(msprintf(gettext("%s : Wrong number of input arguments : %d to %d expected.\n"), "eigs", 1, 6));
15     end
16
17     if(rhs >= 1)
18         if((typeof(varargin(1)) <> "constant")  & typeof(varargin(1)) <> "function" & (typeof(varargin(1)) <> "sparse") | varargin(1) == [])
19             error(msprintf(gettext("%s: Wrong type for input argument #%d: A full or sparse square matrix or a function expected"), "eigs", 1));
20         end
21     end
22
23     if(rhs >= 1 & typeof(varargin(1)) <> "function")
24         if(isreal(varargin(1)))
25             resid = rand(size(varargin(1), "r"), 1);
26         else
27             resid = complex(rand(size(varargin(1), "r"), 1), rand(size(varargin(1), "r"), 1));
28         end
29         A = varargin(1);
30     end
31
32     if(rhs > 1 & typeof(varargin(1)) ==  "function")
33         if(size(varargin(2)) <> 1)
34             error(msprintf(gettext("%s: Wrong type for input argument #%d: A positive integer expected if the first input argument is a function."), "eigs",2));
35         end
36         a_real = 1;
37         a_sym = 0;
38         resid = rand(varargin(2),1);
39         info = 0;
40         Af = varargin(1);
41         Asize = varargin(2);
42     end
43
44     maxiter = 300;
45     tol = %eps;
46     ncv = [];
47     cholB = 0;
48     info = 0;
49     B = [];
50     sigma = "LM";
51     if(rhs == 1)
52         if(~issparse(varargin(1)))
53             info = int32(0);
54         end
55     else
56         if(~issparse(varargin(1)) & ~issparse(varargin(2)))
57             info = int32(0);
58         end
59     end
60
61     if(typeof(varargin(1)) <> "function")
62         select rhs
63         case 1
64             nev =  min(size(A, "r"), 6);
65         case 2
66             nev = min(size(A, "r"), 6);
67             B = varargin(2);
68         case 3
69             B = varargin(2);
70             nev = varargin(3);
71         case 4
72             B = varargin(2);
73             nev = varargin(3);
74             sigma = varargin(4);
75         case 5
76             B = varargin(2);
77             nev = varargin(3);
78             sigma = varargin(4);
79             opts = varargin(5);
80             if(~isstruct(opts))
81                 error(msprintf(gettext("%s: Wrong type for input argument #%d: A structure expected"), "eigs", 5));
82             end
83             if(size(intersect(fieldnames(opts), ["tol", "maxiter", "ncv", "resid", "cholB"]), "*") < size(fieldnames(opts),"*"))
84                 error(msprintf(gettext("%s: Wrong type for input argument: If A is a matrix, use opts with tol, maxiter, ncv, resid, cholB"), "eigs"));
85             end
86             if(isfield(opts, "tol"))
87                 tol = opts.tol;
88             end
89             if(isfield(opts, "maxiter"))
90                 maxiter = opts.maxiter;
91             end
92             if(isfield(opts, "ncv"))
93                 ncv = opts.ncv;
94             end
95             if(isfield(opts, "resid"))
96                 resid = opts.resid;
97                 if(issparse(varargin(1)) | issparse(varargin(2)))
98                     info = 1;
99                 else
100                     info = int32(1);
101                 end
102                 if(and(resid==0))
103                     if(issparse(varargin(1)) | issparse(varargin(2)))
104                         info = 0;
105                     else
106                         info = int32(0);
107                     end
108                 end
109             end
110             if(isfield(opts,"cholB"))
111                 cholB = opts.cholB;
112             end
113         end
114
115         select lhs
116         case 1
117             if(issparse(A) | issparse(B))
118                 d = speigs(A, B, nev, sigma, maxiter, tol, ncv, cholB, resid, info);
119             else
120                 d = %_eigs(A, B, nev, sigma, maxiter, tol, ncv, cholB, resid, info);
121             end
122         case 2
123             if(issparse(A) | issparse(B))
124                 [d, v] = speigs(A, B, nev, "LM", maxiter, tol, ncv, cholB, resid, info);
125             else
126                 [d, v] = %_eigs(A, B, nev, "LM", maxiter, tol, ncv, cholB, resid, info);
127             end
128         end
129     else
130         select rhs
131         case 2
132             nev = min(Asize, 6)
133
134         case 3
135             nev = min(Asize, 6);
136             B = varargin(3);
137
138         case 4
139             B = varagin(3);
140             nev = varargin(4);
141
142         case 5
143             B = varagin(3);
144             nev = varargin(4);
145             sigma = varargin(5);
146
147         case 6
148             B = varargin(3);
149             nev = varargin(4);
150             sigma = varargin(5);
151             opts = varargin(6);
152             if(~isstruct(opts)) then
153                 error(msprintf(gettext("%s: Wrong type for input argument #%d: A structure expected"), "eigs",5));
154             end
155             if(size(intersect(fieldnames(opts), ["tol", "maxiter", "ncv", "resid", "cholB", "issym", "isreal"]), "*") < size(fieldnames(opts),"*"))
156                 error(msprintf(gettext("%s: Wrong type for input argument: If A is a matrix, use opts with tol, maxiter, ncv, resid, cholB"), "eigs"));
157             end
158             if(isfield(opts,"tol"))
159                 tol = opts.tol;
160             end
161             if(isfield(opts,"maxiter"))
162                 maxiter = opts.maxiter;
163             end
164             if(isfield(opts, "ncv"))
165                 ncv = opts.ncv;
166             end
167             if(isfield(opts,"resid"))
168                 resid = opts.resid;
169                 info = 1;
170                 if(and(resid==0))
171                     info = 0;
172                 end
173             end
174             if(isfield(opts,"cholB"))
175                 cholB = opts.cholB;
176             end
177             if(isfield(opts,"issym"))
178                 a_sym = opts.issym;
179             end
180             if(isfield(opts,"isreal"))
181                 a_real = opts.isreal;
182                 if(~a_real & ~isfield(opts,"resid"))
183                     resid = complex(rand(Asize, 1), rand(Asize, 1));
184                 end
185             end
186         end
187         select lhs
188         case 1
189             d = feigs(Af, Asize, B, nev, sigma, maxiter, tol, ncv, cholB, resid, info, a_real, a_sym);
190         case 2
191             [d, v] = feigs(Af, Asize, B, nev, sigma, maxiter, tol, ncv, cholB, resid, info, a_real, a_sym);
192         end
193     end
194 endfunction
195
196 function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, resid, info)
197     lhs = argn(1);
198     rvec = 0;
199     if(lhs > 1)
200         rvec = 1;
201     end
202
203     //**************************
204     //First variable A :
205     //**************************
206     [mA, nA] = size(A);
207
208     //check if A is a square matrix
209     if(mA * nA < 2 | mA <> nA)
210         error(msprintf(gettext("%s: Wrong type for input argument #%d: A square matrix expected.\n"), "eigs", 1));
211     end
212
213     //check if A is complex
214     Areal = isreal(A);
215
216     //check if A is symetric
217     Asym = norm(A-A') == 0;
218
219     //*************************
220     //Second variable B :
221     //*************************
222     if((typeof(%_B) <> "constant") & (typeof(%_B) <> "sparse"))
223         error(msprintf(gettext("%s: Wrong type for input argument #%d: An empty matrix or full or sparse square matrix expected.\n"), "eigs", 2));
224     end
225     [mB, nB] = size(%_B);
226
227     //Check if B is a square matrix
228     if(mB * nB <> 0 & (mB <> mA | nB <> nA))
229         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: B must have the same size as A.\n"), "eigs", 2));
230     end
231
232     //check if B is complex
233     Breal = isreal(%_B);
234     matB = mB * nB;
235
236     //*************************
237     //NEV :
238     //*************************
239     //verification du type de nev
240     //check if nev is complex?
241     if(typeof(nev) <> "constant") | (~isreal(nev)) | (size(nev,"*") <> 1)
242         error(msprintf(gettext("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 3));
243     end
244
245     if(nev <> floor(nev) | (nev<=0))
246         error(msprintf(gettext("%s: Wrong type for input argument #%d: k must be a positive integer.\n"), "eigs", 3));
247     end
248
249     if(Asym & Areal & Breal)
250         if(nev >= nA)
251             error(msprintf(gettext("%s: Wrong value for input argument #%d: For real symmetric problems, k must be an integer in the range 1 to N - 1.\n"), "eigs", 3));
252         end
253     else
254         if(nev >= nA - 1)
255             error(msprintf(gettext("%s: Wrong value for input argument #%d: For real non symmetric or complex problems, k must be an integer in the range 1 to N - 2.\n"), "eigs", 3));
256         end
257     end
258
259     //*************************
260     //SIGMA AND WHICH :
261     //*************************
262     //Check type
263     select type(which)
264     case 1 then
265         if(typeof(which) <> "constant" | which == [] | isnan(which))
266             error(msprintf(gettext("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 4));
267         end
268         sigma = which;
269         which = "LM";
270
271     case 10 then
272         [mWHICH, nWHICH] = size(which);
273         if(mWHICH * nWHICH <> 1)
274             error(msprintf(gettext("%s: Wrong type for input argument #%d: A string expected.\n"), "eigs", 4));
275         end
276         if(strcmp(which,"LM") ~= 0 & strcmp(which,"SM") ~= 0  & strcmp(which,"LR") ~= 0 & strcmp(which,"SR") ~= 0 & strcmp(which,"LI") ~= 0 & strcmp(which,"SI") ~= 0 & strcmp(which,"LA") ~= 0 & strcmp(which,"SA") ~= 0 & strcmp(which,"BE") ~= 0)
277             if(Areal & Breal & Asym)
278                 error(msprintf(gettext("%s: Wrong value for input argument #%d: Unrecognized sigma value.\n Sigma must be one of ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n"), "eigs", 4, "LM", "SM", "LA", "SA", "BE"));
279             else
280                 error(msprintf(gettext("%s: Wrong value for input argument #%d: Unrecognized sigma value.\n Sigma must be one of ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n"), "eigs", 4, "LM", "SM", "LR", "SR", "LI", "SI"));
281             end
282         end
283         if((~Areal | ~Breal | ~Asym) & (strcmp(which,"LA") == 0 | strcmp(which,"SA") == 0 | strcmp(which,"BE") == 0))
284             error(msprintf(gettext("%s: Invalid sigma value for complex or non symmetric problem.\n"), "eigs"));
285         end
286         if(Areal & Breal & Asym & (strcmp(which,"LR") == 0 | strcmp(which,"SR") == 0 | strcmp(which,"LI") == 0 | strcmp(which,"SI") == 0))
287             error(msprintf(gettext("%s: Invalid sigma value for real symmetric problem.\n"), "eigs"));
288         end
289         sigma = 0;
290
291     else
292         error(msprintf(gettext("%s: Wrong type for input argument #%d: a real scalar or a string expected.\n"), "eigs", 4));
293     end
294
295     if(~Areal | ~Breal)
296         sigma = complex(sigma);
297     end
298
299     //*************************
300     //MAXITER :
301     //*************************
302     if(typeof(maxiter) <> "constant" | ~isreal(maxiter) | size(maxiter, "*") <> 1)
303         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be a scalar.\n"), "eigs", 5, "opts.maxiter"));
304     end
305
306     if((maxiter <> floor(maxiter)) | (maxiter <= 0) )
307         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer positive value.\n"), "eigs", 5, "opts.maxiter"));
308     end
309
310     //*************************
311     //TOL :
312     //*************************
313     if(typeof(tol) <> "constant" | ~isreal(tol) | isnan(tol))
314         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be a real scalar.\n"), "eigs", 6, "opts.tol"));
315     end
316
317     if(size(tol, "*") <> 1)
318         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: %s must be 1 by 1 size.\n"), "eigs", 6, "opts.tol"));
319     end
320
321     //*************************
322     //NCV :
323     //*************************
324     if(typeof(ncv) <> "constant" | ~isreal(ncv))
325         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 7, "opts.ncv"));
326     end
327
328     if(size(ncv, "*") > 1 | ncv <> floor(ncv) | (ncv <> [] & ncv <= 0))
329         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 7, "opts.ncv"));
330     end
331
332     if(isempty(ncv))
333         if(~Asym & Areal & Breal)
334             ncv = min(max(2*nev+1, 20), nA);
335         else
336             ncv = min(max(2*nev, 20), nA);
337         end
338     else
339         if(ncv <= nev | ncv > nA)
340             if(Asym & Areal & Breal)
341                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For real symmetric problems, NCV must be k < NCV <= N.\n"), "eigs", 7));
342             elseif(~Asym & Areal & Breal)
343                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For real non symmetric problems, NCV must be k + 2 < NCV < N.\n"), "eigs", 7));
344             else
345                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For complex problems, NCV must be k + 1 < NCV <= N.\n"), "eigs", 7));
346             end
347         end
348     end
349
350     //*************************
351     //CHOL :
352     //*************************
353     select typeof(cholB)
354     case "boolean"
355         if and(opts.cholB <> [%f %t]) then
356             error(msprintf(gettext("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB","%f", "%t"));
357         end
358     case "constant"
359         //check if chol is complex?
360         if(~isreal(cholB) | size(cholB, "*") <> 1)
361             error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB"));
362         end
363
364         if(and(cholB <> [0 1]))
365             error(msprintf(gettext("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB","%f", "%t"));
366         end
367     else
368         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB"));
369     end
370
371     //*************************
372     //RESID :
373     //*************************
374     if(typeof(resid) <> "constant")
375         error(msprintf(gettext("%s: Wrong type for input argument #%d: A real or complex matrix expected.\n"), "eigs", 9));
376     end
377
378     if(size(resid, "*") ~= nA)
379         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: Start vector %s must be N by 1.\n"), "eigs", 9, "opts.resid"));
380     end
381
382     if(Areal & Breal)
383         //resid complexe ?
384         if(~isreal(resid))
385             error(msprintf(gettext("%s: Wrong type for input argument #%d: Start vector %s must be real for real problems.\n"), "eigs", 9, "opts.resid"));
386         end
387     else
388         if(isreal(resid))
389             error(msprintf(gettext("%s: Wrong type for input argument #%d: Start vector %s must be complex for complex problems.\n"), "eigs", 9, "opts.resid"));
390         end
391     end
392
393     iparam = zeros(11,1);
394     iparam(1) = 1;
395     iparam(3) = maxiter;
396     iparam(7) = 1;
397
398     ipntr = zeros(14,1);
399
400     //MODE 1, 2, 3, 4, 5
401     if(~strcmp(which,"SM") | sigma <> 0)
402         iparam(7) = 3;
403         which = "LM";
404     end
405
406     //bmat initialization
407     if(matB == 0 | iparam(7) == 1)
408         bmat = "I";
409     else
410         bmat = "G";
411     end
412
413     if(cholB)
414         if(or(triu(%_B) <> %_B))
415             error(msprintf(gettext("%s: Wrong type for input argument #%d: if opts.cholB is true, B must be upper triangular.\n"), "eigs", 2));
416         end
417         if(issparse(%_B)) //sparse cholesky decomposition is reversed...
418             Rprime = %_B;
419             R = Rprime';
420         else
421             R = %_B;
422             Rprime = R';
423         end
424     end
425
426     if(~cholB & matB & iparam(7) == 1)
427         if(issparse(%_B) & ~Breal)
428             error(msprintf(gettext("%s: Impossible to use the Cholesky factorization with complex sparse matrices.\n"), "eigs"));
429         else
430             if(issparse(%_B))
431                 [R, P] = spchol(%_B);
432                 perm = spget(P);
433                 perm = perm(:,2);
434                 iperm = spget(P');
435                 iperm = iperm(:,2);
436             else
437                 R = chol(%_B);
438                 Rprime = R';
439             end
440         end
441     end
442     if(matB & issparse(%_B) & iparam(7) ==1)
443         Rfact = umf_lufact(R);
444         Rprimefact = umf_lufact(R');
445     end
446
447     //Main
448     howmny = "A";
449     ido = 0;
450     info_eupd = 0;
451     _select = zeros(ncv,1);
452     if(iparam(7) == 3)
453         if(matB == 0)
454             AMSB = A - sigma * speye(nA, nA);
455         else
456             if(cholB)
457                 AMSB = A - (sigma * Rprime * R);
458             else
459                 AMSB = A - sigma * %_B;
460             end
461         end
462         Lup = umf_lufact(AMSB);
463     end
464
465     if(Areal)
466         if(Asym)
467             lworkl = ncv * ncv + 8 * ncv;
468             v = zeros(nA, ncv);
469             workl = zeros(lworkl, 1);
470             workd = zeros(3 * nA, 1);
471             d = zeros(nev, 1);
472             z = zeros(nA, nev);
473         else
474             lworkl = 3 * ncv * (ncv + 2);
475             v = zeros(nA, ncv);
476             workl = zeros(lworkl, 1);
477             workd = zeros(3 * nA, 1);
478             dr = zeros(nev+1, 1);
479             di = zeros(nev+1, 1);
480             z = zeros(nA, nev + 1);
481             workev = zeros(3 * ncv, 1);
482         end
483     else
484         lworkl = 3 * ncv * ncv + 5 * ncv;
485         v = zeros(nA, ncv) + 0 * %i;
486         workl = zeros(lworkl, 1) + 0 * %i;
487         workd = zeros(3 * nA, 1) + 0 * %i;
488         rwork = zeros(ncv, 1);
489         d = zeros(nev + 1, 1) + 0 * %i;
490         z = zeros(nA, nev) + 0 * %i;
491         workev = zeros(2 * ncv, 1) + 0 * %i;
492     end
493
494     while(ido <> 99)
495         if(Areal & Breal)
496             if(Asym)
497                 [ido, resid, v, iparam, ipntr, workd, workl, info] = dsaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info);
498                 if(info < 0)
499                     error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DSAUPD", info));
500                 end
501             else
502                 [ido, resid, v, iparam, ipntr, workd, workl, info] = dnaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info);
503                 if(info < 0)
504                     error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNAUPD", info));
505                 end
506             end
507         else
508             [ido, resid, v, iparam, ipntr, workd, workl, rwork, info] = znaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, rwork, info);
509             if(info < 0)
510                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "ZNAUPD", info));
511             end
512         end
513
514         if(ido == -1 | ido == 1 | ido == 2)
515             if(iparam(7) == 1)
516                 if(ido == 2)
517                     workd(ipntr(2):ipntr(2)+nA-1) = workd(ipntr(1):ipntr(1)+nA-1);
518                 else
519                     if(matB == 0)
520                         workd(ipntr(2):ipntr(2)+nA-1) = A * workd(ipntr(1):ipntr(1)+nA-1);
521                     else
522                         if(issparse(%_B))
523                             y = umf_lusolve(Rprimefact, workd(ipntr(1):ipntr(1)+nA-1));
524                             if(~cholB)
525                                 y = A * y(perm);
526                                 y = y(iperm);
527                             else
528                                 y = A * y;
529                             end
530                             workd(ipntr(2):ipntr(2)+nA-1) = umf_lusolve(Rfact, y);
531                         else
532                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime \ (A * (R \ workd(ipntr(1):ipntr(1)+nA-1)));
533                         end
534                     end
535                 end
536             elseif(iparam(7) == 3)
537                 if(matB == 0)
538                     if(ido == 2)
539                         workd(ipntr(2):ipntr(2)+nA-1) = workd(ipntr(1):ipntr(1)+nA-1);
540                     else
541                         workd(ipntr(2):ipntr(2)+nA-1) = umf_lusolve(Lup, workd(ipntr(1):ipntr(1)+nA-1));
542                     end
543                 else
544                     if(ido == 2)
545                         if(cholB)
546                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime * (R * workd(ipntr(1):ipntr(1)+nA-1));
547                         else
548                             workd(ipntr(2):ipntr(2)+nA-1) = %_B * workd(ipntr(1):ipntr(1)+nA-1);
549                         end
550                     elseif(ido == -1)
551                         if(cholB)
552                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime * (R * workd(ipntr(1):ipntr(1)+nA-1));
553                         else
554                             workd(ipntr(2):ipntr(2)+nA-1) = %_B * workd(ipntr(1):ipntr(1)+nA-1);
555                         end
556                         workd(ipntr(2):ipntr(2)+nA-1) = umf_lusolve(Lup, workd(ipntr(2):ipntr(2)+nA-1));
557                     else
558                         workd(ipntr(2):ipntr(2)+nA-1) = umf_lusolve(Lup, workd(ipntr(3):ipntr(3)+nA-1));
559                     end
560                 end
561             else
562                 if(Areal & Breal)
563                     if(Asym)
564                         error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "DSAUPD"));
565                     else
566                         error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "DNAUPD"));
567                     end
568                 else
569                     error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "ZNAUPD"));
570                 end
571             end
572         end
573     end
574     if(iparam(7) == 3)
575         umf_ludel(Lup);
576     end
577     if(Areal & Breal)
578         if(Asym)
579             [d, z, resid, v, iparam, iptnr, workd, workl, info_eupd] = dseupd(rvec, howmny, _select, d, z, sigma, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_eupd);
580             if(info_eupd <> 0)
581                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DSEUPD", info_eupd));
582             else
583                 res_d = d;
584                 if(rvec)
585                     res_d = diag(res_d);
586                     res_v = z;
587                 end
588             end
589         else
590             sigmar = real(sigma);
591             sigmai = imag(sigma);
592             computevec = rvec;
593             if iparam(7) == 3 & sigmai then
594                 computevec = 1;
595             end
596             [dr, di, z, resid, v, iparam, ipntr, workd, workl, info_eupd] = dneupd(computevec, howmny, _select, dr, di, z, sigmar, sigmai, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_eupd);
597             if(info_eupd <> 0)
598                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNEUPD", info_eupd));
599             else
600                 if iparam(7) == 3 & sigmai then
601                     res_d = complex(zeros(nev + 1,1));
602                     i = 1;
603                     while i <= nev
604                         if(~di(i))
605                             res_d(i) = complex(z(:,i)'*A*z(:,i), 0);
606                             i = i + 1;
607                         else
608                             real_part = z(:,i)' * A * z(:,i) + z(:,i+1)' * A * z(:,i+1);
609                             imag_part = z(:,i)' * A * z(:,i+1) - z(:,i+1)' * A * z(:,i)
610                             res_d(i) = complex(real_part, imag_part);
611                             res_d(i+1) = complex(real_part, -imag_part);
612                             i = i + 2;
613                         end
614                     end
615                 else
616                     res_d = complex(dr, di);
617                 end
618                 res_d = res_d(1:nev);
619                 if(rvec)
620                     index = find(di~=0);
621                     index = index(1:2:$);
622                     res_v = z;
623                     if ~isempty(index) then
624                         res_v(:,[index index+1]) = [complex(res_v(:,index),res_v(:,index+1)), complex(res_v(:,index),-res_v(:,index+1))];
625                     end
626                     res_d = diag(res_d);
627                     res_v = res_v(:,1:nev);
628                 end
629             end
630         end
631     else
632         [d, z, resid, iparam, ipntr, workd, workl, rwork, info_eupd] = zneupd(rvec, howmny, _select, d, z, sigma, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, rwork, info_eupd);
633         if(info_eupd <> 0)
634             error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "ZNEUPD", info_eupd));
635         else
636             d(nev+1) = []
637             res_d = d;
638             if(rvec)
639                 res_d = diag(d);
640                 res_v = z;
641             end
642         end
643     end
644     if(rvec & iparam(7) == 1 & matB)
645         if issparse(%_B)
646             res_v = umf_lusolve(Rprimefact, res_v);
647             if(~cholB)
648                 res_v = res_v(perm, :);
649             end
650         else
651             res_v = R \ res_v;
652         end
653     end
654 endfunction
655
656
657 function [res_d, res_v] = feigs(A_fun, nA, %_B, nev, which, maxiter, tol, ncv, cholB, resid, info, a_real, a_sym)
658     lhs = argn(1);
659     rvec = 0;
660     if(lhs > 1)
661         rvec = 1;
662     end
663
664     //**************************
665     //Second variable nA :
666     //**************************
667     if(size(nA,"*") <> 1 | ~isreal(nA) | typeof(nA) <> "constant")
668         error(msprintf(gettext("%s: Wrong type for input argument #%d: n must be a positive integer.\n"), "eigs", 2));
669     end
670
671     if(floor(nA) <> nA | nA <= 0)
672         error(msprintf(gettext("%s: Wrong type for input argument #%d: n must be a positive integer.\n"), "eigs", 2));
673     end
674
675     //*************************
676     //Third variable B :
677     //*************************
678     if((typeof(%_B) <> "constant") & (typeof(%_B) <> "sparse"))
679         error(msprintf(gettext("%s: Wrong type for input argument #%d: An empty matrix or full or sparse square matrix expected.\n"), "eigs", 3));
680     end
681     [mB, nB] = size(%_B);
682
683     matB = mB * nB;
684     //Check if B is a square matrix
685     if(matB & (mB <> nA |nB <> nA))
686         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: B must have the same size as A.\n"), "eigs", 3));
687     end
688
689     //check if B is complex
690     Breal = isreal(%_B);
691
692     //*************************
693     //NEV :
694     //*************************
695     //verification du type de nev
696     //check if nev is complex?
697     if(typeof(nev) <> "constant") | (~isreal(nev)) | (size(nev,"*") <> 1)
698         error(msprintf(gettext("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 3));
699     end
700
701     if(nev <> floor(nev) | (nev<=0))
702         error(msprintf(gettext("%s: Wrong type for input argument #%d: k must be a positive integer.\n"), "eigs", 4));
703     end
704
705     if(a_sym & a_real & Breal)
706         if(nev >= nA)
707             error(msprintf(gettext("%s: Wrong value for input argument #%d: For real symmetric problems, k must be in the range 1 to N - 1.\n"), "eigs", 4));
708         end
709     else
710         if(nev >= nA - 1)
711             error(msprintf(gettext("%s: Wrong value for input argument #%d: For real non symmetric or complex problems, k must be in the range 1 to N - 2.\n"), "eigs", 4));
712         end
713     end
714
715     //*************************
716     //SIGMA AND WHICH :
717     //*************************
718     //Check type
719     select type(which)
720     case 1 then
721         if(typeof(which) <> "constant" | which == [] | isnan(which))
722             error(msprintf(gettext("%s: Wrong type for input argument #%d: a scalar expected.\n"), "eigs", 5));
723         end
724         sigma = which;
725         which = "LM";
726     case 10 then
727         [mWHICH, nWHICH] = size(which);
728         if(mWHICH * nWHICH <> 1)
729             error(msprintf(gettext("%s: Wrong type for input argument #%d: a string expected.\n"), "eigs", 5));
730         end
731         if(strcmp(which,"LM") ~= 0 & strcmp(which,"SM") ~= 0  & strcmp(which,"LR") ~= 0 & strcmp(which,"SR") ~= 0 & strcmp(which,"LI") ~= 0 & strcmp(which,"SI") ~= 0 & strcmp(which,"LA") ~= 0 & strcmp(which,"SA") ~= 0 & strcmp(which,"BE") ~= 0)
732             if(a_real & Breal & a_sym)
733                 error(msprintf(gettext("%s: Wrong value for input argument #%d: Unrecognized sigma value.\n Sigma must be one of %s, %s, %s, %s or %s.\n"), "eigs", 5, "LM", "SM", "LA", "SA", "BE"));
734             else
735                 error(msprintf(gettext("%s: Wrong value for input argument #%d: Unrecognized sigma value.\n Sigma must be one of %s, %s, %s, %s, %s or %s.\n"), "eigs", 5, "LM", "SM", "LR", "SR", "LI", "SI"));
736             end
737         end
738         if((~a_real | ~Breal | ~a_sym) & (strcmp(which,"LA") == 0 | strcmp(which,"SA") == 0 | strcmp(which,"BE") == 0))
739             error(msprintf(gettext("%s: Invalid sigma value for complex or non symmetric problem.\n"), "eigs"));
740         end
741         if(a_real & Breal & a_sym & (strcmp(which,"LR") == 0 | strcmp(which,"SR") == 0 | strcmp(which,"LI") == 0 | strcmp(which,"SI") == 0))
742             error(msprintf(gettext("%s: Invalid sigma value for real symmetric problem.\n"), "eigs"));
743         end
744         sigma = 0;
745     else
746         error(msprintf(gettext("%s: Wrong type for input argument #%d: a real scalar or a string expected.\n"), "eigs", 5));
747     end
748
749     if(~a_real | ~Breal)
750         sigma = complex(sigma);
751     end
752
753     //*************************
754     //MAXITER :
755     //*************************
756     if(typeof(maxiter) <> "constant" | ~isreal(maxiter) | size(maxiter,"*") <> 1)
757         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be a scalar.\n"), "eigs", 6, "opts.maxiter"));
758     end
759
760     if((maxiter <> floor(maxiter)) | (maxiter <= 0) )
761         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer positive value.\n"), "eigs", 6, "opts.maxiter"));
762     end
763
764     //*************************
765     //TOL :
766     //*************************
767     if(typeof(tol) <> "constant" | ~isreal(tol) | isnan(tol))
768         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be a real scalar.\n"), "eigs", 7, "opts.tol"));
769     end
770
771     if(size(tol,"*") <> 1)
772         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: %s must be 1 by 1 size.\n"), "eigs", 7, "opts.tol"));
773     end
774
775     //*************************
776     //NCV :
777     //*************************
778     if(typeof(ncv) <> "constant" | ~isreal(ncv))
779         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 8, "opts.ncv"));
780     end
781
782     if(size(ncv,"*") > 1 | ncv <> floor(ncv) | (ncv <> [] & ncv <= 0))
783         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: %s must be an integer scalar.\n"), "eigs", 8, "opts.ncv"));
784     end
785
786     if(isempty(ncv))
787         if(~a_sym & a_real & Breal)
788             ncv = min(max(2*nev+1, 20), nA);
789         else
790             ncv = min(max(2*nev, 20), nA);
791         end
792     else
793         if(ncv <= nev | ncv > nA)
794             if(a_sym & a_real & Breal)
795                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For real symmetric problems, NCV must be k < NCV <= N.\n"), "eigs", 8));
796             elseif(~a_sym & a_real & Breal)
797                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For real non symmetric problems, NCV must be k + 2 < NCV < N.\n"), "eigs", 8));
798
799             else
800                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For complex problems, NCV must be k + 1 < NCV <= N.\n"), "eigs", 8));
801             end
802         end
803     end
804
805     //*************************
806     //CHOL :
807     //*************************
808     select typeof(cholB)
809     case "boolean"
810         if and(cholB <> [%f %t]) then
811             error(msprintf(gettext("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB","%f", "%t"));
812         end
813     case "constant"
814         //check if chol is complex?
815         if(~isreal(cholB) | size(cholB, "*") <> 1)
816             error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB"));
817         end
818
819         if(and(cholB <> [0 1]))
820             error(msprintf(gettext("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB","%f", "%t"));
821         end
822     else
823         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB"));
824     end
825
826     //*************************
827     //RESID :
828     //*************************
829     if(typeof(resid) <> "constant")
830         error(msprintf(gettext("%s: Wrong type for input argument #%d: A real or complex matrix expected.\n"), "eigs", 10));
831     end
832
833     if(size(resid,"*") ~= nA)
834         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: Start vector opts.resid must be N by 1.\n"), "eigs", 10));
835     end
836
837     if(a_real & Breal)
838         //resid complexe ?
839         if(~isreal(resid))
840             error(msprintf(gettext("%s: Wrong type for input argument #%d: Start vector opts.resid must be real for real problems.\n"), "eigs", 10));
841         end
842     else
843         if(isreal(resid))
844             error(msprintf(gettext("%s: Wrong type for input argument #%d: Start vector opts.resid must be complex for complex problems.\n"), "eigs", 10));
845         end
846     end
847
848     iparam = zeros(11,1);
849     iparam(1) = 1;
850     iparam(3) = maxiter;
851     iparam(7) = 1;
852
853     ipntr = zeros(14,1);
854
855     //MODE 1, 2, 3, 4, 5
856     if(~strcmp(which,"SM") | sigma <> 0)
857         iparam(7) = 3;
858         which = "LM";
859     end
860
861     //bmat initialization
862     if(matB == 0 | iparam(7) == 1)
863         bmat = "I";
864     else
865         bmat = "G";
866     end
867
868     if(cholB)
869         if(or(triu(%_B) <> %_B))
870             error(msprintf(gettext("%s: Wrong type for input argument #%d: if opts.cholB is true, B must be upper triangular.\n"), "eigs", 2));
871         end
872         if(issparse(%_B)) //sparse cholesky decomposition is reversed...
873             Rprime = %_B;
874             R = Rprime;
875         else
876             R = %_B;
877             Rprime = R';
878         end
879     end
880     if(~cholB & matB & iparam(7) == 1)
881         if(issparse(%_B) & ~Breal)
882             error(msprintf(gettext("%s: Impossible to use the Cholesky factorization with complex sparse matrices.\n"), "eigs"));
883         else
884             if(issparse(%_B))
885                 [R,P] = spchol(%_B);
886                 perm = spget(P);
887                 perm = perm(:,2);
888                 iperm = spget(P');
889                 iperm = iperm(:,2);
890             else
891                 R = chol(%_B);
892                 Rprime = R';
893             end
894         end
895     end
896     if(matB & issparse(%_B) & iparam(7) == 1)
897         Rfact = umf_lufact(R);
898         Rprimefact = umf_lufact(R');
899     end
900
901     //Main
902     howmny = "A";
903     ido = 0;
904     info_aupd = 0;
905     _select = zeros(ncv,1);
906
907     if(a_real)
908         if(a_sym)
909             lworkl = ncv * ncv + 8 * ncv;
910             v = zeros(nA, ncv);
911             workl = zeros(lworkl, 1);
912             workd = zeros(3 * nA, 1);
913             d = zeros(nev, 1);
914             z = zeros(nA, nev);
915         else
916             lworkl = 3 * ncv * (ncv + 2);
917             v = zeros(nA, ncv);
918             workl = zeros(lworkl, 1);
919             workd = zeros(3 * nA, 1);
920             dr = zeros(nev+1, 1);
921             di = zeros(nev+1, 1);
922             z = zeros(nA, nev + 1);
923             workev = zeros(3 * ncv, 1);
924         end
925     else
926         lworkl = 3 * ncv * ncv + 5 * ncv;
927         v = zeros(nA, ncv) + 0 * %i;
928         workl = zeros(lworkl, 1) + 0 * %i;
929         workd = zeros(3 * nA, 1) + 0 * %i;
930         rwork = zeros(ncv, 1);
931         d = zeros(nev + 1, 1) + 0 * %i;
932         z = zeros(nA, nev) + 0 * %i;
933         workev = zeros(2 * ncv, 1) + 0 * %i;
934     end
935     while(ido <> 99)
936         if(a_real & Breal)
937             if(a_sym)
938                 [ido, resid, v, iparam, ipntr, workd, workl, info] = dsaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_aupd);
939                 if(info_aupd <0)
940                     error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DSAUPD", info_aupd));
941                 end
942             else
943                 [ido, resid, v, iparam, ipntr, workd, workl, info] = dnaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_aupd);
944                 if(info_aupd <0)
945                     error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNAUPD", info_aupd));
946                 end
947             end
948         else
949             [ido, resid, v, iparam, ipntr, workd, workl, rwork, info] = znaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, rwork, info_aupd);
950             if(info_aupd <0)
951                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "ZNAUPD", info_aupd));
952             end
953         end
954
955         if(ido == -1 | ido == 1 | ido == 2)
956             if(iparam(7) == 1)
957                 if(ido == 2)
958                     workd(ipntr(2):ipntr(2)+nA-1) = workd(ipntr(1):ipntr(1)+nA-1);
959                 else
960                     if(matB == 0)
961                         ierr = execstr("workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(1):ipntr(1)+nA-1))", "errcatch");
962                         if(ierr <> 0)
963                             break;
964                         end
965                     else
966                         if(issparse(%_B))
967                             y = umf_lusolve(Rprimefact, workd(ipntr(1):ipntr(1)+nA-1));
968                             if(~cholB)
969                                 ierr = execstr("workd(ipntr(2):ipntr(2)+nA-1) = A_fun( y(perm) )", "errcatch");
970                                 if(ierr <> 0)
971                                     break;
972                                 end
973                                 y = y(iperm);
974                             else
975                                 ierr = execstr("workd(ipntr(2):ipntr(2)+nA-1) = A_fun(y)", "errcatch");
976                                 if(ierr <> 0)
977                                     break;
978                                 end
979                             end
980                             workd(ipntr(2):ipntr(2)+nA-1) = umf_lusolve(Rfact, y);
981                         else
982                             ierr = execstr("workd(ipntr(2):ipntr(2)+nA-1) = A_fun( R \ workd(ipntr(1):ipntr(1)+nA-1) )", "errcatch");
983                             if(ierr <> 0)
984                                 break;
985                             end
986                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime \ workd(ipntr(2):ipntr(2)+nA-1);
987                         end
988                     end
989                 end
990             elseif(iparam(7) == 3)
991                 if(matB == 0)
992                     if(ido == 2)
993                         workd(ipntr(2):ipntr(2)+nA-1) = workd(ipntr(1):ipntr(1)+nA-1);
994                     else
995                         ierr = execstr("workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(1):ipntr(1)+nA-1))", "errcatch");
996                         if(ierr <> 0)
997                             break;
998                         end
999                     end
1000                 else
1001                     if(ido == 2)
1002                         if(cholB)
1003                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime * (R * workd(ipntr(1):ipntr(1)+nA-1));
1004                         else
1005                             workd(ipntr(2):ipntr(2)+nA-1) = %_B * workd(ipntr(1):ipntr(1)+nA-1);
1006                         end
1007                     elseif(ido == -1)
1008                         if(cholB)
1009                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime * (R * workd(ipntr(1):ipntr(1)+nA-1));
1010                         else
1011                             workd(ipntr(2):ipntr(2)+nA-1) = %_B * workd(ipntr(1):ipntr(1)+nA-1);
1012                         end
1013                         ierr = execstr("workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(2):ipntr(2)+nA-1))", "errcatch");
1014                         if(ierr <> 0)
1015                             break;
1016                         end
1017                     else
1018                         ierr = execstr("workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(3):ipntr(3)+nA-1))", "errcatch");
1019                         if(ierr <> 0)
1020                             break;
1021                         end
1022                     end
1023                 end
1024             else
1025                 if(a_real & Breal)
1026                     if(a_sym)
1027                         error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "DSAUPD"));
1028                     else
1029                         error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "DNAUPD"));
1030                     end
1031                 else
1032                     error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "ZNAUPD"));
1033                 end
1034             end
1035         end
1036     end
1037
1038     if(ierr <> 0)
1039         if(ierr == 10)
1040             error(msprintf(gettext("%s: Wrong value for input argument #%d: n does not match rows number of matrix A.\n"), "eigs", 2));
1041         end
1042         error(msprintf(gettext("%s: Wrong type or value for input arguments.\n"), "eigs"));
1043     end
1044
1045     if(a_real & Breal)
1046         if(a_sym)
1047             [d, z, resid, v, iparam, iptnr, workd, workl, info_eupd] = dseupd(rvec, howmny, _select, d, z, sigma, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info);
1048             if(info <> 0)
1049                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DSEUPD", info));
1050             else
1051                 res_d = d;
1052                 if(rvec)
1053                     res_d = diag(res_d);
1054                     res_v = z;
1055                 end
1056             end
1057         else
1058             sigmar = real(sigma);
1059             sigmai = imag(sigma);
1060             computevec = rvec;
1061             if iparam(7) == 3 & sigmai then
1062                 computevec = 1;
1063             end
1064             [dr, di, z, resid, v, iparam, ipntr, workd, workl, info_eupd] = dneupd(computevec, howmny, _select, dr, di, z, sigmar, sigmai, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info);
1065             if(info <> 0)
1066                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNEUPD", info));
1067             else
1068                 if iparam(7) == 3 & sigmai then
1069                     res_d = complex(zeros(nev + 1,1));
1070                     i = 1;
1071                     while i <= nev
1072                         if(~di(i))
1073                             res_d(i) = complex(z(:,i)'*A*z(:,i), 0);
1074                             i = i + 1;
1075                         else
1076                             real_part = z(:,i)' * A * z(:,i) + z(:,i+1)' * A * z(:,i+1);
1077                             imag_part = z(:,i)' * A * z(:,i+1) - z(:,i+1)' * A * z(:,i)
1078                             res_d(i) = complex(real_part, imag_part);
1079                             res_d(i+1) = complex(real_part, -imag_part);
1080                             i = i + 2;
1081                         end
1082                     end
1083                 else
1084                     res_d = complex(dr,di);
1085                 end
1086                 res_d = res_d(1:nev);
1087                 if(rvec)
1088                     index = find(di~=0);
1089                     index = index(1:2:$);
1090                     res_v = z;
1091                     res_v(:,[index index+1]) = [complex(res_v(:,index), res_v(:,index+1)), complex(res_v(:,index), -res_v(:,index+1))];
1092                     res_d = diag(res_d);
1093                     res_v = res_v(:,1:nev);
1094                 end
1095             end
1096         end
1097     else
1098         [d, z, resid, iparam, ipntr, workd, workl, rwork, info_eupd] = zneupd(rvec, howmny, _select, d, z, sigma, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, rwork, info);
1099         if(info <> 0)
1100             error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "ZNEUPD", info));
1101         else
1102             d(nev+1) = []
1103             res_d = d;
1104             if(rvec)
1105                 res_d = diag(d);
1106                 res_v = z;
1107             end
1108         end
1109     end
1110     if(rvec & iparam(7) == 1 & matB)
1111         if issparse(%_B)
1112             res_v = umf_lusolve(Rprimefact, res_v);
1113             if(~cholB)
1114                 res_v = res_v(perm, :);
1115             end
1116         else
1117             res_v = R \ res_v;
1118         end
1119     end
1120 endfunction