[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: eigs and ARPACK
|
From: |
David Bateman |
|
Subject: |
Re: eigs and ARPACK |
|
Date: |
Fri, 05 Dec 2008 00:09:42 +0100 |
|
User-agent: |
Mozilla-Thunderbird 2.0.0.17 (X11/20081018) |
John W. Eaton wrote:
On 3-Dec-2008, Søren Hauberg wrote:
| It was just pointed out on the Octave-Forge list that ARPACK seems to
| have been released under a standard BSD license [1]. Doesn't this mean
| that we can move the code from the 'arpack' package into core Octave,
| which would give us the 'eigs' and 'svds' functions?
|
| Søren
|
| [1] http://www.caam.rice.edu/software/ARPACK/RiceBSD.txt
Yes, this is great news. Would someone like to prepare a patch?
jwe
Well attached is an untested port of the octave-forge code. I don't know
if I'll have the time to test this code before a few days and so I'm
sending it even if it is untested. It also doesn't support single
precision matrices and probably should..
Regards
David
--
David Bateman address@hidden
35 rue Gambetta +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE +33 6 72 01 06 33 (Mob)
# HG changeset patch
# User David Bateman <address@hidden>
# Date 1228431849 -3600
# Node ID b066459427388ea53edb1531e8e4bde36fa2cb7e
# Parent b5f0f1c09661795d26a76e54ec77ce1709721ddf
Add eigs and svds functions
diff --git a/Makeconf.in b/Makeconf.in
--- a/Makeconf.in
+++ b/Makeconf.in
@@ -237,6 +237,7 @@
CHOLMOD_LIBS = @CHOLMOD_LIBS@
CXSPARSE_LIBS = @CXSPARSE_LIBS@
OPENGL_LIBS = @OPENGL_LIBS@
+ARPACK_LIBS = @ARPACK_LIBS@
LIBS = @LIBS@
USE_64_BIT_IDX_T = @USE_64_BIT_IDX_T@
diff --git a/configure.in b/configure.in
--- a/configure.in
+++ b/configure.in
@@ -1042,6 +1042,27 @@
fi
if test -n "$warn_cxsparse"; then
AC_MSG_WARN($warn_cxsparse)
+fi
+
+ARPACK_LIBS=
+AC_SUBST(ARPACK_LIBS)
+
+AC_ARG_WITH(arpack,
+ [AS_HELP_STRING([--without-arpack],
+ [don't use ARPACK, disable some sparse functionality])],
+ with_arpack=$withval, with_arpack=yes)
+
+warn_arpack="arpack not found. This will result in a lack of the eigs
function."
+if test "$with_arpack" = yes; then
+ with_arpack=no
+ AC_CHECK_LIB(arpack, F77_FUNC(dseupd,DSEUPD), [ARPACK_LIBS="-larpack";
with_arpack=yes])
+ if test "$with_arpack" = yes; then
+ AC_DEFINE(HAVE_ARPACK, 1, [Define if the ARPACK library is used.])
+ warn_cxsparse=
+ fi
+fi
+if test -n "$warn_arpack"; then
+ AC_MSG_WARN($warn_arpack)
fi
### Handle shared library options.
@@ -2038,6 +2059,7 @@
CCOLAMD libraries: $CCOLAMD_LIBS
CHOLMOD libraries: $CHOLMOD_LIBS
CXSPARSE libraries: $CXSPARSE_LIBS
+ ARPACK libraries: $ARPACK_LIBS
HDF5 libraries: $HDF5_LIBS
CURL libraries: $CURL_LIBS
REGEX libraries: $REGEX_LIBS
@@ -2146,6 +2168,11 @@
if test -n "$warn_cxsparse"; then
AC_MSG_WARN($warn_cxsparse)
+ warn_msg_printed=true
+fi
+
+if test -n "$warn_arpack"; then
+ AC_MSG_WARN($warn_arpack)
warn_msg_printed=true
fi
diff --git a/doc/interpreter/sparse.txi b/doc/interpreter/sparse.txi
--- a/doc/interpreter/sparse.txi
+++ b/doc/interpreter/sparse.txi
@@ -473,9 +473,8 @@
@dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, @dfn{symrcm}
@item Linear algebra:
- @dfn{matrix_type}, @dfn{normest}, @dfn{condest}, @dfn{sprank}
- @dfn{spaugment}
address@hidden @dfn{eigs}, @dfn{svds} but these are in octave-forge for now
+ @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, @dfn{normest}, @dfn{sprank},
+ @dfn{spaugment}, @dfn{svds}
@item Iterative techniques:
@dfn{luinc}, @dfn{pcg}, @dfn{pcr}
@@ -831,6 +830,15 @@
function to find a least squares solution to a linear equation.
@DOCSTRING(spaugment)
+
+Finally, the function @code{eigs} can be used to calculate a limited
+number of eigenvalues and eigenvectors based on a selection criteria
+and likewise for @code{svds} which calculates a limited number of
+singular values and vectors.
+
address@hidden(eigs)
+
address@hidden(svds)
@node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse
Matrices
@section Iterative Techniques applied to sparse matrices
diff --git a/liboctave/Makefile.in b/liboctave/Makefile.in
--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -105,8 +105,8 @@
TEMPLATE_SRC := Array.cc ArrayN.cc DiagArray2.cc \
MArray.cc MArray2.cc MArrayN.cc MDiagArray2.cc \
- base-lu.cc oct-sort.cc sparse-base-lu.cc sparse-base-chol.cc \
- sparse-dmsolve.cc
+ base-lu.cc eigs-base.cc oct-sort.cc sparse-base-lu.cc \
+ sparse-base-chol.cc sparse-dmsolve.cc
TI_SRC := Array-C.cc Array-b.cc Array-ch.cc Array-i.cc Array-d.cc \
Array-f.cc Array-fC.cc Array-s.cc Array-str.cc \
diff --git a/liboctave/eigs-base.cc b/liboctave/eigs-base.cc
new file mode 100755
--- /dev/null
+++ b/liboctave/eigs-base.cc
@@ -0,0 +1,3762 @@
+/*
+
+Copyright (C) 2005 David Bateman
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING. If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cfloat>
+#include <cmath>
+#include <vector>
+#include <iostream>
+
+#include "f77-fcn.h"
+#include "quit.h"
+#include "SparsedbleLU.h"
+#include "SparseCmplxLU.h"
+#include "dSparse.h"
+#include "CSparse.h"
+#include "MatrixType.h"
+#include "SparsedbleCHOL.h"
+#include "SparseCmplxCHOL.h"
+#include "oct-rand.h"
+#include "dbleCHOL.h"
+#include "CmplxCHOL.h"
+#include "dbleLU.h"
+#include "CmplxLU.h"
+
+#ifdef HAVE_ARPACK
+typedef ColumnVector (*EigsFunc) (const ColumnVector &x, int &eigs_error);
+typedef ComplexColumnVector (*EigsComplexFunc)
+ (const ComplexColumnVector &x, int &eigs_error);
+
+// Arpack and blas fortran functions we call.
+extern "C"
+{
+ F77_RET_T
+ F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, const double&,
+ double*, const octave_idx_type&, double*,
+ const octave_idx_type&, octave_idx_type*,
+ octave_idx_type*, double*, double*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (dseupd, DSEUPD) (const int&, F77_CONST_CHAR_ARG_DECL,
+ int*, double*, double*,
+ const octave_idx_type&, const double&,
+ F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ const double&, double*, const octave_idx_type&,
+ double*, const octave_idx_type&, octave_idx_type*,
+ octave_idx_type*, double*, double*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (dnaupd, DNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ octave_idx_type&, const double&,
+ double*, const octave_idx_type&, double*,
+ const octave_idx_type&, octave_idx_type*,
+ octave_idx_type*, double*, double*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (dneupd, DNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL,
+ int*, double*, double*,
+ double*, const octave_idx_type&, const double&,
+ const double&, double*, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ octave_idx_type&, const double&, double*,
+ const octave_idx_type&, double*,
+ const octave_idx_type&, octave_idx_type*,
+ octave_idx_type*, double*, double*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (znaupd, ZNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, const double&,
+ Complex*, const octave_idx_type&, Complex*,
+ const octave_idx_type&, octave_idx_type*,
+ octave_idx_type*, Complex*, Complex*,
+ const octave_idx_type&, double *, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (zneupd, ZNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL,
+ int*, Complex*, Complex*,
+ const octave_idx_type&, const Complex&,
+ Complex*, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, const double&,
+ Complex*, const octave_idx_type&, Complex*,
+ const octave_idx_type&, octave_idx_type*,
+ octave_idx_type*, Complex*, Complex*,
+ const octave_idx_type&, double *, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, const octave_idx_type&,
const double&,
+ const double*, const octave_idx_type&, const double*,
+ const octave_idx_type&, const double&, double*,
+ const octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL);
+
+
+ F77_RET_T
+ F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
+ const octave_idx_type&, const octave_idx_type&,
const Complex&,
+ const Complex*, const octave_idx_type&, const
Complex*,
+ const octave_idx_type&, const Complex&, Complex*,
const octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL);
+
+}
+
+
+#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
+static octave_idx_type
+lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
+
+static octave_idx_type
+lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
+ ComplexMatrix&);
+
+static octave_idx_type
+lusolve (const Matrix&, const Matrix&, Matrix&);
+
+static octave_idx_type
+lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
+
+static ComplexMatrix
+ltsolve (const SparseComplexMatrix&, const ColumnVector&,
+ const ComplexMatrix&);
+
+static Matrix
+ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,);
+
+static ComplexMatrix
+ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+static Matrix
+ltsolve (const Matrix&, const ColumnVector&, const Matrix&,);
+
+static ComplexMatrix
+utsolve (const SparseComplexMatrix&, const ColumnVector&, const
ComplexMatrix&);
+
+static Matrix
+utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
+
+static ComplexMatrix
+utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+static Matrix
+utsolve (const Matrix&, const ColumnVector&, const Matrix&);
+
+#endif
+
+template <class M, class SM>
+static octave_idx_type
+lusolve (const SM& L, const SM& U, M& m)
+{
+ octave_idx_type err = 0;
+ double rcond;
+ MatrixType utyp (MatrixType::Upper);
+
+ // Sparse L is lower triangular, Dense L is permuted lower triangular!!!
+ m = L.solve (m, err, rcond, 0);
+ if (err)
+ return err;
+
+ m = U.solve (utyp, m, err, rcond, 0);
+
+ return err;
+}
+
+template <class SM, class M>
+static M
+ltsolve (const SM& L, const ColumnVector& Q, const M& m)
+{
+ octave_idx_type n = L.cols();
+ octave_idx_type b_nc = m.cols();
+ octave_idx_type err = 0;
+ double rcond;
+ MatrixType ltyp (MatrixType::Lower);
+ M tmp = L.solve (ltyp, m, err, rcond, 0);
+ M retval;
+ const double* qv = Q.fortran_vec();
+
+ if (!err)
+ {
+ retval.resize (n, b_nc);
+ for (octave_idx_type j = 0; j < b_nc; j++)
+ {
+ for (octave_idx_type i = 0; i < n; i++)
+ retval.elem(static_cast<octave_idx_type>(qv[i]), j) =
+ tmp.elem(i,j);
+ }
+ }
+
+ return retval;
+}
+
+template <class SM, class M>
+static M
+utsolve (const SM& U, const ColumnVector& Q, const M& m)
+{
+ octave_idx_type n = U.cols();
+ octave_idx_type b_nc = m.cols();
+ octave_idx_type err = 0;
+ double rcond;
+ MatrixType utyp (MatrixType::Upper);
+
+ M retval (n, b_nc);
+ const double* qv = Q.fortran_vec();
+ for (octave_idx_type j = 0; j < b_nc; j++)
+ {
+ for (octave_idx_type i = 0; i < n; i++)
+ retval.elem(i,j) = m.elem(static_cast<octave_idx_type>(qv[i]), j);
+ }
+ return U.solve (utyp, retval, err, rcond, 0);
+}
+
+static bool
+vector_product (const SparseMatrix& m, const double* x, double* y)
+{
+ octave_idx_type nc = m.cols ();
+
+ for (octave_idx_type j = 0; j < nc; j++)
+ y[j] = 0.;
+
+ for (octave_idx_type j = 0; j < nc; j++)
+ for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
+ y[m.ridx(i)] += m.data(i) * x[j];
+
+ return true;
+}
+
+static bool
+vector_product (const Matrix& m, const double *x, double *y)
+{
+ octave_idx_type nr = m.rows ();
+ octave_idx_type nc = m.cols ();
+
+ F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
+ nr, nc, 1.0, m.data (), nr,
+ x, 1, 0.0, y, 1
+ F77_CHAR_ARG_LEN (1)));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable error in dgemv");
+ return false;
+ }
+ else
+ return true;
+}
+
+static bool
+vector_product (const SparseComplexMatrix& m, const Complex* x,
+ Complex* y)
+{
+ octave_idx_type nc = m.cols ();
+
+ for (octave_idx_type j = 0; j < nc; j++)
+ y[j] = 0.;
+
+ for (octave_idx_type j = 0; j < nc; j++)
+ for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
+ y[m.ridx(i)] += m.data(i) * x[j];
+
+ return true;
+}
+
+static bool
+vector_product (const ComplexMatrix& m, const Complex *x, Complex *y)
+{
+ octave_idx_type nr = m.rows ();
+ octave_idx_type nc = m.cols ();
+
+ F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
+ nr, nc, 1.0, m.data (), nr,
+ x, 1, 0.0, y, 1
+ F77_CHAR_ARG_LEN (1)));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable error in zgemv");
+ return false;
+ }
+ else
+ return true;
+}
+
+static bool
+make_cholb (Matrix& b, Matrix& bt, ColumnVector& permB)
+{
+ octave_idx_type info;
+ CHOL fact (b, info);
+ octave_idx_type n = b.cols();
+
+ if (info != 0)
+ return false;
+ else
+ {
+ bt = fact.chol_matrix ();
+ b = bt.transpose();
+ permB = ColumnVector(n);
+ for (octave_idx_type i = 0; i < n; i++)
+ permB(i) = i;
+ return true;
+ }
+}
+
+static bool
+make_cholb (SparseMatrix& b, SparseMatrix& bt, ColumnVector& permB)
+{
+ octave_idx_type info;
+ SparseCHOL fact (b, info, false);
+
+ if (fact.P() != 0)
+ return false;
+ else
+ {
+ b = fact.L();
+ bt = b.transpose();
+ permB = fact.perm() - 1.0;
+ return true;
+ }
+}
+
+static bool
+make_cholb (ComplexMatrix& b, ComplexMatrix& bt, ColumnVector& permB)
+{
+ octave_idx_type info;
+ ComplexCHOL fact (b, info);
+ octave_idx_type n = b.cols();
+
+ if (info != 0)
+ return false;
+ else
+ {
+ bt = fact.chol_matrix ();
+ b = bt.hermitian();
+ permB = ColumnVector(n);
+ for (octave_idx_type i = 0; i < n; i++)
+ permB(i) = i;
+ return true;
+ }
+}
+
+static bool
+make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt,
+ ColumnVector& permB)
+{
+ octave_idx_type info;
+ SparseComplexCHOL fact (b, info, false);
+
+ if (fact.P() != 0)
+ return false;
+ else
+ {
+ b = fact.L();
+ bt = b.hermitian();
+ permB = fact.perm() - 1.0;
+ return true;
+ }
+}
+
+static bool
+LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b,
+ bool cholB, const ColumnVector& permB, double sigma,
+ SparseMatrix &L, SparseMatrix &U, octave_idx_type *P,
+ octave_idx_type *Q)
+{
+ bool have_b = ! b.is_empty ();
+ octave_idx_type n = m.rows();
+
+ // Caclulate LU decomposition of 'A - sigma * B'
+ SparseMatrix AminusSigmaB (m);
+
+ if (have_b)
+ {
+ if (cholB)
+ {
+ if (permB.length())
+ {
+ SparseMatrix tmp(n,n,n);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ tmp.xcidx(i) = i;
+ tmp.xridx(i) =
+ static_cast<octave_idx_type>(permB(i));
+ tmp.xdata(i) = 1;
+ }
+ tmp.xcidx(n) = n;
+
+ AminusSigmaB = AminusSigmaB - sigma * tmp *
+ b.transpose() * b * tmp.transpose();
+ }
+ else
+ AminusSigmaB = AminusSigmaB - sigma *
+ b.transpose() * b;
+ }
+ else
+ AminusSigmaB = AminusSigmaB - sigma * b;
+ }
+ else
+ {
+ SparseMatrix sigmat (n, n, n);
+
+ // Create sigma * speye(n,n)
+ sigmat.xcidx (0) = 0;
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ sigmat.xdata(i) = sigma;
+ sigmat.xridx(i) = i;
+ sigmat.xcidx(i+1) = i + 1;
+ }
+
+ AminusSigmaB = AminusSigmaB - sigmat;
+ }
+
+ SparseLU fact (AminusSigmaB);
+
+ L = fact.L ();
+ U = fact.U ();
+ const octave_idx_type *P2 = fact.row_perm ();
+ const octave_idx_type *Q2 = fact.col_perm ();
+
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ P[j] = P2[j];
+ Q[j] = Q2[j];
+ }
+
+ // Test condition number of LU decomposition
+ double minU = octave_NaN;
+ double maxU = octave_NaN;
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ double d = 0.;
+ if (U.xcidx(j+1) > U.xcidx(j) &&
+ U.xridx (U.xcidx(j+1)-1) == j)
+ d = std::abs (U.xdata (U.xcidx(j+1)-1));
+
+ if (xisnan (minU) || d < minU)
+ minU = d;
+
+ if (xisnan (maxU) || d > maxU)
+ maxU = d;
+ }
+
+ double rcond = (minU / maxU);
+ volatile double rcond_plus_one = rcond + 1.0;
+
+ if (rcond_plus_one == 1.0 || xisnan (rcond))
+ {
+ (*current_liboctave_warning_handler)
+ ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+ (*current_liboctave_warning_handler)
+ (" an eigenvalue. Convergence is not guaranteed");
+ }
+
+ return true;
+}
+
+static bool
+LuAminusSigmaB (const Matrix &m, const Matrix &b,
+ bool cholB, const ColumnVector& permB, double sigma,
+ Matrix &L, Matrix &U, octave_idx_type *P,
+ octave_idx_type *Q)
+{
+ bool have_b = ! b.is_empty ();
+ octave_idx_type n = m.cols();
+
+ // Caclulate LU decomposition of 'A - sigma * B'
+ Matrix AminusSigmaB (m);
+
+ if (have_b)
+ {
+ if (cholB)
+ {
+ Matrix tmp = sigma * b.transpose() * b;
+ const double *pB = permB.fortran_vec();
+ double *p = AminusSigmaB.fortran_vec();
+
+ if (permB.length())
+ {
+ for (octave_idx_type j = 0;
+ j < b.cols(); j++)
+ for (octave_idx_type i = 0;
+ i < b.rows(); i++)
+ *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]),
+ static_cast<octave_idx_type>(pB[j]));
+ }
+ else
+ AminusSigmaB = AminusSigmaB - tmp;
+ }
+ else
+ AminusSigmaB = AminusSigmaB - sigma * b;
+ }
+ else
+ {
+ double *p = AminusSigmaB.fortran_vec();
+
+ for (octave_idx_type i = 0; i < n; i++)
+ p[i*(n+1)] -= sigma;
+ }
+
+ LU fact (AminusSigmaB);
+
+ L = fact.P().transpose() * fact.L ();
+ U = fact.U ();
+ for (octave_idx_type j = 0; j < n; j++)
+ P[j] = Q[j] = j;
+
+ // Test condition number of LU decomposition
+ double minU = octave_NaN;
+ double maxU = octave_NaN;
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ double d = std::abs (U.xelem(j,j));
+ if (xisnan (minU) || d < minU)
+ minU = d;
+
+ if (xisnan (maxU) || d > maxU)
+ maxU = d;
+ }
+
+ double rcond = (minU / maxU);
+ volatile double rcond_plus_one = rcond + 1.0;
+
+ if (rcond_plus_one == 1.0 || xisnan (rcond))
+ {
+ (*current_liboctave_warning_handler)
+ ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+ (*current_liboctave_warning_handler)
+ (" an eigenvalue. Convergence is not guaranteed");
+ }
+
+ return true;
+}
+
+static bool
+LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b,
+ bool cholB, const ColumnVector& permB, Complex sigma,
+ SparseComplexMatrix &L, SparseComplexMatrix &U,
+ octave_idx_type *P, octave_idx_type *Q)
+{
+ bool have_b = ! b.is_empty ();
+ octave_idx_type n = m.rows();
+
+ // Caclulate LU decomposition of 'A - sigma * B'
+ SparseComplexMatrix AminusSigmaB (m);
+
+ if (have_b)
+ {
+ if (cholB)
+ {
+ if (permB.length())
+ {
+ SparseMatrix tmp(n,n,n);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ tmp.xcidx(i) = i;
+ tmp.xridx(i) =
+ static_cast<octave_idx_type>(permB(i));
+ tmp.xdata(i) = 1;
+ }
+ tmp.xcidx(n) = n;
+
+ AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b *
+ tmp.transpose() * sigma;
+ }
+ else
+ AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b;
+ }
+ else
+ AminusSigmaB = AminusSigmaB - sigma * b;
+ }
+ else
+ {
+ SparseComplexMatrix sigmat (n, n, n);
+
+ // Create sigma * speye(n,n)
+ sigmat.xcidx (0) = 0;
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ sigmat.xdata(i) = sigma;
+ sigmat.xridx(i) = i;
+ sigmat.xcidx(i+1) = i + 1;
+ }
+
+ AminusSigmaB = AminusSigmaB - sigmat;
+ }
+
+ SparseComplexLU fact (AminusSigmaB);
+
+ L = fact.L ();
+ U = fact.U ();
+ const octave_idx_type *P2 = fact.row_perm ();
+ const octave_idx_type *Q2 = fact.col_perm ();
+
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ P[j] = P2[j];
+ Q[j] = Q2[j];
+ }
+
+ // Test condition number of LU decomposition
+ double minU = octave_NaN;
+ double maxU = octave_NaN;
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ double d = 0.;
+ if (U.xcidx(j+1) > U.xcidx(j) &&
+ U.xridx (U.xcidx(j+1)-1) == j)
+ d = std::abs (U.xdata (U.xcidx(j+1)-1));
+
+ if (xisnan (minU) || d < minU)
+ minU = d;
+
+ if (xisnan (maxU) || d > maxU)
+ maxU = d;
+ }
+
+ double rcond = (minU / maxU);
+ volatile double rcond_plus_one = rcond + 1.0;
+
+ if (rcond_plus_one == 1.0 || xisnan (rcond))
+ {
+ (*current_liboctave_warning_handler)
+ ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+ (*current_liboctave_warning_handler)
+ (" an eigenvalue. Convergence is not guaranteed");
+ }
+
+ return true;
+}
+
+static bool
+LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b,
+ bool cholB, const ColumnVector& permB, Complex sigma,
+ ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P,
+ octave_idx_type *Q)
+{
+ bool have_b = ! b.is_empty ();
+ octave_idx_type n = m.cols();
+
+ // Caclulate LU decomposition of 'A - sigma * B'
+ ComplexMatrix AminusSigmaB (m);
+
+ if (have_b)
+ {
+ if (cholB)
+ {
+ ComplexMatrix tmp = sigma * b.hermitian() * b;
+ const double *pB = permB.fortran_vec();
+ Complex *p = AminusSigmaB.fortran_vec();
+
+ if (permB.length())
+ {
+ for (octave_idx_type j = 0;
+ j < b.cols(); j++)
+ for (octave_idx_type i = 0;
+ i < b.rows(); i++)
+ *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]),
+ static_cast<octave_idx_type>(pB[j]));
+ }
+ else
+ AminusSigmaB = AminusSigmaB - tmp;
+ }
+ else
+ AminusSigmaB = AminusSigmaB - sigma * b;
+ }
+ else
+ {
+ Complex *p = AminusSigmaB.fortran_vec();
+
+ for (octave_idx_type i = 0; i < n; i++)
+ p[i*(n+1)] -= sigma;
+ }
+
+ ComplexLU fact (AminusSigmaB);
+
+ L = fact.P().transpose() * fact.L ();
+ U = fact.U ();
+ for (octave_idx_type j = 0; j < n; j++)
+ P[j] = Q[j] = j;
+
+ // Test condition number of LU decomposition
+ double minU = octave_NaN;
+ double maxU = octave_NaN;
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ double d = std::abs (U.xelem(j,j));
+ if (xisnan (minU) || d < minU)
+ minU = d;
+
+ if (xisnan (maxU) || d > maxU)
+ maxU = d;
+ }
+
+ double rcond = (minU / maxU);
+ volatile double rcond_plus_one = rcond + 1.0;
+
+ if (rcond_plus_one == 1.0 || xisnan (rcond))
+ {
+ (*current_liboctave_warning_handler)
+ ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
+ (*current_liboctave_warning_handler)
+ (" an eigenvalue. Convergence is not guaranteed");
+ }
+
+ return true;
+}
+
+template <class M>
+octave_idx_type
+EigsRealSymmetricMatrix (const M& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, Matrix &eig_vec,
+ ColumnVector &eig_val, const M& _b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream& os, double tol, int rvec,
+ bool cholB, int disp, int maxit)
+{
+ M b(_b);
+ octave_idx_type n = m.cols ();
+ octave_idx_type mode = 1;
+ bool have_b = ! b.is_empty();
+ bool note3 = false;
+ char bmat = 'I';
+ double sigma = 0.;
+ M bt;
+
+ if (m.rows() != m.cols())
+ {
+ (*current_liboctave_error_handler) ("eigs: A must be square");
+ return -1;
+ }
+ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: B must be square and the same size as A");
+ return -1;
+ }
+
+ if (resid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ resid = ColumnVector (octave_rand::vector(n));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p <= k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k and n");
+ return -1;
+ }
+
+ if (k > n )
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (have_b && cholB && permB.length() != 0)
+ {
+ // Check the we really have a permutation vector
+ if (permB.length() != n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ else
+ {
+ Array<bool> checked(n,false);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ octave_idx_type bidx =
+ static_cast<octave_idx_type> (permB(i));
+ if (checked(bidx) || bidx < 0 ||
+ bidx >= n || D_NINT (bidx) != bidx)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ }
+ }
+ }
+
+ if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+ typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+ typ != "SI")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecognized sigma value");
+ return -1;
+ }
+
+ if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: invalid sigma value for real symmetric problem");
+ return -1;
+ }
+
+ if (have_b)
+ {
+ // See Note 3 dsaupd
+ note3 = true;
+ if (cholB)
+ {
+ bt = b;
+ b = b.transpose();
+ if (permB.length() == 0)
+ {
+ permB = ColumnVector(n);
+ for (octave_idx_type i = 0; i < n; i++)
+ permB(i) = i;
+ }
+ }
+ else
+ {
+ if (! make_cholb(b, bt, permB))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: The matrix B is not positive definite");
+ return -1;
+ }
+ }
+ }
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ octave_idx_type lwork = p * (p + 8);
+
+ OCTAVE_LOCAL_BUFFER (double, v, n * p);
+ OCTAVE_LOCAL_BUFFER (double, workl, lwork);
+ OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
+ double *presid = resid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (dsaupd, DSAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dsaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ if (have_b)
+ {
+ Matrix mtmp (n,1);
+ for (octave_idx_type i = 0; i < n; i++)
+ mtmp(i,0) = workd[i + iptr(0) - 1];
+
+ mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
+
+ for (octave_idx_type i = 0; i < n; i++)
+ workd[i+iptr(1)-1] = mtmp(i,0);
+ }
+ else if (!vector_product (m, workd + iptr(0) - 1,
+ workd + iptr(1) - 1))
+ break;
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dsaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ eig_vec.resize (n, k);
+ double *z = eig_vec.fortran_vec ();
+
+ eig_val.resize (k);
+ double *d = eig_val.fortran_vec ();
+
+ F77_FUNC (dseupd, DSEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
+ F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
+ F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dseupd");
+ return -1;
+ }
+ else
+ {
+ if (info2 == 0)
+ {
+ octave_idx_type k2 = k / 2;
+ if (typ != "SM" && typ != "BE")
+ {
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ double dtmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = dtmp;
+ }
+ }
+
+ if (rvec)
+ {
+ if (typ != "SM" && typ != "BE")
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ dtmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = dtmp[j];
+ }
+ }
+
+ if (note3)
+ eig_vec = ltsolve(b, permB, eig_vec);
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dseupd", info2);
+ return -1;
+ }
+ }
+
+ return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsRealSymmetricMatrixShift (const M& m, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, Matrix &eig_vec,
+ ColumnVector &eig_val, const M& _b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream& os, double tol, int rvec,
+ bool cholB, int disp, int maxit)
+{
+ M b(_b);
+ octave_idx_type n = m.cols ();
+ octave_idx_type mode = 3;
+ bool have_b = ! b.is_empty();
+ std::string typ = "LM";
+
+ if (m.rows() != m.cols())
+ {
+ (*current_liboctave_error_handler) ("eigs: A must be square");
+ return -1;
+ }
+ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: B must be square and the same size as A");
+ return -1;
+ }
+
+ // FIXME: The "SM" type for mode 1 seems unstable though faster!!
+ //if (! std::abs (sigma))
+ // return EigsRealSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
+ // _b, permB, resid, os, tol, rvec, cholB,
+ // disp, maxit);
+
+ if (resid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ resid = ColumnVector (octave_rand::vector(n));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p <= k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k and n");
+ return -1;
+ }
+
+ if (k > n )
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (have_b && cholB && permB.length() != 0)
+ {
+ // Check the we really have a permutation vector
+ if (permB.length() != n)
+ {
+ (*current_liboctave_error_handler) ("eigs: permB vector invalid");
+ return -1;
+ }
+ else
+ {
+ Array<bool> checked(n,false);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ octave_idx_type bidx =
+ static_cast<octave_idx_type> (permB(i));
+ if (checked(bidx) || bidx < 0 ||
+ bidx >= n || D_NINT (bidx) != bidx)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ }
+ }
+ }
+
+ char bmat = 'I';
+ if (have_b)
+ bmat = 'G';
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ M L, U;
+
+ OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
+ OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
+
+ if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q))
+ return -1;
+
+ octave_idx_type lwork = p * (p + 8);
+
+ OCTAVE_LOCAL_BUFFER (double, v, n * p);
+ OCTAVE_LOCAL_BUFFER (double, workl, lwork);
+ OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
+ double *presid = resid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (dsaupd, DSAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dsaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ if (have_b)
+ {
+ if (ido == -1)
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ vector_product (m, workd+iptr(0)-1, dtmp);
+
+ Matrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = dtmp[P[i]];
+
+ lusolve (L, U, tmp);
+
+ double *ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ else if (ido == 2)
+ vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
+ else
+ {
+ double *ip2 = workd+iptr(2)-1;
+ Matrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = ip2[P[i]];
+
+ lusolve (L, U, tmp);
+
+ ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ }
+ else
+ {
+ if (ido == 2)
+ {
+ for (octave_idx_type i = 0; i < n; i++)
+ workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1];
+ }
+ else
+ {
+ double *ip2 = workd+iptr(0)-1;
+ Matrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = ip2[P[i]];
+
+ lusolve (L, U, tmp);
+
+ ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ }
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dsaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ eig_vec.resize (n, k);
+ double *z = eig_vec.fortran_vec ();
+
+ eig_val.resize (k);
+ double *d = eig_val.fortran_vec ();
+
+ F77_FUNC (dseupd, DSEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
+ F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
+ k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dseupd");
+ return -1;
+ }
+ else
+ {
+ if (info2 == 0)
+ {
+ octave_idx_type k2 = k / 2;
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ double dtmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = dtmp;
+ }
+
+ if (rvec)
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ dtmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = dtmp[j];
+ }
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dseupd", info2);
+ return -1;
+ }
+ }
+
+ return ip(4);
+}
+
+octave_idx_type
+EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
+ const std::string &_typ, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, Matrix &eig_vec,
+ ColumnVector &eig_val, ColumnVector &resid,
+ std::ostream& os, double tol, int rvec, bool cholB,
+ int disp, int maxit)
+{
+ std::string typ (_typ);
+ bool have_sigma = (sigma ? true : false);
+ char bmat = 'I';
+ octave_idx_type mode = 1;
+ int err = 0;
+
+ if (resid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ resid = ColumnVector (octave_rand::vector(n));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p <= k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k and n");
+ return -1;
+ }
+
+ if (k > n )
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (! have_sigma)
+ {
+ if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+ typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+ typ != "SI")
+ (*current_liboctave_error_handler)
+ ("eigs: unrecognized sigma value");
+
+ if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: invalid sigma value for real symmetric problem");
+ return -1;
+ }
+
+ if (typ == "SM")
+ {
+ typ = "LM";
+ sigma = 0.;
+ mode = 3;
+ }
+ }
+ else if (! std::abs (sigma))
+ typ = "SM";
+ else
+ {
+ typ = "LM";
+ mode = 3;
+ }
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ octave_idx_type lwork = p * (p + 8);
+
+ OCTAVE_LOCAL_BUFFER (double, v, n * p);
+ OCTAVE_LOCAL_BUFFER (double, workl, lwork);
+ OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
+ double *presid = resid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (dsaupd, DSAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dsaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ double *ip2 = workd + iptr(0) - 1;
+ ColumnVector x(n);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ x(i) = *ip2++;
+
+ ColumnVector y = fun (x, err);
+
+ if (err)
+ return false;
+
+ ip2 = workd + iptr(1) - 1;
+ for (octave_idx_type i = 0; i < n; i++)
+ *ip2++ = y(i);
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dsaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ eig_vec.resize (n, k);
+ double *z = eig_vec.fortran_vec ();
+
+ eig_val.resize (k);
+ double *d = eig_val.fortran_vec ();
+
+ F77_FUNC (dseupd, DSEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
+ F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
+ k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dseupd");
+ return -1;
+ }
+ else
+ {
+ if (info2 == 0)
+ {
+ octave_idx_type k2 = k / 2;
+ if (typ != "SM" && typ != "BE")
+ {
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ double dtmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = dtmp;
+ }
+ }
+
+ if (rvec)
+ {
+ if (typ != "SM" && typ != "BE")
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ dtmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = dtmp[j];
+ }
+ }
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dseupd", info2);
+ return -1;
+ }
+ }
+
+ return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsRealNonSymmetricMatrix (const M& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val, const M& _b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream& os, double tol, int rvec,
+ bool cholB, int disp, int maxit)
+{
+ M b(_b);
+ octave_idx_type n = m.cols ();
+ octave_idx_type mode = 1;
+ bool have_b = ! b.is_empty();
+ bool note3 = false;
+ char bmat = 'I';
+ double sigmar = 0.;
+ double sigmai = 0.;
+ M bt;
+
+ if (m.rows() != m.cols())
+ {
+ (*current_liboctave_error_handler) ("eigs: A must be square");
+ return -1;
+ }
+ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: B must be square and the same size as A");
+ return -1;
+ }
+
+ if (resid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ resid = ColumnVector (octave_rand::vector(n));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2 + 1;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p < k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k+1 and n");
+ return -1;
+ }
+
+ if (k > n - 1)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (have_b && cholB && permB.length() != 0)
+ {
+ // Check the we really have a permutation vector
+ if (permB.length() != n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ else
+ {
+ Array<bool> checked(n,false);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ octave_idx_type bidx =
+ static_cast<octave_idx_type> (permB(i));
+ if (checked(bidx) || bidx < 0 ||
+ bidx >= n || D_NINT (bidx) != bidx)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ }
+ }
+ }
+
+ if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+ typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+ typ != "SI")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecognized sigma value");
+ return -1;
+ }
+
+ if (typ == "LA" || typ == "SA" || typ == "BE")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: invalid sigma value for unsymmetric problem");
+ return -1;
+ }
+
+ if (have_b)
+ {
+ // See Note 3 dsaupd
+ note3 = true;
+ if (cholB)
+ {
+ bt = b;
+ b = b.transpose();
+ if (permB.length() == 0)
+ {
+ permB = ColumnVector(n);
+ for (octave_idx_type i = 0; i < n; i++)
+ permB(i) = i;
+ }
+ }
+ else
+ {
+ if (! make_cholb(b, bt, permB))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: The matrix B is not positive definite");
+ return -1;
+ }
+ }
+ }
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ octave_idx_type lwork = 3 * p * (p + 2);
+
+ OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
+ OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
+ OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
+ double *presid = resid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (dnaupd, DNAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dnaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ if (have_b)
+ {
+ Matrix mtmp (n,1);
+ for (octave_idx_type i = 0; i < n; i++)
+ mtmp(i,0) = workd[i + iptr(0) - 1];
+
+ mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
+
+ for (octave_idx_type i = 0; i < n; i++)
+ workd[i+iptr(1)-1] = mtmp(i,0);
+ }
+ else if (!vector_product (m, workd + iptr(0) - 1,
+ workd + iptr(1) - 1))
+ break;
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dnaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ Matrix eig_vec2 (n, k + 1);
+ double *z = eig_vec2.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
+ OCTAVE_LOCAL_BUFFER (double, di, k + 1);
+ OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
+
+ F77_FUNC (dneupd, DNEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
+ sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
+ F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dneupd");
+ return -1;
+ }
+ else
+ {
+ eig_val.resize (k+1);
+ Complex *d = eig_val.fortran_vec ();
+
+ if (info2 == 0)
+ {
+ octave_idx_type jj = 0;
+ for (octave_idx_type i = 0; i < k+1; i++)
+ {
+ if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
+ jj++;
+ else
+ d [i-jj] = Complex (dr[i], di[i]);
+ }
+ if (jj == 0 && !rvec)
+ for (octave_idx_type i = 0; i < k; i++)
+ d[i] = d[i+1];
+
+ octave_idx_type k2 = k / 2;
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ Complex dtmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = dtmp;
+ }
+ eig_val.resize(k);
+
+ if (rvec)
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ dtmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = dtmp[j];
+ }
+
+ eig_vec.resize (n, k);
+ octave_idx_type i = 0;
+ while (i < k)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (i+1) * n;
+ if (std::imag(eig_val(i)) == 0)
+ {
+ for (octave_idx_type j = 0; j < n; j++)
+ eig_vec(j,i) =
+ Complex(z[j+off1],0.);
+ i++;
+ }
+ else
+ {
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ eig_vec(j,i) =
+ Complex(z[j+off1],z[j+off2]);
+ if (i < k - 1)
+ eig_vec(j,i+1) =
+ Complex(z[j+off1],-z[j+off2]);
+ }
+ i+=2;
+ }
+ }
+
+ if (note3)
+ eig_vec = ltsolve(M (b), permB, eig_vec);
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dneupd", info2);
+ return -1;
+ }
+ }
+
+ return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsRealNonSymmetricMatrixShift (const M& m, double sigmar,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info,
+ ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val, const M& _b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream& os, double tol, int rvec,
+ bool cholB, int disp, int maxit)
+{
+ M b(_b);
+ octave_idx_type n = m.cols ();
+ octave_idx_type mode = 3;
+ bool have_b = ! b.is_empty();
+ std::string typ = "LM";
+ double sigmai = 0.;
+
+ if (m.rows() != m.cols())
+ {
+ (*current_liboctave_error_handler) ("eigs: A must be square");
+ return -1;
+ }
+ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: B must be square and the same size as A");
+ return -1;
+ }
+
+ // FIXME: The "SM" type for mode 1 seems unstable though faster!!
+ //if (! std::abs (sigmar))
+ // return EigsRealNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
+ // _b, permB, resid, os, tol, rvec, cholB,
+ // disp, maxit);
+
+ if (resid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ resid = ColumnVector (octave_rand::vector(n));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2 + 1;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p < k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k+1 and n");
+ return -1;
+ }
+
+ if (k > n - 1)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (have_b && cholB && permB.length() != 0)
+ {
+ // Check that we really have a permutation vector
+ if (permB.length() != n)
+ {
+ (*current_liboctave_error_handler) ("eigs: permB vector invalid");
+ return -1;
+ }
+ else
+ {
+ Array<bool> checked(n,false);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ octave_idx_type bidx =
+ static_cast<octave_idx_type> (permB(i));
+ if (checked(bidx) || bidx < 0 ||
+ bidx >= n || D_NINT (bidx) != bidx)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ }
+ }
+ }
+
+ char bmat = 'I';
+ if (have_b)
+ bmat = 'G';
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ M L, U;
+
+ OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
+ OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
+
+ if (! LuAminusSigmaB(m, b, cholB, permB, sigmar, L, U, P, Q))
+ return -1;
+
+ octave_idx_type lwork = 3 * p * (p + 2);
+
+ OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
+ OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
+ OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
+ double *presid = resid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (dnaupd, DNAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dsaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ if (have_b)
+ {
+ if (ido == -1)
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ vector_product (m, workd+iptr(0)-1, dtmp);
+
+ Matrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = dtmp[P[i]];
+
+ lusolve (L, U, tmp);
+
+ double *ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ else if (ido == 2)
+ vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
+ else
+ {
+ double *ip2 = workd+iptr(2)-1;
+ Matrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = ip2[P[i]];
+
+ lusolve (L, U, tmp);
+
+ ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ }
+ else
+ {
+ if (ido == 2)
+ {
+ for (octave_idx_type i = 0; i < n; i++)
+ workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1];
+ }
+ else
+ {
+ double *ip2 = workd+iptr(0)-1;
+ Matrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = ip2[P[i]];
+
+ lusolve (L, U, tmp);
+
+ ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ }
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dsaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ Matrix eig_vec2 (n, k + 1);
+ double *z = eig_vec2.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
+ OCTAVE_LOCAL_BUFFER (double, di, k + 1);
+ OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
+
+ F77_FUNC (dneupd, DNEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
+ sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
+ F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dneupd");
+ return -1;
+ }
+ else
+ {
+ eig_val.resize (k+1);
+ Complex *d = eig_val.fortran_vec ();
+
+ if (info2 == 0)
+ {
+ octave_idx_type jj = 0;
+ for (octave_idx_type i = 0; i < k+1; i++)
+ {
+ if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
+ jj++;
+ else
+ d [i-jj] = Complex (dr[i], di[i]);
+ }
+ if (jj == 0 && !rvec)
+ for (octave_idx_type i = 0; i < k; i++)
+ d[i] = d[i+1];
+
+ octave_idx_type k2 = k / 2;
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ Complex dtmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = dtmp;
+ }
+ eig_val.resize(k);
+
+ if (rvec)
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ dtmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = dtmp[j];
+ }
+
+ eig_vec.resize (n, k);
+ octave_idx_type i = 0;
+ while (i < k)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (i+1) * n;
+ if (std::imag(eig_val(i)) == 0)
+ {
+ for (octave_idx_type j = 0; j < n; j++)
+ eig_vec(j,i) =
+ Complex(z[j+off1],0.);
+ i++;
+ }
+ else
+ {
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ eig_vec(j,i) =
+ Complex(z[j+off1],z[j+off2]);
+ if (i < k - 1)
+ eig_vec(j,i+1) =
+ Complex(z[j+off1],-z[j+off2]);
+ }
+ i+=2;
+ }
+ }
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dneupd", info2);
+ return -1;
+ }
+ }
+
+ return ip(4);
+}
+
+octave_idx_type
+EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
+ const std::string &_typ, double sigmar,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val, ColumnVector &resid,
+ std::ostream& os, double tol, int rvec, bool cholB,
+ int disp, int maxit)
+{
+ std::string typ (_typ);
+ bool have_sigma = (sigmar ? true : false);
+ char bmat = 'I';
+ double sigmai = 0.;
+ octave_idx_type mode = 1;
+ int err = 0;
+
+ if (resid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ resid = ColumnVector (octave_rand::vector(n));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2 + 1;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p < k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k+1 and n");
+ return -1;
+ }
+
+ if (k > n - 1)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+
+ if (! have_sigma)
+ {
+ if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+ typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+ typ != "SI")
+ (*current_liboctave_error_handler)
+ ("eigs: unrecognized sigma value");
+
+ if (typ == "LA" || typ == "SA" || typ == "BE")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: invalid sigma value for unsymmetric problem");
+ return -1;
+ }
+
+ if (typ == "SM")
+ {
+ typ = "LM";
+ sigmar = 0.;
+ mode = 3;
+ }
+ }
+ else if (! std::abs (sigmar))
+ typ = "SM";
+ else
+ {
+ typ = "LM";
+ mode = 3;
+ }
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ octave_idx_type lwork = 3 * p * (p + 2);
+
+ OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
+ OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
+ OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
+ double *presid = resid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (dnaupd, DNAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dnaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ double *ip2 = workd + iptr(0) - 1;
+ ColumnVector x(n);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ x(i) = *ip2++;
+
+ ColumnVector y = fun (x, err);
+
+ if (err)
+ return false;
+
+ ip2 = workd + iptr(1) - 1;
+ for (octave_idx_type i = 0; i < n; i++)
+ *ip2++ = y(i);
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dsaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ Matrix eig_vec2 (n, k + 1);
+ double *z = eig_vec2.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
+ OCTAVE_LOCAL_BUFFER (double, di, k + 1);
+ OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
+
+ F77_FUNC (dneupd, DNEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
+ sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
+ F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in dneupd");
+ return -1;
+ }
+ else
+ {
+ eig_val.resize (k+1);
+ Complex *d = eig_val.fortran_vec ();
+
+ if (info2 == 0)
+ {
+ octave_idx_type jj = 0;
+ for (octave_idx_type i = 0; i < k+1; i++)
+ {
+ if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
+ jj++;
+ else
+ d [i-jj] = Complex (dr[i], di[i]);
+ }
+ if (jj == 0 && !rvec)
+ for (octave_idx_type i = 0; i < k; i++)
+ d[i] = d[i+1];
+
+ octave_idx_type k2 = k / 2;
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ Complex dtmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = dtmp;
+ }
+ eig_val.resize(k);
+
+ if (rvec)
+ {
+ OCTAVE_LOCAL_BUFFER (double, dtmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ dtmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = dtmp[j];
+ }
+
+ eig_vec.resize (n, k);
+ octave_idx_type i = 0;
+ while (i < k)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (i+1) * n;
+ if (std::imag(eig_val(i)) == 0)
+ {
+ for (octave_idx_type j = 0; j < n; j++)
+ eig_vec(j,i) =
+ Complex(z[j+off1],0.);
+ i++;
+ }
+ else
+ {
+ for (octave_idx_type j = 0; j < n; j++)
+ {
+ eig_vec(j,i) =
+ Complex(z[j+off1],z[j+off2]);
+ if (i < k - 1)
+ eig_vec(j,i+1) =
+ Complex(z[j+off1],-z[j+off2]);
+ }
+ i+=2;
+ }
+ }
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dneupd", info2);
+ return -1;
+ }
+ }
+
+ return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsComplexNonSymmetricMatrix (const M& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val, const M& _b,
+ ColumnVector &permB,
+ ComplexColumnVector &cresid,
+ std::ostream& os, double tol, int rvec,
+ bool cholB, int disp, int maxit)
+{
+ M b(_b);
+ octave_idx_type n = m.cols ();
+ octave_idx_type mode = 1;
+ bool have_b = ! b.is_empty();
+ bool note3 = false;
+ char bmat = 'I';
+ Complex sigma = 0.;
+ M bt;
+
+ if (m.rows() != m.cols())
+ {
+ (*current_liboctave_error_handler) ("eigs: A must be square");
+ return -1;
+ }
+ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: B must be square and the same size as A");
+ return -1;
+ }
+
+ if (cresid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ Array<double> rr (octave_rand::vector(n));
+ Array<double> ri (octave_rand::vector(n));
+ cresid = ComplexColumnVector (n);
+ for (octave_idx_type i = 0; i < n; i++)
+ cresid(i) = Complex(rr(i),ri(i));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2 + 1;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p < k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k+1 and n");
+ return -1;
+ }
+
+ if (k > n - 1)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (have_b && cholB && permB.length() != 0)
+ {
+ // Check the we really have a permutation vector
+ if (permB.length() != n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ else
+ {
+ Array<bool> checked(n,false);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ octave_idx_type bidx =
+ static_cast<octave_idx_type> (permB(i));
+ if (checked(bidx) || bidx < 0 ||
+ bidx >= n || D_NINT (bidx) != bidx)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ }
+ }
+ }
+
+ if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+ typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+ typ != "SI")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecognized sigma value");
+ return -1;
+ }
+
+ if (typ == "LA" || typ == "SA" || typ == "BE")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: invalid sigma value for complex problem");
+ return -1;
+ }
+
+ if (have_b)
+ {
+ // See Note 3 dsaupd
+ note3 = true;
+ if (cholB)
+ {
+ bt = b;
+ b = b.hermitian();
+ if (permB.length() == 0)
+ {
+ permB = ColumnVector(n);
+ for (octave_idx_type i = 0; i < n; i++)
+ permB(i) = i;
+ }
+ }
+ else
+ {
+ if (! make_cholb(b, bt, permB))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: The matrix B is not positive definite");
+ return -1;
+ }
+ }
+ }
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ octave_idx_type lwork = p * (3 * p + 5);
+
+ OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
+ OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
+ OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
+ OCTAVE_LOCAL_BUFFER (double, rwork, p);
+ Complex *presid = cresid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (znaupd, ZNAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, rwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in znaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ if (have_b)
+ {
+ ComplexMatrix mtmp (n,1);
+ for (octave_idx_type i = 0; i < n; i++)
+ mtmp(i,0) = workd[i + iptr(0) - 1];
+ mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
+ for (octave_idx_type i = 0; i < n; i++)
+ workd[i+iptr(1)-1] = mtmp(i,0);
+
+ }
+ else if (!vector_product (m, workd + iptr(0) - 1,
+ workd + iptr(1) - 1))
+ break;
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in znaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ eig_vec.resize (n, k);
+ Complex *z = eig_vec.fortran_vec ();
+
+ eig_val.resize (k+1);
+ Complex *d = eig_val.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
+
+ F77_FUNC (zneupd, ZNEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
+ F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
+ k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in zneupd");
+ return -1;
+ }
+
+ if (info2 == 0)
+ {
+ octave_idx_type k2 = k / 2;
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ Complex ctmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = ctmp;
+ }
+ eig_val.resize(k);
+
+ if (rvec)
+ {
+ OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ ctmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = ctmp[j];
+ }
+
+ if (note3)
+ eig_vec = ltsolve(b, permB, eig_vec);
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in zneupd", info2);
+ return -1;
+ }
+
+ return ip(4);
+}
+
+template <class M>
+octave_idx_type
+EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info,
+ ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val, const M& _b,
+ ColumnVector &permB,
+ ComplexColumnVector &cresid,
+ std::ostream& os, double tol, int rvec,
+ bool cholB, int disp, int maxit)
+{
+ M b(_b);
+ octave_idx_type n = m.cols ();
+ octave_idx_type mode = 3;
+ bool have_b = ! b.is_empty();
+ std::string typ = "LM";
+
+ if (m.rows() != m.cols())
+ {
+ (*current_liboctave_error_handler) ("eigs: A must be square");
+ return -1;
+ }
+ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: B must be square and the same size as A");
+ return -1;
+ }
+
+ // FIXME: The "SM" type for mode 1 seems unstable though faster!!
+ //if (! std::abs (sigma))
+ // return EigsComplexNonSymmetricMatrix (m, "SM", k, p, info, eig_vec,
+ // eig_val, _b, permB, cresid, os, tol,
+ // rvec, cholB, disp, maxit);
+
+ if (cresid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ Array<double> rr (octave_rand::vector(n));
+ Array<double> ri (octave_rand::vector(n));
+ cresid = ComplexColumnVector (n);
+ for (octave_idx_type i = 0; i < n; i++)
+ cresid(i) = Complex(rr(i),ri(i));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2 + 1;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p < k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k+1 and n");
+ return -1;
+ }
+
+ if (k > n - 1)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (have_b && cholB && permB.length() != 0)
+ {
+ // Check that we really have a permutation vector
+ if (permB.length() != n)
+ {
+ (*current_liboctave_error_handler) ("eigs: permB vector invalid");
+ return -1;
+ }
+ else
+ {
+ Array<bool> checked(n,false);
+ for (octave_idx_type i = 0; i < n; i++)
+ {
+ octave_idx_type bidx =
+ static_cast<octave_idx_type> (permB(i));
+ if (checked(bidx) || bidx < 0 ||
+ bidx >= n || D_NINT (bidx) != bidx)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: permB vector invalid");
+ return -1;
+ }
+ }
+ }
+ }
+
+ char bmat = 'I';
+ if (have_b)
+ bmat = 'G';
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ M L, U;
+
+ OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
+ OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
+
+ if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q))
+ return -1;
+
+ octave_idx_type lwork = p * (3 * p + 5);
+
+ OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
+ OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
+ OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
+ OCTAVE_LOCAL_BUFFER (double, rwork, p);
+ Complex *presid = cresid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (znaupd, ZNAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, rwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in znaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ if (have_b)
+ {
+ if (ido == -1)
+ {
+ OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
+
+ vector_product (m, workd+iptr(0)-1, ctmp);
+
+ ComplexMatrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = ctmp[P[i]];
+
+ lusolve (L, U, tmp);
+
+ Complex *ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ else if (ido == 2)
+ vector_product (b, workd + iptr(0) - 1, workd + iptr(1) - 1);
+ else
+ {
+ Complex *ip2 = workd+iptr(2)-1;
+ ComplexMatrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = ip2[P[i]];
+
+ lusolve (L, U, tmp);
+
+ ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ }
+ else
+ {
+ if (ido == 2)
+ {
+ for (octave_idx_type i = 0; i < n; i++)
+ workd[iptr(0) + i - 1] =
+ workd[iptr(1) + i - 1];
+ }
+ else
+ {
+ Complex *ip2 = workd+iptr(0)-1;
+ ComplexMatrix tmp(n, 1);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp(i,0) = ip2[P[i]];
+
+ lusolve (L, U, tmp);
+
+ ip2 = workd+iptr(1)-1;
+ for (octave_idx_type i = 0; i < n; i++)
+ ip2[Q[i]] = tmp(i,0);
+ }
+ }
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dsaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ eig_vec.resize (n, k);
+ Complex *z = eig_vec.fortran_vec ();
+
+ eig_val.resize (k+1);
+ Complex *d = eig_val.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
+
+ F77_FUNC (zneupd, ZNEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
+ F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
+ k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in zneupd");
+ return -1;
+ }
+
+ if (info2 == 0)
+ {
+ octave_idx_type k2 = k / 2;
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ Complex ctmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = ctmp;
+ }
+ eig_val.resize(k);
+
+ if (rvec)
+ {
+ OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ ctmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = ctmp[j];
+ }
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in zneupd", info2);
+ return -1;
+ }
+
+ return ip(4);
+}
+
+octave_idx_type
+EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
+ const std::string &_typ, Complex sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ ComplexColumnVector &cresid, std::ostream& os,
+ double tol, int rvec, bool cholB, int disp,
+ int maxit)
+{
+ std::string typ (_typ);
+ bool have_sigma = (std::abs(sigma) ? true : false);
+ char bmat = 'I';
+ octave_idx_type mode = 1;
+ int err = 0;
+
+ if (cresid.is_empty())
+ {
+ std::string rand_dist = octave_rand::distribution();
+ octave_rand::distribution("uniform");
+ Array<double> rr (octave_rand::vector(n));
+ Array<double> ri (octave_rand::vector(n));
+ cresid = ComplexColumnVector (n);
+ for (octave_idx_type i = 0; i < n; i++)
+ cresid(i) = Complex(rr(i),ri(i));
+ octave_rand::distribution(rand_dist);
+ }
+
+ if (p < 0)
+ {
+ p = k * 2 + 1;
+
+ if (p < 20)
+ p = 20;
+
+ if (p > n - 1)
+ p = n - 1 ;
+ }
+ else if (p < k || p > n)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: opts.p must be between k+1 and n");
+ return -1;
+ }
+
+ if (k > n - 1)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: Too many eigenvalues to extract (k >= n-1).\n"
+ " Use 'eig(full(A))' instead");
+ return -1;
+ }
+
+ if (! have_sigma)
+ {
+ if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
+ typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
+ typ != "SI")
+ (*current_liboctave_error_handler)
+ ("eigs: unrecognized sigma value");
+
+ if (typ == "LA" || typ == "SA" || typ == "BE")
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: invalid sigma value for complex problem");
+ return -1;
+ }
+
+ if (typ == "SM")
+ {
+ typ = "LM";
+ sigma = 0.;
+ mode = 3;
+ }
+ }
+ else if (! std::abs (sigma))
+ typ = "SM";
+ else
+ {
+ typ = "LM";
+ mode = 3;
+ }
+
+ Array<octave_idx_type> ip (11);
+ octave_idx_type *iparam = ip.fortran_vec ();
+
+ ip(0) = 1; //ishift
+ ip(1) = 0; // ip(1) not referenced
+ ip(2) = maxit; // mxiter, maximum number of iterations
+ ip(3) = 1; // NB blocksize in recurrence
+ ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
+ ip(5) = 0; //ip(5) not referenced
+ ip(6) = mode; // mode
+ ip(7) = 0;
+ ip(8) = 0;
+ ip(9) = 0;
+ ip(10) = 0;
+ // ip(7) to ip(10) return values
+
+ Array<octave_idx_type> iptr (14);
+ octave_idx_type *ipntr = iptr.fortran_vec ();
+
+ octave_idx_type ido = 0;
+ int iter = 0;
+ octave_idx_type lwork = p * (3 * p + 5);
+
+ OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
+ OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
+ OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
+ OCTAVE_LOCAL_BUFFER (double, rwork, p);
+ Complex *presid = cresid.fortran_vec ();
+
+ do
+ {
+ F77_FUNC (znaupd, ZNAUPD)
+ (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
+ k, tol, presid, p, v, n, iparam,
+ ipntr, workd, workl, lwork, rwork, info
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in znaupd");
+ return -1;
+ }
+
+ if (disp > 0 && !xisnan(workl[iptr(5)-1]))
+ {
+ if (iter++)
+ {
+ os << "Iteration " << iter - 1 <<
+ ": a few Ritz values of the " << p << "-by-" <<
+ p << " matrix\n";
+ for (int i = 0 ; i < k; i++)
+ os << " " << workl[iptr(5)+i-1] << "\n";
+ }
+
+ // This is a kludge, as ARPACK doesn't give its
+ // iteration pointer. But as workl[iptr(5)-1] is
+ // an output value updated at each iteration, setting
+ // a value in this array to NaN and testing for it
+ // is a way of obtaining the iteration counter.
+ if (ido != 99)
+ workl[iptr(5)-1] = octave_NaN;
+ }
+
+ if (ido == -1 || ido == 1 || ido == 2)
+ {
+ Complex *ip2 = workd + iptr(0) - 1;
+ ComplexColumnVector x(n);
+
+ for (octave_idx_type i = 0; i < n; i++)
+ x(i) = *ip2++;
+
+ ComplexColumnVector y = fun (x, err);
+
+ if (err)
+ return false;
+
+ ip2 = workd + iptr(1) - 1;
+ for (octave_idx_type i = 0; i < n; i++)
+ *ip2++ = y(i);
+ }
+ else
+ {
+ if (info < 0)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in dsaupd", info);
+ return -1;
+ }
+ break;
+ }
+ }
+ while (1);
+
+ octave_idx_type info2;
+
+ // We have a problem in that the size of the C++ bool
+ // type relative to the fortran logical type. It appears
+ // that fortran uses 4-bytes per logical and C++ 1-byte
+ // per bool, though this might be system dependent. As
+ // long as the HOWMNY arg is not "S", the logical array
+ // is just workspace for ARPACK, so use int type to
+ // avoid problems.
+ Array<int> s (p);
+ int *sel = s.fortran_vec ();
+
+ eig_vec.resize (n, k);
+ Complex *z = eig_vec.fortran_vec ();
+
+ eig_val.resize (k+1);
+ Complex *d = eig_val.fortran_vec ();
+
+ OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
+
+ F77_FUNC (zneupd, ZNEUPD)
+ (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
+ F77_CONST_CHAR_ARG2 (&bmat, 1), n,
+ F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
+ k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
+ F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
+
+ if (f77_exception_encountered)
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: unrecoverable exception encountered in zneupd");
+ return -1;
+ }
+
+ if (info2 == 0)
+ {
+ octave_idx_type k2 = k / 2;
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ Complex ctmp = d[i];
+ d[i] = d[k - i - 1];
+ d[k - i - 1] = ctmp;
+ }
+ eig_val.resize(k);
+
+ if (rvec)
+ {
+ OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
+
+ for (octave_idx_type i = 0; i < k2; i++)
+ {
+ octave_idx_type off1 = i * n;
+ octave_idx_type off2 = (k - i - 1) * n;
+
+ if (off1 == off2)
+ continue;
+
+ for (octave_idx_type j = 0; j < n; j++)
+ ctmp[j] = z[off1 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off1 + j] = z[off2 + j];
+
+ for (octave_idx_type j = 0; j < n; j++)
+ z[off2 + j] = ctmp[j];
+ }
+ }
+ }
+ else
+ {
+ (*current_liboctave_error_handler)
+ ("eigs: error %d in zneupd", info2);
+ return -1;
+ }
+
+ return ip(4);
+}
+
+#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
+extern octave_idx_type
+EigsRealSymmetricMatrix (const Matrix& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, Matrix &eig_vec,
+ ColumnVector &eig_val, const Matrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, Matrix &eig_vec,
+ ColumnVector &eig_val, const SparseMatrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream& os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricMatrixShift (const Matrix& m, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, Matrix &eig_vec,
+ ColumnVector &eig_val, const Matrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, Matrix &eig_vec,
+ ColumnVector &eig_val, const SparseMatrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
+ const std::string &typ, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info,
+ Matrix &eig_vec, ColumnVector &eig_val,
+ ColumnVector &resid, std::ostream &os,
+ double tol = DBL_EPSILON, int rvec = 0,
+ bool cholB = 0, int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val, const Matrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ const SparseMatrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info,
+ ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val, const Matrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info,
+ ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ const SparseMatrix& b,
+ ColumnVector &permB, ColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
+ const std::string &_typ, double sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ ColumnVector &resid, std::ostream& os,
+ double tol = DBL_EPSILON, int rvec = 0,
+ bool cholB = 0, int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ const ComplexMatrix& b, ColumnVector &permB,
+ ComplexColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m,
+ const std::string typ, octave_idx_type k,
+ octave_idx_type p, octave_idx_type &info,
+ ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ const SparseComplexMatrix& b,
+ ColumnVector &permB,
+ ComplexColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0, int disp = 0,
+ int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info,
+ ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ const ComplexMatrix& b,
+ ColumnVector &permB,
+ ComplexColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0,
+ int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m,
+ Complex sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info,
+ ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ const SparseComplexMatrix& b,
+ ColumnVector &permB,
+ ComplexColumnVector &resid,
+ std::ostream &os, double tol = DBL_EPSILON,
+ int rvec = 0, bool cholB = 0,
+ int disp = 0, int maxit = 300);
+
+extern octave_idx_type
+EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
+ const std::string &_typ, Complex sigma,
+ octave_idx_type k, octave_idx_type p,
+ octave_idx_type &info, ComplexMatrix &eig_vec,
+ ComplexColumnVector &eig_val,
+ ComplexColumnVector &resid, std::ostream& os,
+ double tol = DBL_EPSILON, int rvec = 0,
+ bool cholB = 0, int disp = 0, int maxit = 300);
+#endif
+
+#ifndef _MSC_VER
+template static octave_idx_type
+lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
+
+template static octave_idx_type
+lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
+ ComplexMatrix&);
+
+template static octave_idx_type
+lusolve (const Matrix&, const Matrix&, Matrix&);
+
+template static octave_idx_type
+lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
+
+template static ComplexMatrix
+ltsolve (const SparseComplexMatrix&, const ColumnVector&,
+ const ComplexMatrix&);
+
+template static Matrix
+ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
+
+template static ComplexMatrix
+ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+template static Matrix
+ltsolve (const Matrix&, const ColumnVector&, const Matrix&);
+
+template static ComplexMatrix
+utsolve (const SparseComplexMatrix&, const ColumnVector&,
+ const ComplexMatrix&);
+
+template static Matrix
+utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
+
+template static ComplexMatrix
+utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
+
+template static Matrix
+utsolve (const Matrix&, const ColumnVector&, const Matrix&);
+#endif
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/scripts/sparse/Makefile.in b/scripts/sparse/Makefile.in
--- a/scripts/sparse/Makefile.in
+++ b/scripts/sparse/Makefile.in
@@ -35,7 +35,7 @@
SOURCES = cgs.m colperm.m etreeplot.m gplot.m nonzeros.m normest.m \
pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \
spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \
- spvcat.m spy.m treelayout.m treeplot.m
+ spvcat.m spy.m svds.m treelayout.m treeplot.m
DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
diff --git a/scripts/sparse/svds.m b/scripts/sparse/svds.m
new file mode 100755
--- /dev/null
+++ b/scripts/sparse/svds.m
@@ -0,0 +1,231 @@
+## Copyright (C) 2006 David Bateman
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} address@hidden =} svds (@var{a})
+## @deftypefnx {Function File} address@hidden =} svds (@var{a}, @var{k})
+## @deftypefnx {Function File} address@hidden =} svds (@var{a}, @var{k},
@var{sigma})
+## @deftypefnx {Function File} address@hidden =} svds (@var{a}, @var{k},
@var{sigma}, @var{opts})
+## @deftypefnx {Function File} address@hidden, @var{s}, @var{v}, @var{flag}]
=} svds (@dots{})
+##
+## Find a few singular values of the matrix @var{a}. The singular values
+## are calculated using
+##
+## @example
+## @group
+## address@hidden, @var{n}] = size(@var{a})
+## @var{s} = eigs([sparse(@var{m}, @var{m}), @var{a}; @dots{}
+## @var{a}', sparse(@var{n}, @var{n})])
+## @end group
+## @end example
+##
+## The eigenvalues returned by @code{eigs} correspond to the singular
+## values of @var{a}. The number of singular values to calculate is given
+## by @var{k}, whose default value is 6.
+##
+## The argument @var{sigma} can be used to specify which singular values
+## to find. @var{sigma} can be either the string 'L', the default, in
+## which case the largest singular values of @var{a} are found. Otherwise
+## @var{sigma} should be a real scalar, in which case the singular values
+## closest to @var{sigma} are found. Note that for relatively small values
+## of @var{sigma}, there is the chance that the requested number of singular
+## values are not returned. In that case @var{sigma} should be increased.
+##
+## If @var{opts} is given, then it is a structure that defines options
+## that @code{svds} will pass to @var{eigs}. The possible fields of this
+## structure are therefore determined by @code{eigs}. By default three
+## fields of this structure are set by @code{svds}.
+##
+## @table @code
+## @item tol
+## The required convergence tolerance for the singular values. @code{eigs}
+## is passed @var{tol} divided by @code{sqrt(2)}. The default value is
+## 1e-10.
+##
+## @item maxit
+## The maximum number of iterations. The defaut is 300.
+##
+## @item disp
+## The level of diagnostic printout. If @code{disp} is 0 then there is no
+## printout. The default value is 0.
+## @end table
+##
+## If more than one output argument is given, then @code{svds} also
+## calculates the left and right singular vectors of @var{a}. @var{flag}
+## is used to signal the convergence of @code{svds}. If @code{svds}
+## converges to the desired tolerance, then @var{flag} given by
+##
+## @example
+## @group
+## norm (@var{a} * @var{v} - @var{u} * @var{s}, 1) <= @dots{}
+## @var{tol} * norm (@var{a}, 1)
+## @end group
+## @end example
+##
+## will be zero.
+## @end deftypefn
+## @seealso{eigs}
+
+function [u, s, v, flag] = svds (a, k, sigma, opts)
+
+ if (nargin < 1 || nargin > 4)
+ error ("Incorrect number of arguments");
+ endif
+
+ if (nargin < 4)
+ opts.tol = 1e-10 / sqrt(2);
+ opts.disp = 0;
+ opts.maxit = 300;
+ else
+ if (!isstruct(opts))
+ error("opts must be a structure");
+ endif
+ if (!isfield(opts,"tol"))
+ opts.tol = 1e-10 / sqrt(2);
+ endif
+ endif
+
+ if (nargin < 3 || strcmp(sigma,"L"))
+ if (isreal(a))
+ sigma = "LA";
+ else
+ sigma = "LR";
+ endif
+ elseif (isscalar(sigma) && isreal(sigma))
+ if ((sigma < 0))
+ error ("sigma must be a positive real value");
+ endif
+ else
+ error ("sigma must be a positive real value or the string 'L'");
+ endif
+
+ maxA = max(max(abs(a)));
+ if (maxA == 0)
+ u = eye(m, k);
+ s = zeros(k, k);
+ v = eye(n, k);
+ else
+ [m, n] = size(a);
+ if (nargin < 2)
+ k = min([6, m, n]);
+ else
+ k = min([k, m, n]);
+ endif
+
+ ## Scale everything by the 1-norm to make things more stable.
+ B = a / maxA;
+ Bopts = opts;
+ Bopts.tol = opts.tol / maxA;
+ Bsigma = sigma;
+ if (!ischar(Bsigma))
+ Bsigma = Bsigma / maxA;
+ endif
+
+ if (!ischar(Bsigma) && Bsigma == 0)
+ ## The eigenvalues returns by eigs are symmetric about 0. As we
+ ## are only interested in the positive eigenvalues, we have to
+ ## double k. If sigma is smaller than the smallest singular value
+ ## this can also be an issue. However, we'd like to avoid double
+ ## k for all scalar value of sigma...
+ [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)],
+ 2 * k, Bsigma, Bopts);
+ else
+ [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)],
+ k, Bsigma, Bopts);
+ endif
+ s = diag(s);
+
+ if (ischar(sigma))
+ norma = max(s);
+ else
+ norma = normest(a);
+ endif
+ V = sqrt(2) * V;
+ u = V(1:m,:);
+ v = V(m+1:end,:);
+
+ ## We wish to exclude all eigenvalues that are less than zero as these
+ ## are artifacts of the way the matrix passed to eigs is formed. There
+ ## is also the possibility that the value of sigma chosen is exactly
+ ## a singular value, and in that case we're dead!! So have to rely on
+ ## the warning from eigs. We exclude the singular values which are
+ ## less than or equal to zero to within some tolerance scaled by the
+ ## norm since if we don't we might end up with too many singular
+ ## values. What is appropriate for the tolerance?
+ tol = norma * opts.tol;
+ ind = find(s > tol);
+ if (length(ind) < k)
+ ## Find the zero eigenvalues of B, Ignore the eigenvalues that are
+ ## nominally negative.
+ zind = find(abs(s) <= tol);
+ p = min(length(zind), k-length(ind));
+ ind = [ind;zind(1:p)];
+ elseif (length(ind) > k)
+ ind = ind(1:k);
+ endif
+ u = u(:,ind);
+ s = s(ind);
+ v = v(:,ind);
+
+ if (length(s) < k)
+ warning("returning fewer singular values than requested.");
+ if (!ischar(sigma))
+ warning("try increasing the value of sigma");
+ endif
+ endif
+
+ s = s * maxA;
+ endif
+
+ if (nargout < 2)
+ u = s;
+ else
+ s = diag(s);
+ if (nargout > 3)
+ flag = norm(a*v - u*s, 1) > sqrt(2) * opts.tol * norm(a, 1);
+ endif
+ endif
+endfunction
+
+%!shared n, k, a, u, s, v, opts
+%! n = 100;
+%! k = 7;
+%! a =
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),0.4*n*ones(1,n),ones(1,n-2)]);
+%! %%a =
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]);
+%! [u,s,v] = svd(full(a));
+%! s = diag(s);
+%! [dum, idx] = sort(abs(s));
+%! s = s(idx);
+%! u = u(:,idx);
+%! v = v(:,idx);
+%! randn('state',42)
+%!test
+%! [u2,s2,v2,flag] = svds(a,k);
+%! s2 = diag(s2);
+%! assert(flag,!1);
+%! assert(s(end:-1:end-k+1), s2, 1e-10);
+%!test
+%! [u2,s2,v2,flag] = svds(a,k,0);
+%! s2 = diag(s2);
+%! assert(flag,!1);
+%! assert(s(k:-1:1), s2, 1e-10);
+%!test
+%! idx = floor(n/2);
+%! % Don't put sigma right on a singular value or there are convergence
+%! sigma = 0.99*s(idx) + 0.01*s(idx+1);
+%! [u2,s2,v2,flag] = svds(a,k,sigma);
+%! s2 = diag(s2);
+%! assert(flag,!1);
+%! assert(s((idx+floor(k/2)):-1:(idx-floor(k/2))), s2, 1e-10);
diff --git a/src/DLD-FUNCTIONS/eigs.cc b/src/DLD-FUNCTIONS/eigs.cc
new file mode 100755
--- /dev/null
+++ b/src/DLD-FUNCTIONS/eigs.cc
@@ -0,0 +1,1465 @@
+/*
+
+Copyright (C) 2005 David Bateman
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING. If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "ov.h"
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "quit.h"
+#include "variables.h"
+#include "ov-re-sparse.h"
+#include "ov-cx-sparse.h"
+#include "oct-map.h"
+#include "pager.h"
+#include "unwind-prot.h"
+
+#include "eigs-base.cc"
+
+// Global pointer for user defined function.
+static octave_function *eigs_fcn = 0;
+
+// Have we warned about imaginary values returned from user function?
+static bool warned_imaginary = false;
+
+// Is this a recursive call?
+static int call_depth = 0;
+
+ColumnVector
+eigs_func (const ColumnVector &x, int &eigs_error)
+{
+ ColumnVector retval;
+ octave_value_list args;
+ args(0) = x;
+
+ if (eigs_fcn)
+ {
+ octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args);
+
+ if (error_state)
+ {
+ eigs_error = 1;
+ gripe_user_supplied_eval ("eigs");
+ return retval;
+ }
+
+ if (tmp.length () && tmp(0).is_defined ())
+ {
+ if (! warned_imaginary && tmp(0).is_complex_type ())
+ {
+ warning ("eigs: ignoring imaginary part returned from
user-supplied function");
+ warned_imaginary = true;
+ }
+
+ retval = ColumnVector (tmp(0).vector_value ());
+
+ if (error_state)
+ {
+ eigs_error = 1;
+ gripe_user_supplied_eval ("eigs");
+ }
+ }
+ else
+ {
+ eigs_error = 1;
+ gripe_user_supplied_eval ("eigs");
+ }
+ }
+
+ return retval;
+}
+
+ComplexColumnVector
+eigs_complex_func (const ComplexColumnVector &x, int &eigs_error)
+{
+ ComplexColumnVector retval;
+ octave_value_list args;
+ args(0) = x;
+
+ if (eigs_fcn)
+ {
+ octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args);
+
+ if (error_state)
+ {
+ eigs_error = 1;
+ gripe_user_supplied_eval ("eigs");
+ return retval;
+ }
+
+ if (tmp.length () && tmp(0).is_defined ())
+ {
+ retval = ComplexColumnVector (tmp(0).complex_vector_value ());
+
+ if (error_state)
+ {
+ eigs_error = 1;
+ gripe_user_supplied_eval ("eigs");
+ }
+ }
+ else
+ {
+ eigs_error = 1;
+ gripe_user_supplied_eval ("eigs");
+ }
+ }
+
+ return retval;
+}
+
+DEFUN_DLD (eigs, args, nargout,
+ "-*- texinfo -*-\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a}, @var{k})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a}, @var{k},
@var{sigma})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a}, @var{k},
@var{sigma},@var{opts})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a}, @var{b})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a}, @var{b},
@var{k})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a}, @var{b},
@var{k}, @var{sigma})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{a}, @var{b},
@var{k}, @var{sigma}, @var{opts})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{k})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b}, @var{k})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{k}, @var{sigma})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b}, @var{k}, @var{sigma})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{k}, @var{sigma}, @var{opts})\n\
address@hidden {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b}, @var{k}, @var{sigma}, @var{opts})\n\
address@hidden {Loadable Function} address@hidden, @var{d}]} = eigs (@var{a},
@dots{})\n\
address@hidden {Loadable Function} address@hidden, @var{d}]} = eigs (@var{af},
@var{n}, @dots{})\n\
address@hidden {Loadable Function} address@hidden, @var{d}, @var{flag}]} = eigs
(@var{a}, @dots{})\n\
address@hidden {Loadable Function} address@hidden, @var{d}, @var{flag}]} = eigs
(@var{af}, @var{n}, @dots{})\n\
+Calculate a limited number of eigenvalues and eigenvectors of @var{a},\n\
+based on a selection criteria. The number eigenvalues and eigenvectors to\n\
+calculate is given by @var{k} whose default value is 6.\n\
+\n\
+By default @code{eigs} solve the equation\n\
address@hidden
address@hidden
+$A \\nu = \\lambda \\nu$\n\
address@hidden tex\n\
address@hidden iftex\n\
address@hidden
address@hidden * v = lambda * v}\n\
address@hidden ifinfo\n\
+, where\n\
address@hidden
address@hidden
+$\\lambda$ is a scalar representing one of the eigenvalues, and $\\nu$\n\
address@hidden tex\n\
address@hidden iftex\n\
address@hidden
address@hidden is a scalar representing one of the eigenvalues, and @code{v}\n\
address@hidden ifinfo\n\
+is the corresponding eigenvector. If given the positive definite matrix\n\
address@hidden then @code{eigs} solves the general eigenvalue equation\n\
address@hidden
address@hidden
+$A \\nu = \\lambda B \\nu$\n\
address@hidden tex\n\
address@hidden iftex\n\
address@hidden
address@hidden * v = lambda * B * v}\n\
address@hidden ifinfo\n\
+.\n\
+\n\
+The argument @var{sigma} determines which eigenvalues are returned.\n\
address@hidden can be either a scalar or a string. When @var{sigma} is a
scalar,\n\
+the @var{k} eigenvalues closest to @var{sigma} are returned. If @var{sigma}\n\
+is a string, it must have one of the values\n\
+\n\
address@hidden @asis\n\
address@hidden 'lm'\n\
+Largest magnitude (default).\n\
+\n\
address@hidden 'sm'\n\
+Smallest magnitude.\n\
+\n\
address@hidden 'la'\n\
+Largest Algebraic (valid only for real symmetric problems).\n\
+\n\
address@hidden 'sa'\n\
+Smallest Algebraic (valid only for real symmetric problems).\n\
+\n\
address@hidden 'be'\n\
+Both ends, with one more from the high-end if @var{k} is odd (valid only for\n\
+real symmetric problems).\n\
+\n\
address@hidden 'lr'\n\
+Largest real part (valid only for complex or unsymmetric problems).\n\
+\n\
address@hidden 'sr'\n\
+Smallest real part (valid only for complex or unsymmetric problems).\n\
+\n\
address@hidden 'li'\n\
+Largest imaginary part (valid only for complex or unsymmetric problems).\n\
+\n\
address@hidden 'si'\n\
+Smallest imaginary part (valid only for complex or unsymmetric problems).\n\
address@hidden table\n\
+\n\
+If @var{opts} is given, it is a structure defining some of the options that\n\
address@hidden should use. The fields of the structure @var{opts} are\n\
+\n\
address@hidden @code\n\
address@hidden issym\n\
+If @var{af} is given, then flags whether the function @var{af} defines a\n\
+symmetric problem. It is ignored if @var{a} is given. The default is false.\n\
+\n\
address@hidden isreal\n\
+If @var{af} is given, then flags whether the function @var{af} defines a\n\
+real problem. It is ignored if @var{a} is given. The default is true.\n\
+\n\
address@hidden tol\n\
+Defines the required convergence tolerance, given as @code{tol * norm (A)}.\n\
+The default is @code{eps}.\n\
+\n\
address@hidden maxit\n\
+The maximum number of iterations. The default is 300.\n\
+\n\
address@hidden p\n\
+The number of Lanzcos basis vectors to use. More vectors will result in\n\
+faster convergence, but a larger amount of memory. The optimal value of 'p'\n\
+is problem dependent and should be in the range @var{k} to @var{n}. The\n\
+default value is @code{2 * @var{k}}.\n\
+\n\
address@hidden v0\n\
+The starting vector for the computation. The default is to have @sc{Arpack}\n\
+randomly generate a starting vector.\n\
+\n\
address@hidden disp\n\
+The level of diagnostic printout. If @code{disp} is 0 then there is no\n\
+printout. The default value is 1.\n\
+\n\
address@hidden cholB\n\
+Flag if @code{chol (@var{b})} is passed rather than @var{b}. The default is\n\
+false.\n\
+\n\
address@hidden permB\n\
+The permutation vector of the Cholesky factorization of @var{b} if\n\
address@hidden is true. That is @code{chol ( @var{b} (permB, permB))}. The\n\
+default is @code{1:@var{n}}.\n\
+\n\
address@hidden table\n\
+\n\
+It is also possible to represent @var{a} by a function denoted @var{af}.\n\
address@hidden must be followed by a scalar argument @var{n} defining the
length\n\
+of the vector argument accepted by @var{af}. @var{af} can be passed either\n\
+as an inline function, function handle or as a string. In the case where\n\
address@hidden is passed as a string, the name of the string defines the
function\n\
+to use.\n\
+\n\
address@hidden is a function of the form @code{function y = af (x), y =
@dots{};\n\
+endfunction}, where the required return value of @var{af} is determined by\n\
+the value of @var{sigma}, and are\n\
+\n\
address@hidden @code\n\
address@hidden A * x\n\
+If @var{sigma} is not given or is a string other than 'sm'.\n\
+\n\
address@hidden A \\ x\n\
+If @var{sigma} is 'sm'.\n\
+\n\
address@hidden (A - sigma * I) \\ x\n\
+for standard eigenvalue problem, where @code{I} is the identity matrix of\n\
+the same size as @code{A}. If @var{sigma} is zero, this reduces the\n\
address@hidden \\ x}.\n\
+\n\
address@hidden (A - sigma * B) \\ x\n\
+for the general eigenvalue problem.\n\
address@hidden table\n\
+\n\
+The return arguments of @code{eigs} depends on the number of return\n\
+arguments. With a single return argument, a vector @var{d} of length @var{k}\n\
+is returned, represent the @var{k} eigenvalues that have been found. With
two\n\
+return arguments, @var{v} is a @address@hidden matrix whose columns are\n\
+the @var{k} eigenvectors corresponding to the returned eigenvalues. The\n\
+eigenvalues themselves are then returned in @var{d} in the form of a\n\
address@hidden@var{k} matrix, where the elements on the diagonal are the\n\
+eigenvalues.\n\
+\n\
+Given a third return argument @var{flag}, @code{eigs} also returns the
status\n\
+of the convergence. If @var{flag} is 0, then all eigenvalues have converged,\n\
+otherwise not.\n\
+\n\
+This function is based on the @sc{Arpack} package, written by R Lehoucq,\n\
+K Maschhoff, D Sorensen and C Yang. For more information see\n\
address@hidden://www.caam.rice.edu/software/ARPACK/}.\n\
+\n\
address@hidden deftypefn\n\
address@hidden, svds}")
+{
+ octave_value_list retval;
+#ifdef HAVE_ARPACK
+ int nargin = args.length ();
+ std::string fcn_name;
+ octave_idx_type n = 0;
+ octave_idx_type k = 6;
+ Complex sigma = 0.;
+ double sigmar, sigmai;
+ bool have_sigma = false;
+ std::string typ = "LM";
+ Matrix amm, bmm, bmt;
+ ComplexMatrix acm, bcm, bct;
+ SparseMatrix asmm, bsmm, bsmt;
+ SparseComplexMatrix ascm, bscm, bsct;
+ int b_arg = 0;
+ bool have_b = false;
+ bool have_a_fun = false;
+ bool a_is_complex = false;
+ bool b_is_complex = false;
+ bool symmetric = false;
+ bool cholB = false;
+ bool a_is_sparse = false;
+ ColumnVector permB;
+ int arg_offset = 0;
+ double tol = DBL_EPSILON;
+ int maxit = 300;
+ int disp = 0;
+ octave_idx_type p = -1;
+ ColumnVector resid;
+ ComplexColumnVector cresid;
+ octave_idx_type info = 1;
+ char bmat = 'I';
+
+ warned_imaginary = false;
+
+ unwind_protect::begin_frame ("Feigs");
+
+ unwind_protect_int (call_depth);
+ call_depth++;
+
+ if (call_depth > 1)
+ {
+ error ("eigs: invalid recursive call");
+ if (fcn_name.length())
+ clear_function (fcn_name);
+ unwind_protect::run_frame ("Feigs");
+ return retval;
+ }
+
+ if (nargin == 0)
+ print_usage ();
+ else if (args(0).is_function_handle () || args(0).is_inline_function () ||
+ args(0).is_string ())
+ {
+ if (args(0).is_string ())
+ {
+ std::string name = args(0).string_value ();
+ std::string fname = "function y = ";
+ fcn_name = unique_symbol_name ("__eigs_fcn_");
+ fname.append (fcn_name);
+ fname.append ("(x) y = ");
+ eigs_fcn = extract_function (args(0), "eigs", fcn_name, fname,
+ "; endfunction");
+ }
+ else
+ eigs_fcn = args(0).function_value ();
+
+ if (!eigs_fcn)
+ error ("eigs: unknown function");
+
+ if (nargin < 2)
+ error ("eigs: incorrect number of arguments");
+ else
+ {
+ n = args(1).nint_value ();
+ arg_offset = 1;
+ have_a_fun = true;
+ }
+ }
+ else
+ {
+ if (args(0).is_complex_type ())
+ {
+ if (args(0).is_sparse_type ())
+ {
+ ascm = (args(0).sparse_complex_matrix_value());
+ a_is_sparse = true;
+ }
+ else
+ acm = (args(0).complex_matrix_value());
+ a_is_complex = true;
+ symmetric = false; // ARAPACK doesn't special case complex symmetric
+ }
+ else
+ {
+ if (args(0).is_sparse_type ())
+ {
+ asmm = (args(0).sparse_matrix_value());
+ a_is_sparse = true;
+ symmetric = asmm.is_symmetric();
+ }
+ else
+ {
+ amm = (args(0).matrix_value());
+ symmetric = amm.is_symmetric();
+ }
+ }
+
+ // Note hold off reading B till later to avoid issues of double
+ // copies of the matrix if B is full/real while A is complex..
+ if (!error_state && nargin > 1 &&
+ !(args(1).is_real_scalar ()))
+ if (args(1).is_complex_type ())
+ {
+ b_arg = 1+arg_offset;
+ have_b = true;
+ bmat = 'G';
+ b_is_complex = true;
+ arg_offset++;
+ }
+ else
+ {
+ b_arg = 1+arg_offset;
+ have_b = true;
+ bmat = 'G';
+ arg_offset++;
+ }
+ }
+
+ if (!error_state && nargin > (1+arg_offset))
+ k = args(1+arg_offset).nint_value ();
+
+ if (!error_state && nargin > (2+arg_offset))
+ if (args(2+arg_offset).is_real_scalar () ||
+ args(2+arg_offset).is_complex_scalar ())
+ {
+ sigma = args(2+arg_offset).complex_value ();
+ have_sigma = true;
+ }
+ else if (args(2+arg_offset).is_string ())
+ {
+ typ = args(2+arg_offset).string_value ();
+
+ // Use STL function to convert to upper case
+ transform (typ.begin (), typ.end (), typ.begin (), toupper);
+
+ sigma = 0.;
+ }
+ else
+ error ("eigs: sigma must be a scalar or a string");
+
+ sigmar = std::real (sigma);
+ sigmai = std::imag (sigma);
+
+ if (!error_state && nargin > (3+arg_offset))
+ if (args(3+arg_offset).is_map ())
+ {
+ Octave_map map = args(3+arg_offset).map_value ();
+
+ // issym is ignored if A is not a function
+ if (map.contains ("issym"))
+ if (have_a_fun)
+ symmetric = map.contents ("issym")(0).double_value () != 0.;
+
+ // isreal is ignored if A is not a function
+ if (map.contains ("isreal"))
+ if (have_a_fun)
+ a_is_complex = ! (map.contents ("isreal")(0).double_value () != 0.);
+
+ if (map.contains ("tol"))
+ tol = map.contents ("tol")(0).double_value ();
+
+ if (map.contains ("maxit"))
+ maxit = map.contents ("maxit")(0).nint_value ();
+
+ if (map.contains ("p"))
+ p = map.contents ("p")(0).nint_value ();
+
+ if (map.contains ("v0"))
+ {
+ if (a_is_complex || b_is_complex)
+ cresid = ComplexColumnVector
+ (map.contents ("v0")(0).complex_vector_value());
+ else
+ resid = ColumnVector (map.contents ("v0")(0).vector_value());
+ }
+
+ if (map.contains ("disp"))
+ disp = map.contents ("disp")(0).nint_value ();
+
+ if (map.contains ("cholB"))
+ cholB = map.contents ("cholB")(0).double_value () != 0.;
+
+ if (map.contains ("permB"))
+ permB = ColumnVector (map.contents ("permB")(0).vector_value ())
+ - 1.0;
+ }
+ else
+ error ("eigs: options argument must be a structure");
+
+ if (nargin > (4+arg_offset))
+ error ("eigs: incorrect number of arguments");
+
+ if (have_b)
+ {
+ if (a_is_complex || b_is_complex)
+ {
+ if (a_is_sparse)
+ bscm = args(b_arg).sparse_complex_matrix_value ();
+ else
+ bcm = args(b_arg).complex_matrix_value ();
+ }
+ else
+ {
+ if (a_is_sparse)
+ bsmm = args(b_arg).sparse_matrix_value ();
+ else
+ bmm = args(b_arg).matrix_value ();
+ }
+ }
+
+ // Mode 1 for SM mode seems unstable for some reason.
+ // Use Mode 3 instead, with sigma = 0.
+ if (!error_state && !have_sigma && typ == "SM")
+ have_sigma = true;
+
+ if (!error_state)
+ {
+ octave_idx_type nconv;
+ if (a_is_complex || b_is_complex)
+ {
+ ComplexMatrix eig_vec;
+ ComplexColumnVector eig_val;
+
+
+ if (have_a_fun)
+ nconv = EigsComplexNonSymmetricFunc
+ (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, eig_val,
+ cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+ else if (have_sigma)
+ if (a_is_sparse)
+ nconv = EigsComplexNonSymmetricMatrixShift
+ (ascm, sigma, k, p, info, eig_vec, eig_val, bscm, permB,
+ cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else
+ nconv = EigsComplexNonSymmetricMatrixShift
+ (acm, sigma, k, p, info, eig_vec, eig_val, bcm, permB, cresid,
+ octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+ else
+ if (a_is_sparse)
+ nconv = EigsComplexNonSymmetricMatrix
+ (ascm, typ, k, p, info, eig_vec, eig_val, bscm, permB, cresid,
+ octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+ else
+ nconv = EigsComplexNonSymmetricMatrix
+ (acm, typ, k, p, info, eig_vec, eig_val, bcm, permB, cresid,
+ octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+
+ if (nargout < 2)
+ retval (0) = eig_val;
+ else
+ {
+ retval(2) = double (info);
+ retval(1) = ComplexDiagMatrix (eig_val);
+ retval(0) = eig_vec;
+ }
+ }
+ else if (sigmai != 0.)
+ {
+ // Promote real problem to a complex one.
+ ComplexMatrix eig_vec;
+ ComplexColumnVector eig_val;
+
+ if (have_a_fun)
+ nconv = EigsComplexNonSymmetricFunc
+ (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, eig_val,
+ cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+ else
+ if (a_is_sparse)
+ nconv = EigsComplexNonSymmetricMatrixShift
+ (SparseComplexMatrix (asmm), sigma, k, p, info, eig_vec,
+ eig_val, SparseComplexMatrix (bsmm), permB, cresid,
+ octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+ else
+ nconv = EigsComplexNonSymmetricMatrixShift
+ (ComplexMatrix (amm), sigma, k, p, info, eig_vec,
+ eig_val, ComplexMatrix (bmm), permB, cresid,
+ octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+
+ if (nargout < 2)
+ retval (0) = eig_val;
+ else
+ {
+ retval(2) = double (info);
+ retval(1) = ComplexDiagMatrix (eig_val);
+ retval(0) = eig_vec;
+ }
+ }
+ else
+ {
+ if (symmetric)
+ {
+ Matrix eig_vec;
+ ColumnVector eig_val;
+
+ if (have_a_fun)
+ nconv = EigsRealSymmetricFunc
+ (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else if (have_sigma)
+ if (a_is_sparse)
+ nconv = EigsRealSymmetricMatrixShift
+ (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else
+ nconv = EigsRealSymmetricMatrixShift
+ (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else
+ if (a_is_sparse)
+ nconv = EigsRealSymmetricMatrix
+ (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else
+ nconv = EigsRealSymmetricMatrix
+ (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+
+ if (nargout < 2)
+ retval (0) = eig_val;
+ else
+ {
+ retval(2) = double (info);
+ retval(1) = DiagMatrix (eig_val);
+ retval(0) = eig_vec;
+ }
+ }
+ else
+ {
+ ComplexMatrix eig_vec;
+ ComplexColumnVector eig_val;
+
+ if (have_a_fun)
+ nconv = EigsRealNonSymmetricFunc
+ (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else if (have_sigma)
+ if (a_is_sparse)
+ nconv = EigsRealNonSymmetricMatrixShift
+ (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else
+ nconv = EigsRealNonSymmetricMatrixShift
+ (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else
+ if (a_is_sparse)
+ nconv = EigsRealNonSymmetricMatrix
+ (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+ else
+ nconv = EigsRealNonSymmetricMatrix
+ (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
+ resid, octave_stdout, tol, (nargout > 1), cholB, disp,
+ maxit);
+
+ if (nargout < 2)
+ retval (0) = eig_val;
+ else
+ {
+ retval(2) = double (info);
+ retval(1) = ComplexDiagMatrix (eig_val);
+ retval(0) = eig_vec;
+ }
+ }
+ }
+
+ if (nconv <= 0)
+ warning ("eigs: None of the %d requested eigenvalues converged", k);
+ else if (nconv < k)
+ warning ("eigs: Only %d of the %d requested eigenvalues converged",
+ nconv, k);
+ }
+
+ if (! fcn_name.empty ())
+ clear_function (fcn_name);
+
+ unwind_protect::run_frame ("Feigs");
+#else
+ error ("eigs: not available in this version of Octave");
+#endif
+
+ return retval;
+}
+
+/* #### SPARSE MATRIX VERSIONS #### */
+
+/*
+
+%% Real positive definite tests, n must be even
+%!shared n, k, A, d0, d2
+%! n = 20;
+%! k = 4;
+%! A =
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)]);
+%! d0 = eig (A);
+%! d2 = sort(d0);
+%! [dum, idx] = sort( abs(d0));
+%! d0 = d0(idx);
+%!test
+%! d1 = eigs (A, k);
+%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
+%!test
+%! d1 = eigs (A,k+1);
+%! assert (d1, d0(end:-1:(end-k)),1e-12);
+%!test
+%! d1 = eigs (A, k, 'lm');
+%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sm');
+%! assert (d1, d0(k:-1:1), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'la');
+%! assert (d1, d2(end:-1:(end-k+1)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sa');
+%! assert (d1, d2(1:k), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'be');
+%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-12);
+%!test
+%! d1 = eigs (A, k+1, 'be');
+%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-12);
+%!test
+%! d1 = eigs (A, k, 4.1);
+%! [dum,idx0] = sort (abs(d0 - 4.1));
+%! [dum,idx1] = sort (abs(d1 - 4.1));
+%! assert (d1(idx1), d0(idx0(1:k)), 1e-12);
+%!test
+%! d1 = eigs(A, speye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12);
+%!test
+%! fn = @(x) A * x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'lm', opts);
+%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
+%!test
+%! fn = @(x) A \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (d1, d0(k:-1:1), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (d1, eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'la');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sa');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'be');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+
+%% Real unsymmetric tests
+%!shared n, k, A, d0
+%! n = 20;
+%! k = 4;
+%! A =
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]);
+%! d0 = eig (A);
+%! [dum, idx] = sort(abs(d0));
+%! d0 = d0(idx);
+%!test
+%! d1 = eigs (A, k);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A,k+1);
+%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12);
+%!test
+%! d1 = eigs (A, k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sm');
+%! assert (abs(d1), abs(d0(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'lr');
+%! [dum, idx] = sort (real(d0));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sr');
+%! [dum, idx] = sort (real(abs(d0)));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'li');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'si');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 4.1);
+%! [dum,idx0] = sort (abs(d0 - 4.1));
+%! [dum,idx1] = sort (abs(d1 - 4.1));
+%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12);
+%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12);
+%!test
+%! d1 = eigs(A, speye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12);
+%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))),
1e-12);
+%!test
+%! fn = @(x) A * x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'lr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'li');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'si');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+
+%% Complex hermitian tests
+%!shared n, k, A, d0
+%! n = 20;
+%! k = 4;
+%! A =
sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)]);
+%! d0 = eig (A);
+%! [dum, idx] = sort(abs(d0));
+%! d0 = d0(idx);
+%!test
+%! d1 = eigs (A, k);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A,k+1);
+%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12);
+%!test
+%! d1 = eigs (A, k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sm');
+%! assert (abs(d1), abs(d0(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'lr');
+%! [dum, idx] = sort (real(abs(d0)));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sr');
+%! [dum, idx] = sort (real(abs(d0)));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'li');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'si');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 4.1);
+%! [dum,idx0] = sort (abs(d0 - 4.1));
+%! [dum,idx1] = sort (abs(d1 - 4.1));
+%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12);
+%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12);
+%!test
+%! d1 = eigs(A, speye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, speye(n), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
+%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12);
+%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))),
1e-12);
+%!test
+%! fn = @(x) A * x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'lr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'li');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'si');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/* #### FULL MATRIX VERSIONS #### */
+
+/*
+
+%% Real positive definite tests, n must be even
+%!shared n, k, A, d0, d2
+%! n = 20;
+%! k = 4;
+%! A =
full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)]));
+%! d0 = eig (A);
+%! d2 = sort(d0);
+%! [dum, idx] = sort( abs(d0));
+%! d0 = d0(idx);
+%!test
+%! d1 = eigs (A, k);
+%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
+%!test
+%! d1 = eigs (A,k+1);
+%! assert (d1, d0(end:-1:(end-k)),1e-12);
+%!test
+%! d1 = eigs (A, k, 'lm');
+%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sm');
+%! assert (d1, d0(k:-1:1), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'la');
+%! assert (d1, d2(end:-1:(end-k+1)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sa');
+%! assert (d1, d2(1:k), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'be');
+%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-12);
+%!test
+%! d1 = eigs (A, k+1, 'be');
+%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-12);
+%!test
+%! d1 = eigs (A, k, 4.1);
+%! [dum,idx0] = sort (abs(d0 - 4.1));
+%! [dum,idx1] = sort (abs(d1 - 4.1));
+%! assert (d1(idx1), d0(idx0(1:k)), 1e-12);
+%!test
+%! d1 = eigs(A, eye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!assert (eigs(A,k,4.1), eigs(A,eye(n),k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!assert (eigs(A,k,4.1), eigs(A,eye(n),k,4.1), 1e-12);
+%!test
+%! fn = @(x) A * x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'lm', opts);
+%! assert (d1, d0(end:-1:(end-k+1)), 1e-12);
+%!test
+%! fn = @(x) A \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (d1, d0(k:-1:1), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 1; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (d1, eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'la');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sa');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'be');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+
+%% Real unsymmetric tests
+%!shared n, k, A, d0
+%! n = 20;
+%! k = 4;
+%! A =
full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]));
+%! d0 = eig (A);
+%! [dum, idx] = sort(abs(d0));
+%! d0 = d0(idx);
+%!test
+%! d1 = eigs (A, k);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A,k+1);
+%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12);
+%!test
+%! d1 = eigs (A, k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sm');
+%! assert (abs(d1), abs(d0(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'lr');
+%! [dum, idx] = sort (real(d0));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sr');
+%! [dum, idx] = sort (real(abs(d0)));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'li');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'si');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 4.1);
+%! [dum,idx0] = sort (abs(d0 - 4.1));
+%! [dum,idx1] = sort (abs(d1 - 4.1));
+%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12);
+%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12);
+%!test
+%! d1 = eigs(A, eye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12);
+%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(n),k,4.1))), 1e-12);
+%!test
+%! fn = @(x) A * x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 1;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'lr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'li');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'si');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+
+%% Complex hermitian tests
+%!shared n, k, A, d0
+%! n = 20;
+%! k = 4;
+%! A =
full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)]));
+%! d0 = eig (A);
+%! [dum, idx] = sort(abs(d0));
+%! d0 = d0(idx);
+%!test
+%! d1 = eigs (A, k);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A,k+1);
+%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12);
+%!test
+%! d1 = eigs (A, k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sm');
+%! assert (abs(d1), abs(d0(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'lr');
+%! [dum, idx] = sort (real(abs(d0)));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'sr');
+%! [dum, idx] = sort (real(abs(d0)));
+%! d2 = d0(idx);
+%! assert (real(d1), real(d2(1:k)), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'li');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 'si');
+%! [dum, idx] = sort (imag(abs(d0)));
+%! d2 = d0(idx);
+%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12);
+%!test
+%! d1 = eigs (A, k, 4.1);
+%! [dum,idx0] = sort (abs(d0 - 4.1));
+%! [dum,idx1] = sort (abs(d1 - 4.1));
+%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12);
+%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12);
+%!test
+%! d1 = eigs(A, eye(n), k, 'lm');
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! d1 = eigs(A, eye(n), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
+%!test
+%! opts.cholB=true;
+%! q = [2:n,1];
+%! opts.permB=q;
+%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts);
+%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
+%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
+%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12);
+%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(n),k,4.1))), 1e-12);
+%!test
+%! fn = @(x) A * x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 'lm', opts);
+%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
+%!test
+%! fn = @(x) A \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 'sm', opts);
+%! assert (abs(d1), d0(1:k), 1e-12);
+%!test
+%! fn = @(x) (A - 4.1 * eye(n)) \ x;
+%! opts.issym = 0; opts.isreal = 0;
+%! d1 = eigs (fn, n, k, 4.1, opts);
+%! assert (abs(d1), eigs(A,k,4.1), 1e-12);
+%!test
+%! [v1,d1] = eigs(A, k, 'lm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sm');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'lr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'sr');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'li');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+%!test
+%! [v1,d1] = eigs(A, k, 'si');
+%! d1 = diag(d1);
+%! for i=1:k
+%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12)
+%! endfor
+
+*/
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/src/Makefile.in b/src/Makefile.in
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -65,7 +65,7 @@
DLD_XSRC := amd.cc balance.cc besselj.cc betainc.cc bsxfun.cc cellfun.cc \
chol.cc ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
dasrt.cc dassl.cc det.cc dispatch.cc dlmread.cc dmperm.cc eig.cc \
- expm.cc fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc \
+ eigs.cc expm.cc fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc \
fltk_backend.cc \
gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
givens.cc hess.cc hex2num.cc inv.cc kron.cc lookup.cc lsode.cc \
@@ -650,6 +650,7 @@
convhulln.oct: OCT_LINK_DEPS += $(QHULL_LIBS)
__delaunayn__.oct: OCT_LINK_DEPS += $(QHULL_LIBS)
__voronoi__.oct: OCT_LINK_DEPS += $(QHULL_LIBS)
+eigs.oct: OCT_LINK_DEPS += $(ARPACK_LIBS)
regexp.oct: OCT_LINK_DEPS += $(REGEX_LIBS)
urlwrite.oct: OCT_LINK_DEPS += $(CURL_LIBS)
__glpk__.oct: OCT_LINK_DEPS += $(GLPK_LIBS)