Logo Search packages:      
Sourcecode: octave3.2 version File versions

eigs-base.cc

/*

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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }
  
  if (k < 1 || k > n - 2)
    {
      (*current_liboctave_error_handler) 
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n"
       "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler) 
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler) 
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n"
           "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }
  
  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler) 
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }
  
  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler)
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
           "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler)
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2 + 1;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }

  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler) 
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
       "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler) 
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2 + 1;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }

  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler) 
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
           "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler) 
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2 + 1;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }

  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler)
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
           "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler)
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2 + 1;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }

  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler) 
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
       "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler) 
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2 + 1;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }

  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler) 
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
           "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler) 
      ("eigs: opts.p must be greater than k and less than n");
      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 (n < 3)
    {
      (*current_liboctave_error_handler)
      ("eigs: n must be at least 3");
      return -1;
    }

  if (p < 0)
    {
      p = k * 2 + 1;

      if (p < 20)
      p = 20;
      
      if (p > n - 1)
      p = n - 1 ;
    }

  if (k <= 0 || k >= n - 1)
    {
      (*current_liboctave_error_handler)
      ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
           "      Use 'eig(full(A))' instead");
      return -1;
    }

  if (p <= k || p >= n)
    {
      (*current_liboctave_error_handler)
      ("eigs: opts.p must be greater than k and less than n");
      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: ***
*/

Generated by  Doxygen 1.6.0   Back to index