26#ifndef OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
27#define OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
36#include "AnasaziOutputManager.hpp"
37#include "AnasaziBasicOutputManager.hpp"
38#include "AnasaziConfigDefs.hpp"
39#include "AnasaziOperator.hpp"
40#include "AnasaziMultiVec.hpp"
41#include "AnasaziBasicEigenproblem.hpp"
42#include "AnasaziBlockKrylovSchurSolMgr.hpp"
51 class MultiVecTraits<double,
oomph::DoubleMultiVector>
56 static Teuchos::RCP<oomph::DoubleMultiVector>
Clone(
64 static Teuchos::RCP<oomph::DoubleMultiVector>
CloneCopy(
74 static Teuchos::RCP<oomph::DoubleMultiVector>
CloneCopy(
82 static Teuchos::RCP<oomph::DoubleMultiVector>
CloneCopy(
113 static Teuchos::RCP<const oomph::DoubleMultiVector>
CloneView(
124 static Teuchos::RCP<oomph::DoubleMultiVector>
CloneView(
135 static Teuchos::RCP<oomph::DoubleMultiVector>
CloneView(
144 return static_cast<int>(mv.
nrow());
150 return static_cast<int>(mv.
nvector());
157 const Teuchos::SerialDenseMatrix<int, double>&
B,
168 int mv_n_vector = mv.
nvector();
174 for (
int i = 0;
i < n_row_local; ++
i)
176 for (
int j = 0; j < mv_n_vector; j++)
179 for (
int k = 0; k < A_n_vector; k++)
181 res += C(k,
i) *
B(k, j);
195 const unsigned n_vector = A.
nvector();
197 for (
unsigned v = 0; v < n_vector; v++)
199 for (
unsigned i = 0;
i < n_row_local;
i++)
201 mv(v,
i) = alpha * A(v,
i) + beta *
B(v,
i);
217 const std::vector<double>& alpha)
219 const unsigned n_vector = mv.
nvector();
221 for (
unsigned v = 0; v < n_vector; v++)
223 for (
unsigned i = 0;
i < n_row_local;
i++)
225 mv(v,
i) *= alpha[v];
236 Teuchos::SerialDenseMatrix<int, double>&
B)
238 const unsigned A_nvec = A.
nvector();
240 const unsigned mv_nvec = mv.
nvector();
243 for (
unsigned v1 = 0; v1 < A_nvec; v1++)
245 for (
unsigned v2 = 0; v2 < mv_nvec; v2++)
248 for (
unsigned i = 0;
i < A_nrow_local;
i++)
250 B(v1, v2) += A(v1,
i) * mv(v2,
i);
261 const int n_total_val = A_nvec * mv_nvec;
263 double b[n_total_val];
264 double b2[n_total_val];
266 for (
unsigned v1 = 0; v1 < A_nvec; v1++)
268 for (
unsigned v2 = 0; v2 < mv_nvec; v2++)
270 b[count] =
B(v1, v2);
271 b2[count] = b[count];
285 for (
unsigned v1 = 0; v1 < A_nvec; v1++)
287 for (
unsigned v2 = 0; v2 < mv_nvec; v2++)
289 B(v1, v2) = b2[count];
304 std::vector<double>& b)
315 std::vector<double>& normvec)
332 const std::vector<int>& index,
336 const unsigned n_index = index.size();
342 for (
unsigned v = 0; v < n_index; v++)
344 for (
unsigned i = 0;
i < n_row_local;
i++)
346 mv(index[v],
i) = A(v,
i);
364 const Teuchos::Range1D& index,
368 const unsigned n_index = index.size();
374 unsigned range_index = index.lbound();
375 for (
unsigned v = 0; v < n_index; v++)
377 for (
unsigned i = 0;
i < n_row_local;
i++)
379 mv(range_index,
i) = A(v,
i);
398 const unsigned n_vector = mv.
nvector();
400 for (
unsigned v = 0; v < n_vector; v++)
402 for (
unsigned i = 0;
i < n_row_local;
i++)
404 mv(v,
i) = Teuchos::ScalarTraits<double>::random();
413 const double alpha = Teuchos::ScalarTraits<double>::zero())
465 class OperatorTraits<double,
466 oomph::DoubleMultiVector,
506 const double& sigma = 0.0)
526 const unsigned n_vec = x.
nvector();
545 for (
unsigned i = 0;
i < n_row_local;
i++)
551 for (
unsigned v = 1; v < n_vec; ++v)
558 for (
unsigned i = 0;
i < n_row_local;
i++)
574 typedef Teuchos::ScalarTraits<ST>
SCT;
575 typedef SCT::magnitudeType
MT;
576 typedef Anasazi::MultiVec<ST>
MV;
577 typedef Anasazi::Operator<ST>
OP;
578 typedef Anasazi::MultiVecTraits<ST, MV>
MVT;
579 typedef Anasazi::OperatorTraits<ST, MV, OP>
OPT;
622 Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
627 << Anasazi::Anasazi_Version() << std::endl
664 Vector<std::complex<double>>& eigenvalue,
674 Teuchos::RCP<DoubleMultiVector> initial = Teuchos::rcp(
676 Anasazi::MultiVecTraits<double, DoubleMultiVector>::MvRandom(*initial);
679 Teuchos::RCP<DoubleMultiVectorOperator> Op_pt =
684 Teuchos::RCP<Anasazi::BasicEigenproblem<double,
687 anasazi_pt = Teuchos::rcp(
688 new Anasazi::BasicEigenproblem<
double,
694 anasazi_pt->setHermitian(
false);
697 anasazi_pt->setNEV(n_eval);
700 if (anasazi_pt->setProblem() !=
true)
702 oomph_info <<
"Anasazi::BasicEigenproblem::setProblem() had an error."
704 <<
"End Result: TEST FAILED" << std::endl;
712 Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
714 Teuchos::ParameterList MyPL;
715 MyPL.set(
"Which",
"LM");
716 MyPL.set(
"Block Size", 1);
718 MyPL.set(
"Maximum Restarts", 500);
719 MyPL.set(
"Orthogonalization",
"DGKS");
720 MyPL.set(
"Verbosity", verbosity);
721 MyPL.set(
"Convergence Tolerance", tol);
722 Anasazi::BlockKrylovSchurSolMgr<double,
725 BKS(anasazi_pt, MyPL);
728 Anasazi::ReturnType ret = BKS.solve();
729 bool testFailed =
false;
730 if (ret != Anasazi::Converged)
741 Anasazi::Eigensolution<double, DoubleMultiVector> sol =
742 anasazi_pt->getSolution();
743 std::vector<Anasazi::Value<double>> evals = sol.Evals;
744 Teuchos::RCP<DoubleMultiVector> evecs = sol.Evecs;
746 eigenvalue.resize(evals.size());
747 eigenvector.resize(evals.size());
749 for (
unsigned i = 0;
i < evals.size();
i++)
752 double a = evals[
i].realpart;
753 double b = evals[
i].imagpart;
754 double det = a * a + b * b;
755 eigenvalue[
i] = std::complex<double>(a / det +
Sigma, -b / det);
758 eigenvector[
i].build(evecs->distribution_pt());
763 eigenvector[
i][n] = (*evecs)(
i, n);
static Teuchos::RCP< oomph::DoubleMultiVector > CloneViewNonConst(oomph::DoubleMultiVector &mv, const Teuchos::Range1D &index)
Creates a new oomph::DoubleMultiVector that contains shallow copies of selected entries of the oomph:...
static int GetVecLength(const oomph::DoubleMultiVector &mv)
Obtain the global length of the vector.
static Teuchos::RCP< oomph::DoubleMultiVector > CloneCopy(const oomph::DoubleMultiVector &mv)
Creates a deep copy of the DoubleMultiVector mv return Reference-counted pointer to the new DoubleMul...
static Teuchos::RCP< oomph::DoubleMultiVector > CloneCopy(const oomph::DoubleMultiVector &mv, const std::vector< int > &index)
Creates a new oomph::DoubleMultiVector and (deep) copies the contents of the vectors in index into th...
static Teuchos::RCP< const oomph::DoubleMultiVector > CloneView(const oomph::DoubleMultiVector &mv, const std::vector< int > &index)
Creates a new oomph::DoubleMultiVector that contains shallow copies of selected entries of the oomph:...
static Teuchos::RCP< oomph::DoubleMultiVector > CloneView(oomph::DoubleMultiVector &mv, const std::vector< int > &index)
Creates a new oomph::DoubleMultiVector that contains shallow copies of selected entries of the oomph:...
static void MvTransMv(const double alpha, const oomph::DoubleMultiVector &A, const oomph::DoubleMultiVector &mv, Teuchos::SerialDenseMatrix< int, double > &B)
Compute a dense matrix B through the matrix-matrix multiply .
static void MvTimesMatAddMv(const double alpha, const oomph::DoubleMultiVector &A, const Teuchos::SerialDenseMatrix< int, double > &B, const double beta, oomph::DoubleMultiVector &mv)
Update mv with .
static void MvPrint(const oomph::DoubleMultiVector &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
static Teuchos::RCP< oomph::DoubleMultiVector > CloneCopy(const oomph::DoubleMultiVector &mv, const Teuchos::Range1D &index)
Deep copy of specified columns of oomph::DoubleMultiVector return Reference-counted pointer to the ne...
static void SetBlock(const oomph::DoubleMultiVector &A, const Teuchos::Range1D &index, oomph::DoubleMultiVector &mv)
Deep copy of A into specified columns of mv.
static void MvScale(oomph::DoubleMultiVector &mv, const std::vector< double > &alpha)
Scale each element of the i-th vector in mv with alpha[i].
static Teuchos::RCP< oomph::DoubleMultiVector > CloneView(oomph::DoubleMultiVector &mv, const Teuchos::Range1D &index)
Creates a new oomph::DoubleMultiVector that contains shallow copies of selected entries of the oomph:...
static void MvScale(oomph::DoubleMultiVector &mv, const double alpha)
Scale each element of the vectors in mv with alpha.
static void MvAddMv(const double alpha, const oomph::DoubleMultiVector &A, const double beta, const oomph::DoubleMultiVector &B, oomph::DoubleMultiVector &mv)
Replace mv with .
static int GetNumberVecs(const oomph::DoubleMultiVector &mv)
Obtain the number of vectors in the multivector.
static void SetBlock(const oomph::DoubleMultiVector &A, const std::vector< int > &index, oomph::DoubleMultiVector &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index.
static void MvNorm(const oomph::DoubleMultiVector &mv, std::vector< double > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ,...
static void MvDot(const oomph::DoubleMultiVector &mv, const oomph::DoubleMultiVector &A, std::vector< double > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void MvRandom(oomph::DoubleMultiVector &mv)
Replace the vectors in mv with random vectors.
static Teuchos::RCP< oomph::DoubleMultiVector > CloneViewNonConst(oomph::DoubleMultiVector &mv, const std::vector< int > &index)
Creates a new oomph::DoubleMultiVector that contains shallow copies of selected entries of the oomph:...
static Teuchos::RCP< oomph::DoubleMultiVector > Clone(const oomph::DoubleMultiVector &mv, const int numvecs)
Creates a new empty DoubleMultiVector containing numvecs columns. Return a reference-counted pointer ...
static void MvInit(oomph::DoubleMultiVector &mv, const double alpha=Teuchos::ScalarTraits< double >::zero())
Replace each element of the vectors in mv with alpha.
static void Assign(const oomph::DoubleMultiVector &A, oomph::DoubleMultiVector &mv)
mv := A
static void Apply(const oomph::DoubleMultiVectorOperator &Op, const oomph::DoubleMultiVector &x, oomph::DoubleMultiVector &y)
Matrix operator application method.
Class for the Anasazi eigensolver.
const int & narnoldi() const
Access function for the number of Arnoldi vectors (const version)
void track_eigenvalue_magnitude()
Set the magnitude to be the quantity of interest.
LinearSolver * Default_linear_solver_pt
Pointer to a default linear solver.
Anasazi::OutputManager< ST > * Output_manager_pt
Pointer to output manager.
Teuchos::ScalarTraits< ST > SCT
void get_eigenvalues_right_of_shift()
Set the desired eigenvalues to be right of the shift.
int Spectrum
Integer to set whether the real, imaginary or magnitude is required to be small or large.
Anasazi::OperatorTraits< ST, MV, OP > OPT
virtual ~ANASAZI()
Destructor, delete the linear solver.
bool Compute_eigenvectors
Boolean to indicate whether or not to compute the eigenvectors.
int NArnoldi
Number of Arnoldi vectors to compute.
LinearSolver * Linear_solver_pt
Pointer to a linear solver.
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
int & narnoldi()
Access function for the number of Arnoldi vectors.
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
double Sigma
Set the shifted value.
void disable_compute_eigenvectors()
Set to disable the computation of the eigenvectors.
void track_eigenvalue_real_part()
Set the real part to be the quantity of interest (default)
void get_eigenvalues_left_of_shift()
Set the desired eigenvalues to be left of the shift.
Anasazi::MultiVec< ST > MV
void enable_compute_eigenvectors()
Set to enable the computation of the eigenvectors (default)
void track_eigenvalue_imaginary_part()
Set the imaginary part fo the quantity of interest.
Anasazi::MultiVecTraits< ST, MV > MVT
bool Small
Boolean to set which part of the spectrum left (default) or right of the shifted value.
ANASAZI(const ANASAZI &)
Empty copy constructor.
LinearSolver *const & linear_solver_pt() const
Return a pointer to the linear solver object (const version)
Anasazi::Operator< ST > OP
A class for compressed row matrices. This is a distributable object.
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
bool distributed() const
distribution is serial or distributed
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
Base class for Oomph-lib's Vector Operator classes that will be used with the DoubleMultiVector.
DoubleMultiVectorOperator()
Empty constructor.
virtual ~DoubleMultiVectorOperator()
virtual destructor
virtual void apply(const DoubleMultiVector &x, DoubleMultiVector &y) const =0
The apply interface.
A multi vector in the mathematical sense, initially developed for linear algebra type applications....
void output(std::ostream &outfile) const
output the contents of the vector
void initialise(const double &initial_value)
initialise the whole vector with value v
void dot(const DoubleMultiVector &vec, std::vector< double > &result) const
compute the 2 norm of this vector
void norm(std::vector< double > &result) const
compute the 2 norm of this vector
unsigned nvector() const
Return the number of vectors.
DoubleVector & doublevector(const unsigned &i)
access to the DoubleVector representatoin
A vector in the mathematical sense, initially developed for linear algebra type applications....
Base class for all EigenProblem solves. This simply defines standard interfaces so that different sol...
double Sigma_real
Double value that represents the real part of the shift in shifted eigensolvers.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
Base class for all linear solvers. This merely defines standard interfaces for linear solvers,...
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the rhs vector. Solution is returned in...
void disable_doc_time()
Disable documentation of solve times.
Class for the shift invert operation.
void apply(const DoubleMultiVector &x, DoubleMultiVector &y) const
The apply interface.
ProblemBasedShiftInvertOperator(Problem *const &problem_pt, LinearSolver *const &linear_solver_pt, const double &sigma=0.0)
CRDoubleMatrix * AsigmaM_pt
LinearSolver * Linear_solver_pt
////////////////////////////////////////////////////////////////// //////////////////////////////////...
virtual void get_eigenproblem_matrices(CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
Get the matrices required by a eigensolver. If the shift parameter is non-zero the second matrix will...
LinearAlgebraDistribution *const & dof_distribution_pt() const
Return the pointer to the dof distribution (read-only)
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...