--- a/ChangeLog
+++ b/ChangeLog
@@ -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
+++ b/Makeconf.in
@@ -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
+++ b/NEWS
@@ -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
+++ b/configure.in
@@ -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
+++ b/doc/ChangeLog
@@ -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
+++ b/doc/interpreter/sparse.txi
@@ -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
+++ b/liboctave/ChangeLog
@@ -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
+++ b/liboctave/Makefile.in
@@ -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

new file mode 100755
--- /dev/null
+++ b/liboctave/eigs-base.cc
@@ -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: ***
+*/