Mercurial > octave-nkf

--- a/ChangeLog	Tue Dec 23 01:08:59 2008 +0100
+++ b/ChangeLog	Tue Dec 23 08:28:23 2008 +0100
@@ -1,3 +1,10 @@
+2008-12-23  David Bateman  <dbateman@free.fr>
+
+	* configure.in: Add configuration test for ARPACK.
+	* Makeconf.in (ARPACK_LIBS): Add variable with location of ARPACK
+	library.
+	* NEWS: Document that eigs and svds were moved to Octaave.
+
 2008-10-29  Jaroslav Hajek  <highegg@gmail.com>

 	* configure.in: Remove the OCTAVE_LOCAL_BUFFER stuff (moved to
--- a/Makeconf.in	Tue Dec 23 01:08:59 2008 +0100
+++ b/Makeconf.in	Tue Dec 23 08:28:23 2008 +0100
@@ -237,6 +237,7 @@
 CHOLMOD_LIBS = @CHOLMOD_LIBS@
 CXSPARSE_LIBS = @CXSPARSE_LIBS@
 OPENGL_LIBS = @OPENGL_LIBS@
+ARPACK_LIBS = @ARPACK_LIBS@
 LIBS = @LIBS@

 USE_64_BIT_IDX_T = @USE_64_BIT_IDX_T@
--- a/NEWS	Tue Dec 23 01:08:59 2008 +0100
+++ b/NEWS	Tue Dec 23 08:28:23 2008 +0100
@@ -84,6 +84,9 @@
  ** The imwrite and imread function have been included in Octave based
     on the GraphicsMagick library.

+ ** The eigs and svds functions have been included in Octave based on
+    the ARPACK library.
+
  ** Special treatment in the parser of expressions like "a' * b". In
     these cases the transpose is no longer explicitly formed and BLAS
     libraries are called with the transpose flagged. This significantly
--- a/configure.in	Tue Dec 23 01:08:59 2008 +0100
+++ b/configure.in	Tue Dec 23 08:28:23 2008 +0100
@@ -1044,6 +1044,27 @@
   AC_MSG_WARN($warn_cxsparse)
 fi

+ARPACK_LIBS=
+AC_SUBST(ARPACK_LIBS)
+
+AC_ARG_WITH(arpack,
+  [AS_HELP_STRING([--without-arpack],
+     [don't use ARPACK, disable some sparse functionality])],
+  with_arpack=$withval, with_arpack=yes)
+
+warn_arpack="arpack not found. This will result in a lack of the eigs function."
+if test "$with_arpack" = yes; then
+  with_arpack=no
+  AC_CHECK_LIB(arpack, F77_FUNC(dseupd,DSEUPD), [ARPACK_LIBS="-larpack"; with_arpack=yes], , $LAPACK_LIBS $FLIBS)
+  if test "$with_arpack" = yes; then
+    AC_DEFINE(HAVE_ARPACK, 1, [Define if the ARPACK library is used.])
+    warn_arpack=
+  fi
+fi
+if test -n "$warn_arpack"; then
+  AC_MSG_WARN($warn_arpack)
+fi
+
 ### Handle shared library options.

 ### Enable creation of static libraries.
@@ -2012,6 +2033,7 @@
   CCOLAMD libraries:    $CCOLAMD_LIBS
   CHOLMOD libraries:    $CHOLMOD_LIBS
   CXSPARSE libraries:   $CXSPARSE_LIBS
+  ARPACK libraries:     $ARPACK_LIBS
   HDF5 libraries:       $HDF5_LIBS
   CURL libraries:       $CURL_LIBS
   REGEX libraries:      $REGEX_LIBS
@@ -2123,6 +2145,11 @@
   warn_msg_printed=true
 fi

+if test -n "$warn_arpack"; then
+  AC_MSG_WARN($warn_arpack)
+  warn_msg_printed=true
+fi
+
 if test -n "$warn_curl"; then
   AC_MSG_WARN($warn_curl)
   warn_msg_printed=true
--- a/doc/ChangeLog	Tue Dec 23 01:08:59 2008 +0100
+++ b/doc/ChangeLog	Tue Dec 23 08:28:23 2008 +0100
@@ -1,3 +1,7 @@
+2008-12-23  David Bateman  <dbateman@free.fr>
+
+	* interpreter/sparse.txi: Document the eigs and svds functions.
+
 2008-12-02  Thorsten Meyer  <thorsten.meyier@gmx.de>

         * interpreter/container.txi, interpreter/strings.txi:
--- a/doc/interpreter/sparse.txi	Tue Dec 23 01:08:59 2008 +0100
+++ b/doc/interpreter/sparse.txi	Tue Dec 23 08:28:23 2008 +0100
@@ -473,9 +473,8 @@
   @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, @dfn{symrcm}

 @item Linear algebra:
-  @dfn{matrix_type}, @dfn{normest}, @dfn{condest}, @dfn{sprank}
-  @dfn{spaugment}
-@c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now
+  @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, @dfn{normest}, @dfn{sprank},
+  @dfn{spaugment}, @dfn{svds}

 @item Iterative techniques:
   @dfn{luinc}, @dfn{pcg}, @dfn{pcr}
@@ -832,6 +831,15 @@

 @DOCSTRING(spaugment)

+Finally, the function @code{eigs} can be used to calculate a limited
+number of eigenvalues and eigenvectors based on a selection criteria
+and likewise for @code{svds} which calculates a limited number of
+singular values and vectors.
+
+@DOCSTRING(eigs)
+
+@DOCSTRING(svds)
+
 @node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse Matrices
 @section Iterative Techniques applied to sparse matrices
--- a/liboctave/ChangeLog	Tue Dec 23 01:08:59 2008 +0100
+++ b/liboctave/ChangeLog	Tue Dec 23 08:28:23 2008 +0100
@@ -1,3 +1,8 @@
+2008-12-23  David Bateman  <dbateman@free.fr>
+
+	* eigs-base.cc: New file with template wrapper for ARPACK.
+	* Makefile.in (TEMPLATE_SRC): Add it here.
+
 2008-12-16  Jaroslav Hajek  <highegg@gmail.com>

 	* Array.cc (rec_permute_helper): New class.
--- a/liboctave/Makefile.in	Tue Dec 23 01:08:59 2008 +0100
+++ b/liboctave/Makefile.in	Tue Dec 23 08:28:23 2008 +0100
@@ -103,7 +103,7 @@
 	$(VX_OP_INC) \
 	$(SPARSE_MX_OP_INC)

-TEMPLATE_SRC := Array.cc ArrayN.cc DiagArray2.cc \
+TEMPLATE_SRC := Array.cc ArrayN.cc eigs-base.cc DiagArray2.cc \
 	MArray.cc MArray2.cc MArrayN.cc MDiagArray2.cc \
 	base-lu.cc oct-sort.cc sparse-base-lu.cc \
 	sparse-base-chol.cc sparse-dmsolve.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/eigs-base.cc	Tue Dec 23 08:28:23 2008 +0100
@@ -0,0 +1,3768 @@
+/*
+
+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 3 of the License, 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, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cfloat>
+#include <cmath>
+#include <vector>
+#include <iostream>
+
+#include "f77-fcn.h"
+#include "quit.h"
+#include "SparsedbleLU.h"
+#include "SparseCmplxLU.h"
+#include "dSparse.h"
+#include "CSparse.h"
+#include "MatrixType.h"
+#include "SparsedbleCHOL.h"
+#include "SparseCmplxCHOL.h"
+#include "oct-rand.h"
+#include "dbleCHOL.h"
+#include "CmplxCHOL.h"
+#include "dbleLU.h"
+#include "CmplxLU.h"
+
+#ifdef HAVE_ARPACK
+typedef ColumnVector (*EigsFunc) (const ColumnVector &x, int &eigs_error);
+typedef ComplexColumnVector (*EigsComplexFunc)
+  (const ComplexColumnVector &x, int &eigs_error);
+
+// Arpack and blas 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);
+
+  F77_RET_T
+  F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, const octave_idx_type&, const double&,
+			   const double*, const octave_idx_type&, const double*,
+			   const octave_idx_type&, const double&, double*,
+			   const octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL);
+
+
+  F77_RET_T
+  F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
+                           const octave_idx_type&, const octave_idx_type&, const Complex&,
+                           const Complex*, const octave_idx_type&, const Complex*,
+                           const octave_idx_type&, const Complex&, Complex*, const octave_idx_type&
+                           F77_CHAR_ARG_LEN_DECL);
+
+}
+
+
+#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
+static octave_idx_type
+lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
+
+static octave_idx_type
+lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
+	 ComplexMatrix&);
+
+static octave_idx_type
+lusolve (const Matrix&, const Matrix&, Matrix&);
+
+static octave_idx_type
+lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
+
+static ComplexMatrix
+ltsolve (const SparseComplexMatrix&, const ColumnVector&,
+		const ComplexMatrix&);
+
+static Matrix
+ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,);
+
+static ComplexMatrix
+ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+static Matrix
+ltsolve (const Matrix&, const ColumnVector&, const Matrix&,);
+
+static ComplexMatrix
+utsolve (const SparseComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+static Matrix
+utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
+
+static ComplexMatrix
+utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+static Matrix
+utsolve (const Matrix&, const ColumnVector&, const Matrix&);
+
+#endif
+
+template <class M, class SM>
+static octave_idx_type
+lusolve (const SM& L, const SM& U, M& m)
+{
+  octave_idx_type err = 0;
+  double rcond;
+  MatrixType utyp (MatrixType::Upper);
+
+  // Sparse L is lower triangular, Dense L is permuted lower triangular!!!
+  m = L.solve (m, err, rcond, 0);
+  if (err)
+    return err;
+
+  m = U.solve (utyp, m, err, rcond, 0);
+
+  return err;
+}
+
+template <class SM, class M>
+static M
+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>
+static M
+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);
+}
+
+static bool
+vector_product (const SparseMatrix& m, const double* x, double* 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];
+
+  return true;
+}
+
+static bool
+vector_product (const Matrix& m, const double *x, double *y)
+{
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.cols ();
+
+  F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			   nr, nc, 1.0,  m.data (), nr,
+			   x, 1, 0.0, y, 1
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable error in dgemv");
+      return false;
+    }
+  else
+    return true;
+}
+
+static bool
+vector_product (const SparseComplexMatrix& m, const Complex* x,
+			Complex* 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];
+
+  return true;
+}
+
+static bool
+vector_product (const ComplexMatrix& m, const Complex *x, Complex *y)
+{
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.cols ();
+
+  F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			   nr, nc, 1.0,  m.data (), nr,
+			   x, 1, 0.0, y, 1
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (f77_exception_encountered)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable error in zgemv");
+      return false;
+    }
+  else
+    return true;
+}
+
+static bool
+make_cholb (Matrix& b, Matrix& bt, ColumnVector& permB)
+{
+  octave_idx_type info;
+  CHOL fact (b, info);
+  octave_idx_type n = b.cols();
+
+  if (info != 0)
+    return false;
+  else
+    {
+      bt = fact.chol_matrix ();
+      b =  bt.transpose();
+      permB = ColumnVector(n);
+      for (octave_idx_type i = 0; i < n; i++)
+	permB(i) = i;
+      return true;
+    }
+}
+
+static bool
+make_cholb (SparseMatrix& b, SparseMatrix& bt, ColumnVector& permB)
+{
+  octave_idx_type info;
+  SparseCHOL fact (b, info, false);
+
+  if (fact.P() != 0)
+    return false;
+  else
+    {
+      b = fact.L();
+      bt = b.transpose();
+      permB = fact.perm() - 1.0;
+      return true;
+    }
+}
+
+static bool
+make_cholb (ComplexMatrix& b, ComplexMatrix& bt, ColumnVector& permB)
+{
+  octave_idx_type info;
+  ComplexCHOL fact (b, info);
+  octave_idx_type n = b.cols();
+
+  if (info != 0)
+    return false;
+  else
+    {
+      bt = fact.chol_matrix ();
+      b =  bt.hermitian();
+      permB = ColumnVector(n);
+      for (octave_idx_type i = 0; i < n; i++)
+	permB(i) = i;
+      return true;
+    }
+}
+
+static bool
+make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt,
+	    ColumnVector& permB)
+{
+  octave_idx_type info;
+  SparseComplexCHOL fact (b, info, false);
+
+  if (fact.P() != 0)
+    return false;
+  else
+    {
+      b = fact.L();
+      bt = b.hermitian();
+      permB = fact.perm() - 1.0;
+      return true;
+    }
+}
+
+static bool
+LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b,
+		bool cholB, const ColumnVector& permB, double sigma,
+		SparseMatrix &L, SparseMatrix &U, octave_idx_type *P,
+		octave_idx_type *Q)
+{
+  bool have_b = ! b.is_empty ();
+  octave_idx_type n = m.rows();
+
+  // Caclulate LU decomposition of 'A - sigma * B'
+  SparseMatrix AminusSigmaB (m);
+
+  if (have_b)
+    {
+      if (cholB)
+	{
+	  if (permB.length())
+	    {
+	      SparseMatrix tmp(n,n,n);
+	      for (octave_idx_type i = 0; i < n; i++)
+		{
+		  tmp.xcidx(i) = i;
+		  tmp.xridx(i) =
+		    static_cast<octave_idx_type>(permB(i));
+		  tmp.xdata(i) = 1;
+		}
+	      tmp.xcidx(n) = n;
+
+	      AminusSigmaB = AminusSigmaB - sigma * tmp *
+		b.transpose() * b * tmp.transpose();
+	    }
+	  else
+	    AminusSigmaB = AminusSigmaB - sigma *
+	      b.transpose() * b;
+	}
+      else
+	AminusSigmaB = AminusSigmaB - sigma * b;
+    }
+  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) = sigma;
+	      sigmat.xridx(i) = i;
+	      sigmat.xcidx(i+1) = i + 1;
+	    }
+
+	  AminusSigmaB = AminusSigmaB - sigmat;
+	}
+
+  SparseLU fact (AminusSigmaB);
+
+  L = fact.L ();
+  U = fact.U ();
+  const octave_idx_type *P2 = fact.row_perm ();
+  const octave_idx_type *Q2 = fact.col_perm ();
+
+  for (octave_idx_type j = 0; j < n; j++)
+    {
+      P[j] = P2[j];
+      Q[j] = Q2[j];
+    }
+
+  // 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;
+    }
+
+  double rcond = (minU / maxU);
+  volatile double rcond_plus_one = rcond + 1.0;
+
+  if (rcond_plus_one == 1.0 || xisnan (rcond))
+    {
+      (*current_liboctave_warning_handler)
+	("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+      (*current_liboctave_warning_handler)
+	("       an eigenvalue. Convergence is not guaranteed");
+    }
+
+  return true;
+}
+
+static bool
+LuAminusSigmaB (const Matrix &m, const Matrix &b,
+		bool cholB, const ColumnVector& permB, double sigma,
+		Matrix &L, Matrix &U, octave_idx_type *P,
+		octave_idx_type *Q)
+{
+  bool have_b = ! b.is_empty ();
+  octave_idx_type n = m.cols();
+
+  // Caclulate LU decomposition of 'A - sigma * B'
+  Matrix AminusSigmaB (m);
+
+  if (have_b)
+    {
+      if (cholB)
+	{
+	  Matrix tmp = sigma * b.transpose() * b;
+	  const double *pB = permB.fortran_vec();
+	  double *p = AminusSigmaB.fortran_vec();
+
+	  if (permB.length())
+	    {
+	      for (octave_idx_type j = 0;
+		   j < b.cols(); j++)
+		for (octave_idx_type i = 0;
+		     i < b.rows(); i++)
+		  *p++ -=  tmp.xelem (static_cast<octave_idx_type>(pB[i]),
+				      static_cast<octave_idx_type>(pB[j]));
+	    }
+	  else
+	    AminusSigmaB = AminusSigmaB - tmp;
+	}
+      else
+	AminusSigmaB = AminusSigmaB - sigma * b;
+    }
+  else
+    {
+      double *p = AminusSigmaB.fortran_vec();
+
+      for (octave_idx_type i = 0; i < n; i++)
+	p[i*(n+1)] -= sigma;
+    }
+
+  LU fact (AminusSigmaB);
+
+  L = fact.P().transpose() * fact.L ();
+  U = fact.U ();
+  for (octave_idx_type j = 0; j < n; j++)
+    P[j] = Q[j] = j;
+
+  // 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 = std::abs (U.xelem(j,j));
+      if (xisnan (minU) || d < minU)
+	minU = d;
+
+      if (xisnan (maxU) || d > maxU)
+	maxU = d;
+    }
+
+  double rcond = (minU / maxU);
+  volatile double rcond_plus_one = rcond + 1.0;
+
+  if (rcond_plus_one == 1.0 || xisnan (rcond))
+    {
+      (*current_liboctave_warning_handler)
+	("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+      (*current_liboctave_warning_handler)
+	("       an eigenvalue. Convergence is not guaranteed");
+    }
+
+  return true;
+}
+
+static bool
+LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b,
+		bool cholB, const ColumnVector& permB, Complex sigma,
+		SparseComplexMatrix &L, SparseComplexMatrix &U,
+		octave_idx_type *P, octave_idx_type *Q)
+{
+  bool have_b = ! b.is_empty ();
+  octave_idx_type n = m.rows();
+
+  // Caclulate LU decomposition of 'A - sigma * B'
+  SparseComplexMatrix AminusSigmaB (m);
+
+  if (have_b)
+    {
+      if (cholB)
+	{
+	  if (permB.length())
+	    {
+	      SparseMatrix tmp(n,n,n);
+	      for (octave_idx_type i = 0; i < n; i++)
+		{
+		  tmp.xcidx(i) = i;
+		  tmp.xridx(i) =
+		    static_cast<octave_idx_type>(permB(i));
+		  tmp.xdata(i) = 1;
+		}
+	      tmp.xcidx(n) = n;
+
+	      AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b *
+		tmp.transpose() * sigma;
+	    }
+	  else
+	    AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b;
+	}
+      else
+	AminusSigmaB = AminusSigmaB - sigma * b;
+    }
+  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;
+    }
+
+  SparseComplexLU fact (AminusSigmaB);
+
+  L = fact.L ();
+  U = fact.U ();
+  const octave_idx_type *P2 = fact.row_perm ();
+  const octave_idx_type *Q2 = fact.col_perm ();
+
+  for (octave_idx_type j = 0; j < n; j++)
+    {
+      P[j] = P2[j];
+      Q[j] = Q2[j];
+    }
+
+  // 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;
+    }
+
+  double rcond = (minU / maxU);
+  volatile double rcond_plus_one = rcond + 1.0;
+
+  if (rcond_plus_one == 1.0 || xisnan (rcond))
+    {
+      (*current_liboctave_warning_handler)
+	("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+      (*current_liboctave_warning_handler)
+	("       an eigenvalue. Convergence is not guaranteed");
+    }
+
+  return true;
+}
+
+static bool
+LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b,
+		bool cholB, const ColumnVector& permB, Complex sigma,
+		ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P,
+		octave_idx_type *Q)
+{
+  bool have_b = ! b.is_empty ();
+  octave_idx_type n = m.cols();
+
+  // Caclulate LU decomposition of 'A - sigma * B'
+  ComplexMatrix AminusSigmaB (m);
+
+  if (have_b)
+    {
+      if (cholB)
+	{
+	  ComplexMatrix tmp = sigma * b.hermitian() * b;
+	  const double *pB = permB.fortran_vec();
+	  Complex *p = AminusSigmaB.fortran_vec();
+
+	  if (permB.length())
+	    {
+	      for (octave_idx_type j = 0;
+		   j < b.cols(); j++)
+		for (octave_idx_type i = 0;
+		     i < b.rows(); i++)
+		  *p++ -=  tmp.xelem (static_cast<octave_idx_type>(pB[i]),
+				      static_cast<octave_idx_type>(pB[j]));
+	    }
+	  else
+	    AminusSigmaB = AminusSigmaB - tmp;
+	}
+      else
+	AminusSigmaB = AminusSigmaB - sigma * b;
+    }
+  else
+    {
+      Complex *p = AminusSigmaB.fortran_vec();
+
+      for (octave_idx_type i = 0; i < n; i++)
+	p[i*(n+1)] -= sigma;
+    }
+
+  ComplexLU fact (AminusSigmaB);
+
+  L = fact.P().transpose() * fact.L ();
+  U = fact.U ();
+  for (octave_idx_type j = 0; j < n; j++)
+    P[j] = Q[j] = j;
+
+  // 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 = std::abs (U.xelem(j,j));
+      if (xisnan (minU) || d < minU)
+	minU = d;
+
+      if (xisnan (maxU) || d > maxU)
+	maxU = d;
+    }
+
+  double rcond = (minU / maxU);
+  volatile double rcond_plus_one = rcond + 1.0;
+
+  if (rcond_plus_one == 1.0 || xisnan (rcond))
+    {
+      (*current_liboctave_warning_handler)
+	("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+      (*current_liboctave_warning_handler)
+	("       an eigenvalue. Convergence is not guaranteed");
+    }
+
+  return true;
+}
+
+template <class M>
+octave_idx_type
+EigsRealSymmetricMatrix (const M& m, const std::string typ,
+			 octave_idx_type k, octave_idx_type p,
+			 octave_idx_type &info, Matrix &eig_vec,
+			 ColumnVector &eig_val, const M& _b,
+			 ColumnVector &permB, ColumnVector &resid,
+			 std::ostream& os, double tol, int rvec,
+			 bool cholB, int disp, int maxit)
+{
+  M b(_b);
+  octave_idx_type n = m.cols ();
+  octave_idx_type mode = 1;
+  bool have_b = ! b.is_empty();
+  bool note3 = false;
+  char bmat = 'I';
+  double sigma = 0.;
+  M bt;
+
+  if (m.rows() != m.cols())
+    {
+      (*current_liboctave_error_handler) ("eigs: A must be square");
+      return -1;
+    }
+  if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+    {
+      (*current_liboctave_error_handler)
+	("eigs: B must be square and the same size as A");
+      return -1;
+    }
+
+  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)
+    {
+      p = k * 2;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p <= k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k and n");
+      return -1;
+    }
+
+  if (k > n )
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n).\n"
+	 "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (have_b && cholB && permB.length() != 0)
+    {
+      // Check the we really have a permutation vector
+      if (permB.length() != n)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: permB vector invalid");
+	  return -1;
+	}
+      else
+	{
+	  Array<bool> checked(n,false);
+	  for (octave_idx_type i = 0; i < n; i++)
+	    {
+	      octave_idx_type bidx =
+		static_cast<octave_idx_type> (permB(i));
+	      if (checked(bidx) || bidx < 0 ||
+		  bidx >= n || D_NINT (bidx) != bidx)
+		{
+		  (*current_liboctave_error_handler)
+		    ("eigs: permB vector invalid");
+		  return -1;
+		}
+	    }
+	}
+    }
+
+  if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+      typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+      typ != "SI")
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecognized sigma value");
+      return -1;
+    }
+
+  if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
+    {
+      (*current_liboctave_error_handler)
+	("eigs: invalid sigma value for real symmetric problem");
+      return -1;
+    }
+
+  if (have_b)
+    {
+      // See Note 3 dsaupd
+      note3 = true;
+      if (cholB)
+	{
+	  bt = b;
+	  b = b.transpose();
+	  if (permB.length() == 0)
+	    {
+	      permB = ColumnVector(n);
+	      for (octave_idx_type i = 0; i < n; i++)
+		permB(i) = i;
+	    }
+	}
+      else
+	{
+	  if (! make_cholb(b, bt, permB))
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: The matrix B is not positive definite");
+	      return -1;
+	    }
+	}
+    }
+
+  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;
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in dsaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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_b)
+	    {
+	      Matrix mtmp (n,1);
+	      for (octave_idx_type i = 0; i < n; i++)
+		mtmp(i,0) = workd[i + iptr(0) - 1];
+
+	      mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
+
+	      for (octave_idx_type i = 0; i < n; i++)
+		workd[i+iptr(1)-1] = mtmp(i,0);
+	    }
+	  else if (!vector_product (m, workd + iptr(0) - 1,
+				    workd + iptr(1) - 1))
+	    break;
+	}
+      else
+	{
+	  if (info < 0)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dsaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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 ();
+
+  eig_vec.resize (n, k);
+  double *z = eig_vec.fortran_vec ();
+
+  eig_val.resize (k);
+  double *d = eig_val.fortran_vec ();
+
+  F77_FUNC (dseupd, DSEUPD)
+    (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
+     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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in dseupd");
+      return -1;
+    }
+  else
+    {
+      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];
+		    }
+		}
+
+	      if (note3)
+		eig_vec = ltsolve(b, permB, eig_vec);
+	    }
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: error %d in dseupd", info2);
+	  return -1;
+	}
+    }
+
+  return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsRealSymmetricMatrixShift (const M& m, double sigma,
+			      octave_idx_type k, octave_idx_type p,
+			      octave_idx_type &info, Matrix &eig_vec,
+			      ColumnVector &eig_val, const M& _b,
+			      ColumnVector &permB, ColumnVector &resid,
+			      std::ostream& os, double tol, int rvec,
+			      bool cholB, int disp, int maxit)
+{
+  M b(_b);
+  octave_idx_type n = m.cols ();
+  octave_idx_type mode = 3;
+  bool have_b = ! b.is_empty();
+  std::string typ = "LM";
+
+  if (m.rows() != m.cols())
+    {
+      (*current_liboctave_error_handler) ("eigs: A must be square");
+      return -1;
+    }
+  if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+    {
+      (*current_liboctave_error_handler)
+	("eigs: B must be square and the same size as A");
+      return -1;
+    }
+
+  // FIXME: The "SM" type for mode 1 seems unstable though faster!!
+  //if (! std::abs (sigma))
+  //  return EigsRealSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
+  //				    _b, permB, resid, os, tol, rvec, cholB,
+  //				    disp, maxit);
+
+  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)
+    {
+      p = k * 2;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p <= k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k and n");
+      return -1;
+    }
+
+  if (k > n )
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n).\n"
+	     "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (have_b && cholB && permB.length() != 0)
+    {
+      // Check the we really have a permutation vector
+      if (permB.length() != n)
+	{
+	  (*current_liboctave_error_handler) ("eigs: permB vector invalid");
+	  return -1;
+	}
+      else
+	{
+	  Array<bool> checked(n,false);
+	  for (octave_idx_type i = 0; i < n; i++)
+	    {
+	      octave_idx_type bidx =
+		static_cast<octave_idx_type> (permB(i));
+	      if (checked(bidx) || bidx < 0 ||
+		  bidx >= n || D_NINT (bidx) != bidx)
+		{
+		  (*current_liboctave_error_handler)
+		    ("eigs: permB vector invalid");
+		  return -1;
+		}
+	    }
+	}
+    }
+
+  char bmat = 'I';
+  if (have_b)
+    bmat = 'G';
+
+  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;
+  M L, U;
+
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
+
+  if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q))
+    return -1;
+
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in dsaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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_b)
+	    {
+	      if (ido == -1)
+		{
+		  OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+		  vector_product (m, workd+iptr(0)-1, dtmp);
+
+		  Matrix tmp(n, 1);
+
+		  for (octave_idx_type i = 0; i < n; i++)
+		    tmp(i,0) = dtmp[P[i]];
+
+		  lusolve (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)
+		vector_product (b, 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]];
+
+		  lusolve (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]];
+
+		  lusolve (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)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dsaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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 ();
+
+  eig_vec.resize (n, k);
+  double *z = eig_vec.fortran_vec ();
+
+  eig_val.resize (k);
+  double *d = eig_val.fortran_vec ();
+
+  F77_FUNC (dseupd, DSEUPD)
+    (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
+     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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in dseupd");
+      return -1;
+    }
+  else
+    {
+      if (info2 == 0)
+	{
+	  octave_idx_type k2 = k / 2;
+	  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)
+	    {
+	      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];
+		}
+	    }
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: error %d in dseupd", info2);
+	  return -1;
+	}
+    }
+
+  return ip(4);
+}
+
+octave_idx_type
+EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
+		       const std::string &_typ, double sigma,
+		       octave_idx_type k, octave_idx_type p,
+		       octave_idx_type &info, Matrix &eig_vec,
+		       ColumnVector &eig_val, ColumnVector &resid,
+		       std::ostream& os, double tol, int rvec, bool cholB,
+		       int disp, int maxit)
+{
+  std::string typ (_typ);
+  bool have_sigma = (sigma ? true : false);
+  char bmat = 'I';
+  octave_idx_type mode = 1;
+  int err = 0;
+
+  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)
+    {
+      p = k * 2;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p <= k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k and n");
+      return -1;
+    }
+
+  if (k > n )
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n).\n"
+	     "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (! have_sigma)
+    {
+      if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+	  typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+	  typ != "SI")
+	(*current_liboctave_error_handler)
+	  ("eigs: unrecognized sigma value");
+
+      if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: invalid sigma value for real symmetric problem");
+	  return -1;
+	}
+
+      if (typ == "SM")
+	{
+	  typ = "LM";
+	  sigma = 0.;
+	  mode = 3;
+	}
+    }
+  else if (! std::abs (sigma))
+    typ = "SM";
+  else
+    {
+      typ = "LM";
+      mode = 3;
+    }
+
+  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;
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in dsaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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)
+	{
+	  double *ip2 = workd + iptr(0) - 1;
+	  ColumnVector x(n);
+
+	  for (octave_idx_type i = 0; i < n; i++)
+	    x(i) = *ip2++;
+
+	  ColumnVector y = fun (x, err);
+
+	  if (err)
+	    return false;
+
+	  ip2 = workd + iptr(1) - 1;
+	  for (octave_idx_type i = 0; i < n; i++)
+	    *ip2++ = y(i);
+	}
+      else
+	{
+	  if (info < 0)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dsaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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 ();
+
+  eig_vec.resize (n, k);
+  double *z = eig_vec.fortran_vec ();
+
+  eig_val.resize (k);
+  double *d = eig_val.fortran_vec ();
+
+  F77_FUNC (dseupd, DSEUPD)
+    (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
+     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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in dseupd");
+      return -1;
+    }
+  else
+    {
+      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];
+		    }
+		}
+	    }
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: error %d in dseupd", info2);
+	  return -1;
+	}
+    }
+
+  return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsRealNonSymmetricMatrix (const M& m, const std::string typ,
+			    octave_idx_type k, octave_idx_type p,
+			    octave_idx_type &info, ComplexMatrix &eig_vec,
+			    ComplexColumnVector &eig_val, const M& _b,
+			    ColumnVector &permB, ColumnVector &resid,
+			    std::ostream& os, double tol, int rvec,
+			    bool cholB, int disp, int maxit)
+{
+  M b(_b);
+  octave_idx_type n = m.cols ();
+  octave_idx_type mode = 1;
+  bool have_b = ! b.is_empty();
+  bool note3 = false;
+  char bmat = 'I';
+  double sigmar = 0.;
+  double sigmai = 0.;
+  M bt;
+
+  if (m.rows() != m.cols())
+    {
+      (*current_liboctave_error_handler) ("eigs: A must be square");
+      return -1;
+    }
+  if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+    {
+      (*current_liboctave_error_handler)
+	("eigs: B must be square and the same size as A");
+      return -1;
+    }
+
+  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)
+    {
+      p = k * 2 + 1;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p < k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k+1 and n");
+      return -1;
+    }
+
+  if (k > n - 1)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+	 "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (have_b && cholB && permB.length() != 0)
+    {
+      // Check the we really have a permutation vector
+      if (permB.length() != n)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: permB vector invalid");
+	  return -1;
+	}
+      else
+	{
+	  Array<bool> checked(n,false);
+	  for (octave_idx_type i = 0; i < n; i++)
+	    {
+	      octave_idx_type bidx =
+		static_cast<octave_idx_type> (permB(i));
+	      if (checked(bidx) || bidx < 0 ||
+		  bidx >= n || D_NINT (bidx) != bidx)
+		{
+		  (*current_liboctave_error_handler)
+		    ("eigs: permB vector invalid");
+		  return -1;
+		}
+	    }
+	}
+    }
+
+  if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+      typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+      typ != "SI")
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecognized sigma value");
+      return -1;
+    }
+
+  if (typ == "LA" || typ == "SA" || typ == "BE")
+    {
+      (*current_liboctave_error_handler)
+	("eigs: invalid sigma value for unsymmetric problem");
+      return -1;
+    }
+
+  if (have_b)
+    {
+      // See Note 3 dsaupd
+      note3 = true;
+      if (cholB)
+	{
+	  bt = b;
+	  b = b.transpose();
+	  if (permB.length() == 0)
+	    {
+	      permB = ColumnVector(n);
+	      for (octave_idx_type i = 0; i < n; i++)
+		permB(i) = i;
+	    }
+	}
+      else
+	{
+	  if (! make_cholb(b, bt, permB))
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: The matrix B is not positive definite");
+	      return -1;
+	    }
+	}
+    }
+
+  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;
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in dnaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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_b)
+	    {
+	      Matrix mtmp (n,1);
+	      for (octave_idx_type i = 0; i < n; i++)
+		mtmp(i,0) = workd[i + iptr(0) - 1];
+
+	      mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
+
+	      for (octave_idx_type i = 0; i < n; i++)
+		workd[i+iptr(1)-1] = mtmp(i,0);
+	    }
+	  else if (!vector_product (m, workd + iptr(0) - 1,
+				    workd + iptr(1) - 1))
+	    break;
+	}
+      else
+	{
+	  if (info < 0)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dnaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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_vec2 (n, k + 1);
+  double *z = eig_vec2.fortran_vec ();
+
+  OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
+  OCTAVE_LOCAL_BUFFER (double, di, k + 1);
+  OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
+  for (octave_idx_type i = 0; i < k+1; i++)
+    dr[i] = di[i] = 0.;
+
+  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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in dneupd");
+      return -1;
+    }
+  else
+    {
+      eig_val.resize (k+1);
+      Complex *d = eig_val.fortran_vec ();
+
+      if (info2 == 0)
+	{
+	  octave_idx_type jj = 0;
+	  for (octave_idx_type i = 0; i < k+1; i++)
+	    {
+	      if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
+		jj++;
+	      else
+		d [i-jj] = Complex (dr[i], di[i]);
+	    }
+	  if (jj == 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 (rvec)
+	    {
+	      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];
+		}
+
+	      eig_vec.resize (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_vec(j,i) =
+			  Complex(z[j+off1],0.);
+		      i++;
+		    }
+		  else
+		    {
+		      for (octave_idx_type j = 0; j < n; j++)
+			{
+			  eig_vec(j,i) =
+			    Complex(z[j+off1],z[j+off2]);
+			  if (i < k - 1)
+			    eig_vec(j,i+1) =
+			      Complex(z[j+off1],-z[j+off2]);
+			}
+		      i+=2;
+		    }
+		}
+
+	      if (note3)
+		eig_vec = ltsolve(M (b), permB, eig_vec);
+	    }
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: error %d in dneupd", info2);
+	  return -1;
+	}
+    }
+
+  return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsRealNonSymmetricMatrixShift (const M& m, double sigmar,
+				 octave_idx_type k, octave_idx_type p,
+				 octave_idx_type &info,
+				 ComplexMatrix &eig_vec,
+				 ComplexColumnVector &eig_val, const M& _b,
+				 ColumnVector &permB, ColumnVector &resid,
+				 std::ostream& os, double tol, int rvec,
+				 bool cholB, int disp, int maxit)
+{
+  M b(_b);
+  octave_idx_type n = m.cols ();
+  octave_idx_type mode = 3;
+  bool have_b = ! b.is_empty();
+  std::string typ = "LM";
+  double sigmai = 0.;
+
+  if (m.rows() != m.cols())
+    {
+      (*current_liboctave_error_handler) ("eigs: A must be square");
+      return -1;
+    }
+  if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+    {
+      (*current_liboctave_error_handler)
+	("eigs: B must be square and the same size as A");
+      return -1;
+    }
+
+  // FIXME: The "SM" type for mode 1 seems unstable though faster!!
+  //if (! std::abs (sigmar))
+  //  return EigsRealNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
+  //				       _b, permB, resid, os, tol, rvec, cholB,
+  //				       disp, maxit);
+
+  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)
+    {
+      p = k * 2 + 1;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p < k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k+1 and n");
+      return -1;
+    }
+
+  if (k > n - 1)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+	     "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (have_b && cholB && permB.length() != 0)
+    {
+      // Check that we really have a permutation vector
+      if (permB.length() != n)
+	{
+	  (*current_liboctave_error_handler) ("eigs: permB vector invalid");
+	  return -1;
+	}
+      else
+	{
+	  Array<bool> checked(n,false);
+	  for (octave_idx_type i = 0; i < n; i++)
+	    {
+	      octave_idx_type bidx =
+		static_cast<octave_idx_type> (permB(i));
+	      if (checked(bidx) || bidx < 0 ||
+		  bidx >= n || D_NINT (bidx) != bidx)
+		{
+		  (*current_liboctave_error_handler)
+		    ("eigs: permB vector invalid");
+		  return -1;
+		}
+	    }
+	}
+    }
+
+  char bmat = 'I';
+  if (have_b)
+    bmat = 'G';
+
+  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;
+  M L, U;
+
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
+
+  if (! LuAminusSigmaB(m, b, cholB, permB, sigmar, L, U, P, Q))
+    return -1;
+
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in dsaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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_b)
+	    {
+	      if (ido == -1)
+		{
+		  OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+		  vector_product (m, workd+iptr(0)-1, dtmp);
+
+		  Matrix tmp(n, 1);
+
+		  for (octave_idx_type i = 0; i < n; i++)
+		    tmp(i,0) = dtmp[P[i]];
+
+		  lusolve (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)
+		vector_product (b, 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]];
+
+		  lusolve (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]];
+
+		  lusolve (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)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dsaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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_vec2 (n, k + 1);
+  double *z = eig_vec2.fortran_vec ();
+
+  OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
+  OCTAVE_LOCAL_BUFFER (double, di, k + 1);
+  OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
+  for (octave_idx_type i = 0; i < k+1; i++)
+    dr[i] = di[i] = 0.;
+
+  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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in dneupd");
+      return -1;
+    }
+  else
+    {
+      eig_val.resize (k+1);
+      Complex *d = eig_val.fortran_vec ();
+
+      if (info2 == 0)
+	{
+	  octave_idx_type jj = 0;
+	  for (octave_idx_type i = 0; i < k+1; i++)
+	    {
+	      if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
+		jj++;
+	      else
+		d [i-jj] = Complex (dr[i], di[i]);
+	    }
+	  if (jj == 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 (rvec)
+	    {
+	      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];
+		}
+
+	      eig_vec.resize (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_vec(j,i) =
+			  Complex(z[j+off1],0.);
+		      i++;
+		    }
+		  else
+		    {
+		      for (octave_idx_type j = 0; j < n; j++)
+			{
+			  eig_vec(j,i) =
+			    Complex(z[j+off1],z[j+off2]);
+			  if (i < k - 1)
+			    eig_vec(j,i+1) =
+			      Complex(z[j+off1],-z[j+off2]);
+			}
+		      i+=2;
+		    }
+		}
+	    }
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: error %d in dneupd", info2);
+	  return -1;
+	}
+    }
+
+  return ip(4);
+}
+
+octave_idx_type
+EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
+			  const std::string &_typ, double sigmar,
+			  octave_idx_type k, octave_idx_type p,
+			  octave_idx_type &info, ComplexMatrix &eig_vec,
+			  ComplexColumnVector &eig_val, ColumnVector &resid,
+			  std::ostream& os, double tol, int rvec, bool cholB,
+			  int disp, int maxit)
+{
+  std::string typ (_typ);
+  bool have_sigma = (sigmar ? true : false);
+  char bmat = 'I';
+  double sigmai = 0.;
+  octave_idx_type mode = 1;
+  int err = 0;
+
+  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)
+    {
+      p = k * 2 + 1;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p < k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k+1 and n");
+      return -1;
+    }
+
+  if (k > n - 1)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+	     "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+
+  if (! have_sigma)
+    {
+      if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+	  typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+	  typ != "SI")
+	(*current_liboctave_error_handler)
+	  ("eigs: unrecognized sigma value");
+
+      if (typ == "LA" || typ == "SA" || typ == "BE")
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: invalid sigma value for unsymmetric problem");
+	  return -1;
+	}
+
+      if (typ == "SM")
+	{
+	  typ = "LM";
+	  sigmar = 0.;
+	  mode = 3;
+	}
+    }
+  else if (! std::abs (sigmar))
+    typ = "SM";
+  else
+    {
+      typ = "LM";
+      mode = 3;
+    }
+
+  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;
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in dnaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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)
+	{
+	  double *ip2 = workd + iptr(0) - 1;
+	  ColumnVector x(n);
+
+	  for (octave_idx_type i = 0; i < n; i++)
+	    x(i) = *ip2++;
+
+	  ColumnVector y = fun (x, err);
+
+	  if (err)
+	    return false;
+
+	  ip2 = workd + iptr(1) - 1;
+	  for (octave_idx_type i = 0; i < n; i++)
+	    *ip2++ = y(i);
+	}
+      else
+	{
+	  if (info < 0)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dsaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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_vec2 (n, k + 1);
+  double *z = eig_vec2.fortran_vec ();
+
+  OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
+  OCTAVE_LOCAL_BUFFER (double, di, k + 1);
+  OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
+  for (octave_idx_type i = 0; i < k+1; i++)
+    dr[i] = di[i] = 0.;
+
+  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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in dneupd");
+      return -1;
+    }
+  else
+    {
+      eig_val.resize (k+1);
+      Complex *d = eig_val.fortran_vec ();
+
+      if (info2 == 0)
+	{
+	  octave_idx_type jj = 0;
+	  for (octave_idx_type i = 0; i < k+1; i++)
+	    {
+	      if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
+		jj++;
+	      else
+		d [i-jj] = Complex (dr[i], di[i]);
+	    }
+	  if (jj == 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 (rvec)
+	    {
+	      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];
+		}
+
+	      eig_vec.resize (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_vec(j,i) =
+			  Complex(z[j+off1],0.);
+		      i++;
+		    }
+		  else
+		    {
+		      for (octave_idx_type j = 0; j < n; j++)
+			{
+			  eig_vec(j,i) =
+			    Complex(z[j+off1],z[j+off2]);
+			  if (i < k - 1)
+			    eig_vec(j,i+1) =
+			      Complex(z[j+off1],-z[j+off2]);
+			}
+		      i+=2;
+		    }
+		}
+	    }
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: error %d in dneupd", info2);
+	  return -1;
+	}
+    }
+
+  return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsComplexNonSymmetricMatrix (const M& m, const std::string typ,
+			       octave_idx_type k, octave_idx_type p,
+			       octave_idx_type &info, ComplexMatrix &eig_vec,
+			       ComplexColumnVector &eig_val, const M& _b,
+			       ColumnVector &permB,
+			       ComplexColumnVector &cresid,
+			       std::ostream& os, double tol, int rvec,
+			       bool cholB, int disp, int maxit)
+{
+  M b(_b);
+  octave_idx_type n = m.cols ();
+  octave_idx_type mode = 1;
+  bool have_b = ! b.is_empty();
+  bool note3 = false;
+  char bmat = 'I';
+  Complex sigma = 0.;
+  M bt;
+
+  if (m.rows() != m.cols())
+    {
+      (*current_liboctave_error_handler) ("eigs: A must be square");
+      return -1;
+    }
+  if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+    {
+      (*current_liboctave_error_handler)
+	("eigs: B must be square and the same size as A");
+      return -1;
+    }
+
+  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);
+    }
+
+  if (p < 0)
+    {
+      p = k * 2 + 1;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p < k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k+1 and n");
+      return -1;
+    }
+
+  if (k > n - 1)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+	 "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (have_b && cholB && permB.length() != 0)
+    {
+      // Check the we really have a permutation vector
+      if (permB.length() != n)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: permB vector invalid");
+	  return -1;
+	}
+      else
+	{
+	  Array<bool> checked(n,false);
+	  for (octave_idx_type i = 0; i < n; i++)
+	    {
+	      octave_idx_type bidx =
+		static_cast<octave_idx_type> (permB(i));
+	      if (checked(bidx) || bidx < 0 ||
+		  bidx >= n || D_NINT (bidx) != bidx)
+		{
+		  (*current_liboctave_error_handler)
+		    ("eigs: permB vector invalid");
+		  return -1;
+		}
+	    }
+	}
+    }
+
+  if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+      typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+      typ != "SI")
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecognized sigma value");
+      return -1;
+    }
+
+  if (typ == "LA" || typ == "SA" || typ == "BE")
+    {
+      (*current_liboctave_error_handler)
+	("eigs: invalid sigma value for complex problem");
+      return -1;
+    }
+
+  if (have_b)
+    {
+      // See Note 3 dsaupd
+      note3 = true;
+      if (cholB)
+	{
+	  bt = b;
+	  b = b.hermitian();
+	  if (permB.length() == 0)
+	    {
+	      permB = ColumnVector(n);
+	      for (octave_idx_type i = 0; i < n; i++)
+		permB(i) = i;
+	    }
+	}
+      else
+	{
+	  if (! make_cholb(b, bt, permB))
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: The matrix B is not positive definite");
+	      return -1;
+	    }
+	}
+    }
+
+  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;
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in znaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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_b)
+	    {
+	      ComplexMatrix mtmp (n,1);
+	      for (octave_idx_type i = 0; i < n; i++)
+		mtmp(i,0) = workd[i + iptr(0) - 1];
+	      mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
+	      for (octave_idx_type i = 0; i < n; i++)
+		workd[i+iptr(1)-1] = mtmp(i,0);
+
+	    }
+	  else if (!vector_product (m, workd + iptr(0) - 1,
+				    workd + iptr(1) - 1))
+	    break;
+	}
+      else
+	{
+	  if (info < 0)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in znaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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 ();
+
+  eig_vec.resize (n, k);
+  Complex *z = eig_vec.fortran_vec ();
+
+  eig_val.resize (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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in zneupd");
+      return -1;
+    }
+
+  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 (rvec)
+	{
+	  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];
+	    }
+
+	  if (note3)
+	    eig_vec = ltsolve(b, permB, eig_vec);
+	}
+    }
+  else
+    {
+      (*current_liboctave_error_handler)
+	("eigs: error %d in zneupd", info2);
+      return -1;
+    }
+
+  return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma,
+				    octave_idx_type k, octave_idx_type p,
+				    octave_idx_type &info,
+				    ComplexMatrix &eig_vec,
+				    ComplexColumnVector &eig_val, const M& _b,
+				    ColumnVector &permB,
+				    ComplexColumnVector &cresid,
+				    std::ostream& os, double tol, int rvec,
+				    bool cholB, int disp, int maxit)
+{
+  M b(_b);
+  octave_idx_type n = m.cols ();
+  octave_idx_type mode = 3;
+  bool have_b = ! b.is_empty();
+  std::string typ = "LM";
+
+  if (m.rows() != m.cols())
+    {
+      (*current_liboctave_error_handler) ("eigs: A must be square");
+      return -1;
+    }
+  if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+    {
+      (*current_liboctave_error_handler)
+	("eigs: B must be square and the same size as A");
+      return -1;
+    }
+
+  // FIXME: The "SM" type for mode 1 seems unstable though faster!!
+  //if (! std::abs (sigma))
+  //  return EigsComplexNonSymmetricMatrix (m, "SM", k, p, info, eig_vec,
+  //					  eig_val, _b, permB, cresid, os, tol,
+  //					  rvec, cholB, disp, maxit);
+
+  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);
+    }
+
+  if (p < 0)
+    {
+      p = k * 2 + 1;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p < k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k+1 and n");
+      return -1;
+    }
+
+  if (k > n - 1)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+	     "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (have_b && cholB && permB.length() != 0)
+    {
+      // Check that we really have a permutation vector
+      if (permB.length() != n)
+	{
+	  (*current_liboctave_error_handler) ("eigs: permB vector invalid");
+	  return -1;
+	}
+      else
+	{
+	  Array<bool> checked(n,false);
+	  for (octave_idx_type i = 0; i < n; i++)
+	    {
+	      octave_idx_type bidx =
+		static_cast<octave_idx_type> (permB(i));
+	      if (checked(bidx) || bidx < 0 ||
+		  bidx >= n || D_NINT (bidx) != bidx)
+		{
+		  (*current_liboctave_error_handler)
+		    ("eigs: permB vector invalid");
+		  return -1;
+		}
+	    }
+	}
+    }
+
+  char bmat = 'I';
+  if (have_b)
+    bmat = 'G';
+
+  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;
+  M L, U;
+
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
+
+  if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q))
+    return -1;
+
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in znaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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_b)
+	    {
+	      if (ido == -1)
+		{
+		  OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
+
+		  vector_product (m, workd+iptr(0)-1, ctmp);
+
+		  ComplexMatrix tmp(n, 1);
+
+		  for (octave_idx_type i = 0; i < n; i++)
+		    tmp(i,0) = ctmp[P[i]];
+
+		  lusolve (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)
+		vector_product (b, 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]];
+
+		  lusolve (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]];
+
+		  lusolve (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)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dsaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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 ();
+
+  eig_vec.resize (n, k);
+  Complex *z = eig_vec.fortran_vec ();
+
+  eig_val.resize (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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in zneupd");
+      return -1;
+    }
+
+  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 (rvec)
+	{
+	  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];
+	    }
+	}
+    }
+  else
+    {
+      (*current_liboctave_error_handler)
+	("eigs: error %d in zneupd", info2);
+      return -1;
+    }
+
+  return ip(4);
+}
+
+octave_idx_type
+EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
+			     const std::string &_typ, Complex sigma,
+			     octave_idx_type k, octave_idx_type p,
+			     octave_idx_type &info, ComplexMatrix &eig_vec,
+			     ComplexColumnVector &eig_val,
+			     ComplexColumnVector &cresid, std::ostream& os,
+			     double tol, int rvec, bool cholB, int disp,
+			     int maxit)
+{
+  std::string typ (_typ);
+  bool have_sigma = (std::abs(sigma) ? true : false);
+  char bmat = 'I';
+  octave_idx_type mode = 1;
+  int err = 0;
+
+  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);
+    }
+
+  if (p < 0)
+    {
+      p = k * 2 + 1;
+
+      if (p < 20)
+	p = 20;
+
+      if (p > n - 1)
+	p = n - 1 ;
+    }
+  else if (p < k || p > n)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: opts.p must be between k+1 and n");
+      return -1;
+    }
+
+  if (k > n - 1)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+	     "      Use 'eig(full(A))' instead");
+      return -1;
+    }
+
+  if (! have_sigma)
+    {
+      if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+	  typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+	  typ != "SI")
+	(*current_liboctave_error_handler)
+	  ("eigs: unrecognized sigma value");
+
+      if (typ == "LA" || typ == "SA" || typ == "BE")
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: invalid sigma value for complex problem");
+	  return -1;
+	}
+
+      if (typ == "SM")
+	{
+	  typ = "LM";
+	  sigma = 0.;
+	  mode = 3;
+	}
+    }
+  else if (! std::abs (sigma))
+    typ = "SM";
+  else
+    {
+      typ = "LM";
+      mode = 3;
+    }
+
+  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;
+  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)
+	{
+	  (*current_liboctave_error_handler)
+	    ("eigs: unrecoverable exception encountered in znaupd");
+	  return -1;
+	}
+
+      if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+	{
+	  if (iter++)
+	    {
+	      os << "Iteration " << iter - 1 <<
+		": a few Ritz values of the " << p << "-by-" <<
+		p << " matrix\n";
+	      for (int i = 0 ; i < k; i++)
+		os << "    " << 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)
+	{
+	  Complex *ip2 = workd + iptr(0) - 1;
+	  ComplexColumnVector x(n);
+
+	  for (octave_idx_type i = 0; i < n; i++)
+	    x(i) = *ip2++;
+
+	  ComplexColumnVector y = fun (x, err);
+
+	  if (err)
+	    return false;
+
+	  ip2 = workd + iptr(1) - 1;
+	  for (octave_idx_type i = 0; i < n; i++)
+	    *ip2++ = y(i);
+	}
+      else
+	{
+	  if (info < 0)
+	    {
+	      (*current_liboctave_error_handler)
+		("eigs: error %d in dsaupd", info);
+	      return -1;
+	    }
+	  break;
+	}
+    }
+  while (1);
+
+  octave_idx_type info2;
+
+  // 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 ();
+
+  eig_vec.resize (n, k);
+  Complex *z = eig_vec.fortran_vec ();
+
+  eig_val.resize (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)
+    {
+      (*current_liboctave_error_handler)
+	("eigs: unrecoverable exception encountered in zneupd");
+      return -1;
+    }
+
+  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 (rvec)
+	{
+	  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];
+	    }
+	}
+    }
+  else
+    {
+      (*current_liboctave_error_handler)
+	("eigs: error %d in zneupd", info2);
+      return -1;
+    }
+
+  return ip(4);
+}
+
+#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
+extern octave_idx_type
+EigsRealSymmetricMatrix (const Matrix& m, const std::string typ,
+			 octave_idx_type k, octave_idx_type p,
+			 octave_idx_type &info, Matrix &eig_vec,
+			 ColumnVector &eig_val, const Matrix& b,
+			 ColumnVector &permB, ColumnVector &resid,
+			 std::ostream &os, double tol = DBL_EPSILON,
+			 int rvec = 0, bool cholB = 0, int disp = 0,
+			 int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ,
+			 octave_idx_type k, octave_idx_type p,
+			 octave_idx_type &info, Matrix &eig_vec,
+			 ColumnVector &eig_val, const SparseMatrix& b,
+			 ColumnVector &permB, ColumnVector &resid,
+			 std::ostream& os, double tol = DBL_EPSILON,
+			 int rvec = 0, bool cholB = 0, int disp = 0,
+			 int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricMatrixShift (const Matrix& m, double sigma,
+			      octave_idx_type k, octave_idx_type p,
+			      octave_idx_type &info, Matrix &eig_vec,
+			      ColumnVector &eig_val, const Matrix& b,
+			      ColumnVector &permB, ColumnVector &resid,
+			      std::ostream &os, double tol = DBL_EPSILON,
+			      int rvec = 0, bool cholB = 0, int disp = 0,
+			      int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma,
+			      octave_idx_type k, octave_idx_type p,
+			      octave_idx_type &info, Matrix &eig_vec,
+			      ColumnVector &eig_val, const SparseMatrix& b,
+			      ColumnVector &permB, ColumnVector &resid,
+			      std::ostream &os, double tol = DBL_EPSILON,
+			      int rvec = 0, bool cholB = 0, int disp = 0,
+			      int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
+		       const std::string &typ, double sigma,
+		       octave_idx_type k, octave_idx_type p,
+		       octave_idx_type &info,
+		       Matrix &eig_vec, ColumnVector &eig_val,
+		       ColumnVector &resid, std::ostream &os,
+		       double tol = DBL_EPSILON, int rvec = 0,
+		       bool cholB = 0, int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ,
+			    octave_idx_type k, octave_idx_type p,
+			    octave_idx_type &info, ComplexMatrix &eig_vec,
+			    ComplexColumnVector &eig_val, const Matrix& b,
+			    ColumnVector &permB, ColumnVector &resid,
+			    std::ostream &os, double tol = DBL_EPSILON,
+			    int rvec = 0, bool cholB = 0, int disp = 0,
+			    int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ,
+			    octave_idx_type k, octave_idx_type p,
+			    octave_idx_type &info, ComplexMatrix &eig_vec,
+			    ComplexColumnVector &eig_val,
+			    const SparseMatrix& b,
+			    ColumnVector &permB, ColumnVector &resid,
+			    std::ostream &os, double tol = DBL_EPSILON,
+			    int rvec = 0, bool cholB = 0, int disp = 0,
+			    int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma,
+				 octave_idx_type k, octave_idx_type p,
+				 octave_idx_type &info,
+				 ComplexMatrix &eig_vec,
+				 ComplexColumnVector &eig_val, const Matrix& b,
+				 ColumnVector &permB, ColumnVector &resid,
+				 std::ostream &os, double tol = DBL_EPSILON,
+				 int rvec = 0, bool cholB = 0, int disp = 0,
+				 int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma,
+				 octave_idx_type k, octave_idx_type p,
+				 octave_idx_type &info,
+				 ComplexMatrix &eig_vec,
+				 ComplexColumnVector &eig_val,
+				 const SparseMatrix& b,
+				 ColumnVector &permB, ColumnVector &resid,
+				 std::ostream &os, double tol = DBL_EPSILON,
+				 int rvec = 0, bool cholB = 0, int disp = 0,
+				 int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
+			  const std::string &_typ, double sigma,
+			  octave_idx_type k, octave_idx_type p,
+			  octave_idx_type &info, ComplexMatrix &eig_vec,
+			  ComplexColumnVector &eig_val,
+			  ColumnVector &resid, std::ostream& os,
+			  double tol = DBL_EPSILON, int rvec = 0,
+			  bool cholB = 0, int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ,
+			       octave_idx_type k, octave_idx_type p,
+			       octave_idx_type &info, ComplexMatrix &eig_vec,
+			       ComplexColumnVector &eig_val,
+			       const ComplexMatrix& b, ColumnVector &permB,
+			       ComplexColumnVector &resid,
+			       std::ostream &os, double tol = DBL_EPSILON,
+			       int rvec = 0, bool cholB = 0, int disp = 0,
+			       int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m,
+			       const std::string typ, octave_idx_type k,
+			       octave_idx_type p, octave_idx_type &info,
+			       ComplexMatrix &eig_vec,
+			       ComplexColumnVector &eig_val,
+			       const SparseComplexMatrix& b,
+			       ColumnVector &permB,
+			       ComplexColumnVector &resid,
+			       std::ostream &os, double tol = DBL_EPSILON,
+			       int rvec = 0, bool cholB = 0, int disp = 0,
+			       int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma,
+				    octave_idx_type k, octave_idx_type p,
+				    octave_idx_type &info,
+				    ComplexMatrix &eig_vec,
+				    ComplexColumnVector &eig_val,
+				    const ComplexMatrix& b,
+				    ColumnVector &permB,
+				    ComplexColumnVector &resid,
+				    std::ostream &os, double tol = DBL_EPSILON,
+				    int rvec = 0, bool cholB = 0,
+				    int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m,
+				    Complex sigma,
+				    octave_idx_type k, octave_idx_type p,
+				    octave_idx_type &info,
+				    ComplexMatrix &eig_vec,
+				    ComplexColumnVector &eig_val,
+				    const SparseComplexMatrix& b,
+				    ColumnVector &permB,
+				    ComplexColumnVector &resid,
+				    std::ostream &os, double tol = DBL_EPSILON,
+				    int rvec = 0, bool cholB = 0,
+				    int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
+			     const std::string &_typ, Complex sigma,
+			     octave_idx_type k, octave_idx_type p,
+			     octave_idx_type &info, ComplexMatrix &eig_vec,
+			     ComplexColumnVector &eig_val,
+			     ComplexColumnVector &resid, std::ostream& os,
+			     double tol = DBL_EPSILON, int rvec = 0,
+			     bool cholB = 0, int disp = 0, int maxit = 300);
+#endif
+
+#ifndef _MSC_VER
+template static octave_idx_type
+lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
+
+template static octave_idx_type
+lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
+	 ComplexMatrix&);
+
+template static octave_idx_type
+lusolve (const Matrix&, const Matrix&, Matrix&);
+
+template static octave_idx_type
+lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
+
+template static ComplexMatrix
+ltsolve (const SparseComplexMatrix&, const ColumnVector&,
+	 const ComplexMatrix&);
+
+template static Matrix
+ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
+
+template static ComplexMatrix
+ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+template static Matrix
+ltsolve (const Matrix&, const ColumnVector&, const Matrix&);
+
+template static ComplexMatrix
+utsolve (const SparseComplexMatrix&, const ColumnVector&,
+	 const ComplexMatrix&);
+
+template static Matrix
+utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
+
+template static ComplexMatrix
+utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+template static Matrix
+utsolve (const Matrix&, const ColumnVector&, const Matrix&);
+#endif
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- a/scripts/ChangeLog	Tue Dec 23 01:08:59 2008 +0100
+++ b/scripts/ChangeLog	Tue Dec 23 08:28:23 2008 +0100
@@ -1,3 +1,8 @@
+2008-12-23  David Bateman  <dbateman@free.fr>
+
+	* sparse/svds.m: New function.
+	* sparse/Makefile.in (SOURCES): Add it here.
+
 2008-11-21  Radek Salac  <salac.r@gmail.com>

 	* sparse/bicgstab.m: New function.
--- a/scripts/sparse/Makefile.in	Tue Dec 23 01:08:59 2008 +0100
+++ b/scripts/sparse/Makefile.in	Tue Dec 23 08:28:23 2008 +0100
@@ -35,7 +35,7 @@
 SOURCES = bicgstab.m cgs.m colperm.m etreeplot.m gplot.m nonzeros.m normest.m \
   pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \
   spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \
-  spvcat.m spy.m treelayout.m treeplot.m
+  spvcat.m spy.m svds.m treelayout.m treeplot.m

 DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/svds.m	Tue Dec 23 08:28:23 2008 +0100
@@ -0,0 +1,231 @@
+## Copyright (C) 2006 David Bateman
+##
+## This program 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 of the License, or
+## (at your option) any later version.
+##
+## This program 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 this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{s} =} svds (@var{a})
+## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k})
+## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma})
+## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma}, @var{opts})
+## @deftypefnx {Function File} {[@var{u}, @var{s}, @var{v}, @var{flag}] =} svds (@dots{})
+##
+## Find a few singular values of the matrix @var{a}. The singular values
+## are calculated using
+##
+## @example
+## @group
+## [@var{m}, @var{n}] = size(@var{a})
+## @var{s} = eigs([sparse(@var{m}, @var{m}), @var{a}; @dots{}
+##                 @var{a}', sparse(@var{n}, @var{n})])
+## @end group
+## @end example
+##
+## The eigenvalues returned by @code{eigs} correspond to the singular
+## values of @var{a}. The number of singular values to calculate is given
+## by @var{k}, whose default value is 6.
+##
+## The argument @var{sigma} can be used to specify which singular values
+## to find. @var{sigma} can be either the string 'L', the default, in
+## which case the largest singular values of @var{a} are found. Otherwise
+## @var{sigma} should be a real scalar, in which case the singular values
+## closest to @var{sigma} are found. Note that for relatively small values
+## of @var{sigma}, there is the chance that the requested number of singular
+## values are not returned. In that case @var{sigma} should be increased.
+##
+## If @var{opts} is given, then it is a structure that defines options
+## that @code{svds} will pass to @var{eigs}. The possible fields of this
+## structure are therefore determined by @code{eigs}. By default three
+## fields of this structure are set by @code{svds}.
+##
+## @table @code
+## @item tol
+## The required convergence tolerance for the singular values. @code{eigs}
+## is passed @var{tol} divided by @code{sqrt(2)}. The default value is
+## 1e-10.
+##
+## @item maxit
+## The maximum number of iterations. The defaut is 300.
+##
+## @item disp
+## The level of diagnostic printout. If @code{disp} is 0 then there is no
+## printout. The default value is 0.
+## @end table
+##
+## If more than one output argument is given, then @code{svds} also
+## calculates the left and right singular vectors of @var{a}. @var{flag}
+## is used to signal the convergence of @code{svds}. If @code{svds}
+## converges to the desired tolerance, then @var{flag} given by
+##
+## @example
+## @group
+## norm (@var{a} * @var{v} - @var{u} * @var{s}, 1) <= @dots{}
+##         @var{tol} * norm (@var{a}, 1)
+## @end group
+## @end example
+##
+## will be zero.
+## @end deftypefn
+## @seealso{eigs}
+
+function [u, s, v, flag] = svds (a, k, sigma, opts)
+
+  if (nargin < 1 || nargin > 4)
+    error ("Incorrect number of arguments");
+  endif
+
+  if (nargin < 4)
+    opts.tol = 1e-10 / sqrt(2);
+    opts.disp = 0;
+    opts.maxit = 300;
+  else
+    if (!isstruct(opts))
+      error("opts must be a structure");
+    endif
+    if (!isfield(opts,"tol"))
+      opts.tol = 1e-10 / sqrt(2);
+    endif
+  endif
+
+  if (nargin < 3 || strcmp(sigma,"L"))
+    if (isreal(a))
+      sigma = "LA";
+    else
+      sigma = "LR";
+    endif
+  elseif (isscalar(sigma) && isreal(sigma))
+    if ((sigma < 0))
+      error ("sigma must be a positive real value");
+    endif
+  else
+    error ("sigma must be a positive real value or the string 'L'");
+  endif
+
+  maxA = max(max(abs(a)));
+  if (maxA == 0)
+    u = eye(m, k);
+    s = zeros(k, k);
+    v = eye(n, k);
+  else
+    [m, n] = size(a);
+    if (nargin < 2)
+      k = min([6, m, n]);
+    else
+      k = min([k, m, n]);
+    endif
+
+    ## Scale everything by the 1-norm to make things more stable.
+    B = a / maxA;
+    Bopts = opts;
+    Bopts.tol = opts.tol / maxA;
+    Bsigma = sigma;
+    if (!ischar(Bsigma))
+      Bsigma = Bsigma / maxA;
+    endif
+
+    if (!ischar(Bsigma) && Bsigma == 0)
+      ## The eigenvalues returns by eigs are symmetric about 0. As we
+      ## are only interested in the positive eigenvalues, we have to
+      ## double k. If sigma is smaller than the smallest singular value
+      ## this can also be an issue. However, we'd like to avoid double
+      ## k for all scalar value of sigma...
+      [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)],
+			   2 * k, Bsigma, Bopts);
+    else
+      [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)],
+			   k, Bsigma, Bopts);
+    endif
+    s = diag(s);
+
+    if (ischar(sigma))
+      norma = max(s);
+    else
+      norma = normest(a);
+    endif
+    V = sqrt(2) * V;
+    u = V(1:m,:);
+    v = V(m+1:end,:);
+
+    ## We wish to exclude all eigenvalues that are less than zero as these
+    ## are artifacts of the way the matrix passed to eigs is formed. There
+    ## is also the possibility that the value of sigma chosen is exactly
+    ## a singular value, and in that case we're dead!! So have to rely on
+    ## the warning from eigs. We exclude the singular values which are
+    ## less than or equal to zero to within some tolerance scaled by the
+    ## norm since if we don't we might end up with too many singular
+    ## values. What is appropriate for the tolerance?
+    tol = norma * opts.tol;
+    ind = find(s > tol);
+    if (length(ind) < k)
+      ## Find the zero eigenvalues of B, Ignore the eigenvalues that are
+      ## nominally negative.
+      zind = find(abs(s) <= tol);
+      p = min(length(zind), k-length(ind));
+      ind = [ind;zind(1:p)];
+    elseif (length(ind) > k)
+      ind = ind(1:k);
+    endif
+    u = u(:,ind);
+    s = s(ind);
+    v = v(:,ind);
+
+    if (length(s) < k)
+      warning("returning fewer singular values than requested.");
+      if (!ischar(sigma))
+	warning("try increasing the value of sigma");
+      endif
+    endif
+
+    s = s * maxA;
+  endif
+
+  if (nargout < 2)
+    u = s;
+  else
+    s = diag(s);
+    if (nargout > 3)
+      flag = norm(a*v - u*s, 1) > sqrt(2) * opts.tol * norm(a, 1);
+    endif
+  endif
+endfunction
+
+%!shared n, k, a, u, s, v, opts
+%! n = 100;
+%! k = 7;
+%! a = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),0.4*n*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)]);
+%! [u,s,v] = svd(full(a));
+%! s = diag(s);
+%! [dum, idx] = sort(abs(s));
+%! s = s(idx);
+%! u = u(:,idx);
+%! v = v(:,idx);
+%! randn('state',42)
+%!test
+%! [u2,s2,v2,flag] = svds(a,k);
+%! s2 = diag(s2);
+%! assert(flag,!1);
+%! assert(s(end:-1:end-k+1), s2, 1e-10);
+%!test
+%! [u2,s2,v2,flag] = svds(a,k,0);
+%! s2 = diag(s2);
+%! assert(flag,!1);
+%! assert(s(k:-1:1), s2, 1e-10);
+%!test
+%! idx = floor(n/2);
+%! % Don't put sigma right on a singular value or there are convergence
+%! sigma = 0.99*s(idx) + 0.01*s(idx+1);
+%! [u2,s2,v2,flag] = svds(a,k,sigma);
+%! s2 = diag(s2);
+%! assert(flag,!1);
+%! assert(s((idx+floor(k/2)):-1:(idx-floor(k/2))), s2, 1e-10);
--- a/src/ChangeLog	Tue Dec 23 01:08:59 2008 +0100
+++ b/src/ChangeLog	Tue Dec 23 08:28:23 2008 +0100
@@ -1,3 +1,8 @@
+2008-12-23  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/eigs.cc: New file.
+	* Makefile.in (DLD_XSRC): Add it here.
+
 2008-12-22  David Bateman  <dbateman@free.fr>

 	* DLD-FUNCTIONS/__voronoi__.cc (F__voronoi__): Resize AtInf array
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/eigs.cc	Tue Dec 23 08:28:23 2008 +0100
@@ -0,0 +1,1465 @@
+/*
+
+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 3 of the License, 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, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "ov.h"
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "quit.h"
+#include "variables.h"
+#include "ov-re-sparse.h"
+#include "ov-cx-sparse.h"
+#include "oct-map.h"
+#include "pager.h"
+#include "unwind-prot.h"
+
+#include "eigs-base.cc"
+
+// Global pointer for user defined function.
+static octave_function *eigs_fcn = 0;
+
+// Have we warned about imaginary values returned from user function?
+static bool warned_imaginary = false;
+
+// Is this a recursive call?
+static int call_depth = 0;
+
+ColumnVector
+eigs_func (const ColumnVector &x, int &eigs_error)
+{
+  ColumnVector retval;
+  octave_value_list args;
+  args(0) = x;
+
+  if (eigs_fcn)
+    {
+      octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  eigs_error = 1;
+	  gripe_user_supplied_eval ("eigs");
+	  return retval;
+	}
+
+      if (tmp.length () && tmp(0).is_defined ())
+	{
+	  if (! warned_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("eigs: ignoring imaginary part returned from user-supplied function");
+	      warned_imaginary = true;
+	    }
+
+	  retval = ColumnVector (tmp(0).vector_value ());
+
+	  if (error_state)
+	    {
+	      eigs_error = 1;
+	      gripe_user_supplied_eval ("eigs");
+	    }
+	}
+      else
+	{
+	  eigs_error = 1;
+	  gripe_user_supplied_eval ("eigs");
+	}
+    }
+
+  return retval;
+}
+
+ComplexColumnVector
+eigs_complex_func (const ComplexColumnVector &x, int &eigs_error)
+{
+  ComplexColumnVector retval;
+  octave_value_list args;
+  args(0) = x;
+
+  if (eigs_fcn)
+    {
+      octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  eigs_error = 1;
+	  gripe_user_supplied_eval ("eigs");
+	  return retval;
+	}
+
+      if (tmp.length () && tmp(0).is_defined ())
+	{
+	  retval = ComplexColumnVector (tmp(0).complex_vector_value ());
+
+	  if (error_state)
+	    {
+	      eigs_error = 1;
+	      gripe_user_supplied_eval ("eigs");
+	    }
+	}
+      else
+	{
+	  eigs_error = 1;
+	  gripe_user_supplied_eval ("eigs");
+	}
+    }
+
+  return retval;
+}
+
+DEFUN_DLD (eigs, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{d}} = eigs (@var{a})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k}, @var{sigma})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k}, @var{sigma},@var{opts})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k}, @var{sigma})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k}, @var{sigma}, @var{opts})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k}, @var{sigma})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k}, @var{sigma})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts})\n\
+@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k}, @var{sigma}, @var{opts})\n\
+@deftypefnx {Loadable Function} {[@var{v}, @var{d}]} = eigs (@var{a}, @dots{})\n\
+@deftypefnx {Loadable Function} {[@var{v}, @var{d}]} = eigs (@var{af}, @var{n}, @dots{})\n\
+@deftypefnx {Loadable Function} {[@var{v}, @var{d}, @var{flag}]} = eigs (@var{a}, @dots{})\n\
+@deftypefnx {Loadable Function} {[@var{v}, @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 problems).\n\
+\n\
+@item 'sa'\n\
+Smallest Algebraic (valid only for real symmetric problems).\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 problems).\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 @var{n}-by-@var{k} 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\
+@var{n}-by-@var{k} 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 ();
+  std::string fcn_name;
+  octave_idx_type n = 0;
+  octave_idx_type k = 6;
+  Complex sigma = 0.;
+  double sigmar, sigmai;
+  bool have_sigma = false;
+  std::string typ = "LM";
+  Matrix amm, bmm, bmt;
+  ComplexMatrix acm, bcm, bct;
+  SparseMatrix asmm, bsmm, bsmt;
+  SparseComplexMatrix ascm, bscm, bsct;
+  int b_arg = 0;
+  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;
+  bool a_is_sparse = 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';
+
+  warned_imaginary = false;
+
+  unwind_protect::begin_frame ("Feigs");
+
+  unwind_protect_int (call_depth);
+  call_depth++;
+
+  if (call_depth > 1)
+    {
+      error ("eigs: invalid recursive call");
+      if (fcn_name.length())
+	clear_function (fcn_name);
+      unwind_protect::run_frame ("Feigs");
+      return retval;
+    }
+
+  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 = ");
+	  eigs_fcn = extract_function (args(0), "eigs", fcn_name, fname,
+				       "; endfunction");
+	}
+      else
+	eigs_fcn = args(0).function_value ();
+
+      if (!eigs_fcn)
+	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
+    {
+      if (args(0).is_complex_type ())
+	{
+	  if (args(0).is_sparse_type ())
+	    {
+	      ascm = (args(0).sparse_complex_matrix_value());
+	      a_is_sparse = true;
+	    }
+	  else
+	    acm = (args(0).complex_matrix_value());
+	  a_is_complex = true;
+	  symmetric = false; // ARAPACK doesn't special case complex symmetric
+	}
+      else
+	{
+	  if (args(0).is_sparse_type ())
+	    {
+	      asmm = (args(0).sparse_matrix_value());
+	      a_is_sparse = true;
+	      symmetric = asmm.is_symmetric();
+	    }
+	  else
+	    {
+	      amm = (args(0).matrix_value());
+	      symmetric = amm.is_symmetric();
+	    }
+	}
+
+      // 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 &&
+	  !(args(1).is_real_scalar ()))
+	if (args(1).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 upper 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)
+	{
+	  if (a_is_sparse)
+	    bscm = args(b_arg).sparse_complex_matrix_value ();
+	  else
+	    bcm = args(b_arg).complex_matrix_value ();
+	}
+      else
+	{
+	  if (a_is_sparse)
+	    bsmm = args(b_arg).sparse_matrix_value ();
+	  else
+	    bmm = args(b_arg).matrix_value ();
+	}
+    }
+
+  // Mode 1 for SM mode seems unstable for some reason.
+  // Use Mode 3 instead, with sigma = 0.
+  if (!error_state && !have_sigma && typ == "SM")
+    have_sigma = true;
+
+  if (!error_state)
+    {
+      octave_idx_type nconv;
+      if (a_is_complex || b_is_complex)
+	{
+	  ComplexMatrix eig_vec;
+	  ComplexColumnVector eig_val;
+
+
+	  if (have_a_fun)
+	    nconv = EigsComplexNonSymmetricFunc
+	      (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, eig_val,
+	       cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+	  else if (have_sigma)
+	    if (a_is_sparse)
+	      nconv = EigsComplexNonSymmetricMatrixShift
+		(ascm, sigma, k, p, info, eig_vec, eig_val, bscm, permB,
+		 cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		 maxit);
+	    else
+	      nconv = EigsComplexNonSymmetricMatrixShift
+		(acm, sigma, k, p, info, eig_vec, eig_val, bcm, permB, cresid,
+		 octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+	  else
+	    if (a_is_sparse)
+	      nconv = EigsComplexNonSymmetricMatrix
+		(ascm, typ, k, p, info, eig_vec, eig_val, bscm, permB, cresid,
+		 octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+	    else
+	      nconv = EigsComplexNonSymmetricMatrix
+		(acm, typ, k, p, info, eig_vec, eig_val, bcm, permB, cresid,
+		 octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+
+	  if (nargout < 2)
+	    retval (0) = eig_val;
+	  else
+	    {
+	      retval(2) = double (info);
+	      retval(1) = ComplexDiagMatrix (eig_val);
+	      retval(0) = eig_vec;
+	    }
+	}
+      else if (sigmai != 0.)
+	{
+	  // Promote real problem to a complex one.
+	  ComplexMatrix eig_vec;
+	  ComplexColumnVector eig_val;
+
+	  if (have_a_fun)
+	    nconv = EigsComplexNonSymmetricFunc
+	      (eigs_complex_func, n, typ,  sigma, k, p, info, eig_vec, eig_val,
+	       cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+	  else
+	    if (a_is_sparse)
+	      nconv = EigsComplexNonSymmetricMatrixShift
+		(SparseComplexMatrix (asmm), sigma, k, p, info, eig_vec,
+		 eig_val, SparseComplexMatrix (bsmm), permB, cresid,
+		 octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+	    else
+	      nconv = EigsComplexNonSymmetricMatrixShift
+		(ComplexMatrix (amm), sigma, k, p, info, eig_vec,
+		 eig_val, ComplexMatrix (bmm), permB, cresid,
+		 octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+
+	  if (nargout < 2)
+	    retval (0) = eig_val;
+	  else
+	    {
+	      retval(2) = double (info);
+	      retval(1) = ComplexDiagMatrix (eig_val);
+	      retval(0) = eig_vec;
+	    }
+	}
+      else
+	{
+	  if (symmetric)
+	    {
+	      Matrix eig_vec;
+	      ColumnVector eig_val;
+
+	      if (have_a_fun)
+		nconv = EigsRealSymmetricFunc
+		  (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val,
+		   resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		   maxit);
+	      else if (have_sigma)
+		if (a_is_sparse)
+		  nconv = EigsRealSymmetricMatrixShift
+		    (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+		else
+		  nconv = EigsRealSymmetricMatrixShift
+		    (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+	      else
+		if (a_is_sparse)
+		  nconv = EigsRealSymmetricMatrix
+		    (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+		else
+		  nconv = EigsRealSymmetricMatrix
+		    (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+
+	      if (nargout < 2)
+		retval (0) = eig_val;
+	      else
+		{
+		  retval(2) = double (info);
+		  retval(1) = DiagMatrix (eig_val);
+		  retval(0) = eig_vec;
+		}
+	    }
+	  else
+	    {
+	      ComplexMatrix eig_vec;
+	      ComplexColumnVector eig_val;
+
+	      if (have_a_fun)
+		nconv = EigsRealNonSymmetricFunc
+		  (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val,
+		   resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		   maxit);
+	      else if (have_sigma)
+		if (a_is_sparse)
+		  nconv = EigsRealNonSymmetricMatrixShift
+		    (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+		else
+		  nconv = EigsRealNonSymmetricMatrixShift
+		    (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+	      else
+		if (a_is_sparse)
+		  nconv = EigsRealNonSymmetricMatrix
+		    (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+		else
+		  nconv = EigsRealNonSymmetricMatrix
+		    (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
+		     resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+		     maxit);
+
+	      if (nargout < 2)
+		retval (0) = eig_val;
+	      else
+		{
+		  retval(2) = double (info);
+		  retval(1) = ComplexDiagMatrix (eig_val);
+		  retval(0) = eig_vec;
+		}
+	    }
+	}
+
+      if (nconv <= 0)
+	warning ("eigs: None of the %d requested eigenvalues converged", k);
+      else if (nconv < k)
+	warning ("eigs: Only %d of the %d requested eigenvalues converged",
+		 nconv, k);
+    }
+
+  if (! fcn_name.empty ())
+    clear_function (fcn_name);
+
+  unwind_protect::run_frame ("Feigs");
+#else
+  error ("eigs: not available in this version of Octave");
+#endif
+
+  return retval;
+}
+
+/* #### SPARSE MATRIX VERSIONS #### */
+
+/*
+
+%% 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);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.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
+%! fn = @(x) A \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (d1, d0(k:-1:1), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (d1, eigs(A,k,4.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
+%!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],[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);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.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
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.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);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.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
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.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
+
+*/
+
+/* #### FULL MATRIX VERSIONS #### */
+
+/*
+
+%% Real positive definite tests, n must be even
+%!shared n, k, A, d0, d2
+%! n = 20;
+%! k = 4;
+%! A = full(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, eye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!assert (eigs(A,k,4.1), eigs(A,eye(n),k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!assert (eigs(A,k,4.1), eigs(A,eye(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
+%! fn = @(x) A \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (d1, d0(k:-1:1), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (d1, eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%!  assert(max(abs((A - d1(i)*eye(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)*eye(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)*eye(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)*eye(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)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+
+%% Real unsymmetric tests
+%!shared n, k, A, d0
+%! n = 20;
+%! k = 4;
+%! A =  full(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, eye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12);
+%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(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
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%!  assert(max(abs((A - d1(i)*eye(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)*eye(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)*eye(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)*eye(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)*eye(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)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+
+%% Complex hermitian tests
+%!shared n, k, A, d0
+%! n = 20;
+%! k = 4;
+%! A = full(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, eye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
+%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12);
+%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(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
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%!  assert(max(abs((A - d1(i)*eye(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)*eye(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)*eye(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)*eye(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)*eye(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)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- a/src/Makefile.in	Tue Dec 23 01:08:59 2008 +0100
+++ b/src/Makefile.in	Tue Dec 23 08:28:23 2008 +0100
@@ -300,7 +300,7 @@
   -L../libcruft $(LIBCRUFT) -L../liboctave $(LIBOCTAVE) \
   -L. $(LIBOCTINTERP) $(CHOLMOD_LIBS) $(UMFPACK_LIBS) $(AMD_LIBS) \
    $(CAMD_LIBS) $(COLAMD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \
-   $(FFTW_LIBS) $(LIBS) $(FLIBS)
+   $(FFTW_LIBS) $(ARPACK_LIBS) $(LIBS) $(FLIBS)

 BUILT_DISTFILES = DOCSTRINGS oct-gperf.h parse.cc lex.cc y.tab.h \
 	$(OPT_HANDLERS) $(BUILT_EXTRAS)
@@ -371,7 +371,7 @@
 	$(OCTAVE_LIBS) \
 	$(LEXLIB) $(UMFPACK_LIBS) $(AMD_LIBS) $(CAMD_LIBS) $(COLAMD_LIBS) \
 	$(CHOLMOD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \
-	$(FFTW_LIBS) $(OPENGL_LIBS) $(LIBS) $(FLIBS)
+	$(FFTW_LIBS) $(ARPACK_LIBS) $(OPENGL_LIBS) $(LIBS) $(FLIBS)

 stmp-pic: pic
 	@if [ -f stmp-pic ]; then \
@@ -650,6 +650,7 @@
 convhulln.oct: OCT_LINK_DEPS += $(QHULL_LIBS)
 __delaunayn__.oct: OCT_LINK_DEPS += $(QHULL_LIBS)
 __voronoi__.oct: OCT_LINK_DEPS += $(QHULL_LIBS)
+eigs.oct: OCT_LINK_DEPS += $(ARPACK_LIBS)
 regexp.oct: OCT_LINK_DEPS += $(REGEX_LIBS)
 urlwrite.oct: OCT_LINK_DEPS += $(CURL_LIBS)
 __glpk__.oct: OCT_LINK_DEPS += $(GLPK_LIBS)