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