/*

Copyright (C) 2005, 2008, 2009 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: ***
*/