trilinos_eigen_solver.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
31 #define OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
32 
33 //oomph-lib includes
34 #include "double_multi_vector.h"
35 #include "problem.h"
36 #include "linear_solver.h"
37 #include "eigen_solver.h"
38 
39 //Trilinos includes
40 #include "AnasaziOutputManager.hpp"
41 #include "AnasaziBasicOutputManager.hpp"
42 #include "AnasaziConfigDefs.hpp"
43 #include "AnasaziOperator.hpp"
44 #include "AnasaziMultiVec.hpp"
45 #include "AnasaziBasicEigenproblem.hpp"
46 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
47 
48 namespace Anasazi
49 {
50  //========================================================================
51  ///Specialize the Anasazi traits class for the oomph-lib DoubleMultiVector.
52  ///This provides the interfaces required by the Anasazi eigensolvers.
53  //========================================================================
54  template<>
55  class MultiVecTraits<double,oomph::DoubleMultiVector>
56  {
57  public:
58 
59  /// \brief Creates a new empty DoubleMultiVector containing numvecs columns.
60  /// Return a reference-counted pointer to the new multivector
61  static Teuchos::RCP<oomph::DoubleMultiVector> Clone(
62  const oomph::DoubleMultiVector& mv, const int numvecs )
63  {
64  return Teuchos::rcp(new oomph::DoubleMultiVector(numvecs,mv));
65  }
66 
67  /// \brief Creates a deep copy of the DoubleMultiVector mv
68  /// return Reference-counted pointer to the new DoubleMultiVector
69  static Teuchos::RCP<oomph::DoubleMultiVector>
71  {
72  return Teuchos::rcp(new oomph::DoubleMultiVector(mv));
73  }
74 
75  /// \brief Creates a new oomph::DoubleMultiVector and (deep)
76  /// copies the contents of
77  /// the vectors in index into the new vector
78  ///return Reference-counted pointer to the new oomph::DoubleMultiVector
79  static Teuchos::RCP<oomph::DoubleMultiVector>
81  const std::vector<int>& index )
82  {
83  return Teuchos::rcp(new oomph::DoubleMultiVector(mv,index));
84  }
85 
86  /// \brief Deep copy of specified columns of oomph::DoubleMultiVector
87  /// return Reference-counted pointer to the new oomph::DoubleMultiVector
88  static Teuchos::RCP<oomph::DoubleMultiVector> CloneCopy(
89  const oomph::DoubleMultiVector& mv, const Teuchos::Range1D& index )
90  {
91  return Teuchos::rcp(new oomph::DoubleMultiVector(mv,index));
92  }
93 
94  /// \brief Creates a new oomph::DoubleMultiVector that contains
95  /// shallow copies
96  /// of selected entries of the oomph::DoubleMultiVector mv
97  /// return Reference-counted pointer to the new oomph::DoubleMultiVector
98  static Teuchos::RCP<oomph::DoubleMultiVector>
100  const std::vector<int>& index )
101  {
102  return Teuchos::rcp(new oomph::DoubleMultiVector(mv,index,false));
103  }
104 
105  /// \brief Creates a new oomph::DoubleMultiVector that
106  /// contains shallow copies
107  /// of selected entries of the oomph::DoubleMultiVector mv
108  /// return Reference-counted pointer to the new oomph::DoubleMultiVector
109  static Teuchos::RCP<oomph::DoubleMultiVector>
111  const Teuchos::Range1D& index )
112  {
113  return Teuchos::rcp(new oomph::DoubleMultiVector(mv,index,false));
114  }
115 
116  /// \brief Creates a new oomph::DoubleMultiVector that contains
117  /// shallow copies
118  /// of selected entries of the oomph::DoubleMultiVector mv
119  /// return Reference-counted pointer to the new oomph::DoubleMultiVector
120  /// (const version)
121  static Teuchos::RCP<const oomph::DoubleMultiVector>
123  const std::vector<int>& index )
124  {
125  return Teuchos::rcp(new oomph::DoubleMultiVector(mv,index,false));
126  }
127 
128  /// \brief Creates a new oomph::DoubleMultiVector that contains
129  /// shallow copies
130  /// of selected entries of the oomph::DoubleMultiVector mv
131  /// return Reference-counted pointer to the new oomph::DoubleMultiVector
132  /// (Non-const version for Trilinos 9 interface)
133  static Teuchos::RCP<oomph::DoubleMultiVector>
135  const std::vector<int>& index )
136  {
137  return Teuchos::rcp(new oomph::DoubleMultiVector(mv,index,false));
138  }
139 
140  /// \brief Creates a new oomph::DoubleMultiVector that
141  /// contains shallow copies
142  /// of selected entries of the oomph::DoubleMultiVector mv
143  /// return Reference-counted pointer to the new oomph::DoubleMultiVector
144  /// (const version)
145  static Teuchos::RCP<oomph::DoubleMultiVector>
146  CloneView( oomph::DoubleMultiVector & mv, const Teuchos::Range1D& index )
147  {
148  return Teuchos::rcp(new oomph::DoubleMultiVector(mv,index,false));
149  }
150 
151  ///Obtain the global length of the vector
152  static int GetVecLength( const oomph::DoubleMultiVector& mv )
153  { return static_cast<int>(mv.nrow());}
154 
155  ///Obtain the number of vectors in the multivector
156  static int GetNumberVecs( const oomph::DoubleMultiVector& mv )
157  { return static_cast<int>(mv.nvector()); }
158 
159  /// \brief Update \c mv with \f$ \alpha AB + \beta mv \f$.
160  static void MvTimesMatAddMv(
161  const double alpha,
162  const oomph::DoubleMultiVector& A,
163  const Teuchos::SerialDenseMatrix<int,double>& B,
164  const double beta, oomph::DoubleMultiVector& mv )
165  {
166  //For safety let's (deep) copy A
168  //First pre-multiply by the scalars
169  mv *= beta;
170  C *= alpha;
171 
172  //Find the number of vectors in each multivector
173  int mv_n_vector = mv.nvector();
174  int A_n_vector = A.nvector();
175  //Find the number of rows
176  int n_row_local = A.nrow_local();
177  //Now need to multiply by the matrix and add the result
178  //to the appropriate entry in the multivector
179  for(int i=0;i<n_row_local;++i)
180  {
181  for(int j=0;j<mv_n_vector;j++)
182  {
183  double res = 0.0;
184  for(int k=0;k<A_n_vector;k++) {res += C(k,i)*B(k,j);}
185  mv(j,i) += res;
186  }
187  }
188  }
189 
190  /// \brief Replace \c mv with \f$\alpha A + \beta B\f$.
191  static void MvAddMv(const double alpha, const oomph::DoubleMultiVector& A,
192  const double beta, const oomph::DoubleMultiVector& B,
194  {
195  const unsigned n_vector = A.nvector();
196  const unsigned n_row_local = A.nrow_local();
197  for(unsigned v=0;v<n_vector;v++)
198  {
199  for(unsigned i=0;i<n_row_local;i++)
200  {
201  mv(v,i) = alpha*A(v,i) + beta*B(v,i);
202  }
203  }
204  }
205 
206  /*! \brief Scale each element of the vectors in \c mv with \c alpha.
207  */
208  static void MvScale ( oomph::DoubleMultiVector& mv, const double alpha )
209  {mv *= alpha;}
210 
211  /*! \brief Scale each element of the \c i-th vector in \c mv with \c alpha[i].
212  */
214  const std::vector<double>& alpha )
215  {
216  const unsigned n_vector = mv.nvector();
217  const unsigned n_row_local = mv.nrow_local();
218  for(unsigned v=0;v<n_vector;v++)
219  {
220  for(unsigned i=0;i<n_row_local;i++)
221  {
222  mv(v,i) *= alpha[v];
223  }
224  }
225  }
226 
227  /*! \brief Compute a dense matrix \c B through the matrix-matrix multiply \f$ \alpha A^Hmv \f$.
228  */
229  static void MvTransMv( const double alpha, const oomph::DoubleMultiVector& A,
230  const oomph::DoubleMultiVector& mv,
231  Teuchos::SerialDenseMatrix<int,double>& B)
232  {
233  const unsigned A_nvec = A.nvector();
234  const unsigned A_nrow_local = A.nrow_local();
235  const unsigned mv_nvec = mv.nvector();
236  //const unsigned mv_nrow_local = mv.nrow_local();
237 
238  for(unsigned v1=0;v1<A_nvec;v1++)
239  {
240  for(unsigned v2=0;v2<mv_nvec;v2++)
241  {
242  B(v1,v2) = 0.0;
243  for(unsigned i=0;i<A_nrow_local;i++)
244  {
245  B(v1,v2) += A(v1,i)*mv(v2,i);
246  }
247  B(v1,v2) *= alpha;
248  }
249  }
250 
251  //This will need to be communicated if we're in parallel and distributed
252 #ifdef OOMPH_HAS_MPI
253  if (A.distributed() &&
254  A.distribution_pt()->communicator_pt()->nproc() > 1)
255  {
256  const int n_total_val = A_nvec*mv_nvec;
257  //Pack entries into a vector
258  double b[n_total_val]; double b2[n_total_val];
259  unsigned count=0;
260  for(unsigned v1=0;v1<A_nvec;v1++)
261  {
262  for(unsigned v2=0;v2<mv_nvec;v2++)
263  {
264  b[count] = B(v1,v2);
265  b2[count] = b[count];
266  ++count;
267  }
268  }
269 
270 
271  MPI_Allreduce(&b,&b2,n_total_val,MPI_DOUBLE,MPI_SUM,
272  A.distribution_pt()->communicator_pt()->mpi_comm());
273 
274  count=0;
275  for(unsigned v1=0;v1<A_nvec;v1++)
276  {
277  for(unsigned v2=0;v2<mv_nvec;v2++)
278  {
279  B(v1,v2) = b2[count];
280  ++count;
281  }
282  }
283  }
284 #endif
285  }
286 
287 
288  /*! \brief Compute a vector \c b where the components are the individual dot-products of the \c i-th columns of \c A and \c mv, i.e.\f$b[i] = A[i]^Hmv[i]\f$.
289  */
290  static void MvDot ( const oomph::DoubleMultiVector& mv,
291  const oomph::DoubleMultiVector& A,
292  std::vector<double> &b)
293  {mv.dot(A,b);}
294 
295 
296  /*! \brief Compute the 2-norm of each individual vector of \c mv.
297  Upon return, \c normvec[i] holds the value of \f$||mv_i||_2\f$, the \c i-th column of \c mv.
298  */
299  static void MvNorm( const oomph::DoubleMultiVector& mv,
300  std::vector<double> &normvec )
301  { mv.norm(normvec); }
302 
303  //@}
304 
305  //! @name Initialization methods
306  //@{
307  /*! \brief Copy the vectors in \c A to a set of vectors in \c mv indicated by the indices given in \c index.
308 
309  The \c numvecs vectors in \c A are copied to a subset of vectors in \c mv indicated by the indices given in \c index,
310  i.e.<tt> mv[index[i]] = A[i]</tt>.
311  */
312  static void SetBlock( const oomph::DoubleMultiVector& A,
313  const std::vector<int>& index,
315  {
316  //Check some stuff
317  const unsigned n_index = index.size();
318  if(n_index==0) {return;}
319  const unsigned n_row_local = mv.nrow_local();
320  for(unsigned v=0;v<n_index;v++)
321  {
322  for(unsigned i=0;i<n_row_local;i++)
323  {
324  mv(index[v],i) = A(v,i);
325  }
326  }
327  }
328 
329  /// \brief Deep copy of A into specified columns of mv
330  ///
331  /// (Deeply) copy the first <tt>index.size()</tt> columns of \c A
332  /// into the columns of \c mv specified by the given index range.
333  ///
334  /// Postcondition: <tt>mv[i] = A[i - index.lbound()]</tt>
335  /// for all <tt>i</tt> in <tt>[index.lbound(), index.ubound()]</tt>
336  ///
337  /// \param A [in] Source multivector
338  /// \param index [in] Inclusive index range of columns of mv;
339  /// index set of the target
340  /// \param mv [out] Target multivector
341  static void SetBlock( const oomph::DoubleMultiVector& A,
342  const Teuchos::Range1D& index,
344  {
345  //Check some stuff
346  const unsigned n_index = index.size();
347  if(n_index==0) {return;}
348  const unsigned n_row_local = mv.nrow_local();
349  unsigned range_index = index.lbound();
350  for(unsigned v=0;v<n_index;v++)
351  {
352  for(unsigned i=0;i<n_row_local;i++)
353  {
354  mv(range_index,i) = A(v,i);
355  }
356  ++range_index;
357  }
358  }
359 
360  /// \brief mv := A
361  ///
362  /// Assign (deep copy) A into mv.
363  static void Assign( const oomph::DoubleMultiVector& A,
365  { mv = A; }
366 
367  /*! \brief Replace the vectors in \c mv with random vectors.
368  */
370  {
371  const unsigned n_vector = mv.nvector();
372  const unsigned n_row_local = mv.nrow_local();
373  for(unsigned v=0;v<n_vector;v++)
374  {
375  for(unsigned i=0;i<n_row_local;i++)
376  {
377  mv(v,i) = Teuchos::ScalarTraits<double>::random();
378  }
379  }
380  }
381 
382  /*! \brief Replace each element of the vectors in \c mv with \c alpha.
383  */
385  const double alpha =
386  Teuchos::ScalarTraits<double>::zero() )
387  { mv.initialise(alpha); }
388 
389  //@}
390 
391  //! @name Print method
392  //@{
393 
394  /*! \brief Print the \c mv multi-vector to the \c os output stream.
395  */
396  static void MvPrint( const oomph::DoubleMultiVector& mv, std::ostream& os )
397  { mv.output(os); }
398 
399  //@}
400  };
401 
402 
403 }
404 
405 namespace oomph
406 {
407 
408 //===================================================================
409 /// Base class for Oomph-lib's Vector Operator classes that will be
410 /// used with the DoubleMultiVector
411 //==================================================================
413 {
414 public:
415 
416  ///Empty constructor
418 
419  ///virtual destructor
421 
422  ///The apply interface
423  virtual void apply(const DoubleMultiVector &x,
425  const = 0;
426 };
427 
428 }
429 
430 namespace Anasazi
431 {
432 //======================================================================
433 /// Anasazi Traits class that specialises the apply operator based on
434 /// oomph-lib's DoubleVector class and double as the primitive type
435 //======================================================================
436 template <>
437  class OperatorTraits<double,oomph::DoubleMultiVector,
439  {
440  public:
441 
442  /// Matrix operator application method
443  static void Apply ( const oomph::DoubleMultiVectorOperator& Op,
444  const oomph::DoubleMultiVector &x,
446  {
447  Op.apply(x,y);
448  }
449 
450  }; // End of the specialised traits class
451 
452 }
453 
454 
455 namespace oomph
456 {
457 
458 //====================================================================
459 /// Class for the shift invert operation
460 //====================================================================
462 {
463 private:
464  //Pointer to the problem
466 
467  //Pointer to the linear solver used
469 
470  //Storage for the shift
471  double Sigma;
472 
473  //Storage for the matrices
475 
476 public:
477 
479  LinearSolver* const &linear_solver_pt,
480  const double &sigma=0.0) :
481  Problem_pt(problem_pt), Linear_solver_pt(linear_solver_pt), Sigma(sigma)
482  {
483  //Before we get into the Arnoldi loop solve the shifted matrix problem
484  //Allocated Row compressed matrices for the mass matrix and shifted main
485  //matrix
486  M_pt = new CRDoubleMatrix(problem_pt->dof_distribution_pt());
487  AsigmaM_pt = new CRDoubleMatrix(problem_pt->dof_distribution_pt());
488 
489  //Assemble the matrices
490  problem_pt->get_eigenproblem_matrices(*M_pt,*AsigmaM_pt,Sigma);
491 
492  //Do not report the time taken
494  }
495 
496 
497  //Now specify how to apply the operator
498  void apply (const DoubleMultiVector &x,
499  DoubleMultiVector &y) const
500  {
501  const unsigned n_vec = x.nvector();
502  const unsigned n_row_local = x.nrow_local();
503  if(n_vec > 1) {Linear_solver_pt->enable_resolve();}
504  //Solve the first vector
505 
507  //Premultiply by mass matrix
508  M_pt->multiply(x.doublevector(0),X);
509  //For some reason I need to create a new vector and copy here
510  //This is odd and not expected, examine carefully
513  //Need to synchronise
514 //#ifdef OOMPH_HAS_MPI
515 // Problem_pt->synchronise_all_dofs();
516 //#endif
517  for(unsigned i=0;i<n_row_local;i++) {y(0,i) = Y[i];}
518 
519  //Now loop over the others and resolve
520  for(unsigned v=1;v<n_vec;++v)
521  {
522  M_pt->multiply(x.doublevector(v),X);
524 //#ifdef OOMPH_HAS_MPI
525 // Problem_pt->synchronise_all_dofs();
526 //#endif
527  for(unsigned i=0;i<n_row_local;i++) {y(v,i) = Y[i];}
528  }
529  }
530 
531 
532 };
533 
534 
535 //=====================================================================
536 /// Class for the Anasazi eigensolver
537 //=====================================================================
538 class ANASAZI : public EigenSolver
539 {
540  private:
541 
542  typedef double ST;
543  typedef Teuchos::ScalarTraits<ST> SCT;
544  typedef SCT::magnitudeType MT;
545  typedef Anasazi::MultiVec<ST> MV;
546  typedef Anasazi::Operator<ST> OP;
547  typedef Anasazi::MultiVecTraits<ST,MV> MVT;
548  typedef Anasazi::OperatorTraits<ST,MV,OP> OPT;
549 
550  /// \short Pointer to output manager
551  Anasazi::OutputManager<ST> *Output_manager_pt;
552 
553  /// \short Pointer to a linear solver
555 
556  /// \short Pointer to a default linear solver
558 
559  /// \short Integer to set whether the real, imaginary or magnitude is required
560  ///to be small or large.
561  int Spectrum;
562 
563  /// \short Number of Arnoldi vectors to compute
564  int NArnoldi;
565 
566  /// \short Set the shifted value
567  double Sigma;
568 
569  ///\short Boolean to set which part of the spectrum left (default) or right
570  /// of the shifted value.
571  bool Small;
572 
573  /// \short Boolean to indicate whether or not to compute the eigenvectors
575 
576 
577  public:
578 
579  ///Constructor
581  NArnoldi(10), Sigma(0.0), Compute_eigenvectors(true)
582  {
583  Output_manager_pt = new Anasazi::BasicOutputManager<ST>();
584  // Set verbosity level
585  int verbosity =
586  Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
587  Output_manager_pt->setVerbosity(verbosity);
588 
589  // print greeting
590  Output_manager_pt->stream(Anasazi::Warnings)
591  << Anasazi::Anasazi_Version() << std::endl << std::endl;
592 
593  }
594 
595  ///Empty copy constructor
596  ANASAZI(const ANASAZI&): Sigma(0.0) {}
597 
598  ///Destructor, delete the linear solver
599  virtual ~ANASAZI()
600  {
601 
602  }
603 
604  /// Access function for the number of Arnoldi vectors
605  int &narnoldi() {return NArnoldi;}
606 
607  /// Access function for the number of Arnoldi vectors (const version)
608  const int &narnoldi() const {return NArnoldi;}
609 
610  /// \short Set to enable the computation of the eigenvectors (default)
612 
613  /// \short Set to disable the computation of the eigenvectors
615 
616  /// Solve the eigen problem
617  void solve_eigenproblem(Problem* const &problem_pt,
618  const int &n_eval,
619  Vector<std::complex<double> > &eigenvalue,
620  Vector<DoubleVector> &eigenvector)
621  {
622  //No access to sigma, so set from sigma real
623  Sigma = Sigma_real;
624 
625  //Initially be dumb here
626  Linear_solver_pt = problem_pt->linear_solver_pt();
627 
628  //Let's make the initial vector
629  Teuchos::RCP<DoubleMultiVector> initial =
630  Teuchos::rcp( new DoubleMultiVector(1,problem_pt->dof_distribution_pt()));
631  Anasazi::MultiVecTraits<double,DoubleMultiVector>::MvRandom(*initial);
632 
633  //Make the operator
634  Teuchos::RCP<DoubleMultiVectorOperator> Op_pt =
635  Teuchos::rcp(
636  new ProblemBasedShiftInvertOperator(problem_pt,
637  this->linear_solver_pt(),Sigma));
638 
639  //Create the basic problem
640  Teuchos::RCP<
641  Anasazi::BasicEigenproblem<double,DoubleMultiVector,
643  anasazi_pt =
644  Teuchos::rcp(
645  new Anasazi::BasicEigenproblem<double,DoubleMultiVector,
647  (Op_pt,initial));
648 
649  //Think I have it?
650  anasazi_pt->setHermitian(false);
651 
652  // set the number of eigenvalues requested
653  anasazi_pt->setNEV(n_eval);
654 
655  // Inform the eigenproblem that you are done passing it information
656  if (anasazi_pt->setProblem() != true) {
657  oomph_info
658  << "Anasazi::BasicEigenproblem::setProblem() had an error." << std::endl
659  << "End Result: TEST FAILED" << std::endl;
660  }
661 
662  // Create the solver manager
663  // No need to have ncv specificed, Triliinos has a sensible default
664 // int ncv = 10;
665  MT tol = 1.0e-10;
666  int verbosity =
667  Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
668 
669  Teuchos::ParameterList MyPL;
670  MyPL.set( "Which", "LM" );
671  MyPL.set( "Block Size", 1 );
672 // MyPL.set( "Num Blocks", ncv );
673  MyPL.set( "Maximum Restarts", 500 );
674  MyPL.set( "Orthogonalization", "DGKS" );
675  MyPL.set( "Verbosity", verbosity );
676  MyPL.set( "Convergence Tolerance", tol );
677  Anasazi::BlockKrylovSchurSolMgr<
678  double,DoubleMultiVector,
680  BKS(anasazi_pt, MyPL);
681 
682  // Solve the problem to the specified tolerances or length
683  Anasazi::ReturnType ret = BKS.solve();
684  bool testFailed = false;
685  if (ret != Anasazi::Converged) {
686  testFailed = true;
687  }
688 
689  if(testFailed) {oomph_info << "Eigensolver not converged\n";}
690 
691  // Get the eigenvalues and eigenvectors from the eigenproblem
692  Anasazi::Eigensolution<double,DoubleMultiVector> sol =
693  anasazi_pt->getSolution();
694  std::vector<Anasazi::Value<double> > evals = sol.Evals;
695  Teuchos::RCP<DoubleMultiVector> evecs = sol.Evecs;
696 
697  eigenvalue.resize(evals.size());
698  eigenvector.resize(evals.size());
699 
700  for(unsigned i=0;i<evals.size();i++)
701  {
702  //Undo shift and invert
703  double a = evals[i].realpart; double b = evals[i].imagpart;
704  double det = a*a + b*b;
705  eigenvalue[i] = std::complex<double>(a/det+Sigma,-b/det);
706 
707  //Now set the eigenvectors, I hope
708  eigenvector[i].build(evecs->distribution_pt());
709  unsigned nrow_local= evecs->nrow_local();
710  //Would be faster with pointers, but I'll sort that out later!
711  for(unsigned n=0;n<nrow_local;n++)
712  {
713  eigenvector[i][n] = (*evecs)(i,n);
714  }
715  }
716  }
717 
718  /// Set the desired eigenvalues to be left of the shift
720 
721  /// Set the desired eigenvalues to be right of the shift
723 
724  /// Set the real part to be the quantity of interest (default)
726 
727  /// Set the imaginary part fo the quantity of interest
729 
730  /// Set the magnitude to be the quantity of interest
732 
733  /// Return a pointer to the linear solver object
735 
736  /// Return a pointer to the linear solver object (const version)
738 
739 };
740 
741 }
742 
743 
744 #endif
static int GetVecLength(const oomph::DoubleMultiVector &mv)
Obtain the global length of the vector.
virtual ~DoubleMultiVectorOperator()
virtual destructor
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...
unsigned nrow_local() const
access function for the num of local rows on this processor.
static void MvAddMv(const double alpha, const oomph::DoubleMultiVector &A, const double beta, const oomph::DoubleMultiVector &B, oomph::DoubleMultiVector &mv)
Replace mv with .
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 ...
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
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:...
void get_eigenvalues_left_of_shift()
Set the desired eigenvalues to be left of the shift.
bool Small
Boolean to set which part of the spectrum left (default) or right of the shifted value.
static void MvScale(oomph::DoubleMultiVector &mv, const double alpha)
Scale each element of the vectors in mv with alpha.
virtual void apply(const DoubleMultiVector &x, DoubleMultiVector &y) const =0
The apply interface.
ANASAZI(const ANASAZI &)
Empty copy constructor.
void enable_compute_eigenvectors()
Set to enable the computation of the eigenvectors (default)
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...
Definition: problem.cc:8337
const int & narnoldi() const
Access function for the number of Arnoldi vectors (const version)
static void MvPrint(const oomph::DoubleMultiVector &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
int & narnoldi()
Access function for the number of Arnoldi vectors.
void disable_compute_eigenvectors()
Set to disable the computation of the eigenvectors.
Anasazi::OperatorTraits< ST, MV, OP > OPT
cstr elem_len * i
Definition: cfortran.h:607
unsigned nrow() const
access function to the number of global rows.
bool Compute_eigenvectors
Boolean to indicate whether or not to compute the eigenvectors.
int Spectrum
Integer to set whether the real, imaginary or magnitude is required to be small or large...
The Problem class.
Definition: problem.h:152
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...
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 void SetBlock(const oomph::DoubleMultiVector &A, const Teuchos::Range1D &index, oomph::DoubleMultiVector &mv)
Deep copy of A into specified columns of mv.
static void MvRandom(oomph::DoubleMultiVector &mv)
Replace the vectors in mv with random vectors.
void dot(const DoubleMultiVector &vec, std::vector< double > &result) const
compute the 2 norm of this vector
LinearSolver * Default_linear_solver_pt
Pointer to a default linear solver.
bool distributed() const
distribution is serial or distributed
OomphInfo oomph_info
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:...
Class for the shift invert operation.
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:...
DoubleMultiVectorOperator()
Empty constructor.
static void Apply(const oomph::DoubleMultiVectorOperator &Op, const oomph::DoubleMultiVector &x, oomph::DoubleMultiVector &y)
Matrix operator application method.
void get_eigenvalues_right_of_shift()
Set the desired eigenvalues to be right of the shift.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void output(std::ostream &outfile) const
output the contents of the vector
Class for the Anasazi eigensolver.
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:...
double Sigma_real
Double value that represents the real part of the shift in shifted eigensolvers.
Definition: eigen_solver.h:71
double Sigma
Set the shifted value.
Anasazi::MultiVec< ST > MV
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...
int NArnoldi
Number of Arnoldi vectors to compute.
DoubleVector & doublevector(const unsigned &i)
access to the DoubleVector representatoin
void initialise(const double &initial_value)
initialise the whole vector with value v
static void MvScale(oomph::DoubleMultiVector &mv, const std::vector< double > &alpha)
Scale each element of the i-th vector in mv with alpha[i].
virtual ~ANASAZI()
Destructor, delete the linear solver.
LinearAlgebraDistribution *const & dof_distribution_pt() const
Return the pointer to the dof distribution (read-only)
Definition: problem.h:1570
Anasazi::MultiVecTraits< ST, MV > MVT
static int GetNumberVecs(const oomph::DoubleMultiVector &mv)
Obtain the number of vectors in the multivector.
LinearSolver * Linear_solver_pt
Pointer to a linear solver.
Anasazi::Operator< ST > OP
static void Assign(const oomph::DoubleMultiVector &A, oomph::DoubleMultiVector &mv)
mv := A
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:...
ProblemBasedShiftInvertOperator(Problem *const &problem_pt, LinearSolver *const &linear_solver_pt, const double &sigma=0.0)
Teuchos::ScalarTraits< ST > SCT
void track_eigenvalue_real_part()
Set the real part to be the quantity of interest (default)
void disable_doc_time()
Disable documentation of solve times.
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 .
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
ANASAZI()
Constructor.
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 .
void norm(std::vector< double > &result) const
compute the 2 norm of this vector
SCT::magnitudeType MT
void track_eigenvalue_imaginary_part()
Set the imaginary part fo the quantity of interest.
A multi vector in the mathematical sense, initially developed for linear algebra type applications...
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1424
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...
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1805
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...
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
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 ...
void track_eigenvalue_magnitude()
Set the magnitude to be the quantity of interest.
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 ...
void apply(const DoubleMultiVector &x, DoubleMultiVector &y) const
The apply interface.
static void MvInit(oomph::DoubleMultiVector &mv, const double alpha=Teuchos::ScalarTraits< double >::zero())
Replace each element of the vectors in mv with alpha.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
LinearSolver *const & linear_solver_pt() const
Return a pointer to the linear solver object (const version)
unsigned nvector() const
Return the number of vectors.
Anasazi::OutputManager< ST > * Output_manager_pt
Pointer to output manager.