scilab/modules/arnoldi/macros/eigs.sci

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