Mercurial > hg > octave-terminal
changeset 8417:654bcfb937bf
Add the eigs and svds functions
author | David Bateman <dbateman@free.fr> |
---|---|
date | Tue, 23 Dec 2008 08:28:23 +0100 |
parents | 0a56e5c21c29 |
children | 679c22082ac7 |
files | ChangeLog Makeconf.in NEWS configure.in doc/ChangeLog doc/interpreter/sparse.txi liboctave/ChangeLog liboctave/Makefile.in liboctave/eigs-base.cc scripts/ChangeLog scripts/sparse/Makefile.in scripts/sparse/svds.m src/ChangeLog src/DLD-FUNCTIONS/eigs.cc src/Makefile.in |
diffstat | 15 files changed, 5537 insertions(+), 7 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,10 @@ +2008-12-23 David Bateman <dbateman@free.fr> + + * configure.in: Add configuration test for ARPACK. + * Makeconf.in (ARPACK_LIBS): Add variable with location of ARPACK + library. + * NEWS: Document that eigs and svds were moved to Octaave. + 2008-10-29 Jaroslav Hajek <highegg@gmail.com> * configure.in: Remove the OCTAVE_LOCAL_BUFFER stuff (moved to
--- a/Makeconf.in +++ b/Makeconf.in @@ -237,6 +237,7 @@ CHOLMOD_LIBS = @CHOLMOD_LIBS@ CXSPARSE_LIBS = @CXSPARSE_LIBS@ OPENGL_LIBS = @OPENGL_LIBS@ +ARPACK_LIBS = @ARPACK_LIBS@ LIBS = @LIBS@ USE_64_BIT_IDX_T = @USE_64_BIT_IDX_T@
--- a/NEWS +++ b/NEWS @@ -84,6 +84,9 @@ ** The imwrite and imread function have been included in Octave based on the GraphicsMagick library. + ** The eigs and svds functions have been included in Octave based on + the ARPACK library. + ** Special treatment in the parser of expressions like "a' * b". In these cases the transpose is no longer explicitly formed and BLAS libraries are called with the transpose flagged. This significantly
--- a/configure.in +++ b/configure.in @@ -1044,6 +1044,27 @@ AC_MSG_WARN($warn_cxsparse) fi +ARPACK_LIBS= +AC_SUBST(ARPACK_LIBS) + +AC_ARG_WITH(arpack, + [AS_HELP_STRING([--without-arpack], + [don't use ARPACK, disable some sparse functionality])], + with_arpack=$withval, with_arpack=yes) + +warn_arpack="arpack not found. This will result in a lack of the eigs function." +if test "$with_arpack" = yes; then + with_arpack=no + AC_CHECK_LIB(arpack, F77_FUNC(dseupd,DSEUPD), [ARPACK_LIBS="-larpack"; with_arpack=yes], , $LAPACK_LIBS $FLIBS) + if test "$with_arpack" = yes; then + AC_DEFINE(HAVE_ARPACK, 1, [Define if the ARPACK library is used.]) + warn_arpack= + fi +fi +if test -n "$warn_arpack"; then + AC_MSG_WARN($warn_arpack) +fi + ### Handle shared library options. ### Enable creation of static libraries. @@ -2012,6 +2033,7 @@ CCOLAMD libraries: $CCOLAMD_LIBS CHOLMOD libraries: $CHOLMOD_LIBS CXSPARSE libraries: $CXSPARSE_LIBS + ARPACK libraries: $ARPACK_LIBS HDF5 libraries: $HDF5_LIBS CURL libraries: $CURL_LIBS REGEX libraries: $REGEX_LIBS @@ -2123,6 +2145,11 @@ warn_msg_printed=true fi +if test -n "$warn_arpack"; then + AC_MSG_WARN($warn_arpack) + warn_msg_printed=true +fi + if test -n "$warn_curl"; then AC_MSG_WARN($warn_curl) warn_msg_printed=true
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,7 @@ +2008-12-23 David Bateman <dbateman@free.fr> + + * interpreter/sparse.txi: Document the eigs and svds functions. + 2008-12-02 Thorsten Meyer <thorsten.meyier@gmx.de> * interpreter/container.txi, interpreter/strings.txi:
--- a/doc/interpreter/sparse.txi +++ b/doc/interpreter/sparse.txi @@ -473,9 +473,8 @@ @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, @dfn{symrcm} @item Linear algebra: - @dfn{matrix_type}, @dfn{normest}, @dfn{condest}, @dfn{sprank} - @dfn{spaugment} -@c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now + @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, @dfn{normest}, @dfn{sprank}, + @dfn{spaugment}, @dfn{svds} @item Iterative techniques: @dfn{luinc}, @dfn{pcg}, @dfn{pcr} @@ -832,6 +831,15 @@ @DOCSTRING(spaugment) +Finally, the function @code{eigs} can be used to calculate a limited +number of eigenvalues and eigenvectors based on a selection criteria +and likewise for @code{svds} which calculates a limited number of +singular values and vectors. + +@DOCSTRING(eigs) + +@DOCSTRING(svds) + @node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse Matrices @section Iterative Techniques applied to sparse matrices
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +2008-12-23 David Bateman <dbateman@free.fr> + + * eigs-base.cc: New file with template wrapper for ARPACK. + * Makefile.in (TEMPLATE_SRC): Add it here. + 2008-12-16 Jaroslav Hajek <highegg@gmail.com> * Array.cc (rec_permute_helper): New class.
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -103,7 +103,7 @@ $(VX_OP_INC) \ $(SPARSE_MX_OP_INC) -TEMPLATE_SRC := Array.cc ArrayN.cc DiagArray2.cc \ +TEMPLATE_SRC := Array.cc ArrayN.cc eigs-base.cc DiagArray2.cc \ MArray.cc MArray2.cc MArrayN.cc MDiagArray2.cc \ base-lu.cc oct-sort.cc sparse-base-lu.cc \ sparse-base-chol.cc sparse-dmsolve.cc
new file mode 100755 --- /dev/null +++ b/liboctave/eigs-base.cc @@ -0,0 +1,3768 @@ +/* + +Copyright (C) 2005 David Bateman + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <cfloat> +#include <cmath> +#include <vector> +#include <iostream> + +#include "f77-fcn.h" +#include "quit.h" +#include "SparsedbleLU.h" +#include "SparseCmplxLU.h" +#include "dSparse.h" +#include "CSparse.h" +#include "MatrixType.h" +#include "SparsedbleCHOL.h" +#include "SparseCmplxCHOL.h" +#include "oct-rand.h" +#include "dbleCHOL.h" +#include "CmplxCHOL.h" +#include "dbleLU.h" +#include "CmplxLU.h" + +#ifdef HAVE_ARPACK +typedef ColumnVector (*EigsFunc) (const ColumnVector &x, int &eigs_error); +typedef ComplexColumnVector (*EigsComplexFunc) + (const ComplexColumnVector &x, int &eigs_error); + +// Arpack and blas fortran functions we call. +extern "C" +{ + F77_RET_T + F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const double&, + double*, const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dseupd, DSEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, + int*, double*, double*, + const octave_idx_type&, const double&, + F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + const double&, double*, const octave_idx_type&, + double*, const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dnaupd, DNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + octave_idx_type&, const double&, + double*, const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dneupd, DNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, + int*, double*, double*, + double*, const octave_idx_type&, const double&, + const double&, double*, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + octave_idx_type&, const double&, double*, + const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (znaupd, ZNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const double&, + Complex*, const octave_idx_type&, Complex*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, Complex*, Complex*, + const octave_idx_type&, double *, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zneupd, ZNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, + int*, Complex*, Complex*, + const octave_idx_type&, const Complex&, + Complex*, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const double&, + Complex*, const octave_idx_type&, Complex*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, Complex*, Complex*, + const octave_idx_type&, double *, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const double&, + const double*, const octave_idx_type&, const double*, + const octave_idx_type&, const double&, double*, + const octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + + + F77_RET_T + F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const Complex&, + const Complex*, const octave_idx_type&, const Complex*, + const octave_idx_type&, const Complex&, Complex*, const octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + +} + + +#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) +static octave_idx_type +lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&); + +static octave_idx_type +lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, + ComplexMatrix&); + +static octave_idx_type +lusolve (const Matrix&, const Matrix&, Matrix&); + +static octave_idx_type +lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&); + +static ComplexMatrix +ltsolve (const SparseComplexMatrix&, const ColumnVector&, + const ComplexMatrix&); + +static Matrix +ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,); + +static ComplexMatrix +ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); + +static Matrix +ltsolve (const Matrix&, const ColumnVector&, const Matrix&,); + +static ComplexMatrix +utsolve (const SparseComplexMatrix&, const ColumnVector&, const ComplexMatrix&); + +static Matrix +utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); + +static ComplexMatrix +utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); + +static Matrix +utsolve (const Matrix&, const ColumnVector&, const Matrix&); + +#endif + +template <class M, class SM> +static octave_idx_type +lusolve (const SM& L, const SM& U, M& m) +{ + octave_idx_type err = 0; + double rcond; + MatrixType utyp (MatrixType::Upper); + + // Sparse L is lower triangular, Dense L is permuted lower triangular!!! + m = L.solve (m, err, rcond, 0); + if (err) + return err; + + m = U.solve (utyp, m, err, rcond, 0); + + return err; +} + +template <class SM, class M> +static M +ltsolve (const SM& L, const ColumnVector& Q, const M& m) +{ + octave_idx_type n = L.cols(); + octave_idx_type b_nc = m.cols(); + octave_idx_type err = 0; + double rcond; + MatrixType ltyp (MatrixType::Lower); + M tmp = L.solve (ltyp, m, err, rcond, 0); + M retval; + const double* qv = Q.fortran_vec(); + + if (!err) + { + retval.resize (n, b_nc); + for (octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < n; i++) + retval.elem(static_cast<octave_idx_type>(qv[i]), j) = + tmp.elem(i,j); + } + } + + return retval; +} + +template <class SM, class M> +static M +utsolve (const SM& U, const ColumnVector& Q, const M& m) +{ + octave_idx_type n = U.cols(); + octave_idx_type b_nc = m.cols(); + octave_idx_type err = 0; + double rcond; + MatrixType utyp (MatrixType::Upper); + + M retval (n, b_nc); + const double* qv = Q.fortran_vec(); + for (octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < n; i++) + retval.elem(i,j) = m.elem(static_cast<octave_idx_type>(qv[i]), j); + } + return U.solve (utyp, retval, err, rcond, 0); +} + +static bool +vector_product (const SparseMatrix& m, const double* x, double* y) +{ + octave_idx_type nc = m.cols (); + + for (octave_idx_type j = 0; j < nc; j++) + y[j] = 0.; + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) + y[m.ridx(i)] += m.data(i) * x[j]; + + return true; +} + +static bool +vector_product (const Matrix& m, const double *x, double *y) +{ + octave_idx_type nr = m.rows (); + octave_idx_type nc = m.cols (); + + F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1), + nr, nc, 1.0, m.data (), nr, + x, 1, 0.0, y, 1 + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable error in dgemv"); + return false; + } + else + return true; +} + +static bool +vector_product (const SparseComplexMatrix& m, const Complex* x, + Complex* y) +{ + octave_idx_type nc = m.cols (); + + for (octave_idx_type j = 0; j < nc; j++) + y[j] = 0.; + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) + y[m.ridx(i)] += m.data(i) * x[j]; + + return true; +} + +static bool +vector_product (const ComplexMatrix& m, const Complex *x, Complex *y) +{ + octave_idx_type nr = m.rows (); + octave_idx_type nc = m.cols (); + + F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1), + nr, nc, 1.0, m.data (), nr, + x, 1, 0.0, y, 1 + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable error in zgemv"); + return false; + } + else + return true; +} + +static bool +make_cholb (Matrix& b, Matrix& bt, ColumnVector& permB) +{ + octave_idx_type info; + CHOL fact (b, info); + octave_idx_type n = b.cols(); + + if (info != 0) + return false; + else + { + bt = fact.chol_matrix (); + b = bt.transpose(); + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + return true; + } +} + +static bool +make_cholb (SparseMatrix& b, SparseMatrix& bt, ColumnVector& permB) +{ + octave_idx_type info; + SparseCHOL fact (b, info, false); + + if (fact.P() != 0) + return false; + else + { + b = fact.L(); + bt = b.transpose(); + permB = fact.perm() - 1.0; + return true; + } +} + +static bool +make_cholb (ComplexMatrix& b, ComplexMatrix& bt, ColumnVector& permB) +{ + octave_idx_type info; + ComplexCHOL fact (b, info); + octave_idx_type n = b.cols(); + + if (info != 0) + return false; + else + { + bt = fact.chol_matrix (); + b = bt.hermitian(); + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + return true; + } +} + +static bool +make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt, + ColumnVector& permB) +{ + octave_idx_type info; + SparseComplexCHOL fact (b, info, false); + + if (fact.P() != 0) + return false; + else + { + b = fact.L(); + bt = b.hermitian(); + permB = fact.perm() - 1.0; + return true; + } +} + +static bool +LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b, + bool cholB, const ColumnVector& permB, double sigma, + SparseMatrix &L, SparseMatrix &U, octave_idx_type *P, + octave_idx_type *Q) +{ + bool have_b = ! b.is_empty (); + octave_idx_type n = m.rows(); + + // Caclulate LU decomposition of 'A - sigma * B' + SparseMatrix AminusSigmaB (m); + + if (have_b) + { + if (cholB) + { + if (permB.length()) + { + SparseMatrix tmp(n,n,n); + for (octave_idx_type i = 0; i < n; i++) + { + tmp.xcidx(i) = i; + tmp.xridx(i) = + static_cast<octave_idx_type>(permB(i)); + tmp.xdata(i) = 1; + } + tmp.xcidx(n) = n; + + AminusSigmaB = AminusSigmaB - sigma * tmp * + b.transpose() * b * tmp.transpose(); + } + else + AminusSigmaB = AminusSigmaB - sigma * + b.transpose() * b; + } + else + AminusSigmaB = AminusSigmaB - sigma * b; + } + else + { + SparseMatrix sigmat (n, n, n); + + // Create sigma * speye(n,n) + sigmat.xcidx (0) = 0; + for (octave_idx_type i = 0; i < n; i++) + { + sigmat.xdata(i) = sigma; + sigmat.xridx(i) = i; + sigmat.xcidx(i+1) = i + 1; + } + + AminusSigmaB = AminusSigmaB - sigmat; + } + + SparseLU fact (AminusSigmaB); + + L = fact.L (); + U = fact.U (); + const octave_idx_type *P2 = fact.row_perm (); + const octave_idx_type *Q2 = fact.col_perm (); + + for (octave_idx_type j = 0; j < n; j++) + { + P[j] = P2[j]; + Q[j] = Q2[j]; + } + + // Test condition number of LU decomposition + double minU = octave_NaN; + double maxU = octave_NaN; + for (octave_idx_type j = 0; j < n; j++) + { + double d = 0.; + if (U.xcidx(j+1) > U.xcidx(j) && + U.xridx (U.xcidx(j+1)-1) == j) + d = std::abs (U.xdata (U.xcidx(j+1)-1)); + + if (xisnan (minU) || d < minU) + minU = d; + + if (xisnan (maxU) || d > maxU) + maxU = d; + } + + double rcond = (minU / maxU); + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + (*current_liboctave_warning_handler) + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + (*current_liboctave_warning_handler) + (" an eigenvalue. Convergence is not guaranteed"); + } + + return true; +} + +static bool +LuAminusSigmaB (const Matrix &m, const Matrix &b, + bool cholB, const ColumnVector& permB, double sigma, + Matrix &L, Matrix &U, octave_idx_type *P, + octave_idx_type *Q) +{ + bool have_b = ! b.is_empty (); + octave_idx_type n = m.cols(); + + // Caclulate LU decomposition of 'A - sigma * B' + Matrix AminusSigmaB (m); + + if (have_b) + { + if (cholB) + { + Matrix tmp = sigma * b.transpose() * b; + const double *pB = permB.fortran_vec(); + double *p = AminusSigmaB.fortran_vec(); + + if (permB.length()) + { + for (octave_idx_type j = 0; + j < b.cols(); j++) + for (octave_idx_type i = 0; + i < b.rows(); i++) + *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), + static_cast<octave_idx_type>(pB[j])); + } + else + AminusSigmaB = AminusSigmaB - tmp; + } + else + AminusSigmaB = AminusSigmaB - sigma * b; + } + else + { + double *p = AminusSigmaB.fortran_vec(); + + for (octave_idx_type i = 0; i < n; i++) + p[i*(n+1)] -= sigma; + } + + LU fact (AminusSigmaB); + + L = fact.P().transpose() * fact.L (); + U = fact.U (); + for (octave_idx_type j = 0; j < n; j++) + P[j] = Q[j] = j; + + // Test condition number of LU decomposition + double minU = octave_NaN; + double maxU = octave_NaN; + for (octave_idx_type j = 0; j < n; j++) + { + double d = std::abs (U.xelem(j,j)); + if (xisnan (minU) || d < minU) + minU = d; + + if (xisnan (maxU) || d > maxU) + maxU = d; + } + + double rcond = (minU / maxU); + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + (*current_liboctave_warning_handler) + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + (*current_liboctave_warning_handler) + (" an eigenvalue. Convergence is not guaranteed"); + } + + return true; +} + +static bool +LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b, + bool cholB, const ColumnVector& permB, Complex sigma, + SparseComplexMatrix &L, SparseComplexMatrix &U, + octave_idx_type *P, octave_idx_type *Q) +{ + bool have_b = ! b.is_empty (); + octave_idx_type n = m.rows(); + + // Caclulate LU decomposition of 'A - sigma * B' + SparseComplexMatrix AminusSigmaB (m); + + if (have_b) + { + if (cholB) + { + if (permB.length()) + { + SparseMatrix tmp(n,n,n); + for (octave_idx_type i = 0; i < n; i++) + { + tmp.xcidx(i) = i; + tmp.xridx(i) = + static_cast<octave_idx_type>(permB(i)); + tmp.xdata(i) = 1; + } + tmp.xcidx(n) = n; + + AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b * + tmp.transpose() * sigma; + } + else + AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b; + } + else + AminusSigmaB = AminusSigmaB - sigma * b; + } + else + { + SparseComplexMatrix sigmat (n, n, n); + + // Create sigma * speye(n,n) + sigmat.xcidx (0) = 0; + for (octave_idx_type i = 0; i < n; i++) + { + sigmat.xdata(i) = sigma; + sigmat.xridx(i) = i; + sigmat.xcidx(i+1) = i + 1; + } + + AminusSigmaB = AminusSigmaB - sigmat; + } + + SparseComplexLU fact (AminusSigmaB); + + L = fact.L (); + U = fact.U (); + const octave_idx_type *P2 = fact.row_perm (); + const octave_idx_type *Q2 = fact.col_perm (); + + for (octave_idx_type j = 0; j < n; j++) + { + P[j] = P2[j]; + Q[j] = Q2[j]; + } + + // Test condition number of LU decomposition + double minU = octave_NaN; + double maxU = octave_NaN; + for (octave_idx_type j = 0; j < n; j++) + { + double d = 0.; + if (U.xcidx(j+1) > U.xcidx(j) && + U.xridx (U.xcidx(j+1)-1) == j) + d = std::abs (U.xdata (U.xcidx(j+1)-1)); + + if (xisnan (minU) || d < minU) + minU = d; + + if (xisnan (maxU) || d > maxU) + maxU = d; + } + + double rcond = (minU / maxU); + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + (*current_liboctave_warning_handler) + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + (*current_liboctave_warning_handler) + (" an eigenvalue. Convergence is not guaranteed"); + } + + return true; +} + +static bool +LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b, + bool cholB, const ColumnVector& permB, Complex sigma, + ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P, + octave_idx_type *Q) +{ + bool have_b = ! b.is_empty (); + octave_idx_type n = m.cols(); + + // Caclulate LU decomposition of 'A - sigma * B' + ComplexMatrix AminusSigmaB (m); + + if (have_b) + { + if (cholB) + { + ComplexMatrix tmp = sigma * b.hermitian() * b; + const double *pB = permB.fortran_vec(); + Complex *p = AminusSigmaB.fortran_vec(); + + if (permB.length()) + { + for (octave_idx_type j = 0; + j < b.cols(); j++) + for (octave_idx_type i = 0; + i < b.rows(); i++) + *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), + static_cast<octave_idx_type>(pB[j])); + } + else + AminusSigmaB = AminusSigmaB - tmp; + } + else + AminusSigmaB = AminusSigmaB - sigma * b; + } + else + { + Complex *p = AminusSigmaB.fortran_vec(); + + for (octave_idx_type i = 0; i < n; i++) + p[i*(n+1)] -= sigma; + } + + ComplexLU fact (AminusSigmaB); + + L = fact.P().transpose() * fact.L (); + U = fact.U (); + for (octave_idx_type j = 0; j < n; j++) + P[j] = Q[j] = j; + + // Test condition number of LU decomposition + double minU = octave_NaN; + double maxU = octave_NaN; + for (octave_idx_type j = 0; j < n; j++) + { + double d = std::abs (U.xelem(j,j)); + if (xisnan (minU) || d < minU) + minU = d; + + if (xisnan (maxU) || d > maxU) + maxU = d; + } + + double rcond = (minU / maxU); + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + (*current_liboctave_warning_handler) + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + (*current_liboctave_warning_handler) + (" an eigenvalue. Convergence is not guaranteed"); + } + + return true; +} + +template <class M> +octave_idx_type +EigsRealSymmetricMatrix (const M& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) +{ + M b(_b); + octave_idx_type n = m.cols (); + octave_idx_type mode = 1; + bool have_b = ! b.is_empty(); + bool note3 = false; + char bmat = 'I'; + double sigma = 0.; + M bt; + + if (m.rows() != m.cols()) + { + (*current_liboctave_error_handler) ("eigs: A must be square"); + return -1; + } + if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) + { + (*current_liboctave_error_handler) + ("eigs: B must be square and the same size as A"); + return -1; + } + + if (resid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + resid = ColumnVector (octave_rand::vector(n)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p <= k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k and n"); + return -1; + } + + if (k > n ) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (have_b && cholB && permB.length() != 0) + { + // Check the we really have a permutation vector + if (permB.length() != n) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + else + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } + } + + if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + { + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); + return -1; + } + + if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for real symmetric problem"); + return -1; + } + + if (have_b) + { + // See Note 3 dsaupd + note3 = true; + if (cholB) + { + bt = b; + b = b.transpose(); + if (permB.length() == 0) + { + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + } + } + else + { + if (! make_cholb(b, bt, permB)) + { + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + return -1; + } + } + } + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + octave_idx_type lwork = p * (p + 8); + + OCTAVE_LOCAL_BUFFER (double, v, n * p); + OCTAVE_LOCAL_BUFFER (double, workl, lwork); + OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); + double *presid = resid.fortran_vec (); + + do + { + F77_FUNC (dsaupd, DSAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + if (have_b) + { + Matrix mtmp (n,1); + for (octave_idx_type i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + + mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); + + for (octave_idx_type i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + } + else if (!vector_product (m, workd + iptr(0) - 1, + workd + iptr(1) - 1)) + break; + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + eig_vec.resize (n, k); + double *z = eig_vec.fortran_vec (); + + eig_val.resize (k); + double *d = eig_val.fortran_vec (); + + F77_FUNC (dseupd, DSEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, + F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) + F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dseupd"); + return -1; + } + else + { + if (info2 == 0) + { + octave_idx_type k2 = k / 2; + if (typ != "SM" && typ != "BE") + { + for (octave_idx_type i = 0; i < k2; i++) + { + double dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + } + + if (rvec) + { + if (typ != "SM" && typ != "BE") + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + } + + if (note3) + eig_vec = ltsolve(b, permB, eig_vec); + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in dseupd", info2); + return -1; + } + } + + return ip(4); +} + +template <class M> +octave_idx_type +EigsRealSymmetricMatrixShift (const M& m, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) +{ + M b(_b); + octave_idx_type n = m.cols (); + octave_idx_type mode = 3; + bool have_b = ! b.is_empty(); + std::string typ = "LM"; + + if (m.rows() != m.cols()) + { + (*current_liboctave_error_handler) ("eigs: A must be square"); + return -1; + } + if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) + { + (*current_liboctave_error_handler) + ("eigs: B must be square and the same size as A"); + return -1; + } + + // FIXME: The "SM" type for mode 1 seems unstable though faster!! + //if (! std::abs (sigma)) + // return EigsRealSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val, + // _b, permB, resid, os, tol, rvec, cholB, + // disp, maxit); + + if (resid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + resid = ColumnVector (octave_rand::vector(n)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p <= k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k and n"); + return -1; + } + + if (k > n ) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (have_b && cholB && permB.length() != 0) + { + // Check the we really have a permutation vector + if (permB.length() != n) + { + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + return -1; + } + else + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } + } + + char bmat = 'I'; + if (have_b) + bmat = 'G'; + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + M L, U; + + OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows())); + OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols())); + + if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q)) + return -1; + + octave_idx_type lwork = p * (p + 8); + + OCTAVE_LOCAL_BUFFER (double, v, n * p); + OCTAVE_LOCAL_BUFFER (double, workl, lwork); + OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); + double *presid = resid.fortran_vec (); + + do + { + F77_FUNC (dsaupd, DSAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + if (have_b) + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + vector_product (m, workd+iptr(0)-1, dtmp); + + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = dtmp[P[i]]; + + lusolve (L, U, tmp); + + double *ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + else if (ido == 2) + vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); + else + { + double *ip2 = workd+iptr(2)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + else + { + if (ido == 2) + { + for (octave_idx_type i = 0; i < n; i++) + workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; + } + else + { + double *ip2 = workd+iptr(0)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + eig_vec.resize (n, k); + double *z = eig_vec.fortran_vec (); + + eig_val.resize (k); + double *d = eig_val.fortran_vec (); + + F77_FUNC (dseupd, DSEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, + F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), + k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dseupd"); + return -1; + } + else + { + if (info2 == 0) + { + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + double dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in dseupd", info2); + return -1; + } + } + + return ip(4); +} + +octave_idx_type +EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, + const std::string &_typ, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, ColumnVector &resid, + std::ostream& os, double tol, int rvec, bool cholB, + int disp, int maxit) +{ + std::string typ (_typ); + bool have_sigma = (sigma ? true : false); + char bmat = 'I'; + octave_idx_type mode = 1; + int err = 0; + + if (resid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + resid = ColumnVector (octave_rand::vector(n)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p <= k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k and n"); + return -1; + } + + if (k > n ) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (! have_sigma) + { + if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); + + if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for real symmetric problem"); + return -1; + } + + if (typ == "SM") + { + typ = "LM"; + sigma = 0.; + mode = 3; + } + } + else if (! std::abs (sigma)) + typ = "SM"; + else + { + typ = "LM"; + mode = 3; + } + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + octave_idx_type lwork = p * (p + 8); + + OCTAVE_LOCAL_BUFFER (double, v, n * p); + OCTAVE_LOCAL_BUFFER (double, workl, lwork); + OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); + double *presid = resid.fortran_vec (); + + do + { + F77_FUNC (dsaupd, DSAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + + if (ido == -1 || ido == 1 || ido == 2) + { + double *ip2 = workd + iptr(0) - 1; + ColumnVector x(n); + + for (octave_idx_type i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (octave_idx_type i = 0; i < n; i++) + *ip2++ = y(i); + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + eig_vec.resize (n, k); + double *z = eig_vec.fortran_vec (); + + eig_val.resize (k); + double *d = eig_val.fortran_vec (); + + F77_FUNC (dseupd, DSEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, + F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), + k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dseupd"); + return -1; + } + else + { + if (info2 == 0) + { + octave_idx_type k2 = k / 2; + if (typ != "SM" && typ != "BE") + { + for (octave_idx_type i = 0; i < k2; i++) + { + double dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + } + + if (rvec) + { + if (typ != "SM" && typ != "BE") + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + } + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in dseupd", info2); + return -1; + } + } + + return ip(4); +} + +template <class M> +octave_idx_type +EigsRealNonSymmetricMatrix (const M& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) +{ + M b(_b); + octave_idx_type n = m.cols (); + octave_idx_type mode = 1; + bool have_b = ! b.is_empty(); + bool note3 = false; + char bmat = 'I'; + double sigmar = 0.; + double sigmai = 0.; + M bt; + + if (m.rows() != m.cols()) + { + (*current_liboctave_error_handler) ("eigs: A must be square"); + return -1; + } + if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) + { + (*current_liboctave_error_handler) + ("eigs: B must be square and the same size as A"); + return -1; + } + + if (resid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + resid = ColumnVector (octave_rand::vector(n)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2 + 1; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p < k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k+1 and n"); + return -1; + } + + if (k > n - 1) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n-1).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (have_b && cholB && permB.length() != 0) + { + // Check the we really have a permutation vector + if (permB.length() != n) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + else + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } + } + + if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + { + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); + return -1; + } + + if (typ == "LA" || typ == "SA" || typ == "BE") + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for unsymmetric problem"); + return -1; + } + + if (have_b) + { + // See Note 3 dsaupd + note3 = true; + if (cholB) + { + bt = b; + b = b.transpose(); + if (permB.length() == 0) + { + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + } + } + else + { + if (! make_cholb(b, bt, permB)) + { + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + return -1; + } + } + } + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + octave_idx_type lwork = 3 * p * (p + 2); + + OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); + OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); + OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); + double *presid = resid.fortran_vec (); + + do + { + F77_FUNC (dnaupd, DNAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dnaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + if (have_b) + { + Matrix mtmp (n,1); + for (octave_idx_type i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + + mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); + + for (octave_idx_type i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + } + else if (!vector_product (m, workd + iptr(0) - 1, + workd + iptr(1) - 1)) + break; + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dnaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + Matrix eig_vec2 (n, k + 1); + double *z = eig_vec2.fortran_vec (); + + OCTAVE_LOCAL_BUFFER (double, dr, k + 1); + OCTAVE_LOCAL_BUFFER (double, di, k + 1); + OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); + for (octave_idx_type i = 0; i < k+1; i++) + dr[i] = di[i] = 0.; + + F77_FUNC (dneupd, DNEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, + sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) + F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dneupd"); + return -1; + } + else + { + eig_val.resize (k+1); + Complex *d = eig_val.fortran_vec (); + + if (info2 == 0) + { + octave_idx_type jj = 0; + for (octave_idx_type i = 0; i < k+1; i++) + { + if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) + jj++; + else + d [i-jj] = Complex (dr[i], di[i]); + } + if (jj == 0 && !rvec) + for (octave_idx_type i = 0; i < k; i++) + d[i] = d[i+1]; + + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + + eig_vec.resize (n, k); + octave_idx_type i = 0; + while (i < k) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (i+1) * n; + if (std::imag(eig_val(i)) == 0) + { + for (octave_idx_type j = 0; j < n; j++) + eig_vec(j,i) = + Complex(z[j+off1],0.); + i++; + } + else + { + for (octave_idx_type j = 0; j < n; j++) + { + eig_vec(j,i) = + Complex(z[j+off1],z[j+off2]); + if (i < k - 1) + eig_vec(j,i+1) = + Complex(z[j+off1],-z[j+off2]); + } + i+=2; + } + } + + if (note3) + eig_vec = ltsolve(M (b), permB, eig_vec); + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in dneupd", info2); + return -1; + } + } + + return ip(4); +} + +template <class M> +octave_idx_type +EigsRealNonSymmetricMatrixShift (const M& m, double sigmar, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) +{ + M b(_b); + octave_idx_type n = m.cols (); + octave_idx_type mode = 3; + bool have_b = ! b.is_empty(); + std::string typ = "LM"; + double sigmai = 0.; + + if (m.rows() != m.cols()) + { + (*current_liboctave_error_handler) ("eigs: A must be square"); + return -1; + } + if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) + { + (*current_liboctave_error_handler) + ("eigs: B must be square and the same size as A"); + return -1; + } + + // FIXME: The "SM" type for mode 1 seems unstable though faster!! + //if (! std::abs (sigmar)) + // return EigsRealNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val, + // _b, permB, resid, os, tol, rvec, cholB, + // disp, maxit); + + if (resid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + resid = ColumnVector (octave_rand::vector(n)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2 + 1; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p < k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k+1 and n"); + return -1; + } + + if (k > n - 1) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n-1).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (have_b && cholB && permB.length() != 0) + { + // Check that we really have a permutation vector + if (permB.length() != n) + { + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + return -1; + } + else + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } + } + + char bmat = 'I'; + if (have_b) + bmat = 'G'; + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + M L, U; + + OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows())); + OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols())); + + if (! LuAminusSigmaB(m, b, cholB, permB, sigmar, L, U, P, Q)) + return -1; + + octave_idx_type lwork = 3 * p * (p + 2); + + OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); + OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); + OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); + double *presid = resid.fortran_vec (); + + do + { + F77_FUNC (dnaupd, DNAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + if (have_b) + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + vector_product (m, workd+iptr(0)-1, dtmp); + + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = dtmp[P[i]]; + + lusolve (L, U, tmp); + + double *ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + else if (ido == 2) + vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); + else + { + double *ip2 = workd+iptr(2)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + else + { + if (ido == 2) + { + for (octave_idx_type i = 0; i < n; i++) + workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; + } + else + { + double *ip2 = workd+iptr(0)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + Matrix eig_vec2 (n, k + 1); + double *z = eig_vec2.fortran_vec (); + + OCTAVE_LOCAL_BUFFER (double, dr, k + 1); + OCTAVE_LOCAL_BUFFER (double, di, k + 1); + OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); + for (octave_idx_type i = 0; i < k+1; i++) + dr[i] = di[i] = 0.; + + F77_FUNC (dneupd, DNEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, + sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) + F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dneupd"); + return -1; + } + else + { + eig_val.resize (k+1); + Complex *d = eig_val.fortran_vec (); + + if (info2 == 0) + { + octave_idx_type jj = 0; + for (octave_idx_type i = 0; i < k+1; i++) + { + if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) + jj++; + else + d [i-jj] = Complex (dr[i], di[i]); + } + if (jj == 0 && !rvec) + for (octave_idx_type i = 0; i < k; i++) + d[i] = d[i+1]; + + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + + eig_vec.resize (n, k); + octave_idx_type i = 0; + while (i < k) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (i+1) * n; + if (std::imag(eig_val(i)) == 0) + { + for (octave_idx_type j = 0; j < n; j++) + eig_vec(j,i) = + Complex(z[j+off1],0.); + i++; + } + else + { + for (octave_idx_type j = 0; j < n; j++) + { + eig_vec(j,i) = + Complex(z[j+off1],z[j+off2]); + if (i < k - 1) + eig_vec(j,i+1) = + Complex(z[j+off1],-z[j+off2]); + } + i+=2; + } + } + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in dneupd", info2); + return -1; + } + } + + return ip(4); +} + +octave_idx_type +EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, + const std::string &_typ, double sigmar, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, ColumnVector &resid, + std::ostream& os, double tol, int rvec, bool cholB, + int disp, int maxit) +{ + std::string typ (_typ); + bool have_sigma = (sigmar ? true : false); + char bmat = 'I'; + double sigmai = 0.; + octave_idx_type mode = 1; + int err = 0; + + if (resid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + resid = ColumnVector (octave_rand::vector(n)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2 + 1; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p < k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k+1 and n"); + return -1; + } + + if (k > n - 1) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n-1).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + + if (! have_sigma) + { + if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); + + if (typ == "LA" || typ == "SA" || typ == "BE") + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for unsymmetric problem"); + return -1; + } + + if (typ == "SM") + { + typ = "LM"; + sigmar = 0.; + mode = 3; + } + } + else if (! std::abs (sigmar)) + typ = "SM"; + else + { + typ = "LM"; + mode = 3; + } + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + octave_idx_type lwork = 3 * p * (p + 2); + + OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); + OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); + OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); + double *presid = resid.fortran_vec (); + + do + { + F77_FUNC (dnaupd, DNAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dnaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + double *ip2 = workd + iptr(0) - 1; + ColumnVector x(n); + + for (octave_idx_type i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (octave_idx_type i = 0; i < n; i++) + *ip2++ = y(i); + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + Matrix eig_vec2 (n, k + 1); + double *z = eig_vec2.fortran_vec (); + + OCTAVE_LOCAL_BUFFER (double, dr, k + 1); + OCTAVE_LOCAL_BUFFER (double, di, k + 1); + OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); + for (octave_idx_type i = 0; i < k+1; i++) + dr[i] = di[i] = 0.; + + F77_FUNC (dneupd, DNEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, + sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) + F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dneupd"); + return -1; + } + else + { + eig_val.resize (k+1); + Complex *d = eig_val.fortran_vec (); + + if (info2 == 0) + { + octave_idx_type jj = 0; + for (octave_idx_type i = 0; i < k+1; i++) + { + if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) + jj++; + else + d [i-jj] = Complex (dr[i], di[i]); + } + if (jj == 0 && !rvec) + for (octave_idx_type i = 0; i < k; i++) + d[i] = d[i+1]; + + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + + eig_vec.resize (n, k); + octave_idx_type i = 0; + while (i < k) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (i+1) * n; + if (std::imag(eig_val(i)) == 0) + { + for (octave_idx_type j = 0; j < n; j++) + eig_vec(j,i) = + Complex(z[j+off1],0.); + i++; + } + else + { + for (octave_idx_type j = 0; j < n; j++) + { + eig_vec(j,i) = + Complex(z[j+off1],z[j+off2]); + if (i < k - 1) + eig_vec(j,i+1) = + Complex(z[j+off1],-z[j+off2]); + } + i+=2; + } + } + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in dneupd", info2); + return -1; + } + } + + return ip(4); +} + +template <class M> +octave_idx_type +EigsComplexNonSymmetricMatrix (const M& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, + ComplexColumnVector &cresid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) +{ + M b(_b); + octave_idx_type n = m.cols (); + octave_idx_type mode = 1; + bool have_b = ! b.is_empty(); + bool note3 = false; + char bmat = 'I'; + Complex sigma = 0.; + M bt; + + if (m.rows() != m.cols()) + { + (*current_liboctave_error_handler) ("eigs: A must be square"); + return -1; + } + if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) + { + (*current_liboctave_error_handler) + ("eigs: B must be square and the same size as A"); + return -1; + } + + if (cresid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + Array<double> rr (octave_rand::vector(n)); + Array<double> ri (octave_rand::vector(n)); + cresid = ComplexColumnVector (n); + for (octave_idx_type i = 0; i < n; i++) + cresid(i) = Complex(rr(i),ri(i)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2 + 1; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p < k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k+1 and n"); + return -1; + } + + if (k > n - 1) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n-1).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (have_b && cholB && permB.length() != 0) + { + // Check the we really have a permutation vector + if (permB.length() != n) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + else + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } + } + + if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + { + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); + return -1; + } + + if (typ == "LA" || typ == "SA" || typ == "BE") + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for complex problem"); + return -1; + } + + if (have_b) + { + // See Note 3 dsaupd + note3 = true; + if (cholB) + { + bt = b; + b = b.hermitian(); + if (permB.length() == 0) + { + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + } + } + else + { + if (! make_cholb(b, bt, permB)) + { + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + return -1; + } + } + } + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + octave_idx_type lwork = p * (3 * p + 5); + + OCTAVE_LOCAL_BUFFER (Complex, v, n * p); + OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); + OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); + OCTAVE_LOCAL_BUFFER (double, rwork, p); + Complex *presid = cresid.fortran_vec (); + + do + { + F77_FUNC (znaupd, ZNAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, rwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in znaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + if (have_b) + { + ComplexMatrix mtmp (n,1); + for (octave_idx_type i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); + for (octave_idx_type i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + + } + else if (!vector_product (m, workd + iptr(0) - 1, + workd + iptr(1) - 1)) + break; + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in znaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + eig_vec.resize (n, k); + Complex *z = eig_vec.fortran_vec (); + + eig_val.resize (k+1); + Complex *d = eig_val.fortran_vec (); + + OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); + + F77_FUNC (zneupd, ZNEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, + F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), + k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in zneupd"); + return -1; + } + + if (info2 == 0) + { + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex ctmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = ctmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + ctmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = ctmp[j]; + } + + if (note3) + eig_vec = ltsolve(b, permB, eig_vec); + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in zneupd", info2); + return -1; + } + + return ip(4); +} + +template <class M> +octave_idx_type +EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, + ComplexColumnVector &cresid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) +{ + M b(_b); + octave_idx_type n = m.cols (); + octave_idx_type mode = 3; + bool have_b = ! b.is_empty(); + std::string typ = "LM"; + + if (m.rows() != m.cols()) + { + (*current_liboctave_error_handler) ("eigs: A must be square"); + return -1; + } + if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) + { + (*current_liboctave_error_handler) + ("eigs: B must be square and the same size as A"); + return -1; + } + + // FIXME: The "SM" type for mode 1 seems unstable though faster!! + //if (! std::abs (sigma)) + // return EigsComplexNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, + // eig_val, _b, permB, cresid, os, tol, + // rvec, cholB, disp, maxit); + + if (cresid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + Array<double> rr (octave_rand::vector(n)); + Array<double> ri (octave_rand::vector(n)); + cresid = ComplexColumnVector (n); + for (octave_idx_type i = 0; i < n; i++) + cresid(i) = Complex(rr(i),ri(i)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2 + 1; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p < k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k+1 and n"); + return -1; + } + + if (k > n - 1) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n-1).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (have_b && cholB && permB.length() != 0) + { + // Check that we really have a permutation vector + if (permB.length() != n) + { + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + return -1; + } + else + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } + } + + char bmat = 'I'; + if (have_b) + bmat = 'G'; + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + M L, U; + + OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows())); + OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols())); + + if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q)) + return -1; + + octave_idx_type lwork = p * (3 * p + 5); + + OCTAVE_LOCAL_BUFFER (Complex, v, n * p); + OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); + OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); + OCTAVE_LOCAL_BUFFER (double, rwork, p); + Complex *presid = cresid.fortran_vec (); + + do + { + F77_FUNC (znaupd, ZNAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, rwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in znaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + if (have_b) + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + vector_product (m, workd+iptr(0)-1, ctmp); + + ComplexMatrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ctmp[P[i]]; + + lusolve (L, U, tmp); + + Complex *ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + else if (ido == 2) + vector_product (b, workd + iptr(0) - 1, workd + iptr(1) - 1); + else + { + Complex *ip2 = workd+iptr(2)-1; + ComplexMatrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + else + { + if (ido == 2) + { + for (octave_idx_type i = 0; i < n; i++) + workd[iptr(0) + i - 1] = + workd[iptr(1) + i - 1]; + } + else + { + Complex *ip2 = workd+iptr(0)-1; + ComplexMatrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + eig_vec.resize (n, k); + Complex *z = eig_vec.fortran_vec (); + + eig_val.resize (k+1); + Complex *d = eig_val.fortran_vec (); + + OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); + + F77_FUNC (zneupd, ZNEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, + F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), + k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in zneupd"); + return -1; + } + + if (info2 == 0) + { + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex ctmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = ctmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + ctmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = ctmp[j]; + } + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in zneupd", info2); + return -1; + } + + return ip(4); +} + +octave_idx_type +EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, + const std::string &_typ, Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + ComplexColumnVector &cresid, std::ostream& os, + double tol, int rvec, bool cholB, int disp, + int maxit) +{ + std::string typ (_typ); + bool have_sigma = (std::abs(sigma) ? true : false); + char bmat = 'I'; + octave_idx_type mode = 1; + int err = 0; + + if (cresid.is_empty()) + { + std::string rand_dist = octave_rand::distribution(); + octave_rand::distribution("uniform"); + Array<double> rr (octave_rand::vector(n)); + Array<double> ri (octave_rand::vector(n)); + cresid = ComplexColumnVector (n); + for (octave_idx_type i = 0; i < n; i++) + cresid(i) = Complex(rr(i),ri(i)); + octave_rand::distribution(rand_dist); + } + + if (p < 0) + { + p = k * 2 + 1; + + if (p < 20) + p = 20; + + if (p > n - 1) + p = n - 1 ; + } + else if (p < k || p > n) + { + (*current_liboctave_error_handler) + ("eigs: opts.p must be between k+1 and n"); + return -1; + } + + if (k > n - 1) + { + (*current_liboctave_error_handler) + ("eigs: Too many eigenvalues to extract (k >= n-1).\n" + " Use 'eig(full(A))' instead"); + return -1; + } + + if (! have_sigma) + { + if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); + + if (typ == "LA" || typ == "SA" || typ == "BE") + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for complex problem"); + return -1; + } + + if (typ == "SM") + { + typ = "LM"; + sigma = 0.; + mode = 3; + } + } + else if (! std::abs (sigma)) + typ = "SM"; + else + { + typ = "LM"; + mode = 3; + } + + Array<octave_idx_type> ip (11); + octave_idx_type *iparam = ip.fortran_vec (); + + ip(0) = 1; //ishift + ip(1) = 0; // ip(1) not referenced + ip(2) = maxit; // mxiter, maximum number of iterations + ip(3) = 1; // NB blocksize in recurrence + ip(4) = 0; // nconv, number of Ritz values that satisfy convergence + ip(5) = 0; //ip(5) not referenced + ip(6) = mode; // mode + ip(7) = 0; + ip(8) = 0; + ip(9) = 0; + ip(10) = 0; + // ip(7) to ip(10) return values + + Array<octave_idx_type> iptr (14); + octave_idx_type *ipntr = iptr.fortran_vec (); + + octave_idx_type ido = 0; + int iter = 0; + octave_idx_type lwork = p * (3 * p + 5); + + OCTAVE_LOCAL_BUFFER (Complex, v, n * p); + OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); + OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); + OCTAVE_LOCAL_BUFFER (double, rwork, p); + Complex *presid = cresid.fortran_vec (); + + do + { + F77_FUNC (znaupd, ZNAUPD) + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, rwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in znaupd"); + return -1; + } + + if (disp > 0 && !xisnan(workl[iptr(5)-1])) + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } + + if (ido == -1 || ido == 1 || ido == 2) + { + Complex *ip2 = workd + iptr(0) - 1; + ComplexColumnVector x(n); + + for (octave_idx_type i = 0; i < n; i++) + x(i) = *ip2++; + + ComplexColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (octave_idx_type i = 0; i < n; i++) + *ip2++ = y(i); + } + else + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } + } + while (1); + + octave_idx_type info2; + + // We have a problem in that the size of the C++ bool + // type relative to the fortran logical type. It appears + // that fortran uses 4-bytes per logical and C++ 1-byte + // per bool, though this might be system dependent. As + // long as the HOWMNY arg is not "S", the logical array + // is just workspace for ARPACK, so use int type to + // avoid problems. + Array<int> s (p); + int *sel = s.fortran_vec (); + + eig_vec.resize (n, k); + Complex *z = eig_vec.fortran_vec (); + + eig_val.resize (k+1); + Complex *d = eig_val.fortran_vec (); + + OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); + + F77_FUNC (zneupd, ZNEUPD) + (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, + F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), + k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in zneupd"); + return -1; + } + + if (info2 == 0) + { + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex ctmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = ctmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + ctmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = ctmp[j]; + } + } + } + else + { + (*current_liboctave_error_handler) + ("eigs: error %d in zneupd", info2); + return -1; + } + + return ip(4); +} + +#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) +extern octave_idx_type +EigsRealSymmetricMatrix (const Matrix& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealSymmetricMatrixShift (const Matrix& m, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, + const std::string &typ, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + Matrix &eig_vec, ColumnVector &eig_val, + ColumnVector &resid, std::ostream &os, + double tol = DBL_EPSILON, int rvec = 0, + bool cholB = 0, int disp = 0, int maxit = 300); + +extern octave_idx_type +EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, + const std::string &_typ, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + ColumnVector &resid, std::ostream& os, + double tol = DBL_EPSILON, int rvec = 0, + bool cholB = 0, int disp = 0, int maxit = 300); + +extern octave_idx_type +EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const ComplexMatrix& b, ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m, + const std::string typ, octave_idx_type k, + octave_idx_type p, octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseComplexMatrix& b, + ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); + +extern octave_idx_type +EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const ComplexMatrix& b, + ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, + int disp = 0, int maxit = 300); + +extern octave_idx_type +EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m, + Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseComplexMatrix& b, + ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, + int disp = 0, int maxit = 300); + +extern octave_idx_type +EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, + const std::string &_typ, Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + ComplexColumnVector &resid, std::ostream& os, + double tol = DBL_EPSILON, int rvec = 0, + bool cholB = 0, int disp = 0, int maxit = 300); +#endif + +#ifndef _MSC_VER +template static octave_idx_type +lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&); + +template static octave_idx_type +lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, + ComplexMatrix&); + +template static octave_idx_type +lusolve (const Matrix&, const Matrix&, Matrix&); + +template static octave_idx_type +lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&); + +template static ComplexMatrix +ltsolve (const SparseComplexMatrix&, const ColumnVector&, + const ComplexMatrix&); + +template static Matrix +ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); + +template static ComplexMatrix +ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); + +template static Matrix +ltsolve (const Matrix&, const ColumnVector&, const Matrix&); + +template static ComplexMatrix +utsolve (const SparseComplexMatrix&, const ColumnVector&, + const ComplexMatrix&); + +template static Matrix +utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); + +template static ComplexMatrix +utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); + +template static Matrix +utsolve (const Matrix&, const ColumnVector&, const Matrix&); +#endif + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2008-12-23 David Bateman <dbateman@free.fr> + + * sparse/svds.m: New function. + * sparse/Makefile.in (SOURCES): Add it here. + 2008-11-21 Radek Salac <salac.r@gmail.com> * sparse/bicgstab.m: New function.
--- a/scripts/sparse/Makefile.in +++ b/scripts/sparse/Makefile.in @@ -35,7 +35,7 @@ SOURCES = bicgstab.m cgs.m colperm.m etreeplot.m gplot.m nonzeros.m normest.m \ pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \ spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \ - spvcat.m spy.m treelayout.m treeplot.m + spvcat.m spy.m svds.m treelayout.m treeplot.m DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
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} {@var{s} =} svds (@var{a}) +## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}) +## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma}) +## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma}, @var{opts}) +## @deftypefnx {Function File} {[@var{u}, @var{s}, @var{v}, @var{flag}] =} svds (@dots{}) +## +## Find a few singular values of the matrix @var{a}. The singular values +## are calculated using +## +## @example +## @group +## [@var{m}, @var{n}] = size(@var{a}) +## @var{s} = eigs([sparse(@var{m}, @var{m}), @var{a}; @dots{} +## @var{a}', sparse(@var{n}, @var{n})]) +## @end group +## @end example +## +## The eigenvalues returned by @code{eigs} correspond to the singular +## values of @var{a}. The number of singular values to calculate is given +## by @var{k}, whose default value is 6. +## +## The argument @var{sigma} can be used to specify which singular values +## to find. @var{sigma} can be either the string 'L', the default, in +## which case the largest singular values of @var{a} are found. Otherwise +## @var{sigma} should be a real scalar, in which case the singular values +## closest to @var{sigma} are found. Note that for relatively small values +## of @var{sigma}, there is the chance that the requested number of singular +## values are not returned. In that case @var{sigma} should be increased. +## +## If @var{opts} is given, then it is a structure that defines options +## that @code{svds} will pass to @var{eigs}. The possible fields of this +## structure are therefore determined by @code{eigs}. By default three +## fields of this structure are set by @code{svds}. +## +## @table @code +## @item tol +## The required convergence tolerance for the singular values. @code{eigs} +## is passed @var{tol} divided by @code{sqrt(2)}. The default value is +## 1e-10. +## +## @item maxit +## The maximum number of iterations. The defaut is 300. +## +## @item disp +## The level of diagnostic printout. If @code{disp} is 0 then there is no +## printout. The default value is 0. +## @end table +## +## If more than one output argument is given, then @code{svds} also +## calculates the left and right singular vectors of @var{a}. @var{flag} +## is used to signal the convergence of @code{svds}. If @code{svds} +## converges to the desired tolerance, then @var{flag} given by +## +## @example +## @group +## norm (@var{a} * @var{v} - @var{u} * @var{s}, 1) <= @dots{} +## @var{tol} * norm (@var{a}, 1) +## @end group +## @end example +## +## will be zero. +## @end deftypefn +## @seealso{eigs} + +function [u, s, v, flag] = svds (a, k, sigma, opts) + + if (nargin < 1 || nargin > 4) + error ("Incorrect number of arguments"); + endif + + if (nargin < 4) + opts.tol = 1e-10 / sqrt(2); + opts.disp = 0; + opts.maxit = 300; + else + if (!isstruct(opts)) + error("opts must be a structure"); + endif + if (!isfield(opts,"tol")) + opts.tol = 1e-10 / sqrt(2); + endif + endif + + if (nargin < 3 || strcmp(sigma,"L")) + if (isreal(a)) + sigma = "LA"; + else + sigma = "LR"; + endif + elseif (isscalar(sigma) && isreal(sigma)) + if ((sigma < 0)) + error ("sigma must be a positive real value"); + endif + else + error ("sigma must be a positive real value or the string 'L'"); + endif + + maxA = max(max(abs(a))); + if (maxA == 0) + u = eye(m, k); + s = zeros(k, k); + v = eye(n, k); + else + [m, n] = size(a); + if (nargin < 2) + k = min([6, m, n]); + else + k = min([k, m, n]); + endif + + ## Scale everything by the 1-norm to make things more stable. + B = a / maxA; + Bopts = opts; + Bopts.tol = opts.tol / maxA; + Bsigma = sigma; + if (!ischar(Bsigma)) + Bsigma = Bsigma / maxA; + endif + + if (!ischar(Bsigma) && Bsigma == 0) + ## The eigenvalues returns by eigs are symmetric about 0. As we + ## are only interested in the positive eigenvalues, we have to + ## double k. If sigma is smaller than the smallest singular value + ## this can also be an issue. However, we'd like to avoid double + ## k for all scalar value of sigma... + [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)], + 2 * k, Bsigma, Bopts); + else + [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)], + k, Bsigma, Bopts); + endif + s = diag(s); + + if (ischar(sigma)) + norma = max(s); + else + norma = normest(a); + endif + V = sqrt(2) * V; + u = V(1:m,:); + v = V(m+1:end,:); + + ## We wish to exclude all eigenvalues that are less than zero as these + ## are artifacts of the way the matrix passed to eigs is formed. There + ## is also the possibility that the value of sigma chosen is exactly + ## a singular value, and in that case we're dead!! So have to rely on + ## the warning from eigs. We exclude the singular values which are + ## less than or equal to zero to within some tolerance scaled by the + ## norm since if we don't we might end up with too many singular + ## values. What is appropriate for the tolerance? + tol = norma * opts.tol; + ind = find(s > tol); + if (length(ind) < k) + ## Find the zero eigenvalues of B, Ignore the eigenvalues that are + ## nominally negative. + zind = find(abs(s) <= tol); + p = min(length(zind), k-length(ind)); + ind = [ind;zind(1:p)]; + elseif (length(ind) > k) + ind = ind(1:k); + endif + u = u(:,ind); + s = s(ind); + v = v(:,ind); + + if (length(s) < k) + warning("returning fewer singular values than requested."); + if (!ischar(sigma)) + warning("try increasing the value of sigma"); + endif + endif + + s = s * maxA; + endif + + if (nargout < 2) + u = s; + else + s = diag(s); + if (nargout > 3) + flag = norm(a*v - u*s, 1) > sqrt(2) * opts.tol * norm(a, 1); + endif + endif +endfunction + +%!shared n, k, a, u, s, v, opts +%! n = 100; +%! k = 7; +%! a = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),0.4*n*ones(1,n),ones(1,n-2)]); +%! %%a = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]); +%! [u,s,v] = svd(full(a)); +%! s = diag(s); +%! [dum, idx] = sort(abs(s)); +%! s = s(idx); +%! u = u(:,idx); +%! v = v(:,idx); +%! randn('state',42) +%!test +%! [u2,s2,v2,flag] = svds(a,k); +%! s2 = diag(s2); +%! assert(flag,!1); +%! assert(s(end:-1:end-k+1), s2, 1e-10); +%!test +%! [u2,s2,v2,flag] = svds(a,k,0); +%! s2 = diag(s2); +%! assert(flag,!1); +%! assert(s(k:-1:1), s2, 1e-10); +%!test +%! idx = floor(n/2); +%! % Don't put sigma right on a singular value or there are convergence +%! sigma = 0.99*s(idx) + 0.01*s(idx+1); +%! [u2,s2,v2,flag] = svds(a,k,sigma); +%! s2 = diag(s2); +%! assert(flag,!1); +%! assert(s((idx+floor(k/2)):-1:(idx-floor(k/2))), s2, 1e-10);
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2008-12-23 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/eigs.cc: New file. + * Makefile.in (DLD_XSRC): Add it here. + 2008-12-22 David Bateman <dbateman@free.fr> * DLD-FUNCTIONS/__voronoi__.cc (F__voronoi__): Resize AtInf array
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\ +@deftypefn {Loadable Function} {@var{d}} = eigs (@var{a})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k}, @var{sigma},@var{opts})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k}, @var{sigma}, @var{opts})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k}, @var{sigma})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts})\n\ +@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k}, @var{sigma}, @var{opts})\n\ +@deftypefnx {Loadable Function} {[@var{v}, @var{d}]} = eigs (@var{a}, @dots{})\n\ +@deftypefnx {Loadable Function} {[@var{v}, @var{d}]} = eigs (@var{af}, @var{n}, @dots{})\n\ +@deftypefnx {Loadable Function} {[@var{v}, @var{d}, @var{flag}]} = eigs (@var{a}, @dots{})\n\ +@deftypefnx {Loadable Function} {[@var{v}, @var{d}, @var{flag}]} = eigs (@var{af}, @var{n}, @dots{})\n\ +Calculate a limited number of eigenvalues and eigenvectors of @var{a},\n\ +based on a selection criteria. The number eigenvalues and eigenvectors to\n\ +calculate is given by @var{k} whose default value is 6.\n\ +\n\ +By default @code{eigs} solve the equation\n\ +@iftex\n\ +@tex\n\ +$A \\nu = \\lambda \\nu$\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +@code{A * v = lambda * v}\n\ +@end ifinfo\n\ +, where\n\ +@iftex\n\ +@tex\n\ +$\\lambda$ is a scalar representing one of the eigenvalues, and $\\nu$\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +@code{lambda} is a scalar representing one of the eigenvalues, and @code{v}\n\ +@end ifinfo\n\ +is the corresponding eigenvector. If given the positive definite matrix\n\ +@var{B} then @code{eigs} solves the general eigenvalue equation\n\ +@iftex\n\ +@tex\n\ +$A \\nu = \\lambda B \\nu$\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +@code{A * v = lambda * B * v}\n\ +@end ifinfo\n\ +.\n\ +\n\ +The argument @var{sigma} determines which eigenvalues are returned.\n\ +@var{sigma} can be either a scalar or a string. When @var{sigma} is a scalar,\n\ +the @var{k} eigenvalues closest to @var{sigma} are returned. If @var{sigma}\n\ +is a string, it must have one of the values\n\ +\n\ +@table @asis\n\ +@item 'lm'\n\ +Largest magnitude (default).\n\ +\n\ +@item 'sm'\n\ +Smallest magnitude.\n\ +\n\ +@item 'la'\n\ +Largest Algebraic (valid only for real symmetric problems).\n\ +\n\ +@item 'sa'\n\ +Smallest Algebraic (valid only for real symmetric problems).\n\ +\n\ +@item 'be'\n\ +Both ends, with one more from the high-end if @var{k} is odd (valid only for\n\ +real symmetric problems).\n\ +\n\ +@item 'lr'\n\ +Largest real part (valid only for complex or unsymmetric problems).\n\ +\n\ +@item 'sr'\n\ +Smallest real part (valid only for complex or unsymmetric problems).\n\ +\n\ +@item 'li'\n\ +Largest imaginary part (valid only for complex or unsymmetric problems).\n\ +\n\ +@item 'si'\n\ +Smallest imaginary part (valid only for complex or unsymmetric problems).\n\ +@end table\n\ +\n\ +If @var{opts} is given, it is a structure defining some of the options that\n\ +@code{eigs} should use. The fields of the structure @var{opts} are\n\ +\n\ +@table @code\n\ +@item issym\n\ +If @var{af} is given, then flags whether the function @var{af} defines a\n\ +symmetric problem. It is ignored if @var{a} is given. The default is false.\n\ +\n\ +@item isreal\n\ +If @var{af} is given, then flags whether the function @var{af} defines a\n\ +real problem. It is ignored if @var{a} is given. The default is true.\n\ +\n\ +@item tol\n\ +Defines the required convergence tolerance, given as @code{tol * norm (A)}.\n\ +The default is @code{eps}.\n\ +\n\ +@item maxit\n\ +The maximum number of iterations. The default is 300.\n\ +\n\ +@item p\n\ +The number of Lanzcos basis vectors to use. More vectors will result in\n\ +faster convergence, but a larger amount of memory. The optimal value of 'p'\n\ +is problem dependent and should be in the range @var{k} to @var{n}. The\n\ +default value is @code{2 * @var{k}}.\n\ +\n\ +@item v0\n\ +The starting vector for the computation. The default is to have @sc{Arpack}\n\ +randomly generate a starting vector.\n\ +\n\ +@item disp\n\ +The level of diagnostic printout. If @code{disp} is 0 then there is no\n\ +printout. The default value is 1.\n\ +\n\ +@item cholB\n\ +Flag if @code{chol (@var{b})} is passed rather than @var{b}. The default is\n\ +false.\n\ +\n\ +@item permB\n\ +The permutation vector of the Cholesky factorization of @var{b} if\n\ +@code{cholB} is true. That is @code{chol ( @var{b} (permB, permB))}. The\n\ +default is @code{1:@var{n}}.\n\ +\n\ +@end table\n\ +\n\ +It is also possible to represent @var{a} by a function denoted @var{af}.\n\ +@var{af} must be followed by a scalar argument @var{n} defining the length\n\ +of the vector argument accepted by @var{af}. @var{af} can be passed either\n\ +as an inline function, function handle or as a string. In the case where\n\ +@var{af} is passed as a string, the name of the string defines the function\n\ +to use.\n\ +\n\ +@var{af} is a function of the form @code{function y = af (x), y = @dots{};\n\ +endfunction}, where the required return value of @var{af} is determined by\n\ +the value of @var{sigma}, and are\n\ +\n\ +@table @code\n\ +@item A * x\n\ +If @var{sigma} is not given or is a string other than 'sm'.\n\ +\n\ +@item A \\ x\n\ +If @var{sigma} is 'sm'.\n\ +\n\ +@item (A - sigma * I) \\ x\n\ +for standard eigenvalue problem, where @code{I} is the identity matrix of\n\ +the same size as @code{A}. If @var{sigma} is zero, this reduces the\n\ +@code{A \\ x}.\n\ +\n\ +@item (A - sigma * B) \\ x\n\ +for the general eigenvalue problem.\n\ +@end table\n\ +\n\ +The return arguments of @code{eigs} depends on the number of return\n\ +arguments. With a single return argument, a vector @var{d} of length @var{k}\n\ +is returned, represent the @var{k} eigenvalues that have been found. With two\n\ +return arguments, @var{v} is a @var{n}-by-@var{k} matrix whose columns are\n\ +the @var{k} eigenvectors corresponding to the returned eigenvalues. The\n\ +eigenvalues themselves are then returned in @var{d} in the form of a\n\ +@var{n}-by-@var{k} matrix, where the elements on the diagonal are the\n\ +eigenvalues.\n\ +\n\ +Given a third return argument @var{flag}, @code{eigs} also returns the status\n\ +of the convergence. If @var{flag} is 0, then all eigenvalues have converged,\n\ +otherwise not.\n\ +\n\ +This function is based on the @sc{Arpack} package, written by R Lehoucq,\n\ +K Maschhoff, D Sorensen and C Yang. For more information see\n\ +@url{http://www.caam.rice.edu/software/ARPACK/}.\n\ +\n\ +@end deftypefn\n\ +@seealso{eig, svds}") +{ + octave_value_list retval; +#ifdef HAVE_ARPACK + int nargin = args.length (); + std::string fcn_name; + octave_idx_type n = 0; + octave_idx_type k = 6; + Complex sigma = 0.; + double sigmar, sigmai; + bool have_sigma = false; + std::string typ = "LM"; + Matrix amm, bmm, bmt; + ComplexMatrix acm, bcm, bct; + SparseMatrix asmm, bsmm, bsmt; + SparseComplexMatrix ascm, bscm, bsct; + int b_arg = 0; + bool have_b = false; + bool have_a_fun = false; + bool a_is_complex = false; + bool b_is_complex = false; + bool symmetric = false; + bool cholB = false; + bool a_is_sparse = false; + ColumnVector permB; + int arg_offset = 0; + double tol = DBL_EPSILON; + int maxit = 300; + int disp = 0; + octave_idx_type p = -1; + ColumnVector resid; + ComplexColumnVector cresid; + octave_idx_type info = 1; + char bmat = 'I'; + + warned_imaginary = false; + + unwind_protect::begin_frame ("Feigs"); + + unwind_protect_int (call_depth); + call_depth++; + + if (call_depth > 1) + { + error ("eigs: invalid recursive call"); + if (fcn_name.length()) + clear_function (fcn_name); + unwind_protect::run_frame ("Feigs"); + return retval; + } + + if (nargin == 0) + print_usage (); + else if (args(0).is_function_handle () || args(0).is_inline_function () || + args(0).is_string ()) + { + if (args(0).is_string ()) + { + std::string name = args(0).string_value (); + std::string fname = "function y = "; + fcn_name = unique_symbol_name ("__eigs_fcn_"); + fname.append (fcn_name); + fname.append ("(x) y = "); + eigs_fcn = extract_function (args(0), "eigs", fcn_name, fname, + "; endfunction"); + } + else + eigs_fcn = args(0).function_value (); + + if (!eigs_fcn) + error ("eigs: unknown function"); + + if (nargin < 2) + error ("eigs: incorrect number of arguments"); + else + { + n = args(1).nint_value (); + arg_offset = 1; + have_a_fun = true; + } + } + else + { + if (args(0).is_complex_type ()) + { + if (args(0).is_sparse_type ()) + { + ascm = (args(0).sparse_complex_matrix_value()); + a_is_sparse = true; + } + else + acm = (args(0).complex_matrix_value()); + a_is_complex = true; + symmetric = false; // ARAPACK doesn't special case complex symmetric + } + else + { + if (args(0).is_sparse_type ()) + { + asmm = (args(0).sparse_matrix_value()); + a_is_sparse = true; + symmetric = asmm.is_symmetric(); + } + else + { + amm = (args(0).matrix_value()); + symmetric = amm.is_symmetric(); + } + } + + // Note hold off reading B till later to avoid issues of double + // copies of the matrix if B is full/real while A is complex.. + if (!error_state && nargin > 1 && + !(args(1).is_real_scalar ())) + if (args(1).is_complex_type ()) + { + b_arg = 1+arg_offset; + have_b = true; + bmat = 'G'; + b_is_complex = true; + arg_offset++; + } + else + { + b_arg = 1+arg_offset; + have_b = true; + bmat = 'G'; + arg_offset++; + } + } + + if (!error_state && nargin > (1+arg_offset)) + k = args(1+arg_offset).nint_value (); + + if (!error_state && nargin > (2+arg_offset)) + if (args(2+arg_offset).is_real_scalar () || + args(2+arg_offset).is_complex_scalar ()) + { + sigma = args(2+arg_offset).complex_value (); + have_sigma = true; + } + else if (args(2+arg_offset).is_string ()) + { + typ = args(2+arg_offset).string_value (); + + // Use STL function to convert to upper case + transform (typ.begin (), typ.end (), typ.begin (), toupper); + + sigma = 0.; + } + else + error ("eigs: sigma must be a scalar or a string"); + + sigmar = std::real (sigma); + sigmai = std::imag (sigma); + + if (!error_state && nargin > (3+arg_offset)) + if (args(3+arg_offset).is_map ()) + { + Octave_map map = args(3+arg_offset).map_value (); + + // issym is ignored if A is not a function + if (map.contains ("issym")) + if (have_a_fun) + symmetric = map.contents ("issym")(0).double_value () != 0.; + + // isreal is ignored if A is not a function + if (map.contains ("isreal")) + if (have_a_fun) + a_is_complex = ! (map.contents ("isreal")(0).double_value () != 0.); + + if (map.contains ("tol")) + tol = map.contents ("tol")(0).double_value (); + + if (map.contains ("maxit")) + maxit = map.contents ("maxit")(0).nint_value (); + + if (map.contains ("p")) + p = map.contents ("p")(0).nint_value (); + + if (map.contains ("v0")) + { + if (a_is_complex || b_is_complex) + cresid = ComplexColumnVector + (map.contents ("v0")(0).complex_vector_value()); + else + resid = ColumnVector (map.contents ("v0")(0).vector_value()); + } + + if (map.contains ("disp")) + disp = map.contents ("disp")(0).nint_value (); + + if (map.contains ("cholB")) + cholB = map.contents ("cholB")(0).double_value () != 0.; + + if (map.contains ("permB")) + permB = ColumnVector (map.contents ("permB")(0).vector_value ()) + - 1.0; + } + else + error ("eigs: options argument must be a structure"); + + if (nargin > (4+arg_offset)) + error ("eigs: incorrect number of arguments"); + + if (have_b) + { + if (a_is_complex || b_is_complex) + { + if (a_is_sparse) + bscm = args(b_arg).sparse_complex_matrix_value (); + else + bcm = args(b_arg).complex_matrix_value (); + } + else + { + if (a_is_sparse) + bsmm = args(b_arg).sparse_matrix_value (); + else + bmm = args(b_arg).matrix_value (); + } + } + + // Mode 1 for SM mode seems unstable for some reason. + // Use Mode 3 instead, with sigma = 0. + if (!error_state && !have_sigma && typ == "SM") + have_sigma = true; + + if (!error_state) + { + octave_idx_type nconv; + if (a_is_complex || b_is_complex) + { + ComplexMatrix eig_vec; + ComplexColumnVector eig_val; + + + if (have_a_fun) + nconv = EigsComplexNonSymmetricFunc + (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, eig_val, + cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit); + else if (have_sigma) + if (a_is_sparse) + nconv = EigsComplexNonSymmetricMatrixShift + (ascm, sigma, k, p, info, eig_vec, eig_val, bscm, permB, + cresid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else + nconv = EigsComplexNonSymmetricMatrixShift + (acm, sigma, k, p, info, eig_vec, eig_val, bcm, permB, cresid, + octave_stdout, tol, (nargout > 1), cholB, disp, maxit); + else + if (a_is_sparse) + nconv = EigsComplexNonSymmetricMatrix + (ascm, typ, k, p, info, eig_vec, eig_val, bscm, permB, cresid, + octave_stdout, tol, (nargout > 1), cholB, disp, maxit); + else + nconv = EigsComplexNonSymmetricMatrix + (acm, typ, k, p, info, eig_vec, eig_val, bcm, permB, cresid, + octave_stdout, tol, (nargout > 1), cholB, disp, maxit); + + if (nargout < 2) + retval (0) = eig_val; + else + { + retval(2) = double (info); + retval(1) = ComplexDiagMatrix (eig_val); + retval(0) = eig_vec; + } + } + else if (sigmai != 0.) + { + // Promote real problem to a complex one. + ComplexMatrix eig_vec; + ComplexColumnVector eig_val; + + if (have_a_fun) + nconv = EigsComplexNonSymmetricFunc + (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, eig_val, + cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit); + else + if (a_is_sparse) + nconv = EigsComplexNonSymmetricMatrixShift + (SparseComplexMatrix (asmm), sigma, k, p, info, eig_vec, + eig_val, SparseComplexMatrix (bsmm), permB, cresid, + octave_stdout, tol, (nargout > 1), cholB, disp, maxit); + else + nconv = EigsComplexNonSymmetricMatrixShift + (ComplexMatrix (amm), sigma, k, p, info, eig_vec, + eig_val, ComplexMatrix (bmm), permB, cresid, + octave_stdout, tol, (nargout > 1), cholB, disp, maxit); + + if (nargout < 2) + retval (0) = eig_val; + else + { + retval(2) = double (info); + retval(1) = ComplexDiagMatrix (eig_val); + retval(0) = eig_vec; + } + } + else + { + if (symmetric) + { + Matrix eig_vec; + ColumnVector eig_val; + + if (have_a_fun) + nconv = EigsRealSymmetricFunc + (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else if (have_sigma) + if (a_is_sparse) + nconv = EigsRealSymmetricMatrixShift + (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else + nconv = EigsRealSymmetricMatrixShift + (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else + if (a_is_sparse) + nconv = EigsRealSymmetricMatrix + (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else + nconv = EigsRealSymmetricMatrix + (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + + if (nargout < 2) + retval (0) = eig_val; + else + { + retval(2) = double (info); + retval(1) = DiagMatrix (eig_val); + retval(0) = eig_vec; + } + } + else + { + ComplexMatrix eig_vec; + ComplexColumnVector eig_val; + + if (have_a_fun) + nconv = EigsRealNonSymmetricFunc + (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else if (have_sigma) + if (a_is_sparse) + nconv = EigsRealNonSymmetricMatrixShift + (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else + nconv = EigsRealNonSymmetricMatrixShift + (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else + if (a_is_sparse) + nconv = EigsRealNonSymmetricMatrix + (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + else + nconv = EigsRealNonSymmetricMatrix + (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB, + resid, octave_stdout, tol, (nargout > 1), cholB, disp, + maxit); + + if (nargout < 2) + retval (0) = eig_val; + else + { + retval(2) = double (info); + retval(1) = ComplexDiagMatrix (eig_val); + retval(0) = eig_vec; + } + } + } + + if (nconv <= 0) + warning ("eigs: None of the %d requested eigenvalues converged", k); + else if (nconv < k) + warning ("eigs: Only %d of the %d requested eigenvalues converged", + nconv, k); + } + + if (! fcn_name.empty ()) + clear_function (fcn_name); + + unwind_protect::run_frame ("Feigs"); +#else + error ("eigs: not available in this version of Octave"); +#endif + + return retval; +} + +/* #### SPARSE MATRIX VERSIONS #### */ + +/* + +%% Real positive definite tests, n must be even +%!shared n, k, A, d0, d2 +%! n = 20; +%! k = 4; +%! A = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)]); +%! d0 = eig (A); +%! d2 = sort(d0); +%! [dum, idx] = sort( abs(d0)); +%! d0 = d0(idx); +%!test +%! d1 = eigs (A, k); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); +%!test +%! d1 = eigs (A,k+1); +%! assert (d1, d0(end:-1:(end-k)),1e-12); +%!test +%! d1 = eigs (A, k, 'lm'); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); +%!test +%! d1 = eigs (A, k, 'sm'); +%! assert (d1, d0(k:-1:1), 1e-12); +%!test +%! d1 = eigs (A, k, 'la'); +%! assert (d1, d2(end:-1:(end-k+1)), 1e-12); +%!test +%! d1 = eigs (A, k, 'sa'); +%! assert (d1, d2(1:k), 1e-12); +%!test +%! d1 = eigs (A, k, 'be'); +%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-12); +%!test +%! d1 = eigs (A, k+1, 'be'); +%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-12); +%!test +%! d1 = eigs (A, k, 4.1); +%! [dum,idx0] = sort (abs(d0 - 4.1)); +%! [dum,idx1] = sort (abs(d1 - 4.1)); +%! assert (d1(idx1), d0(idx0(1:k)), 1e-12); +%!test +%! d1 = eigs(A, speye(n), k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, speye(n), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, speye(n), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12); +%!test +%! fn = @(x) A * x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'lm', opts); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); +%!test +%! fn = @(x) A \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'sm', opts); +%! assert (d1, d0(k:-1:1), 1e-12); +%!test +%! fn = @(x) (A - 4.1 * eye(n)) \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (d1, eigs(A,k,4.1), 1e-12); +%!test +%! [v1,d1] = eigs(A, k, 'lm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'la'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sa'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'be'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor + +*/ + +/* + +%% Real unsymmetric tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]); +%! d0 = eig (A); +%! [dum, idx] = sort(abs(d0)); +%! d0 = d0(idx); +%!test +%! d1 = eigs (A, k); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A,k+1); +%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); +%!test +%! d1 = eigs (A, k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sm'); +%! assert (abs(d1), abs(d0(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'lr'); +%! [dum, idx] = sort (real(d0)); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sr'); +%! [dum, idx] = sort (real(abs(d0))); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'li'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); +%!test +%! d1 = eigs (A, k, 'si'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); +%!test +%! d1 = eigs (A, k, 4.1); +%! [dum,idx0] = sort (abs(d0 - 4.1)); +%! [dum,idx1] = sort (abs(d1 - 4.1)); +%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); +%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); +%!test +%! d1 = eigs(A, speye(n), k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, speye(n), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, speye(n), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12); +%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))), 1e-12); +%!test +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'sm', opts); +%! assert (abs(d1), d0(1:k), 1e-12); +%!test +%! fn = @(x) (A - 4.1 * eye(n)) \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! [v1,d1] = eigs(A, k, 'lm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'lr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'li'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'si'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor + +*/ + +/* + +%% Complex hermitian tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)]); +%! d0 = eig (A); +%! [dum, idx] = sort(abs(d0)); +%! d0 = d0(idx); +%!test +%! d1 = eigs (A, k); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A,k+1); +%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); +%!test +%! d1 = eigs (A, k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sm'); +%! assert (abs(d1), abs(d0(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'lr'); +%! [dum, idx] = sort (real(abs(d0))); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sr'); +%! [dum, idx] = sort (real(abs(d0))); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'li'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); +%!test +%! d1 = eigs (A, k, 'si'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); +%!test +%! d1 = eigs (A, k, 4.1); +%! [dum,idx0] = sort (abs(d0 - 4.1)); +%! [dum,idx1] = sort (abs(d1 - 4.1)); +%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); +%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); +%!test +%! d1 = eigs(A, speye(n), k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, speye(n), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, speye(n), k, 4.1, opts); +%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); +%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts); +%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); +%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); +%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12); +%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))), 1e-12); +%!test +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 'sm', opts); +%! assert (abs(d1), d0(1:k), 1e-12); +%!test +%! fn = @(x) (A - 4.1 * eye(n)) \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! [v1,d1] = eigs(A, k, 'lm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'lr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'li'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'si'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) +%! endfor + +*/ + +/* #### FULL MATRIX VERSIONS #### */ + +/* + +%% Real positive definite tests, n must be even +%!shared n, k, A, d0, d2 +%! n = 20; +%! k = 4; +%! A = full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)])); +%! d0 = eig (A); +%! d2 = sort(d0); +%! [dum, idx] = sort( abs(d0)); +%! d0 = d0(idx); +%!test +%! d1 = eigs (A, k); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); +%!test +%! d1 = eigs (A,k+1); +%! assert (d1, d0(end:-1:(end-k)),1e-12); +%!test +%! d1 = eigs (A, k, 'lm'); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); +%!test +%! d1 = eigs (A, k, 'sm'); +%! assert (d1, d0(k:-1:1), 1e-12); +%!test +%! d1 = eigs (A, k, 'la'); +%! assert (d1, d2(end:-1:(end-k+1)), 1e-12); +%!test +%! d1 = eigs (A, k, 'sa'); +%! assert (d1, d2(1:k), 1e-12); +%!test +%! d1 = eigs (A, k, 'be'); +%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-12); +%!test +%! d1 = eigs (A, k+1, 'be'); +%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-12); +%!test +%! d1 = eigs (A, k, 4.1); +%! [dum,idx0] = sort (abs(d0 - 4.1)); +%! [dum,idx1] = sort (abs(d1 - 4.1)); +%! assert (d1(idx1), d0(idx0(1:k)), 1e-12); +%!test +%! d1 = eigs(A, eye(n), k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!assert (eigs(A,k,4.1), eigs(A,eye(n),k,4.1), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, eye(n), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, eye(n), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!assert (eigs(A,k,4.1), eigs(A,eye(n),k,4.1), 1e-12); +%!test +%! fn = @(x) A * x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'lm', opts); +%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); +%!test +%! fn = @(x) A \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'sm', opts); +%! assert (d1, d0(k:-1:1), 1e-12); +%!test +%! fn = @(x) (A - 4.1 * eye(n)) \ x; +%! opts.issym = 1; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (d1, eigs(A,k,4.1), 1e-12); +%!test +%! [v1,d1] = eigs(A, k, 'lm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'la'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sa'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'be'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor + +*/ + +/* + +%% Real unsymmetric tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)])); +%! d0 = eig (A); +%! [dum, idx] = sort(abs(d0)); +%! d0 = d0(idx); +%!test +%! d1 = eigs (A, k); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A,k+1); +%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); +%!test +%! d1 = eigs (A, k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sm'); +%! assert (abs(d1), abs(d0(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'lr'); +%! [dum, idx] = sort (real(d0)); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sr'); +%! [dum, idx] = sort (real(abs(d0))); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'li'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); +%!test +%! d1 = eigs (A, k, 'si'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); +%!test +%! d1 = eigs (A, k, 4.1); +%! [dum,idx0] = sort (abs(d0 - 4.1)); +%! [dum,idx1] = sort (abs(d1 - 4.1)); +%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); +%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); +%!test +%! d1 = eigs(A, eye(n), k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, eye(n), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, eye(n), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12); +%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(n),k,4.1))), 1e-12); +%!test +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 'sm', opts); +%! assert (abs(d1), d0(1:k), 1e-12); +%!test +%! fn = @(x) (A - 4.1 * eye(n)) \ x; +%! opts.issym = 0; opts.isreal = 1; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! [v1,d1] = eigs(A, k, 'lm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'lr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'li'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'si'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor + +*/ + +/* + +%% Complex hermitian tests +%!shared n, k, A, d0 +%! n = 20; +%! k = 4; +%! A = full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)])); +%! d0 = eig (A); +%! [dum, idx] = sort(abs(d0)); +%! d0 = d0(idx); +%!test +%! d1 = eigs (A, k); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A,k+1); +%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); +%!test +%! d1 = eigs (A, k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sm'); +%! assert (abs(d1), abs(d0(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'lr'); +%! [dum, idx] = sort (real(abs(d0))); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); +%!test +%! d1 = eigs (A, k, 'sr'); +%! [dum, idx] = sort (real(abs(d0))); +%! d2 = d0(idx); +%! assert (real(d1), real(d2(1:k)), 1e-12); +%!test +%! d1 = eigs (A, k, 'li'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); +%!test +%! d1 = eigs (A, k, 'si'); +%! [dum, idx] = sort (imag(abs(d0))); +%! d2 = d0(idx); +%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); +%!test +%! d1 = eigs (A, k, 4.1); +%! [dum,idx0] = sort (abs(d0 - 4.1)); +%! [dum,idx1] = sort (abs(d1 - 4.1)); +%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); +%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); +%!test +%! d1 = eigs(A, eye(n), k, 'lm'); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, eye(n), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! opts.cholB=true; +%! d1 = eigs(A, eye(n), k, 4.1, opts); +%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); +%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); +%!test +%! opts.cholB=true; +%! q = [2:n,1]; +%! opts.permB=q; +%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts); +%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); +%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); +%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12); +%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(n),k,4.1))), 1e-12); +%!test +%! fn = @(x) A * x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 'lm', opts); +%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); +%!test +%! fn = @(x) A \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 'sm', opts); +%! assert (abs(d1), d0(1:k), 1e-12); +%!test +%! fn = @(x) (A - 4.1 * eye(n)) \ x; +%! opts.issym = 0; opts.isreal = 0; +%! d1 = eigs (fn, n, k, 4.1, opts); +%! assert (abs(d1), eigs(A,k,4.1), 1e-12); +%!test +%! [v1,d1] = eigs(A, k, 'lm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sm'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'lr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'sr'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'li'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor +%!test +%! [v1,d1] = eigs(A, k, 'si'); +%! d1 = diag(d1); +%! for i=1:k +%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) +%! endfor + +*/ + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -300,7 +300,7 @@ -L../libcruft $(LIBCRUFT) -L../liboctave $(LIBOCTAVE) \ -L. $(LIBOCTINTERP) $(CHOLMOD_LIBS) $(UMFPACK_LIBS) $(AMD_LIBS) \ $(CAMD_LIBS) $(COLAMD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \ - $(FFTW_LIBS) $(LIBS) $(FLIBS) + $(FFTW_LIBS) $(ARPACK_LIBS) $(LIBS) $(FLIBS) BUILT_DISTFILES = DOCSTRINGS oct-gperf.h parse.cc lex.cc y.tab.h \ $(OPT_HANDLERS) $(BUILT_EXTRAS) @@ -371,7 +371,7 @@ $(OCTAVE_LIBS) \ $(LEXLIB) $(UMFPACK_LIBS) $(AMD_LIBS) $(CAMD_LIBS) $(COLAMD_LIBS) \ $(CHOLMOD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \ - $(FFTW_LIBS) $(OPENGL_LIBS) $(LIBS) $(FLIBS) + $(FFTW_LIBS) $(ARPACK_LIBS) $(OPENGL_LIBS) $(LIBS) $(FLIBS) stmp-pic: pic @if [ -f stmp-pic ]; then \ @@ -650,6 +650,7 @@ convhulln.oct: OCT_LINK_DEPS += $(QHULL_LIBS) __delaunayn__.oct: OCT_LINK_DEPS += $(QHULL_LIBS) __voronoi__.oct: OCT_LINK_DEPS += $(QHULL_LIBS) +eigs.oct: OCT_LINK_DEPS += $(ARPACK_LIBS) regexp.oct: OCT_LINK_DEPS += $(REGEX_LIBS) urlwrite.oct: OCT_LINK_DEPS += $(CURL_LIBS) __glpk__.oct: OCT_LINK_DEPS += $(GLPK_LIBS)