* Bug #12238 fixed - [d v] = eigs(A) was broken for sparse matrices.
[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: if opts.cholB is true, B must be upper triangular.\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             computevec = rvec;
727             if iparam(7) == 3 & sigmai then
728                 computevec = 1;
729             end
730             [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);
731             if(info_eupd <> 0)
732                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNEUPD", info_eupd));
733             else
734                 if iparam(7) == 3 & sigmai then
735                     res_d = complex(zeros(nev + 1,1));
736                     i = 1;
737                     while i <= nev
738                         if(~di(i))
739                             res_d(i) = complex(z(:,i)'*A*z(:,i), 0);
740                             i = i + 1;
741                         else
742                             real_part = z(:,i)' * A * z(:,i) + z(:,i+1)' * A * z(:,i+1);
743                             imag_part = z(:,i)' * A * z(:,i+1) - z(:,i+1)' * A * z(:,i)
744                             res_d(i) = complex(real_part, imag_part);
745                             res_d(i+1) = complex(real_part, -imag_part);
746                             i = i + 2;
747                         end
748                     end
749                 else
750                     res_d = complex(dr, di);
751                 end
752                 res_d = res_d(1:nev);
753                 if(rvec)
754                     index = find(di~=0);
755                     index = index(1:2:$);
756                     res_v = z;
757                     res_v(:,[index index+1]) = [complex(res_v(:,index),res_v(:,index+1)), complex(res_v(:,index),-res_v(:,index+1))];
758                     res_d = diag(res_d);
759                     res_v = res_v(:,1:nev);
760                 end
761             end
762         end
763     else
764         [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);
765         if(info_eupd <> 0)
766             error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "ZNEUPD", info_eupd));
767         else
768             d(nev+1) = []
769             res_d = d;
770             if(rvec)
771                 res_d = diag(d);
772                 res_v = z;
773             end
774         end
775     end
776     if(rvec & iparam(7) == 1 & matB)
777         if issparse(%_B)
778             res_v = umf_lusolve(Rprimefact, res_v);
779             if(~cholB)
780                 res_v = res_v(perm, :);
781             end
782         else
783             res_v = R \ res_v;
784         end
785     end
786 endfunction
787
788
789 function [res_d, res_v] = feigs(A_fun, nA, %_B, nev, which, maxiter, tol, ncv, cholB, resid, info, a_real, a_sym)
790     lhs = argn(1);
791     rvec = 0;
792     if(lhs > 1)
793         rvec = 1;
794     end
795
796     //**************************
797     //Second variable nA :
798     //**************************
799     if(size(nA,'*') <> 1 | ~isreal(nA) | typeof(nA) <> "constant")
800         error(msprintf(gettext("%s: Wrong type for input argument #%d: n must be a positive integer.\n"), "eigs", 2));
801     end
802
803     if(floor(nA) <> nA | nA <= 0)
804         error(msprintf(gettext("%s: Wrong type for input argument #%d: n must be a positive integer.\n"), "eigs", 2));
805     end
806
807     //*************************
808     //Third variable B :
809     //*************************
810     if((typeof(%_B) <> "constant") & (typeof(%_B) <> "sparse"))
811         error(msprintf(gettext("%s: Wrong type for input argument #%d: An empty matrix or full or sparse square matrix expected.\n"), "eigs", 3));
812     end
813     [mB, nB] = size(%_B);
814
815     //Check if B is a square matrix
816     if(mB * nB == 1 | mB <> nB)
817         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: B must have the same size as A.\n"), "eigs", 3));
818     end
819
820     //check if B is complex
821     Breal = isreal(%_B);
822     matB = mB * nB;
823
824     //*************************
825     //NEV :
826     //*************************
827     //verification du type de nev
828     //check if nev is complex?
829     if(typeof(nev) <> "constant") | (~isreal(nev)) | (size(nev,"*") <> 1)
830         error(msprintf(gettext("%s: Wrong type for input argument #%d: A scalar expected.\n"), "eigs", 3));
831     end
832
833     if(nev <> floor(nev) | (nev<=0))
834         error(msprintf(gettext("%s: Wrong type for input argument #%d: k must be a positive integer.\n"), "eigs", 4));
835     end
836
837     if(a_sym & a_real & Breal)
838         if(nev >= nA)
839             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));
840         end
841     else
842         if(nev >= nA - 1)
843             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));
844         end
845     end
846
847     //*************************
848     //SIGMA AND WHICH :
849     //*************************
850     //Check type
851     select type(which)
852     case 1 then
853         if(typeof(which) <> "constant" | which == [] | isnan(which))
854             error(msprintf(gettext("%s: Wrong type for input argument #%d: a scalar expected.\n"), "eigs", 5));
855         end
856         sigma = which;
857         which = 'LM';
858     case 10 then
859         [mWHICH, nWHICH] = size(which);
860         if(mWHICH * nWHICH <> 1)
861             error(msprintf(gettext("%s: Wrong type for input argument #%d: a string expected.\n"), "eigs", 5));
862         end
863         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)
864             if(a_real & Breal & a_sym)
865                 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"));
866             else
867                 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"));
868             end
869         end
870         if((~a_real | ~Breal | ~a_sym) & (strcmp(which,'LA') == 0 | strcmp(which,'SA') == 0 | strcmp(which,'BE') == 0))
871             error(msprintf(gettext("%s: Invalid sigma value for complex or non symmetric problem.\n"), "eigs"));
872         end
873         if(a_real & Breal & a_sym & (strcmp(which,'LR') == 0 | strcmp(which,'SR') == 0 | strcmp(which,'LI') == 0 | strcmp(which,'SI') == 0))
874             error(msprintf(gettext("%s: Invalid sigma value for real symmetric problem.\n"), "eigs"));
875         end
876         sigma = 0;
877     else
878         error(msprintf(gettext("%s: Wrong type for input argument #%d: a real scalar or a string expected.\n"), "eigs", 5));
879     end
880
881     if(~a_real | ~Breal)
882         sigma = complex(sigma);
883     end
884
885     //*************************
886     //MAXITER :
887     //*************************
888     if(typeof(maxiter) <> "constant" | ~isreal(maxiter) | size(maxiter,'*') <> 1)
889         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be a scalar.\n"), "eigs", 6, "opts.maxiter"));
890     end
891
892     if((maxiter <> floor(maxiter)) | (maxiter <= 0) )
893         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer positive value.\n"), "eigs", 6, "opts.maxiter"));
894     end
895
896     //*************************
897     //TOL :
898     //*************************
899     if(typeof(tol) <> "constant" | ~isreal(tol) | isnan(tol))
900         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be a real scalar.\n"), "eigs", 7, "opts.tol"));
901     end
902
903     if(size(tol,'*') <> 1)
904         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: %s must be 1 by 1 size.\n"), "eigs", 7, "opts.tol"));
905     end
906
907     //*************************
908     //NCV :
909     //*************************
910     if(typeof(ncv) <> "constant" | ~isreal(ncv))
911         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar.\n"), "eigs", 8, "opts.ncv"));
912     end
913
914     if(size(ncv,'*') > 1 | ncv <> floor(ncv) | (ncv <> [] & ncv <= 0))
915         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: %s must be an integer scalar.\n"), "eigs", 8, "opts.ncv"));
916     end
917
918     if(isempty(ncv))
919         if(~a_sym & a_real & Breal)
920             ncv = min(max(2*nev+1, 20), nA);
921         else
922             ncv = min(max(2*nev, 20), nA);
923         end
924     else
925         if(ncv <= nev | ncv > nA)
926             if(a_sym & a_real & Breal)
927                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For real symmetric problems, NCV must be k < NCV <= N.\n"), "eigs", 8));
928             elseif(~a_sym & a_real & Breal)
929                 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));
930
931             else
932                 error(msprintf(gettext("%s: Wrong value for input argument #%d: For complex problems, NCV must be k + 1 < NCV <= N.\n"), "eigs", 8));
933             end
934         end
935     end
936
937     //*************************
938     //CHOL :
939     //*************************
940     select typeof(cholB)
941     case "boolean"
942         if and(cholB <> [%f %t]) then
943             error(msprintf(gettext("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB","%f", "%t"));
944         end
945     case "constant"
946         //check if chol is complex?
947         if(~isreal(cholB) | size(cholB, "*") <> 1)
948             error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB"));
949         end
950
951         if(and(cholB <> [0 1]))
952             error(msprintf(gettext("%s: Wrong value for input argument #%d: %s must be %s or %s.\n"), "eigs", 8, "opts.cholB","%f", "%t"));
953         end
954     else
955         error(msprintf(gettext("%s: Wrong type for input argument #%d: %s must be an integer scalar or a boolean.\n"), "eigs", 8, "opts.cholB"));
956     end
957
958     //*************************
959     //RESID :
960     //*************************
961     if(typeof(resid) <> "constant")
962         error(msprintf(gettext("%s: Wrong type for input argument #%d: A real or complex matrix expected.\n"), "eigs", 10));
963     end
964
965     if(size(resid,'*') ~= nA)
966         error(msprintf(gettext("%s: Wrong dimension for input argument #%d: Start vector opts.resid must be N by 1.\n"), "eigs", 10));
967     end
968
969     if(a_real & Breal)
970         //resid complexe ?
971         if(~isreal(resid))
972             error(msprintf(gettext("%s: Wrong type for input argument #%d: Start vector opts.resid must be real for real problems.\n"), "eigs", 10));
973         end
974     else
975         if(isreal(resid))
976             error(msprintf(gettext("%s: Wrong type for input argument #%d: Start vector opts.resid must be complex for complex problems.\n"), "eigs", 10));
977         end
978     end
979
980     iparam = zeros(11,1);
981     iparam(1) = 1;
982     iparam(3) = maxiter;
983     iparam(7) = 1;
984
985     ipntr = zeros(14,1);
986
987     //MODE 1, 2, 3, 4, 5
988     if(~strcmp(which,'SM') | sigma <> 0)
989         iparam(7) = 3;
990         which = 'LM';
991     end
992
993     //bmat initialization
994     if(matB == 0 | iparam(7) == 1)
995         bmat = 'I';
996     else
997         bmat = 'G';
998     end
999
1000     if(cholB)
1001         if(~and(triu(%_B) == %_B))
1002             error(msprintf(gettext("%s: Wrong type for input argument #%d: B must be symmetric or hermitian, definite, semi positive.\n"), "eigs", 2));
1003         end
1004         R = %_B;
1005         Rprime = R';
1006     end
1007
1008     if(~cholB & matB <> 0 & iparam(7) == 1)
1009         if(~Breal)
1010             error(msprintf(gettext("%s: Impossible to use the Cholesky factorisation with complex sparse matrices.\n"), "eigs"));
1011         else
1012             if(issparse(%_B))
1013                 [R,P] = spchol(%_B);
1014             else
1015                 R = chol(%_B);
1016             end
1017             Rprime = R';
1018         end
1019     end
1020
1021     //Main
1022     howmny = 'A';
1023     ido = 0;
1024     info_aupd = 0;
1025     _select = zeros(ncv,1);
1026
1027     if(a_real)
1028         if(a_sym)
1029             lworkl = ncv * ncv + 8 * ncv;
1030             v = zeros(nA, ncv);
1031             workl = zeros(lworkl, 1);
1032             workd = zeros(3 * nA, 1);
1033             d = zeros(nev, 1); 
1034             z = zeros(nA, nev); 
1035         else
1036             lworkl = 3 * ncv * (ncv + 2);
1037             v = zeros(nA, ncv);
1038             workl = zeros(lworkl, 1);
1039             workd = zeros(3 * nA, 1);
1040             dr = zeros(nev+1, 1);
1041             di = zeros(nev+1, 1);
1042             z = zeros(nA, nev + 1);
1043             workev = zeros(3 * ncv, 1);
1044         end
1045     else
1046         lworkl = 3 * ncv * ncv + 5 * ncv;
1047         v = zeros(nA, ncv) + 0 * %i;
1048         workl = zeros(lworkl, 1) + 0 * %i;
1049         workd = zeros(3 * nA, 1) + 0 * %i;
1050         rwork = zeros(ncv, 1);
1051         d = zeros(nev + 1, 1) + 0 * %i;
1052         z = zeros(nA, nev) + 0 * %i;
1053         workev = zeros(2 * ncv, 1) + 0 * %i;
1054     end
1055     while(ido <> 99)
1056         if(a_real & Breal)
1057             if(a_sym)
1058                 [ido, resid, v, iparam, ipntr, workd, workl, info] = dsaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_aupd);
1059                 if(info_aupd <0)
1060                     error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DSAUPD", info_aupd));
1061                 end
1062             else
1063                 [ido, resid, v, iparam, ipntr, workd, workl, info] = dnaupd(ido, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_aupd);
1064                 if(info_aupd <0)
1065                     error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNAUPD", info_aupd));
1066                 end
1067             end
1068         else
1069             [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);
1070             if(info_aupd <0)
1071                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "ZNAUPD", info_aupd));
1072             end
1073         end
1074
1075         if(ido == -1 | ido == 1 | ido == 2)
1076             if(iparam(7) == 1)
1077                 if(matB == 0)
1078                     ierr = execstr('A_fun(workd(ipntr(1):ipntr(1)+nA-1))', 'errcatch');
1079                     if(ierr <> 0)
1080                         break;
1081                     end
1082                     workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(1):ipntr(1)+nA-1));
1083                 else
1084                     ierr = execstr('A_fun(inv(R) * workd(ipntr(1):ipntr(1)+nA-1))', 'errcatch');
1085                     if(ierr <> 0)
1086                         break;
1087                     end
1088                     workd(ipntr(2):ipntr(2)+nA-1) = inv(Rprime) * A_fun(inv(R) * workd(ipntr(1):ipntr(1)+nA-1));
1089                 end
1090             elseif(iparam(7) == 3)
1091                 if(matB == 0)
1092                     if(ido == 2)
1093                         workd(ipntr(2):ipntr(2)+nA-1) = workd(ipntr(1):ipntr(1)+nA-1);
1094                     else
1095                         ierr = execstr('A_fun(workd(ipntr(1):ipntr(1)+nA-1))', 'errcatch');
1096                         if(ierr <> 0)
1097                             break;
1098                         end
1099                         workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(1):ipntr(1)+nA-1));
1100                     end
1101                 else
1102                     if(ido == 2)
1103                         if(cholB)
1104                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime * R * workd(ipntr(1):ipntr(1)+nA-1);
1105                         else
1106                             workd(ipntr(2):ipntr(2)+nA-1) = %_B * workd(ipntr(1):ipntr(1)+nA-1);
1107                         end
1108                     elseif(ido == -1)
1109                         if(cholB)
1110                             workd(ipntr(2):ipntr(2)+nA-1) = Rprime * R * workd(ipntr(1):ipntr(1)+nA-1);
1111                         else
1112                             workd(ipntr(2):ipntr(2)+nA-1) = %_B * workd(ipntr(1):ipntr(1)+nA-1);
1113                         end
1114                         ierr = execstr('A_fun(workd(ipntr(2):ipntr(2)+nA-1))', 'errcatch');
1115                         if(ierr <> 0)
1116                             break;
1117                         end
1118                         workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(2):ipntr(2)+nA-1));
1119                     else
1120                         ierr = execstr('A_fun(workd(ipntr(3):ipntr(3)+nA-1))', 'errcatch');
1121                         if(ierr <> 0)
1122                             break;
1123                         end
1124                         workd(ipntr(2):ipntr(2)+nA-1) = A_fun(workd(ipntr(3):ipntr(3)+nA-1));
1125                     end
1126                 end
1127             else
1128                 if(a_real & Breal)
1129                     if(a_sym)
1130                         error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "DSAUPD"));
1131                     else
1132                         error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "DNAUPD"));
1133                     end
1134                 else
1135                     error(msprintf(gettext("%s: Error with %s: unknown mode returned.\n"), "eigs", "ZNAUPD"));
1136                 end
1137             end
1138         end
1139     end
1140
1141     if(ierr <> 0)
1142         if(ierr == 10)
1143             error(msprintf(gettext("%s: Wrong value for input argument #%d: n does not match rows number of matrix A.\n"), "eigs", 2));
1144         end
1145         error(msprintf(gettext("%s: Wrong type or value for input arguments.\n"), "eigs"));
1146     end
1147
1148     if(a_real & Breal)
1149         if(a_sym)
1150             [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);
1151             if(info <> 0)
1152                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DSEUPD", info));
1153             else
1154                 res_d = d;
1155                 if(rvec)
1156                     res_d = diag(res_d);
1157                     res_v = z;
1158                 end
1159             end
1160         else
1161             sigmar = real(sigma);
1162             sigmai = imag(sigma);
1163             computevec = rvec;
1164             if iparam(7) == 3 & sigmai then
1165                 computevec = 1;
1166             end
1167             [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);
1168             if(info <> 0)
1169                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNEUPD", info));
1170             else
1171                 if iparam(7) == 3 & sigmai then
1172                     res_d = complex(zeros(nev + 1,1));
1173                     i = 1;
1174                     while i <= nev
1175                         if(~di(i))
1176                             res_d(i) = complex(z(:,i)'*A*z(:,i), 0);
1177                             i = i + 1;
1178                         else
1179                             real_part = z(:,i)' * A * z(:,i) + z(:,i+1)' * A * z(:,i+1);
1180                             imag_part = z(:,i)' * A * z(:,i+1) - z(:,i+1)' * A * z(:,i)
1181                             res_d(i) = complex(real_part, imag_part);
1182                             res_d(i+1) = complex(real_part, -imag_part);
1183                             i = i + 2;
1184                         end
1185                     end
1186                 else
1187                     res_d = complex(dr,di);
1188                 end
1189                 res_d = res_d(1:nev);
1190                 if(rvec)
1191                     res_d = diag(res_d)
1192                     res_v = z;
1193                     c1 = 1:2:nev + 1;
1194                     c2 = 2:2:nev + 1;
1195                     if(modulo(nev,2) == 1)
1196                         c1($) = [];
1197                     end
1198                     res_v(:,[c1, c2]) = [res_v(:,c1) + res_v(:,c2) * %i res_v(:,c1) - res_v(:,c2) * %i];
1199                     res_v(:,$) = [];
1200                 end
1201             end
1202         end
1203     else
1204         [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);
1205         if(info <> 0)
1206             error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "ZNEUPD", info));
1207         else
1208             d(nev+1) = []
1209             res_d = d;
1210             if(rvec)
1211                 res_d = diag(d);
1212                 res_v = z;
1213             end
1214         end
1215     end
1216 endfunction