octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Eigs again - Success!!


From: David Bateman
Subject: Eigs again - Success!!
Date: Sat, 04 Nov 2006 21:27:27 +0100
User-agent: Thunderbird 1.5.0.7 (X11/20060921)

Ok, I think I've figured out the seg fqault I was getting. Basically, I
believe my test example turned up a bug in ARPACK itself, and was not an
issue with gfortran/g95 as I'd previously thought..

The problem is that the Ritz vectors are calculated over a N-by-N
sub-matrix that is typically about twice as large as the number of
requested eigenvalues. When this sub-matrix needs to be reordered to
extract the correct eigenvalues and vectors, such as with the "sr" flag,
then the xTRSEN lapack function is used for the non-symmetric and
complex cases, though not for the symmetric real case (which is why the
problem doesn't occur there). xTRSEN is passed the variable NCONV for
its output variable M, where NCONV on input to xTRSEN is the number of
converged eigenvalues (upto the number of eigenvalues requested), and in
fact effectively is the same thing on output. Except for the case where
there are more converged eigenvalues than wanted. In that case NCONV is
increased and the following code that copies the converged eigenvalues
and vectors to the output variables will seg-fault. Enlarging these just
means that the segfault will be moved elsewhere.

The fix is a simple 5 line change to xNEUPD.f in ARPACK, to force the M
output of xTRSEN to not exceed the original NCONV value. I've reported
the error to the maintainers of ARPACK and will see if there is any
follow-up. In any case with the above change I can no longer generate
the segfault even after millions of iterations of the segtest code.. I
attach the ARPACK patch, last version of eigs and the test code for
those that are interested.

Basically, before finally submitting this function I'd like to split
eigs up into a number of classes, and use these template classes to
allow both sparse and full cases to be treated. So I hope to have
something relatively soon. So after 18 months of on again off again work
on eigs, I think I finally have something I can be confident in :-)

Cheers
David

/*

Copyright (C) 2005 David Bateman

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, write to the Free
Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.

*/

// XXX FIXME XXX. Comment out HAVE_CONFIG_H for build outside of octave!!
#if 1
#include <config.h>
#define HAVE_ARPACK
#else
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#endif

#include <cfloat>
#include <cmath>

#include "ov.h"
#include "defun-dld.h"
#include "error.h"
#include "f77-fcn.h"
#include "quit.h"
#include "variables.h"
#include "SparsedbleLU.h"
#include "SparseCmplxLU.h"
#include "ov-re-sparse.h"
#include "ov-cx-sparse.h"
#include "MatrixType.h"
#include "oct-map.h"
#include "pager.h"
#include "SparsedbleCHOL.h"
#include "SparseCmplxCHOL.h"
#include "oct-rand.h"

#ifdef HAVE_ARPACK
// Arpack fortran functions we call.
extern "C"
{
  F77_RET_T
  F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
                             const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, 
                             const octave_idx_type&, const double&,
                             double*, const octave_idx_type&, double*,
                             const octave_idx_type&, octave_idx_type*,
                             octave_idx_type*, double*, double*, 
                             const octave_idx_type&, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (dseupd, DSEUPD) (const int&, F77_CONST_CHAR_ARG_DECL,
                             int*, double*, double*,
                             const octave_idx_type&, const double&,
                             F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, 
                             F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, 
                             const double&, double*, const octave_idx_type&, 
                             double*, const octave_idx_type&, octave_idx_type*,
                             octave_idx_type*, double*, double*, 
                             const octave_idx_type&, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (dnaupd, DNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
                             const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, 
                             octave_idx_type&, const double&,
                             double*, const octave_idx_type&, double*,
                             const octave_idx_type&, octave_idx_type*,
                             octave_idx_type*, double*, double*, 
                             const octave_idx_type&, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (dneupd, DNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL,
                             int*, double*, double*,
                             double*, const octave_idx_type&, const double&,
                             const double&, double*, F77_CONST_CHAR_ARG_DECL, 
                             const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, 
                             octave_idx_type&, const double&, double*, 
                             const octave_idx_type&, double*, 
                             const octave_idx_type&, octave_idx_type*, 
                             octave_idx_type*, double*, double*, 
                             const octave_idx_type&, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (znaupd, ZNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
                             const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, 
                             const octave_idx_type&, const double&,
                             Complex*, const octave_idx_type&, Complex*,
                             const octave_idx_type&, octave_idx_type*,
                             octave_idx_type*, Complex*, Complex*, 
                             const octave_idx_type&, double *, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (zneupd, ZNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL,
                             int*, Complex*, Complex*, 
                             const octave_idx_type&, const Complex&,
                             Complex*, F77_CONST_CHAR_ARG_DECL,
                             const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, 
                             const octave_idx_type&, const double&,
                             Complex*, const octave_idx_type&, Complex*,
                             const octave_idx_type&, octave_idx_type*,
                             octave_idx_type*, Complex*, Complex*, 
                             const octave_idx_type&, double *, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);
}
#endif


#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
extern void
sparse_vector_product (const Sparse<double>&, const double*, double*);

extern void
sparse_vector_product (const Sparse<Complex>&, const Complex*, Complex*);

extern octave_idx_type
sparse_lu_solve (const SparseMatrix&, const SparseMatrix&, Matrix&);

extern octave_idx_type
sparse_lu_solve (const SparseComplexMatrix&, const SparseComplexMatrix&, 
                 ComplexMatrix&);

extern ComplexMatrix
sparse_ltsolve (const SparseComplexMatrix&, const ColumnVector&, 
                const ComplexMatrix&);

extern Matrix
sparse_ltsolve (const SparseMatrix&, const ColumnVector&, 
                const Matrix&,);

extern ComplexMatrix
sparse_utsolve (const SparseComplexMatrix&, const ColumnVector&, 
                const ComplexMatrix&);

extern Matrix
sparse_utsolve (const SparseMatrix&, const ColumnVector&, 
                const Matrix&);
#endif

template <class T>
void
sparse_vector_product (const Sparse<T>& m, const T* x, T* y)
{
  octave_idx_type nc = m.cols ();

  for (octave_idx_type j = 0; j < nc; j++)
    y[j] = 0.;

  for (octave_idx_type j = 0; j < nc; j++)
    for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
      y[m.ridx(i)] += m.data(i) * x[j];
}

template <class M, class SM>
octave_idx_type
sparse_lu_solve (const SM& L, const SM& U, M& m)
{
  octave_idx_type nc = L.cols ();
  octave_idx_type err = 0;
  double rcond;
  MatrixType ltyp (MatrixType::Lower);
  MatrixType utyp (MatrixType::Upper);

  m = L.solve (ltyp, m, err, rcond, 0);
  if (err)
    return err;

  m = U.solve (utyp, m, err, rcond, 0);
  if (err)
    return err;

  return err;
}

template <class SM, class M>
M
sparse_ltsolve (const SM& L, const ColumnVector& Q, const M& m)
{
  octave_idx_type n = L.cols();
  octave_idx_type b_nc = m.cols();
  octave_idx_type err = 0;
  double rcond;
  MatrixType ltyp (MatrixType::Lower);
  M tmp = L.solve (ltyp, m, err, rcond, 0);
  M retval;
  const double* qv = Q.fortran_vec();

  if (!err)
    {
      retval.resize (n, b_nc);
      for (octave_idx_type j = 0; j < b_nc; j++)
        {
          for (octave_idx_type i = 0; i < n; i++)
            retval.elem(static_cast<octave_idx_type>(qv[i]), j)  = 
              tmp.elem(i,j);
        }
    }

  return retval;
}

template <class SM, class M>
M
sparse_utsolve (const SM& U, const ColumnVector& Q, const M& m)
{
  octave_idx_type n = U.cols();
  octave_idx_type b_nc = m.cols();
  octave_idx_type err = 0;
  double rcond;
  MatrixType utyp (MatrixType::Upper);

  M retval (n, b_nc);
  const double* qv = Q.fortran_vec();
  for (octave_idx_type j = 0; j < b_nc; j++)
    {
      for (octave_idx_type i = 0; i < n; i++)
        retval.elem(i,j) = m.elem(static_cast<octave_idx_type>(qv[i]), j);
    }
  return U.solve (utyp, retval, err, rcond, 0);
}

DEFUN_DLD (eigs, args, nargout,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} address@hidden = eigs (@var{a})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{k}, 
@var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{k}, 
@var{sigma},@var{opts})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b}, 
@var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b}, 
@var{k}, @var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b}, 
@var{k}, @var{sigma}, @var{opts})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n}, 
@var{b})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n}, 
@var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n}, 
@var{b}, @var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n}, 
@var{k}, @var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n}, 
@var{b}, @var{k}, @var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n}, 
@var{k}, @var{sigma}, @var{opts})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n}, 
@var{b}, @var{k}, @var{sigma}, @var{opts})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}]} = eigs (@var{a}, 
@dots{})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}]} = eigs (@var{af}, 
@var{n}, @dots{})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}, @var{flag}]} = eigs 
(@var{a}, @dots{})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}, @var{flag}]} = eigs 
(@var{af}, @var{n}, @dots{})\n\
Calculate a limited number of eigenvalues and eigenvectors of @var{a},\n\
based on a selection criteria. The number eigenvalues and eigenvectors to\n\
calculate is given by @var{k} whose default value is 6.\n\
\n\
By default @code{eigs} solve the equation\n\
@iftex\n\
@tex\n\
$A \\nu = \\lambda \\nu$\n\
@end tex\n\
@end iftex\n\
@ifinfo\n\
@code{A * v = lambda * v}\n\
@end ifinfo\n\
, where\n\
@iftex\n\
@tex\n\
$\\lambda$ is a scalar representing one of the eigenvalues, and $\\nu$\n\
@end tex\n\
@end iftex\n\
@ifinfo\n\
@code{lambda} is a scalar representing one of the eigenvalues, and @code{v}\n\
@end ifinfo\n\
is the corresponding eigenvector. If given the positive definite matrix\n\
@var{B} then @code{eigs} solves the general eigenvalue equation\n\
@iftex\n\
@tex\n\
$A \\nu = \\lambda B \\nu$\n\
@end tex\n\
@end iftex\n\
@ifinfo\n\
@code{A * v = lambda * B * v}\n\
@end ifinfo\n\
.\n\
\n\
The argument @var{sigma} determines which eigenvalues are returned.\n\
@var{sigma} can be either a scalar or a string. When @var{sigma} is a scalar,\n\
the @var{k} eigenvalues closest to @var{sigma} are returned. If @var{sigma}\n\
is a string, it must have one of the values\n\
\n\
@table @asis\n\
@item 'lm'\n\
Largest magnitude (default).\n\
\n\
@item 'sm'\n\
Smallest magnitude.\n\
\n\
@item 'la'\n\
Largest Algebraic (valid only for real symmetric problem).\n\
\n\
@item 'sa'\n\
Smallest Algebraic (valid only for real symmetric problem).\n\
\n\
@item 'be'\n\
Both ends, with one more from the high-end if @var{k} is odd (valid only for\n\
real symmetric problem).\n\
\n\
@item 'lr'\n\
Largest real part (valid only for complex or unsymmetric problems).\n\
\n\
@item 'sr'\n\
Smallest real part (valid only for complex or unsymmetric problems).\n\
\n\
@item 'li'\n\
Largest imaginary part (valid only for complex or unsymmetric problems).\n\
\n\
@item 'si'\n\
Smallest imaginary part (valid only for complex or unsymmetric problems).\n\
@end table\n\
\n\
If @var{opts} is given, it is a structure defining some of the options that\n\
@code{eigs} should use. The fields of the structure @var{opts} are\n\
\n\
@table @code\n\
@item issym\n\
If @var{af} is given, then flags whether the function @var{af} defines a\n\
symmetric problem. It is ignored if @var{a} is given. The default is false.\n\
\n\
@item isreal\n\
If @var{af} is given, then flags whether the function @var{af} defines a\n\
real problem. It is ignored if @var{a} is given. The default is true.\n\
\n\
@item tol\n\
Defines the required convergence tolerance, given as @code{tol * norm (A)}.\n\
The default is @code{eps}.\n\
\n\
@item maxit\n\
The maximum number of iterations. The default is 300.\n\
\n\
@item p\n\
The number of Lanzcos basis vectors to use. More vectors will result in\n\
faster convergence, but a larger amount of memory. The optimal value of 'p'\n\
is problem dependent and should be in the range @var{k} to @var{n}. The\n\
default value is @code{2 * @var{k}}.\n\
\n\
@item v0\n\
The starting vector for the computation. The default is to have @sc{Arpack}\n\
randomly generate a starting vector.\n\
\n\
@item disp\n\
The level of diagnostic printout. If @code{disp} is 0 then there is no\n\
printout. The default value is 1.\n\
\n\
@item cholB\n\
Flag if @code{chol (@var{b})} is passed rather than @var{b}. The default is\n\
false.\n\
\n\
@item permB\n\
The permutation vector of the Cholesky factorization of @var{b} if\n\
@code{cholB} is true. That is @code{chol ( @var{b} (permB, permB))}. The\n\
default is @code{1:@var{n}}.\n\
\n\
@end table\n\
\n\
It is also possible to represent @var{a} by a function denoted @var{af}.\n\
@var{af} must be followed by a scalar argument @var{n} defining the length\n\
of the vector argument accepted by @var{af}. @var{af} can be passed either\n\
as an inline function, function handle or as a string. In the case where\n\
@var{af} is passed as a string, the name of the string defines the function\n\
to use.\n\
\n\
@var{af} is a function of the form @code{function y = af (x), y = @dots{};\n\
endfunction}, where the required return value of @var{af} is determined by\n\
the value of @var{sigma}, and are\n\
\n\
@table @code\n\
@item A * x\n\
If @var{sigma} is not given or is a string other than 'sm'.\n\
\n\
@item A \\ x\n\
If @var{sigma} is 'sm'.\n\
\n\
@item (A - sigma * I) \\ x\n\
for standard eigenvalue problem, where @code{I} is the identity matrix of\n\
the same size as @code{A}. If @var{sigma} is zero, this reduces the\n\
@code{A \\ x}.\n\
\n\
@item (A - sigma * B) \\ x\n\
for the general eigenvalue problem.\n\
@end table\n\
\n\
The return arguments of @code{eigs} depends on the number of return\n\
arguments. With a single return argument, a vector @var{d} of length @var{k}\n\
is returned, represent the @var{k} eigenvalues that have been found. With two\n\
return arguments, @var{v} is a @address@hidden matrix whose columns are\n\
the @var{k} eigenvectors corresponding to the returned eigenvalues. The\n\
eigenvalues themselves are then returned in @var{d} in the form of a\n\
@address@hidden matrix, where the elements on the diagonal are the\n\
eigenvalues.\n\
\n\
Given a third return argument @var{flag}, @code{eigs} also returns the status\n\
of the convergence. If @var{flag} is 0, then all eigenvalues have converged,\n\
otherwise not.\n\
\n\
This function is based on the @sc{Arpack} package, written by R Lehoucq,\n\
K Maschhoff, D Sorensen and C Yang. For more information see\n\
@url{http://www.caam.rice.edu/software/ARPACK/}.\n\
\n\
@end deftypefn\n\
@seealso{eig, svds}")
{
  octave_value_list retval;
#ifdef HAVE_ARPACK
  int nargin = args.length ();
  octave_function *func;
  std::string fcn_name;
  octave_idx_type n;
  octave_idx_type k = 6;
  Complex sigma;
  double sigmar, sigmai;
  bool have_sigma = false;
  std::string typ = "LM";
  SparseMatrix asmm, bsmm, bsmt;
  SparseComplexMatrix ascm, bscm, bsct;
  int b_arg;
  bool have_b = false;
  bool have_a_fun = false;
  bool a_is_complex = false;
  bool b_is_complex = false;
  bool symmetric = false;
  bool cholB = false;
  ColumnVector permB;
  int arg_offset = 0;
  double tol = DBL_EPSILON;
  int maxit = 300;
  int disp = 0;
  octave_idx_type p = -1;
  ColumnVector resid;
  ComplexColumnVector cresid;
  octave_idx_type info = 1;
  char bmat = 'I';
  octave_idx_type mode = 1;
  bool note3 = false;

  if (nargin == 0)
    print_usage ();
  else if (args(0).is_function_handle () || args(0).is_inline_function () ||
      args(0).is_string ())
    {
      if (args(0).is_string ())
        {
          std::string name = args(0).string_value ();
          std::string fname = "function y = ";
          fcn_name = unique_symbol_name ("__eigs_fcn_");
          fname.append (fcn_name);
          fname.append ("(x) y = ");
          func = extract_function (args(0), "eigs", fcn_name, fname,
                                       "; endfunction");
        }
      else
        func = args(0).function_value ();

      if (!func)
        error ("eigs: unknown function");

      if (nargin < 2)
        error ("eigs: incorrect number of arguments");
      else
        {
          n = args(1).nint_value ();
          arg_offset = 1;
          have_a_fun = true;
        }
    }
  else
    {
      MatrixType tmp;

      if (args(0).is_complex_type ())
        {
          ascm = (args(0).sparse_complex_matrix_value());
          a_is_complex = true;
          n = ascm.cols ();
          tmp = MatrixType (ascm);
        }
      else
        {
          asmm = (args(0).sparse_matrix_value());
          n = asmm.cols ();
          tmp = MatrixType (asmm);
        }

      if ((tmp.is_diagonal () && !a_is_complex) || tmp.is_hermitian ())
        symmetric = true;
    }


  // Note hold off reading B till later to avoid issues of double copies of 
  // the matrix if B is full/real while A is complex..
  if (!error_state && nargin > (1+arg_offset) && 
      !(args(1+arg_offset).is_real_scalar ()))
    if (args(1+arg_offset).is_complex_type ())
      {
        b_arg = 1+arg_offset;
        have_b = true;
        bmat = 'G';
        b_is_complex = true;
        arg_offset++;
      }
    else
      {
        b_arg = 1+arg_offset;
        have_b = true;
        bmat = 'G';
        arg_offset++;
      }

  if (!error_state && nargin > (1+arg_offset))
    k = args(1+arg_offset).nint_value ();

  if (!error_state && nargin > (2+arg_offset))
    if (args(2+arg_offset).is_real_scalar () ||
        args(2+arg_offset).is_complex_scalar ())
      {
        sigma = args(2+arg_offset).complex_value ();
        have_sigma = true;
      }
    else if (args(2+arg_offset).is_string ())
      {
        typ = args(2+arg_offset).string_value ();

        // Use STL function to convert to lower case
        transform (typ.begin (), typ.end (), typ.begin (), toupper);

        sigma = 0.;
      }
    else
      error ("eigs: sigma must be a scalar or a string");

  sigmar = std::real (sigma);
  sigmai = std::imag (sigma);

  if (!error_state && nargin > (3+arg_offset))
    if (args(3+arg_offset).is_map ())
      {
        Octave_map map = args(3+arg_offset).map_value ();

        // issym is ignored if A is not a function
        if (map.contains ("issym"))
          if (have_a_fun)
            symmetric = map.contents ("issym")(0).double_value () != 0.;

        // isreal is ignored if A is not a function
        if (map.contains ("isreal"))
          if (have_a_fun)
            a_is_complex = ! (map.contents ("isreal")(0).double_value () != 0.);

        if (map.contains ("tol"))
          tol = map.contents ("tol")(0).double_value ();

        if (map.contains ("maxit"))
          maxit = map.contents ("maxit")(0).nint_value ();

        if (map.contains ("p"))
          p = map.contents ("p")(0).nint_value ();

        if (map.contains ("v0"))
          {
            if (a_is_complex || b_is_complex)
              cresid = ComplexColumnVector 
                (map.contents ("v0")(0).complex_vector_value());
            else
              resid = ColumnVector (map.contents ("v0")(0).vector_value());
          }

        if (map.contains ("disp"))
          disp = map.contents ("disp")(0).nint_value ();

        if (map.contains ("cholB"))
          cholB = map.contents ("cholB")(0).double_value () != 0.;

        if (map.contains ("permB"))
          permB = ColumnVector (map.contents ("permB")(0).vector_value ()) 
            - 1.0;
      }
    else
      error ("eigs: options argument must be a structure");

  if (nargin > (4+arg_offset))
    error ("eigs: incorrect number of arguments");

  if (have_b)
    {
      if (a_is_complex || b_is_complex)
        bscm = args(b_arg).sparse_complex_matrix_value ();
      else
        bsmm = args(b_arg).sparse_matrix_value ();
    }

  if (!error_state)
    {
      if (a_is_complex || b_is_complex)
        {
          if (cresid.is_empty ())
            {
              std::string rand_dist = octave_rand::distribution();
              octave_rand::distribution("uniform");
              Array<double> rr (octave_rand::vector(n));
              Array<double> ri (octave_rand::vector(n));
              cresid = ComplexColumnVector (n);
              for (octave_idx_type i = 0; i < n; i++)
                cresid(i) = Complex(rr(i),ri(i));
              octave_rand::distribution(rand_dist);
            }
        }
      else if (resid.is_empty ())
        {
          std::string rand_dist = octave_rand::distribution();
          octave_rand::distribution("uniform");
          resid = ColumnVector (octave_rand::vector(n));
          octave_rand::distribution(rand_dist);
        }

      if (p < 0)
        {
          if (symmetric && !(a_is_complex || b_is_complex))
            p = k * 2;
          else
            p = k * 2 + 1;
          if (p < 20)
            p = 20;

          if (p > n - 1)
            p = n - 1 ;
        }
      else if (symmetric && (p <= k || p > n))
        error ("eigs: opts.p must be between k and n");
      else if (!symmetric && (p < k || p > n))
        error ("eigs: opts.p must be between k+1 and n");

      if (!a_is_complex && !b_is_complex && symmetric)
        {
          if (k > n )
            error ("eigs: Too many eigenvalues to extract (k >= n). Use 
'eig(full(A))' instead");
        }
      else if (k > n - 1)
        error ("eigs: Too many eigenvalues to extract (k >= n-1). Use 
'eig(full(A))' instead");
      
      if (! have_sigma)
        {
          if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && 
              typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
              typ != "SI")
            error ("eigs: unrecognized sigma value");

          if ((a_is_complex || b_is_complex || !symmetric) && 
              (typ == "LA" || typ == "SA" || typ == "BE"))
            error ("eigs: invalid sigma value for complex or unsymmetric 
problem");

          if (!a_is_complex && !b_is_complex && symmetric && 
              (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR"))
            error ("eigs: invalid sigma value for real symmetric problem");

          // Mode 1 for SM mode seems unstable for some reason.
          // Use Mode 3 instead, with sigma=0.
          if (typ == "SM")
            {
              typ = "LM";
              mode = 3;
            }

          if (have_b)
            {
              // See Note 3 dsaupd
              note3 = true;
              bmat = 'I';
              octave_idx_type info;
              if (a_is_complex || b_is_complex)
                {
                  SparseComplexCHOL fact (bscm, info, false);

                  if (fact.P() != 0)
                    error ("eigs: The matrix B is not positive definite");
                  else
                    {
                      bscm = fact.L();
                      bsct = bscm.transpose();
                      permB = fact.perm() - 1.0;

                    }
                }
              else
                {
                  SparseCHOL fact (bsmm, info, false);

                  if (fact.P() != 0)
                    error ("eigs: The matrix B is not positive definite");
                  else
                    {
                      bsmm = fact.L();
                      bsmt = bsmm.transpose();
                      permB = fact.perm() - 1.0;
                    }
                }
            }
        }
      else if (! std::abs (sigma))
        typ = "SM";
      else
        mode = 3;

      if (!error_state)
        {
          Array<octave_idx_type> ip (11);
          octave_idx_type *iparam = ip.fortran_vec ();

          ip(0) = 1; //ishift
          ip(1) = 0;   // ip(1) not referenced
          ip(2) = maxit; // mxiter, maximum number of iterations
          ip(3) = 1; // NB blocksize in recurrence
          ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
          ip(5) = 0; //ip(5) not referenced
          ip(6) = mode; // mode
          ip(7) = 0;
          ip(8) = 0;
          ip(9) = 0;
          ip(10) = 0;
          // ip(7) to ip(10) return values
 
          Array<octave_idx_type> iptr (14);
          octave_idx_type *ipntr = iptr.fortran_vec ();

          octave_idx_type ido = 0;
          int iter = 0;
      
          if (!a_is_complex && !b_is_complex)
            {
              SparseMatrix L, U;
              const octave_idx_type *P, *Q;

              // Need to explicitly declare this here, so that the permutation
              // vectors aren't deleted ahead of time.
              SparseLU fact;  

              if (! have_a_fun && mode == 3)
                {
                  // Caclulate LU decomposition of 'A - sigma * B'
                  SparseMatrix AminusSigmaB (asmm);

                  if (have_b)
                    {
                      if (cholB)
                        {
                          bsmt = bsmm.transpose ();
                          AminusSigmaB = AminusSigmaB - sigmar * bsmt * bsmm;
                        }
                      else
                        AminusSigmaB = AminusSigmaB - sigmar * bsmm;
                    }
                  else
                    {
                      SparseMatrix sigmat (n, n, n);

                      // Create sigma * speye(n,n)
                      sigmat.xcidx (0) = 0;
                      for (octave_idx_type i = 0; i < n; i++)
                        {
                          sigmat.xdata(i) = sigmar;
                          sigmat.xridx(i) = i;
                          sigmat.xcidx(i+1) = i + 1;
                        }

                      AminusSigmaB = AminusSigmaB - sigmat;
                    }

                  fact = SparseLU (AminusSigmaB);
                  
                  L = fact.L ();
                  U = fact.U ();
                  P = fact.row_perm ();
                  Q = fact.col_perm ();

                  // Test condition number of LU decomposition
                  double minU = octave_NaN;
                  double maxU = octave_NaN;
                  for (octave_idx_type j = 0; j < n; j++)
                    {
                      double d = 0.;
                      if (U.xcidx(j+1) > U.xcidx(j) &&
                          U.xridx (U.xcidx(j+1)-1) == j)
                        d = std::abs (U.xdata (U.xcidx(j+1)-1));

                      if (xisnan (minU) || d < minU)
                        minU = d;

                      if (xisnan (maxU) || d < maxU)
                        maxU = d;
                    }

                  if ((minU / maxU + 1.0) == 1.0)
                    {
                      warning ("eigs: 'A - sigma*B' is singular, indicating 
sigma is exactly an eigenvalue.");
                      warning ("      Convergence is not guaranteed");
                    }
                }

              if (symmetric)
                {
                  octave_idx_type lwork = p * (p + 8);

                  OCTAVE_LOCAL_BUFFER (double, v, n * p);
                  OCTAVE_LOCAL_BUFFER (double, workl, lwork);
                  OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
                  double *presid = resid.fortran_vec ();

                  do 
                    {
                      F77_FUNC (dsaupd, DSAUPD) 
                        (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                         F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
                         k, tol, presid, p, v, n, iparam,
                         ipntr, workd, workl, lwork, info
                         F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));

                      if (f77_exception_encountered)
                        {
                          error ("eigs: unrecoverable exception encountered in 
dsaupd");
                          goto eigs_err;
                        }

                      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
                        {
                          if (iter++)
                            {
                              octave_stdout << "Iteration " << iter - 1 << 
                                ": a few Ritz values of the " << p << "-by-" <<
                                p << " matrix\n";
                              for (int i = 0 ; i < k; i++)
                                octave_stdout << "    " <<
                                  workl[iptr(5)+i-1] << "\n";

                            }

                          // This is a kludge, as ARPACK doesn't give its
                          // iteration pointer. But as workl[iptr(5)-1] is
                          // an output value updated at each iteration, setting
                          // a value in this array to NaN and testing for it
                          // is a way of obtaining the iteration counter.
                          if (ido != 99)
                            workl[iptr(5)-1] = octave_NaN; 
                        }

                      if (ido == -1 || ido == 1 || ido == 2)
                        {
                          if (have_a_fun)
                            {
                              double *ip2 = workd + iptr(0) - 1;
                              ColumnVector x(n);

                              for (octave_idx_type i = 0; i < n; i++)
                                x(i) = *ip2++;

                              octave_value_list tmp = 
                                func->do_multi_index_op (1, octave_value (x));

                              ColumnVector y = tmp(0).column_vector_value ();

                              ip2 = workd + iptr(1) - 1;
                              for (octave_idx_type i = 0; i < n; i++)
                                *ip2++ = y(i);
                            }
                          else if (mode == 1)
                            {
                              if (have_b)
                                {
                                  Matrix mtmp (n,1);
                                  for (octave_idx_type i = 0; i < n; i++)
                                    mtmp(i,0) = workd[i + iptr(0) - 1];
                                  mtmp = sparse_utsolve(bsmt, permB, asmm * 
                                    sparse_ltsolve(bsmm, permB, mtmp));
                                  for (octave_idx_type i = 0; i < n; i++)
                                    workd[i+iptr(1)-1] = mtmp(i,0);
                                }
                              else
                                sparse_vector_product (asmm, 
                                                       workd + iptr(0) - 1, 
                                                       workd + iptr(1) - 1);
                            }
                          else if (mode == 2)
                            error ("eigs: impossible");
                          else if (mode == 3)
                            {
                              if (have_b)
                                {
                                  if (ido == -1)
                                    {
                                      if (cholB)
                                        {
                                          error ("fixme!!");
                                        }
                                      else
                                        {
                                          OCTAVE_LOCAL_BUFFER (double, dtmp, n);

                                          sparse_vector_product 
                                            (asmm, workd+iptr(0)-1, dtmp);

                                          Matrix tmp(n, 1);

                                          for (octave_idx_type i = 0; i < n; 
i++)
                                            tmp(i,0) = dtmp[P[i]];
                                  
                                          sparse_lu_solve (L, U, tmp);

                                          double *ip2 = workd+iptr(1)-1;
                                          for (octave_idx_type i = 0; i < n; 
i++)
                                            ip2[Q[i]] = tmp(i,0);
                                        }
                                    }
                                  else if (ido == 2)
                                    sparse_vector_product (bsmm, 
                                                           workd+iptr(0)-1, 
                                                           workd+iptr(1)-1);
                                  else
                                    {
                                      double *ip2 = workd+iptr(2)-1;
                                      Matrix tmp(n, 1);

                                      for (octave_idx_type i = 0; i < n; i++)
                                        tmp(i,0) = ip2[P[i]];
                                  
                                      sparse_lu_solve (L, U, tmp);

                                      ip2 = workd+iptr(1)-1;
                                      for (octave_idx_type i = 0; i < n; i++)
                                        ip2[Q[i]] = tmp(i,0);
                                    }
                                }
                              else
                                {
                                  if (ido == 2)
                                    {
                                      for (octave_idx_type i = 0; i < n; i++)
                                        workd[iptr(0) + i - 1] =
                                          workd[iptr(1) + i - 1];
                                    }
                                  else
                                    {
                                      double *ip2 = workd+iptr(0)-1;
                                      Matrix tmp(n, 1);

                                      for (octave_idx_type i = 0; i < n; i++)
                                        tmp(i,0) = ip2[P[i]];
                                  
                                      sparse_lu_solve (L, U, tmp);

                                      ip2 = workd+iptr(1)-1;
                                      for (octave_idx_type i = 0; i < n; i++)
                                        ip2[Q[i]] = tmp(i,0);
                                    }
                                }
                            }
                        }
                      else
                        {
                          if (info < 0)
                            error ("eigs: error %d in dsaupd", info);
                          break;
                        }
                    } 
                  while (1);
                  
                  if (!error_state)
                    {
                      octave_idx_type info2;

                      int rvec (nargout > 1);

                      // We have a problem in that the size of the C++ bool 
                      // type relative to the fortran logical type. It appears 
                      // that fortran uses 4-bytes per logical and C++ 1-byte 
                      // per bool, though this might be system dependent. As 
                      // long as the HOWMNY arg is not "S", the logical array
                      // is just workspace for ARPACK, so use int type to 
                      // avoid problems.
                      Array<int> s (p);
                      int *sel = s.fortran_vec ();
                        
                      Matrix eig_vec (n, k);
                      double *z = eig_vec.fortran_vec ();

                      ColumnVector eig_val (k);
                      double *d = eig_val.fortran_vec ();

                      F77_FUNC (dseupd, DSEUPD) 
                        (rvec, F77_CONST_CHAR_ARG2 ("A", 1),
                         sel, d, z, n, sigmar, 
                         F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                         F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
                         k, tol, presid, p, v, n, iparam,
                         ipntr, workd, workl, lwork, info2
                         F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) 
                         F77_CHAR_ARG_LEN(2));

                      if (f77_exception_encountered)
                        {
                          error ("eigs: unrecoverable exception encountered in 
dseupd");
                          goto eigs_err;
                        }


                      if (info2 == 0)
                        {
                          octave_idx_type k2 = k / 2;
                          if (typ != "SM" && typ != "BE")
                            {
                              for (octave_idx_type i = 0; i < k2; i++)
                                {
                                  double dtmp = d[i];
                                  d[i] = d[k - i - 1];
                                  d[k - i - 1] = dtmp;
                                }
                            }

                          if (rvec)
                            {
                              if (typ != "SM" && typ != "BE")
                                {
                                  OCTAVE_LOCAL_BUFFER (double, dtmp, n);

                                  for (octave_idx_type i = 0; i < k2; i++)
                                    {
                                      octave_idx_type off1 = i * n;
                                      octave_idx_type off2 = (k - i - 1) * n;

                                      if (off1 == off2)
                                        continue;

                                      for (octave_idx_type j = 0; j < n; j++)
                                        dtmp[j] = z[off1 + j];

                                      for (octave_idx_type j = 0; j < n; j++)
                                        z[off1 + j] = z[off2 + j];

                                      for (octave_idx_type j = 0; j < n; j++)
                                        z[off2 + j] = dtmp[j];
                                    }
                                }
                              retval(2) = double (info);
                              retval(1) = DiagMatrix (eig_val);

                              if (note3)
                                eig_vec = sparse_ltsolve(bsmm, permB, eig_vec);
                              retval(0) = eig_vec;
                            }
                          else
                            retval (0) = eig_val;
                        }
                      else
                        error ("eigs: error %d in dseupd", info2);
                    }
                }
              else
                {
                  octave_idx_type lwork = 3 * p * (p + 2);

                  OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
                  OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
                  OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);

                  double *presid = resid.fortran_vec ();

                  do 
                    {
                      F77_FUNC (dnaupd, DNAUPD) 
                        (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                         F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
                         k, tol, presid, p, v, n, iparam,
                         ipntr, workd, workl, lwork, info
                         F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));

                      if (f77_exception_encountered)
                        {
                          error ("eigs: unrecoverable exception encountered in 
dnaupd");
                          goto eigs_err;
                        }

                      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
                        {
                          if (iter++)
                            {
                              octave_stdout << "Iteration " << iter - 1 << 
                                ": a few Ritz values of the " << p << "-by-" <<
                                p << " matrix\n";
                              for (int i = 0 ; i < k; i++)
                                octave_stdout << "    " <<
                                  workl[iptr(5)+i-1] << "\n";
                            }
                          
                          // This is a kludge, as ARPACK doesn't give its
                          // iteration pointer. But as workl[iptr(5)-1] is
                          // an output value updated at each iteration, setting
                          // a value in this array to NaN and testing for it
                          // is a way of obtaining the iteration counter.
                          if (ido != 99)
                            workl[iptr(5)-1] = octave_NaN; 
                        }

                      if (ido == -1 || ido == 1 || ido == 2)
                        {
                          if (have_a_fun)
                            {
                              double *ip2 = workd + iptr(0) - 1;
                              ColumnVector x(n);

                              for (octave_idx_type i = 0; i < n; i++)
                                x(i) = *ip2++;

                              octave_value_list tmp = 
                                func->do_multi_index_op (1, octave_value (x));

                              ColumnVector y = tmp(0).column_vector_value ();

                              ip2 = workd + iptr(1) - 1;
                              for (octave_idx_type i = 0; i < n; i++)
                                *ip2++ = y(i);
                            }
                          else if (mode == 1)
                            {
                              if (have_b)
                                {
                                  Matrix mtmp (n,1);
                                  for (octave_idx_type i = 0; i < n; i++)
                                    mtmp(i,0) = workd[i + iptr(0) - 1];
                                  mtmp = sparse_utsolve(bsmt, permB, asmm * 
                                    sparse_ltsolve(bsmm, permB, mtmp));
                                  for (octave_idx_type i = 0; i < n; i++)
                                    workd[i+iptr(1)-1] = mtmp(i,0);
                                }
                              else
                                sparse_vector_product (asmm, 
                                                       workd + iptr(0) - 1, 
                                                       workd + iptr(1) - 1);
                            }
                          else if (mode == 2)
                            {
                              error ("eigs: impossible");
                            }
                          else if (mode == 3)
                            {
                              if (have_b)
                                {
                                  if (ido == -1)
                                    {
                                      if (cholB)
                                        {
                                          error ("fixme!!");
                                        }
                                      else
                                        {
                                          OCTAVE_LOCAL_BUFFER (double, dtmp, n);

                                          sparse_vector_product 
                                            (asmm, workd+iptr(0)-1, dtmp);

                                          Matrix tmp(n, 1);

                                          for (octave_idx_type i = 0; i < n; 
i++)
                                            tmp(i,0) = dtmp[P[i]];
                                  
                                          sparse_lu_solve (L, U, tmp);

                                          double *ip2 = workd+iptr(1)-1;
                                          for (octave_idx_type i = 0; i < n; 
i++)
                                            ip2[Q[i]] = tmp(i,0);
                                        }
                                    }
                                  else if (ido == 2)
                                    sparse_vector_product (bsmm, 
                                                           workd+iptr(0)-1, 
                                                           workd+iptr(1)-1);
                                  else
                                    {
                                      double *ip2 = workd+iptr(2)-1;
                                      Matrix tmp(n, 1);

                                      for (octave_idx_type i = 0; i < n; i++)
                                        tmp(i,0) = ip2[P[i]];
                                  
                                      sparse_lu_solve (L, U, tmp);

                                      ip2 = workd+iptr(1)-1;
                                      for (octave_idx_type i = 0; i < n; i++)
                                        ip2[Q[i]] = tmp(i,0);
                                    }
                                }
                              else
                                {
                                  if (ido == 2)
                                    {
                                      for (octave_idx_type i = 0; i < n; i++)
                                        workd[iptr(0) + i - 1] =
                                          workd[iptr(1) + i - 1];
                                    }
                                  else
                                    {
                                      double *ip2 = workd+iptr(0)-1;
                                      Matrix tmp(n, 1);

                                      for (octave_idx_type i = 0; i < n; i++)
                                        tmp(i,0) = ip2[P[i]];
                                  
                                      sparse_lu_solve (L, U, tmp);

                                      ip2 = workd+iptr(1)-1;
                                      for (octave_idx_type i = 0; i < n; i++)
                                        ip2[Q[i]] = tmp(i,0);
                                    }
                                }
                            }
                        }
                      else
                        {
                          if (info < 0)
                            error ("eigs: error %d in dnaupd", info);
                          break;
                        }
                    } 
                  while (1);

                  if (!error_state)
                    {
                      octave_idx_type info2;
                      int rvec (nargout > 1);

                      // We have a problem in that the size of the C++ bool 
                      // type relative to the fortran logical type. It appears 
                      // that fortran uses 4-bytes per logical and C++ 1-byte 
                      // per bool, though this might be system dependent. As 
                      // long as the HOWMNY arg is not "S", the logical array
                      // is just workspace for ARPACK, so use int to avoid
                      // problems!!
                      Array<int> s (p);
                      int *sel = s.fortran_vec ();

                      Matrix eig_vec (n, k + 1);
                      double *z = eig_vec.fortran_vec ();

                      OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
                      OCTAVE_LOCAL_BUFFER (double, di, k + 1);
                      OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);

                      F77_FUNC (dneupd, DNEUPD) 
                        (rvec, F77_CONST_CHAR_ARG2 ("A", 1),
                         sel, dr, di, z, n, sigmar, sigmai, workev, 
                         F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                         F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
                         k, tol, presid, p, v, n, iparam,
                         ipntr, workd, workl, lwork, info2
                         F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) 
                         F77_CHAR_ARG_LEN(2));

                      if (f77_exception_encountered)
                        {
                          error ("eigs: unrecoverable exception encountered in 
dneupd");
                          goto eigs_err;
                        }

                      ComplexColumnVector eig_val (k+1);
                      Complex *d = eig_val.fortran_vec ();

                      if (info2 == 0)
                        {
                          octave_idx_type j = 0;
                          for (octave_idx_type i = 0; i < k+1; i++)
                            {
                              if (dr[i] == 0.0 && di[i] == 0.0 && j == 0)
                                j++;
                              else
                                d [i-j] = Complex (dr[i], di[i]);
                            }
                          if (j == 0 && !rvec)
                            for (octave_idx_type i = 0; i < k; i++)
                              d[i] = d[i+1];

                          octave_idx_type k2 = k / 2;
                          for (octave_idx_type i = 0; i < k2; i++)
                            {
                              Complex dtmp = d[i];
                              d[i] = d[k - i - 1];
                              d[k - i - 1] = dtmp;
                            }
                          eig_val.resize(k);

                          if (nargout < 2)
                            retval (0) = eig_val;
                          else
                            {
                              OCTAVE_LOCAL_BUFFER (double, dtmp, n);

                              for (octave_idx_type i = 0; i < k2; i++)
                                {
                                  octave_idx_type off1 = i * n;
                                  octave_idx_type off2 = (k - i - 1) * n;

                                  if (off1 == off2)
                                    continue;

                                  for (octave_idx_type j = 0; j < n; j++)
                                    dtmp[j] = z[off1 + j];

                                  for (octave_idx_type j = 0; j < n; j++)
                                    z[off1 + j] = z[off2 + j];

                                  for (octave_idx_type j = 0; j < n; j++)
                                    z[off2 + j] = dtmp[j];
                                }

                              ComplexMatrix eig_vec2 (n, k);
                              octave_idx_type i = 0;
                              while (i < k)
                                {
                                  octave_idx_type off1 = i * n;
                                  octave_idx_type off2 = (i+1) * n;
                                  if (std::imag(eig_val(i)) == 0)
                                    {
                                      for (octave_idx_type j = 0; j < n; j++)
                                        eig_vec2(j,i) = 
                                          Complex(z[j+off1],0.);
                                      i++;
                                    }
                                  else
                                    {
                                      for (octave_idx_type j = 0; j < n; j++)
                                        {
                                          eig_vec2(j,i) = 
                                            Complex(z[j+off1],z[j+off2]);
                                          if (i < k - 1)
                                            eig_vec2(j,i+1) = 
                                              Complex(z[j+off1],-z[j+off2]);
                                        }
                                      i+=2;
                                    }
                                }

                              retval(2) = double (info);
                              retval(1) = ComplexDiagMatrix(eig_val);

                              if (note3)
                                eig_vec2 = 
                                  sparse_ltsolve(SparseComplexMatrix(bsmm),
                                                 permB, eig_vec2);
                              retval(0) = eig_vec2;
                            }
                        }
                      else
                        error ("eigs: error %d in dneupd", info2);

                    }
                }
            }
          else
            {
              SparseComplexMatrix L, U;
              const octave_idx_type *P, *Q;

              // Need to explicitly declare this here, so that the permutation
              // vectors aren't deleted ahead of time.
              SparseComplexLU fact;  

              if (! have_a_fun && mode == 3)
                {
                  // Caclulate LU decomposition of 'A - sigma * B'
                  SparseComplexMatrix AminusSigmaB (ascm);

                  if (have_b)
                    {
                      if (cholB)
                        {
                          bsct = bscm.transpose ();
                          AminusSigmaB = AminusSigmaB - sigma * bsct * bscm;
                        }
                      else
                        AminusSigmaB = AminusSigmaB - sigma * bscm;
                    }
                  else
                    {
                      SparseComplexMatrix sigmat (n, n, n);

                      // Create sigma * speye(n,n)
                      sigmat.xcidx (0) = 0;
                      for (octave_idx_type i = 0; i < n; i++)
                        {
                          sigmat.xdata(i) = sigma;
                          sigmat.xridx(i) = i;
                          sigmat.xcidx(i+1) = i + 1;
                        }

                      AminusSigmaB = AminusSigmaB - sigmat;
                    }

                  fact = SparseComplexLU (AminusSigmaB);
                  
                  L = fact.L ();
                  U = fact.U ();
                  P = fact.row_perm ();
                  Q = fact.col_perm ();

                  // Test condition number of LU decomposition
                  double minU = octave_NaN;
                  double maxU = octave_NaN;
                  for (octave_idx_type j = 0; j < n; j++)
                    {
                      double d = 0.;
                      if (U.xcidx(j+1) > U.xcidx(j) &&
                          U.xridx (U.xcidx(j+1)-1) == j)
                        d = std::abs (U.xdata (U.xcidx(j+1)-1));

                      if (xisnan (minU) || d < minU)
                        minU = d;

                      if (xisnan (maxU) || d < maxU)
                        maxU = d;
                    }

                  if ((minU / maxU + 1.0) == 1.0)
                    {
                      warning ("eigs: 'A - sigma*B' is singular, indicating 
sigma is exactly an eigenvalue.");
                      warning ("      Convergence is not guaranteed");
                    }
                }

              octave_idx_type lwork = p * (3 * p + 5);
              
              OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
              OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
              OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
              OCTAVE_LOCAL_BUFFER (double, rwork, p);
              Complex *presid = cresid.fortran_vec ();

              do 
                {
                  F77_FUNC (znaupd, ZNAUPD) 
                    (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                     F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
                     k, tol, presid, p, v, n, iparam,
                     ipntr, workd, workl, lwork, rwork, info
                     F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));

                  if (f77_exception_encountered)
                    {
                      error ("eigs: unrecoverable exception encountered in 
znaupd");
                      goto eigs_err;
                    }

                  if (disp > 0 && !xisnan(workl[iptr(5)-1]))
                    {
                      if (iter++)
                        {
                          octave_stdout << "Iteration " << iter - 1 << 
                            ": a few Ritz values of the " << p << "-by-" <<
                            p << " matrix\n";
                          for (int i = 0 ; i < k; i++)
                            octave_stdout << "    " << 
                              workl[iptr(5)+i-1] << "\n";
                        }
                          
                      // This is a kludge, as ARPACK doesn't give its
                      // iteration pointer. But as workl[iptr(5)-1] is
                      // an output value updated at each iteration, setting
                      // a value in this array to NaN and testing for it
                      // is a way of obtaining the iteration counter.
                      if (ido != 99)
                        workl[iptr(5)-1] = octave_NaN; 
                    }

                  if (ido == -1 || ido == 1 || ido == 2)
                    {

                      if (have_a_fun)
                        {
                          Complex *ip2 = workd + iptr(0) - 1;
                          ComplexColumnVector x(n);

                          for (octave_idx_type i = 0; i < n; i++)
                            x(i) = *ip2++;

                          octave_value_list tmp = 
                            func->do_multi_index_op (1, octave_value (x));

                          ComplexColumnVector y = 
                            tmp(0).complex_column_vector_value ();

                          ip2 = workd + iptr(1) - 1;
                          for (octave_idx_type i = 0; i < n; i++)
                            *ip2++ = y(i);
                        }
                      else if (mode == 1)
                        {
                          if (have_b)
                            {
                                  ComplexMatrix mtmp (n,1);
                                  for (octave_idx_type i = 0; i < n; i++)
                                    mtmp(i,0) = workd[i + iptr(0) - 1];
                                  mtmp = sparse_utsolve(bsct, permB, ascm * 
                                    sparse_ltsolve(bscm, permB, mtmp));
                                  for (octave_idx_type i = 0; i < n; i++)
                                    workd[i+iptr(1)-1] = mtmp(i,0);

                            }
                          else
                            sparse_vector_product (ascm, workd + iptr(0) - 1, 
                                                   workd + iptr(1) - 1);
                        }
                      else if (mode == 2)
                        {
                          error ("eigs: Impossible");
                        }
                      else if (mode == 3)
                        {
                          if (have_b)
                            {
                              if (ido == -1)
                                {
                                  if (cholB)
                                    {
                                      error ("fixme!!");
                                    }
                                  else
                                    {
                                      OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);

                                      sparse_vector_product 
                                        (ascm, workd+iptr(0)-1, ctmp);

                                      ComplexMatrix tmp(n, 1);

                                      for (octave_idx_type i = 0; i < n; i++)
                                        tmp(i,0) = ctmp[P[i]];
                                  
                                      sparse_lu_solve (L, U, tmp);

                                      Complex *ip2 = workd+iptr(1)-1;
                                      for (octave_idx_type i = 0; i < n; i++)
                                        ip2[Q[i]] = tmp(i,0);
                                    }
                                }
                              else if (ido == 2)
                                sparse_vector_product (bscm, 
                                                       workd + iptr(0) - 1, 
                                                       workd + iptr(1) - 1);
                              else
                                {
                                  Complex *ip2 = workd+iptr(2)-1;
                                  ComplexMatrix tmp(n, 1);

                                  for (octave_idx_type i = 0; i < n; i++)
                                    tmp(i,0) = ip2[P[i]];
                                  
                                  sparse_lu_solve (L, U, tmp);

                                  ip2 = workd+iptr(1)-1;
                                  for (octave_idx_type i = 0; i < n; i++)
                                    ip2[Q[i]] = tmp(i,0);
                                }
                            }
                          else
                            {
                              if (ido == 2)
                                {
                                  for (octave_idx_type i = 0; i < n; i++)
                                    workd[iptr(0) + i - 1] =
                                      workd[iptr(1) + i - 1];
                                }
                              else
                                {
                                  Complex *ip2 = workd+iptr(0)-1;
                                  ComplexMatrix tmp(n, 1);

                                  for (octave_idx_type i = 0; i < n; i++)
                                    tmp(i,0) = ip2[P[i]];
                                  
                                  sparse_lu_solve (L, U, tmp);

                                  ip2 = workd+iptr(1)-1;
                                  for (octave_idx_type i = 0; i < n; i++)
                                    ip2[Q[i]] = tmp(i,0);
                                }
                            }
                        }
                    }
                  else
                    {
                      if (info < 0)
                        error ("eigs: error %d in znaupd", info);
                      break;
                    }
                } 
              while (1);

              if (!error_state)
                {
                  octave_idx_type info2;

                  int rvec (nargout > 1);

                  // We have a problem in that the size of the C++ bool 
                  // type relative to the fortran logical type. It appears 
                  // that fortran uses 4-bytes per logical and C++ 1-byte 
                  // per bool, though this might be system dependent. As 
                  // long as the HOWMNY arg is not "S", the logical array
                  // is just workspace for ARPACK, so use int to avoid
                  //  problems!!
                  Array<int> s (p);
                  int *sel = s.fortran_vec ();

                  ComplexMatrix eig_vec (n, k);
                  Complex *z = eig_vec.fortran_vec ();

                  ComplexColumnVector eig_val (k+1);
                  Complex *d = eig_val.fortran_vec ();

                  OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);

                  F77_FUNC (zneupd, ZNEUPD) 
                    (rvec, F77_CONST_CHAR_ARG2 ("A", 1),
                     sel, d, z, n, sigma, workev,
                     F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                     F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
                     k, tol, presid, p, v, n, iparam,
                     ipntr, workd, workl, lwork, rwork, info2
                     F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) 
                     F77_CHAR_ARG_LEN(2));

                  if (f77_exception_encountered)
                    {
                      error ("eigs: unrecoverable exception encountered in 
zneupd");
                      goto eigs_err;
                    }

                  if (info2 == 0)
                    {
                      octave_idx_type k2 = k / 2;
                      for (octave_idx_type i = 0; i < k2; i++)
                        {
                          Complex ctmp = d[i];
                          d[i] = d[k - i - 1];
                          d[k - i - 1] = ctmp;
                        }
                      eig_val.resize(k);

                      if (nargout < 2)
                        retval (0) = eig_val;
                      else
                        {
                          OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);

                          for (octave_idx_type i = 0; i < k2; i++)
                            {
                              octave_idx_type off1 = i * n;
                              octave_idx_type off2 = (k - i - 1) * n;

                              if (off1 == off2)
                                continue;

                              for (octave_idx_type j = 0; j < n; j++)
                                ctmp[j] = z[off1 + j];

                              for (octave_idx_type j = 0; j < n; j++)
                                z[off1 + j] = z[off2 + j];

                              for (octave_idx_type j = 0; j < n; j++)
                                z[off2 + j] = ctmp[j];
                            }

                          retval(2) = double (info);
                          retval(1) = ComplexDiagMatrix (eig_val);
                          if (note3)
                            eig_vec = sparse_ltsolve(bscm, permB, eig_vec);

                          retval(0) = eig_vec;
                        }
                    }
                  else
                    error ("eigs: error %d in zneupd", info2);
                }
            }

          if (ip(4) == 0)
            warning ("eigs: None of the %d requested eigenvalues converged", k);
          else if (ip(4) < k)
            warning ("eigs: Only %d of the %d requested eigenvalues converged", 
                     ip(4), k);

        }
    }

 eigs_err:

  if (! fcn_name.empty ())
    clear_function (fcn_name);
#else
  error ("eigs: not available in this version of Octave");
#endif
  return retval;
}

template void
sparse_vector_product (const Sparse<double>&, const double*, double*);

template void
sparse_vector_product (const Sparse<Complex>&, const Complex*, Complex*);

template octave_idx_type
sparse_lu_solve (const SparseMatrix&, const SparseMatrix&, Matrix&);

template octave_idx_type
sparse_lu_solve (const SparseComplexMatrix&, const SparseComplexMatrix&, 
                 ComplexMatrix&);

template ComplexMatrix
sparse_ltsolve (const SparseComplexMatrix&, const ColumnVector&, 
                const ComplexMatrix&);

template Matrix
sparse_ltsolve (const SparseMatrix&, const ColumnVector&, 
                const Matrix&);

template ComplexMatrix
sparse_utsolve (const SparseComplexMatrix&, const ColumnVector&, 
                const ComplexMatrix&);

template Matrix
sparse_utsolve (const SparseMatrix&, const ColumnVector&,
                const Matrix&);

/*

%% Real positive definite tests, n must be even
%!shared n, k, A, d0, d2
%! n = 20;
%! k = 4;
%! A = 
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)]);
%! d0 = eig (A);
%! d2 = sort(d0);
%! [dum, idx] = sort( abs(d0));
%! d0 = d0(idx);
%!test
%! d1 = eigs (A, k);
%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
%!test
%! d1 = eigs (A,k+1);
%! assert (d1, d0(end:-1:(end-k)),1e-12);
%!test
%! d1 = eigs (A, k, "lm");
%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
%!test
%! d1 = eigs (A, k, "sm");
%! assert (d1, d0(k:-1:1), 1e-12);
%!test
%! d1 = eigs (A, k, "la");
%! assert (d1, d2(end:-1:(end-k+1)), 1e-12);
%!test
%! d1 = eigs (A, k, "sa");
%! assert (d1, d2(1:k), 1e-12);
%!test
%! d1 = eigs (A, k, "be");
%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-12);
%!test
%! d1 = eigs (A, k+1, "be");
%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-12);
%!test
%! d1 = eigs (A, k, 4.1);
%! [dum,idx0] = sort (abs(d0 - 4.1));
%! [dum,idx1] = sort (abs(d1 - 4.1));
%! assert (d1(idx1), d0(idx0(1:k)), 1e-12);
%!test
%! d1 = eigs(A, speye(n), k, 'lm');
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12);
%!#test
%! opts.cholB=true;
%! d1 = eigs(A, speye(n), k, 'lm',opts);
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12);
%!test
%! fn = @(x) A * x;
%! opts.issym = 1; opts.isreal = 1;
%! d1 = eigs (fn, n, k, "lm", opts);
%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
%!test
%! [v1,d1] = eigs(A, k, "lm");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "sm");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "la");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "sa");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "be");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor

*/

/*

%% Real unsymmetric tests, n mist be odd to force eigs to use dnaupd/dneupd
%!shared n, k, A, d0
%! n = 20;
%! k = 4;
%! ## FIXME I can get the matrix below to seg-fault the tests and give
%! ## incorrect results for "SI" and very very ocassionally for "SM"..
%! A = 
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),-ones(1,n-2)]);
%! ## A =  
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]);
%! d0 = eig (A);
%! [dum, idx] = sort(abs(d0));
%! d0 = d0(idx);
%!test
%! d1 = eigs (A, k);
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!test
%! d1 = eigs (A,k+1);
%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12);
%!test
%! d1 = eigs (A, k, "lm");
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!test
%! d1 = eigs (A, k, "sm");
%! assert (abs(d1), abs(d0(1:k)), 1e-12);
%!test
%! d1 = eigs (A, k, "lr");
%! [dum, idx] = sort (real(d0));
%! d2 = d0(idx);
%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12);
%!test
%! d1 = eigs (A, k, "sr");
%! [dum, idx] = sort (real(abs(d0)));
%! d2 = d0(idx);
%! assert (real(d1), real(d2(1:k)), 1e-12);
%!test
%! d1 = eigs (A, k, "li");
%! [dum, idx] = sort (imag(abs(d0)));
%! d2 = d0(idx);
%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12);
%!test
%! d1 = eigs (A, k, "si");
%! [dum, idx] = sort (imag(abs(d0)));
%! d2 = d0(idx);
%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12);
%!test
%! d1 = eigs (A, k, 4.1);
%! [dum,idx0] = sort (abs(d0 - 4.1));
%! [dum,idx1] = sort (abs(d1 - 4.1));
%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12);
%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12);
%!test
%! d1 = eigs(A, speye(n), k, 'lm');
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!#test
%! opts.cholB=true;
%! d1 = eigs(A, speye(n), k, 'lm',opts);
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12);
%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))), 1e-12);
%!test
%! fn = @(x) A * x;
%! opts.issym = 0; opts.isreal = 1;
%! d1 = eigs (fn, n, k, "lm", opts);
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!test
%! [v1,d1] = eigs(A, k, "lm");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "sm");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "lr");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "sr");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "li");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "si");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor

*/

/*

%% Complex hermitian tests
%!shared n, k, A, d0
%! n = 20;
%! k = 4;
%! A = 
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)]);
%! d0 = eig (A);
%! [dum, idx] = sort(abs(d0));
%! d0 = d0(idx);
%!test
%! d1 = eigs (A, k);
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!test
%! d1 = eigs (A,k+1);
%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12);
%!test
%! d1 = eigs (A, k, "lm");
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!test
%! d1 = eigs (A, k, "sm");
%! assert (abs(d1), abs(d0(1:k)), 1e-12);
%!test
%! d1 = eigs (A, k, "lr");
%! [dum, idx] = sort (real(abs(d0)));
%! d2 = d0(idx);
%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12);
%!test
%! d1 = eigs (A, k, "sr");
%! [dum, idx] = sort (real(abs(d0)));
%! d2 = d0(idx);
%! assert (real(d1), real(d2(1:k)), 1e-12);
%!test
%! d1 = eigs (A, k, "li");
%! [dum, idx] = sort (imag(abs(d0)));
%! d2 = d0(idx);
%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12);
%!test
%! d1 = eigs (A, k, "si");
%! [dum, idx] = sort (imag(abs(d0)));
%! d2 = d0(idx);
%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12);
%!test
%! d1 = eigs (A, k, 4.1);
%! [dum,idx0] = sort (abs(d0 - 4.1));
%! [dum,idx1] = sort (abs(d1 - 4.1));
%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12);
%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12);
%!test
%! d1 = eigs(A, speye(n), k, 'lm');
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!#test
%! opts.cholB=true;
%! d1 = eigs(A, speye(n), k, 'lm',opts);
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12);
%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))), 1e-12);
%!test
%! fn = @(x) A * x;
%! opts.issym = 0; opts.isreal = 0;
%! d1 = eigs (fn, n, k, "lm", opts);
%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
%!test
%! [v1,d1] = eigs(A, k, "lm");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "sm");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "lr");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "sr");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "li");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor
%!test
%! [v1,d1] = eigs(A, k, "si");
%! d1 = diag(d1);
%! for i=1:k
%!  assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
%! endfor

*/

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
function aerr = segtest (iter)
  %% This will seg-fault octave consistently, but not matlab.
  n=20;
  k=4; 
  A = 
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),-ones(1,n-2)]);
  opts.disp = 0; opts.p = 19;
  aerr = 0;
  for i=1:iter
    [v1,d1] = eigs(A, k, 'sr', opts);
    d1 = diag(d1);
    merr = 0;
    for i=1:k
      newerr = max(abs((A - d1(i)*speye(n))*v1(:,i)));
      if (newerr > merr)
        merr = newerr;
      end
    end
    fprintf('Max Err: %g\n', merr);
    if (merr > aerr)
      aerr = merr;
    end
  end
end
*** ARPACK.orig/SRC/dneupd.f    2000-09-20 22:58:34.000000000 +0200
--- ARPACK/SRC/dneupd.f 2006-11-03 11:27:21.159925329 +0100
***************
*** 353,359 ****
       &           mode  , msglvl, outncv, ritzr   ,
       &           ritzi , wri   , wrr   , irr     ,
       &           iri   , ibd   , ishift, numcnv  ,
!      &           np    , jj 
        logical    reord
        Double precision 
       &           conds  , rnorm, sep  , temp,
--- 353,359 ----
       &           mode  , msglvl, outncv, ritzr   ,
       &           ritzi , wri   , wrr   , irr     ,
       &           iri   , ibd   , ishift, numcnv  ,
!      &           np    , jj    , nconv2
        logical    reord
        Double precision 
       &           conds  , rnorm, sep  , temp,
***************
*** 661,676 ****
       &                   workl(iuptri), ldh          , 
       &                   workl(invsub), ldq          , 
       &                   workl(iheigr), workl(iheigi), 
!      &                   nconv        , conds        ,
       &                   sep          , workl(ihbds) , 
       &                   ncv          , iwork        ,
       &                   1            , ierr)
  c
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000
              end if
  c
              if (msglvl .gt. 2) then
                  call dvout (logfil, ncv, workl(iheigr), ndigit,
       &           '_neupd: Real part of the eigenvalues of H--reordered')
--- 661,681 ----
       &                   workl(iuptri), ldh          , 
       &                   workl(invsub), ldq          , 
       &                   workl(iheigr), workl(iheigi), 
!      &                   nconv2       , conds        ,
       &                   sep          , workl(ihbds) , 
       &                   ncv          , iwork        ,
       &                   1            , ierr)
  c
+             if (nconv2 .lt. nconv) then
+                nconv = nconv2
+             end if
+ 
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000
              end if
  c
+ 
              if (msglvl .gt. 2) then
                  call dvout (logfil, ncv, workl(iheigr), ndigit,
       &           '_neupd: Real part of the eigenvalues of H--reordered')
*** ARPACK.orig/SRC/sneupd.f    2000-09-20 22:58:34.000000000 +0200
--- ARPACK/SRC/sneupd.f 2006-11-03 11:32:54.691913071 +0100
***************
*** 353,359 ****
       &           mode  , msglvl, outncv, ritzr   ,
       &           ritzi , wri   , wrr   , irr     ,
       &           iri   , ibd   , ishift, numcnv  ,
!      &           np    , jj 
        logical    reord
        Real 
       &           conds  , rnorm, sep  , temp,
--- 353,359 ----
       &           mode  , msglvl, outncv, ritzr   ,
       &           ritzi , wri   , wrr   , irr     ,
       &           iri   , ibd   , ishift, numcnv  ,
!      &           np    , jj    , nconv2
        logical    reord
        Real 
       &           conds  , rnorm, sep  , temp,
***************
*** 661,671 ****
       &                   workl(iuptri), ldh          , 
       &                   workl(invsub), ldq          , 
       &                   workl(iheigr), workl(iheigi), 
!      &                   nconv        , conds        ,
       &                   sep          , workl(ihbds) , 
       &                   ncv          , iwork        ,
       &                   1            , ierr)
  c
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000
--- 661,675 ----
       &                   workl(iuptri), ldh          , 
       &                   workl(invsub), ldq          , 
       &                   workl(iheigr), workl(iheigi), 
!      &                   nconv2       , conds        ,
       &                   sep          , workl(ihbds) , 
       &                   ncv          , iwork        ,
       &                   1            , ierr)
  c
+             if (nconv2 .lt. nconv) then
+                nconv = nconv2
+             end if
+ 
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000
*** ARPACK.orig/SRC/zneupd.f    2002-08-15 07:51:12.000000000 +0200
--- ARPACK/SRC/zneupd.f 2006-11-03 11:33:28.032312460 +0100
***************
*** 301,307 ****
       &           invsub, iuptri, iwev  , j    , ldh   , ldq   ,
       &           mode  , msglvl, ritz  , wr   , k     , irz   ,
       &           ibd   , outncv, iq    , np   , numcnv, jj    ,
!      &           ishift
        Complex*16
       &           rnorm, temp, vl(1)
        Double precision
--- 301,307 ----
       &           invsub, iuptri, iwev  , j    , ldh   , ldq   ,
       &           mode  , msglvl, ritz  , wr   , k     , irz   ,
       &           ibd   , outncv, iq    , np   , numcnv, jj    ,
!      &           ishift, nconv2
        Complex*16
       &           rnorm, temp, vl(1)
        Double precision
***************
*** 592,600 ****
              call ztrsen('None'       , 'V'          , select      ,
       &                   ncv          , workl(iuptri), ldh         ,
       &                   workl(invsub), ldq          , workl(iheig),
!      &                   nconv        , conds        , sep         , 
       &                   workev       , ncv          , ierr)
  c
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000
--- 592,604 ----
              call ztrsen('None'       , 'V'          , select      ,
       &                   ncv          , workl(iuptri), ldh         ,
       &                   workl(invsub), ldq          , workl(iheig),
!      &                   nconv2       , conds        , sep         , 
       &                   workev       , ncv          , ierr)
  c
+             if (nconv2 .lt. nconv) then
+                nconv = nconv2
+             end if
+ 
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000
*** ARPACK.orig/SRC/cneupd.f    2002-08-15 07:51:12.000000000 +0200
--- ARPACK/SRC/cneupd.f 2006-11-03 11:31:54.499802784 +0100
***************
*** 301,307 ****
       &           invsub, iuptri, iwev  , j    , ldh   , ldq   ,
       &           mode  , msglvl, ritz  , wr   , k     , irz   ,
       &           ibd   , outncv, iq    , np   , numcnv, jj    ,
!      &           ishift
        Complex
       &           rnorm, temp, vl(1)
        Real
--- 301,307 ----
       &           invsub, iuptri, iwev  , j    , ldh   , ldq   ,
       &           mode  , msglvl, ritz  , wr   , k     , irz   ,
       &           ibd   , outncv, iq    , np   , numcnv, jj    ,
!      &           ishift, nconv2
        Complex
       &           rnorm, temp, vl(1)
        Real
***************
*** 592,600 ****
              call ctrsen('None'       , 'V'          , select      ,
       &                   ncv          , workl(iuptri), ldh         ,
       &                   workl(invsub), ldq          , workl(iheig),
!      &                   nconv        , conds        , sep         , 
       &                   workev       , ncv          , ierr)
  c
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000
--- 592,604 ----
              call ctrsen('None'       , 'V'          , select      ,
       &                   ncv          , workl(iuptri), ldh         ,
       &                   workl(invsub), ldq          , workl(iheig),
!      &                   nconv2       , conds        , sep         , 
       &                   workev       , ncv          , ierr)
  c
+             if (nconv2 .lt. nconv) then
+                nconv = nconv2
+             end if
+ 
              if (ierr .eq. 1) then
                 info = 1
                 go to 9000

reply via email to

[Prev in Thread] Current Thread [Next in Thread]