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