[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Eigs again - Success!!
From: |
David Bateman |
Subject: |
Eigs again - Success!! |
Date: |
Sat, 04 Nov 2006 21:27:27 +0100 |
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
Ok, I think I've figured out the seg fqault I was getting. Basically, I
believe my test example turned up a bug in ARPACK itself, and was not an
issue with gfortran/g95 as I'd previously thought..
The problem is that the Ritz vectors are calculated over a N-by-N
sub-matrix that is typically about twice as large as the number of
requested eigenvalues. When this sub-matrix needs to be reordered to
extract the correct eigenvalues and vectors, such as with the "sr" flag,
then the xTRSEN lapack function is used for the non-symmetric and
complex cases, though not for the symmetric real case (which is why the
problem doesn't occur there). xTRSEN is passed the variable NCONV for
its output variable M, where NCONV on input to xTRSEN is the number of
converged eigenvalues (upto the number of eigenvalues requested), and in
fact effectively is the same thing on output. Except for the case where
there are more converged eigenvalues than wanted. In that case NCONV is
increased and the following code that copies the converged eigenvalues
and vectors to the output variables will seg-fault. Enlarging these just
means that the segfault will be moved elsewhere.
The fix is a simple 5 line change to xNEUPD.f in ARPACK, to force the M
output of xTRSEN to not exceed the original NCONV value. I've reported
the error to the maintainers of ARPACK and will see if there is any
follow-up. In any case with the above change I can no longer generate
the segfault even after millions of iterations of the segtest code.. I
attach the ARPACK patch, last version of eigs and the test code for
those that are interested.
Basically, before finally submitting this function I'd like to split
eigs up into a number of classes, and use these template classes to
allow both sparse and full cases to be treated. So I hope to have
something relatively soon. So after 18 months of on again off again work
on eigs, I think I finally have something I can be confident in :-)
Cheers
David
/*
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 2, 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, write to the Free
Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.
*/
// XXX FIXME XXX. Comment out HAVE_CONFIG_H for build outside of octave!!
#if 1
#include <config.h>
#define HAVE_ARPACK
#else
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#endif
#include <cfloat>
#include <cmath>
#include "ov.h"
#include "defun-dld.h"
#include "error.h"
#include "f77-fcn.h"
#include "quit.h"
#include "variables.h"
#include "SparsedbleLU.h"
#include "SparseCmplxLU.h"
#include "ov-re-sparse.h"
#include "ov-cx-sparse.h"
#include "MatrixType.h"
#include "oct-map.h"
#include "pager.h"
#include "SparsedbleCHOL.h"
#include "SparseCmplxCHOL.h"
#include "oct-rand.h"
#ifdef HAVE_ARPACK
// Arpack 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);
}
#endif
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
extern void
sparse_vector_product (const Sparse<double>&, const double*, double*);
extern void
sparse_vector_product (const Sparse<Complex>&, const Complex*, Complex*);
extern octave_idx_type
sparse_lu_solve (const SparseMatrix&, const SparseMatrix&, Matrix&);
extern octave_idx_type
sparse_lu_solve (const SparseComplexMatrix&, const SparseComplexMatrix&,
ComplexMatrix&);
extern ComplexMatrix
sparse_ltsolve (const SparseComplexMatrix&, const ColumnVector&,
const ComplexMatrix&);
extern Matrix
sparse_ltsolve (const SparseMatrix&, const ColumnVector&,
const Matrix&,);
extern ComplexMatrix
sparse_utsolve (const SparseComplexMatrix&, const ColumnVector&,
const ComplexMatrix&);
extern Matrix
sparse_utsolve (const SparseMatrix&, const ColumnVector&,
const Matrix&);
#endif
template <class T>
void
sparse_vector_product (const Sparse<T>& m, const T* x, T* 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];
}
template <class M, class SM>
octave_idx_type
sparse_lu_solve (const SM& L, const SM& U, M& m)
{
octave_idx_type nc = L.cols ();
octave_idx_type err = 0;
double rcond;
MatrixType ltyp (MatrixType::Lower);
MatrixType utyp (MatrixType::Upper);
m = L.solve (ltyp, m, err, rcond, 0);
if (err)
return err;
m = U.solve (utyp, m, err, rcond, 0);
if (err)
return err;
return err;
}
template <class SM, class M>
M
sparse_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>
M
sparse_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);
}
DEFUN_DLD (eigs, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} address@hidden = eigs (@var{a})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{k},
@var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{k},
@var{sigma},@var{opts})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b},
@var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b},
@var{k}, @var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{a}, @var{b},
@var{k}, @var{sigma}, @var{opts})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b}, @var{k})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{k}, @var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b}, @var{k}, @var{sigma})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{k}, @var{sigma}, @var{opts})\n\
@deftypefnx {Loadable Function} address@hidden = eigs (@var{af}, @var{n},
@var{b}, @var{k}, @var{sigma}, @var{opts})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}]} = eigs (@var{a},
@dots{})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}]} = eigs (@var{af},
@var{n}, @dots{})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}, @var{flag}]} = eigs
(@var{a}, @dots{})\n\
@deftypefnx {Loadable Function} address@hidden, @var{d}, @var{flag}]} = eigs
(@var{af}, @var{n}, @dots{})\n\
Calculate a limited number of eigenvalues and eigenvectors of @var{a},\n\
based on a selection criteria. The number eigenvalues and eigenvectors to\n\
calculate is given by @var{k} whose default value is 6.\n\
\n\
By default @code{eigs} solve the equation\n\
@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 problem).\n\
\n\
@item 'sa'\n\
Smallest Algebraic (valid only for real symmetric problem).\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 problem).\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 @address@hidden matrix whose columns are\n\
the @var{k} eigenvectors corresponding to the returned eigenvalues. The\n\
eigenvalues themselves are then returned in @var{d} in the form of a\n\
@address@hidden 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 ();
octave_function *func;
std::string fcn_name;
octave_idx_type n;
octave_idx_type k = 6;
Complex sigma;
double sigmar, sigmai;
bool have_sigma = false;
std::string typ = "LM";
SparseMatrix asmm, bsmm, bsmt;
SparseComplexMatrix ascm, bscm, bsct;
int b_arg;
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;
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';
octave_idx_type mode = 1;
bool note3 = false;
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 = ");
func = extract_function (args(0), "eigs", fcn_name, fname,
"; endfunction");
}
else
func = args(0).function_value ();
if (!func)
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
{
MatrixType tmp;
if (args(0).is_complex_type ())
{
ascm = (args(0).sparse_complex_matrix_value());
a_is_complex = true;
n = ascm.cols ();
tmp = MatrixType (ascm);
}
else
{
asmm = (args(0).sparse_matrix_value());
n = asmm.cols ();
tmp = MatrixType (asmm);
}
if ((tmp.is_diagonal () && !a_is_complex) || tmp.is_hermitian ())
symmetric = true;
}
// 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+arg_offset) &&
!(args(1+arg_offset).is_real_scalar ()))
if (args(1+arg_offset).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 lower 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)
bscm = args(b_arg).sparse_complex_matrix_value ();
else
bsmm = args(b_arg).sparse_matrix_value ();
}
if (!error_state)
{
if (a_is_complex || b_is_complex)
{
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);
}
}
else 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)
{
if (symmetric && !(a_is_complex || b_is_complex))
p = k * 2;
else
p = k * 2 + 1;
if (p < 20)
p = 20;
if (p > n - 1)
p = n - 1 ;
}
else if (symmetric && (p <= k || p > n))
error ("eigs: opts.p must be between k and n");
else if (!symmetric && (p < k || p > n))
error ("eigs: opts.p must be between k+1 and n");
if (!a_is_complex && !b_is_complex && symmetric)
{
if (k > n )
error ("eigs: Too many eigenvalues to extract (k >= n). Use
'eig(full(A))' instead");
}
else if (k > n - 1)
error ("eigs: Too many eigenvalues to extract (k >= n-1). Use
'eig(full(A))' instead");
if (! have_sigma)
{
if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
typ != "SI")
error ("eigs: unrecognized sigma value");
if ((a_is_complex || b_is_complex || !symmetric) &&
(typ == "LA" || typ == "SA" || typ == "BE"))
error ("eigs: invalid sigma value for complex or unsymmetric
problem");
if (!a_is_complex && !b_is_complex && symmetric &&
(typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR"))
error ("eigs: invalid sigma value for real symmetric problem");
// Mode 1 for SM mode seems unstable for some reason.
// Use Mode 3 instead, with sigma=0.
if (typ == "SM")
{
typ = "LM";
mode = 3;
}
if (have_b)
{
// See Note 3 dsaupd
note3 = true;
bmat = 'I';
octave_idx_type info;
if (a_is_complex || b_is_complex)
{
SparseComplexCHOL fact (bscm, info, false);
if (fact.P() != 0)
error ("eigs: The matrix B is not positive definite");
else
{
bscm = fact.L();
bsct = bscm.transpose();
permB = fact.perm() - 1.0;
}
}
else
{
SparseCHOL fact (bsmm, info, false);
if (fact.P() != 0)
error ("eigs: The matrix B is not positive definite");
else
{
bsmm = fact.L();
bsmt = bsmm.transpose();
permB = fact.perm() - 1.0;
}
}
}
}
else if (! std::abs (sigma))
typ = "SM";
else
mode = 3;
if (!error_state)
{
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;
if (!a_is_complex && !b_is_complex)
{
SparseMatrix L, U;
const octave_idx_type *P, *Q;
// Need to explicitly declare this here, so that the permutation
// vectors aren't deleted ahead of time.
SparseLU fact;
if (! have_a_fun && mode == 3)
{
// Caclulate LU decomposition of 'A - sigma * B'
SparseMatrix AminusSigmaB (asmm);
if (have_b)
{
if (cholB)
{
bsmt = bsmm.transpose ();
AminusSigmaB = AminusSigmaB - sigmar * bsmt * bsmm;
}
else
AminusSigmaB = AminusSigmaB - sigmar * bsmm;
}
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) = sigmar;
sigmat.xridx(i) = i;
sigmat.xcidx(i+1) = i + 1;
}
AminusSigmaB = AminusSigmaB - sigmat;
}
fact = SparseLU (AminusSigmaB);
L = fact.L ();
U = fact.U ();
P = fact.row_perm ();
Q = fact.col_perm ();
// 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;
}
if ((minU / maxU + 1.0) == 1.0)
{
warning ("eigs: 'A - sigma*B' is singular, indicating
sigma is exactly an eigenvalue.");
warning (" Convergence is not guaranteed");
}
}
if (symmetric)
{
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)
{
error ("eigs: unrecoverable exception encountered in
dsaupd");
goto eigs_err;
}
if (disp > 0 && !xisnan(workl[iptr(5)-1]))
{
if (iter++)
{
octave_stdout << "Iteration " << iter - 1 <<
": a few Ritz values of the " << p << "-by-" <<
p << " matrix\n";
for (int i = 0 ; i < k; i++)
octave_stdout << " " <<
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_a_fun)
{
double *ip2 = workd + iptr(0) - 1;
ColumnVector x(n);
for (octave_idx_type i = 0; i < n; i++)
x(i) = *ip2++;
octave_value_list tmp =
func->do_multi_index_op (1, octave_value (x));
ColumnVector y = tmp(0).column_vector_value ();
ip2 = workd + iptr(1) - 1;
for (octave_idx_type i = 0; i < n; i++)
*ip2++ = y(i);
}
else if (mode == 1)
{
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 = sparse_utsolve(bsmt, permB, asmm *
sparse_ltsolve(bsmm, permB, mtmp));
for (octave_idx_type i = 0; i < n; i++)
workd[i+iptr(1)-1] = mtmp(i,0);
}
else
sparse_vector_product (asmm,
workd + iptr(0) - 1,
workd + iptr(1) - 1);
}
else if (mode == 2)
error ("eigs: impossible");
else if (mode == 3)
{
if (have_b)
{
if (ido == -1)
{
if (cholB)
{
error ("fixme!!");
}
else
{
OCTAVE_LOCAL_BUFFER (double, dtmp, n);
sparse_vector_product
(asmm, workd+iptr(0)-1, dtmp);
Matrix tmp(n, 1);
for (octave_idx_type i = 0; i < n;
i++)
tmp(i,0) = dtmp[P[i]];
sparse_lu_solve (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)
sparse_vector_product (bsmm,
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]];
sparse_lu_solve (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]];
sparse_lu_solve (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)
error ("eigs: error %d in dsaupd", info);
break;
}
}
while (1);
if (!error_state)
{
octave_idx_type info2;
int rvec (nargout > 1);
// 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_vec (n, k);
double *z = eig_vec.fortran_vec ();
ColumnVector eig_val (k);
double *d = eig_val.fortran_vec ();
F77_FUNC (dseupd, DSEUPD)
(rvec, F77_CONST_CHAR_ARG2 ("A", 1),
sel, d, z, n, sigmar,
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)
{
error ("eigs: unrecoverable exception encountered in
dseupd");
goto eigs_err;
}
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];
}
}
retval(2) = double (info);
retval(1) = DiagMatrix (eig_val);
if (note3)
eig_vec = sparse_ltsolve(bsmm, permB, eig_vec);
retval(0) = eig_vec;
}
else
retval (0) = eig_val;
}
else
error ("eigs: error %d in dseupd", info2);
}
}
else
{
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)
{
error ("eigs: unrecoverable exception encountered in
dnaupd");
goto eigs_err;
}
if (disp > 0 && !xisnan(workl[iptr(5)-1]))
{
if (iter++)
{
octave_stdout << "Iteration " << iter - 1 <<
": a few Ritz values of the " << p << "-by-" <<
p << " matrix\n";
for (int i = 0 ; i < k; i++)
octave_stdout << " " <<
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_a_fun)
{
double *ip2 = workd + iptr(0) - 1;
ColumnVector x(n);
for (octave_idx_type i = 0; i < n; i++)
x(i) = *ip2++;
octave_value_list tmp =
func->do_multi_index_op (1, octave_value (x));
ColumnVector y = tmp(0).column_vector_value ();
ip2 = workd + iptr(1) - 1;
for (octave_idx_type i = 0; i < n; i++)
*ip2++ = y(i);
}
else if (mode == 1)
{
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 = sparse_utsolve(bsmt, permB, asmm *
sparse_ltsolve(bsmm, permB, mtmp));
for (octave_idx_type i = 0; i < n; i++)
workd[i+iptr(1)-1] = mtmp(i,0);
}
else
sparse_vector_product (asmm,
workd + iptr(0) - 1,
workd + iptr(1) - 1);
}
else if (mode == 2)
{
error ("eigs: impossible");
}
else if (mode == 3)
{
if (have_b)
{
if (ido == -1)
{
if (cholB)
{
error ("fixme!!");
}
else
{
OCTAVE_LOCAL_BUFFER (double, dtmp, n);
sparse_vector_product
(asmm, workd+iptr(0)-1, dtmp);
Matrix tmp(n, 1);
for (octave_idx_type i = 0; i < n;
i++)
tmp(i,0) = dtmp[P[i]];
sparse_lu_solve (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)
sparse_vector_product (bsmm,
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]];
sparse_lu_solve (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]];
sparse_lu_solve (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)
error ("eigs: error %d in dnaupd", info);
break;
}
}
while (1);
if (!error_state)
{
octave_idx_type info2;
int rvec (nargout > 1);
// 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 to avoid
// problems!!
Array<int> s (p);
int *sel = s.fortran_vec ();
Matrix eig_vec (n, k + 1);
double *z = eig_vec.fortran_vec ();
OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
OCTAVE_LOCAL_BUFFER (double, di, k + 1);
OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
F77_FUNC (dneupd, DNEUPD)
(rvec, F77_CONST_CHAR_ARG2 ("A", 1),
sel, dr, di, z, n, sigmar, sigmai, workev,
F77_CONST_CHAR_ARG2 (&bmat, 1), n,
F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
k, tol, presid, p, v, n, iparam,
ipntr, workd, workl, lwork, info2
F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
F77_CHAR_ARG_LEN(2));
if (f77_exception_encountered)
{
error ("eigs: unrecoverable exception encountered in
dneupd");
goto eigs_err;
}
ComplexColumnVector eig_val (k+1);
Complex *d = eig_val.fortran_vec ();
if (info2 == 0)
{
octave_idx_type j = 0;
for (octave_idx_type i = 0; i < k+1; i++)
{
if (dr[i] == 0.0 && di[i] == 0.0 && j == 0)
j++;
else
d [i-j] = Complex (dr[i], di[i]);
}
if (j == 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 (nargout < 2)
retval (0) = eig_val;
else
{
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];
}
ComplexMatrix eig_vec2 (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_vec2(j,i) =
Complex(z[j+off1],0.);
i++;
}
else
{
for (octave_idx_type j = 0; j < n; j++)
{
eig_vec2(j,i) =
Complex(z[j+off1],z[j+off2]);
if (i < k - 1)
eig_vec2(j,i+1) =
Complex(z[j+off1],-z[j+off2]);
}
i+=2;
}
}
retval(2) = double (info);
retval(1) = ComplexDiagMatrix(eig_val);
if (note3)
eig_vec2 =
sparse_ltsolve(SparseComplexMatrix(bsmm),
permB, eig_vec2);
retval(0) = eig_vec2;
}
}
else
error ("eigs: error %d in dneupd", info2);
}
}
}
else
{
SparseComplexMatrix L, U;
const octave_idx_type *P, *Q;
// Need to explicitly declare this here, so that the permutation
// vectors aren't deleted ahead of time.
SparseComplexLU fact;
if (! have_a_fun && mode == 3)
{
// Caclulate LU decomposition of 'A - sigma * B'
SparseComplexMatrix AminusSigmaB (ascm);
if (have_b)
{
if (cholB)
{
bsct = bscm.transpose ();
AminusSigmaB = AminusSigmaB - sigma * bsct * bscm;
}
else
AminusSigmaB = AminusSigmaB - sigma * bscm;
}
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;
}
fact = SparseComplexLU (AminusSigmaB);
L = fact.L ();
U = fact.U ();
P = fact.row_perm ();
Q = fact.col_perm ();
// 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;
}
if ((minU / maxU + 1.0) == 1.0)
{
warning ("eigs: 'A - sigma*B' is singular, indicating
sigma is exactly an eigenvalue.");
warning (" Convergence is not guaranteed");
}
}
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)
{
error ("eigs: unrecoverable exception encountered in
znaupd");
goto eigs_err;
}
if (disp > 0 && !xisnan(workl[iptr(5)-1]))
{
if (iter++)
{
octave_stdout << "Iteration " << iter - 1 <<
": a few Ritz values of the " << p << "-by-" <<
p << " matrix\n";
for (int i = 0 ; i < k; i++)
octave_stdout << " " <<
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_a_fun)
{
Complex *ip2 = workd + iptr(0) - 1;
ComplexColumnVector x(n);
for (octave_idx_type i = 0; i < n; i++)
x(i) = *ip2++;
octave_value_list tmp =
func->do_multi_index_op (1, octave_value (x));
ComplexColumnVector y =
tmp(0).complex_column_vector_value ();
ip2 = workd + iptr(1) - 1;
for (octave_idx_type i = 0; i < n; i++)
*ip2++ = y(i);
}
else if (mode == 1)
{
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 = sparse_utsolve(bsct, permB, ascm *
sparse_ltsolve(bscm, permB, mtmp));
for (octave_idx_type i = 0; i < n; i++)
workd[i+iptr(1)-1] = mtmp(i,0);
}
else
sparse_vector_product (ascm, workd + iptr(0) - 1,
workd + iptr(1) - 1);
}
else if (mode == 2)
{
error ("eigs: Impossible");
}
else if (mode == 3)
{
if (have_b)
{
if (ido == -1)
{
if (cholB)
{
error ("fixme!!");
}
else
{
OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
sparse_vector_product
(ascm, workd+iptr(0)-1, ctmp);
ComplexMatrix tmp(n, 1);
for (octave_idx_type i = 0; i < n; i++)
tmp(i,0) = ctmp[P[i]];
sparse_lu_solve (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)
sparse_vector_product (bscm,
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]];
sparse_lu_solve (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]];
sparse_lu_solve (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)
error ("eigs: error %d in znaupd", info);
break;
}
}
while (1);
if (!error_state)
{
octave_idx_type info2;
int rvec (nargout > 1);
// 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 to avoid
// problems!!
Array<int> s (p);
int *sel = s.fortran_vec ();
ComplexMatrix eig_vec (n, k);
Complex *z = eig_vec.fortran_vec ();
ComplexColumnVector eig_val (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)
{
error ("eigs: unrecoverable exception encountered in
zneupd");
goto eigs_err;
}
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 (nargout < 2)
retval (0) = eig_val;
else
{
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];
}
retval(2) = double (info);
retval(1) = ComplexDiagMatrix (eig_val);
if (note3)
eig_vec = sparse_ltsolve(bscm, permB, eig_vec);
retval(0) = eig_vec;
}
}
else
error ("eigs: error %d in zneupd", info2);
}
}
if (ip(4) == 0)
warning ("eigs: None of the %d requested eigenvalues converged", k);
else if (ip(4) < k)
warning ("eigs: Only %d of the %d requested eigenvalues converged",
ip(4), k);
}
}
eigs_err:
if (! fcn_name.empty ())
clear_function (fcn_name);
#else
error ("eigs: not available in this version of Octave");
#endif
return retval;
}
template void
sparse_vector_product (const Sparse<double>&, const double*, double*);
template void
sparse_vector_product (const Sparse<Complex>&, const Complex*, Complex*);
template octave_idx_type
sparse_lu_solve (const SparseMatrix&, const SparseMatrix&, Matrix&);
template octave_idx_type
sparse_lu_solve (const SparseComplexMatrix&, const SparseComplexMatrix&,
ComplexMatrix&);
template ComplexMatrix
sparse_ltsolve (const SparseComplexMatrix&, const ColumnVector&,
const ComplexMatrix&);
template Matrix
sparse_ltsolve (const SparseMatrix&, const ColumnVector&,
const Matrix&);
template ComplexMatrix
sparse_utsolve (const SparseComplexMatrix&, const ColumnVector&,
const ComplexMatrix&);
template Matrix
sparse_utsolve (const SparseMatrix&, const ColumnVector&,
const Matrix&);
/*
%% 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);
%!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
%! [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, n mist be odd to force eigs to use dnaupd/dneupd
%!shared n, k, A, d0
%! n = 20;
%! k = 4;
%! ## FIXME I can get the matrix below to seg-fault the tests and give
%! ## incorrect results for "SI" and very very ocassionally for "SM"..
%! 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)]);
%! ## 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);
%!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
%! [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);
%!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
%! [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
*/
/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
function aerr = segtest (iter)
%% This will seg-fault octave consistently, but not matlab.
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)]);
opts.disp = 0; opts.p = 19;
aerr = 0;
for i=1:iter
[v1,d1] = eigs(A, k, 'sr', opts);
d1 = diag(d1);
merr = 0;
for i=1:k
newerr = max(abs((A - d1(i)*speye(n))*v1(:,i)));
if (newerr > merr)
merr = newerr;
end
end
fprintf('Max Err: %g\n', merr);
if (merr > aerr)
aerr = merr;
end
end
end
*** ARPACK.orig/SRC/dneupd.f 2000-09-20 22:58:34.000000000 +0200
--- ARPACK/SRC/dneupd.f 2006-11-03 11:27:21.159925329 +0100
***************
*** 353,359 ****
& mode , msglvl, outncv, ritzr ,
& ritzi , wri , wrr , irr ,
& iri , ibd , ishift, numcnv ,
! & np , jj
logical reord
Double precision
& conds , rnorm, sep , temp,
--- 353,359 ----
& mode , msglvl, outncv, ritzr ,
& ritzi , wri , wrr , irr ,
& iri , ibd , ishift, numcnv ,
! & np , jj , nconv2
logical reord
Double precision
& conds , rnorm, sep , temp,
***************
*** 661,676 ****
& workl(iuptri), ldh ,
& workl(invsub), ldq ,
& workl(iheigr), workl(iheigi),
! & nconv , conds ,
& sep , workl(ihbds) ,
& ncv , iwork ,
& 1 , ierr)
c
if (ierr .eq. 1) then
info = 1
go to 9000
end if
c
if (msglvl .gt. 2) then
call dvout (logfil, ncv, workl(iheigr), ndigit,
& '_neupd: Real part of the eigenvalues of H--reordered')
--- 661,681 ----
& workl(iuptri), ldh ,
& workl(invsub), ldq ,
& workl(iheigr), workl(iheigi),
! & nconv2 , conds ,
& sep , workl(ihbds) ,
& ncv , iwork ,
& 1 , ierr)
c
+ if (nconv2 .lt. nconv) then
+ nconv = nconv2
+ end if
+
if (ierr .eq. 1) then
info = 1
go to 9000
end if
c
+
if (msglvl .gt. 2) then
call dvout (logfil, ncv, workl(iheigr), ndigit,
& '_neupd: Real part of the eigenvalues of H--reordered')
*** ARPACK.orig/SRC/sneupd.f 2000-09-20 22:58:34.000000000 +0200
--- ARPACK/SRC/sneupd.f 2006-11-03 11:32:54.691913071 +0100
***************
*** 353,359 ****
& mode , msglvl, outncv, ritzr ,
& ritzi , wri , wrr , irr ,
& iri , ibd , ishift, numcnv ,
! & np , jj
logical reord
Real
& conds , rnorm, sep , temp,
--- 353,359 ----
& mode , msglvl, outncv, ritzr ,
& ritzi , wri , wrr , irr ,
& iri , ibd , ishift, numcnv ,
! & np , jj , nconv2
logical reord
Real
& conds , rnorm, sep , temp,
***************
*** 661,671 ****
& workl(iuptri), ldh ,
& workl(invsub), ldq ,
& workl(iheigr), workl(iheigi),
! & nconv , conds ,
& sep , workl(ihbds) ,
& ncv , iwork ,
& 1 , ierr)
c
if (ierr .eq. 1) then
info = 1
go to 9000
--- 661,675 ----
& workl(iuptri), ldh ,
& workl(invsub), ldq ,
& workl(iheigr), workl(iheigi),
! & nconv2 , conds ,
& sep , workl(ihbds) ,
& ncv , iwork ,
& 1 , ierr)
c
+ if (nconv2 .lt. nconv) then
+ nconv = nconv2
+ end if
+
if (ierr .eq. 1) then
info = 1
go to 9000
*** ARPACK.orig/SRC/zneupd.f 2002-08-15 07:51:12.000000000 +0200
--- ARPACK/SRC/zneupd.f 2006-11-03 11:33:28.032312460 +0100
***************
*** 301,307 ****
& invsub, iuptri, iwev , j , ldh , ldq ,
& mode , msglvl, ritz , wr , k , irz ,
& ibd , outncv, iq , np , numcnv, jj ,
! & ishift
Complex*16
& rnorm, temp, vl(1)
Double precision
--- 301,307 ----
& invsub, iuptri, iwev , j , ldh , ldq ,
& mode , msglvl, ritz , wr , k , irz ,
& ibd , outncv, iq , np , numcnv, jj ,
! & ishift, nconv2
Complex*16
& rnorm, temp, vl(1)
Double precision
***************
*** 592,600 ****
call ztrsen('None' , 'V' , select ,
& ncv , workl(iuptri), ldh ,
& workl(invsub), ldq , workl(iheig),
! & nconv , conds , sep ,
& workev , ncv , ierr)
c
if (ierr .eq. 1) then
info = 1
go to 9000
--- 592,604 ----
call ztrsen('None' , 'V' , select ,
& ncv , workl(iuptri), ldh ,
& workl(invsub), ldq , workl(iheig),
! & nconv2 , conds , sep ,
& workev , ncv , ierr)
c
+ if (nconv2 .lt. nconv) then
+ nconv = nconv2
+ end if
+
if (ierr .eq. 1) then
info = 1
go to 9000
*** ARPACK.orig/SRC/cneupd.f 2002-08-15 07:51:12.000000000 +0200
--- ARPACK/SRC/cneupd.f 2006-11-03 11:31:54.499802784 +0100
***************
*** 301,307 ****
& invsub, iuptri, iwev , j , ldh , ldq ,
& mode , msglvl, ritz , wr , k , irz ,
& ibd , outncv, iq , np , numcnv, jj ,
! & ishift
Complex
& rnorm, temp, vl(1)
Real
--- 301,307 ----
& invsub, iuptri, iwev , j , ldh , ldq ,
& mode , msglvl, ritz , wr , k , irz ,
& ibd , outncv, iq , np , numcnv, jj ,
! & ishift, nconv2
Complex
& rnorm, temp, vl(1)
Real
***************
*** 592,600 ****
call ctrsen('None' , 'V' , select ,
& ncv , workl(iuptri), ldh ,
& workl(invsub), ldq , workl(iheig),
! & nconv , conds , sep ,
& workev , ncv , ierr)
c
if (ierr .eq. 1) then
info = 1
go to 9000
--- 592,604 ----
call ctrsen('None' , 'V' , select ,
& ncv , workl(iuptri), ldh ,
& workl(invsub), ldq , workl(iheig),
! & nconv2 , conds , sep ,
& workev , ncv , ierr)
c
+ if (nconv2 .lt. nconv) then
+ nconv = nconv2
+ end if
+
if (ierr .eq. 1) then
info = 1
go to 9000
- Eigs again - Success!!,
David Bateman <=