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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
27#define OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
28
29// oomph-lib includes
30#include "double_multi_vector.h"
31#include "problem.h"
32#include "linear_solver.h"
33#include "eigen_solver.h"
34
35// Trilinos includes
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"
43
44namespace Anasazi
45{
46 //========================================================================
47 /// Specialize the Anasazi traits class for the oomph-lib DoubleMultiVector.
48 /// This provides the interfaces required by the Anasazi eigensolvers.
49 //========================================================================
50 template<>
51 class MultiVecTraits<double, oomph::DoubleMultiVector>
52 {
53 public:
54 /// \brief Creates a new empty DoubleMultiVector containing numvecs columns.
55 /// Return a reference-counted pointer to the new multivector
56 static Teuchos::RCP<oomph::DoubleMultiVector> Clone(
57 const oomph::DoubleMultiVector& mv, const int numvecs)
58 {
59 return Teuchos::rcp(new oomph::DoubleMultiVector(numvecs, mv));
60 }
61
62 /// \brief Creates a deep copy of the DoubleMultiVector mv
63 /// return Reference-counted pointer to the new DoubleMultiVector
64 static Teuchos::RCP<oomph::DoubleMultiVector> CloneCopy(
66 {
67 return Teuchos::rcp(new oomph::DoubleMultiVector(mv));
68 }
69
70 /// \brief Creates a new oomph::DoubleMultiVector and (deep)
71 /// copies the contents of
72 /// the vectors in index into the new vector
73 /// return Reference-counted pointer to the new oomph::DoubleMultiVector
74 static Teuchos::RCP<oomph::DoubleMultiVector> CloneCopy(
75 const oomph::DoubleMultiVector& mv, const std::vector<int>& index)
76 {
77 return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index));
78 }
79
80 /// \brief Deep copy of specified columns of oomph::DoubleMultiVector
81 /// return Reference-counted pointer to the new oomph::DoubleMultiVector
82 static Teuchos::RCP<oomph::DoubleMultiVector> CloneCopy(
83 const oomph::DoubleMultiVector& mv, const Teuchos::Range1D& index)
84 {
85 return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index));
86 }
87
88 /// \brief Creates a new oomph::DoubleMultiVector that contains
89 /// shallow copies
90 /// of selected entries of the oomph::DoubleMultiVector mv
91 /// return Reference-counted pointer to the new oomph::DoubleMultiVector
92 static Teuchos::RCP<oomph::DoubleMultiVector> CloneViewNonConst(
93 oomph::DoubleMultiVector& mv, const std::vector<int>& index)
94 {
95 return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
96 }
97
98 /// \brief Creates a new oomph::DoubleMultiVector that
99 /// contains shallow copies
100 /// of selected entries of the oomph::DoubleMultiVector mv
101 /// return Reference-counted pointer to the new oomph::DoubleMultiVector
102 static Teuchos::RCP<oomph::DoubleMultiVector> CloneViewNonConst(
103 oomph::DoubleMultiVector& mv, const Teuchos::Range1D& index)
104 {
105 return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
106 }
107
108 /// \brief Creates a new oomph::DoubleMultiVector that contains
109 /// shallow copies
110 /// of selected entries of the oomph::DoubleMultiVector mv
111 /// return Reference-counted pointer to the new oomph::DoubleMultiVector
112 /// (const version)
113 static Teuchos::RCP<const oomph::DoubleMultiVector> CloneView(
114 const oomph::DoubleMultiVector& mv, const std::vector<int>& index)
115 {
116 return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
117 }
118
119 /// \brief Creates a new oomph::DoubleMultiVector that contains
120 /// shallow copies
121 /// of selected entries of the oomph::DoubleMultiVector mv
122 /// return Reference-counted pointer to the new oomph::DoubleMultiVector
123 /// (Non-const version for Trilinos 9 interface)
124 static Teuchos::RCP<oomph::DoubleMultiVector> CloneView(
125 oomph::DoubleMultiVector& mv, const std::vector<int>& index)
126 {
127 return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
128 }
129
130 /// \brief Creates a new oomph::DoubleMultiVector that
131 /// contains shallow copies
132 /// of selected entries of the oomph::DoubleMultiVector mv
133 /// return Reference-counted pointer to the new oomph::DoubleMultiVector
134 /// (const version)
135 static Teuchos::RCP<oomph::DoubleMultiVector> CloneView(
136 oomph::DoubleMultiVector& mv, const Teuchos::Range1D& index)
137 {
138 return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
139 }
140
141 /// Obtain the global length of the vector
143 {
144 return static_cast<int>(mv.nrow());
145 }
146
147 /// Obtain the number of vectors in the multivector
149 {
150 return static_cast<int>(mv.nvector());
151 }
152
153 /// \brief Update \c mv with \f$ \alpha AB + \beta mv \f$.
154 static void MvTimesMatAddMv(
155 const double alpha,
157 const Teuchos::SerialDenseMatrix<int, double>& B,
158 const double beta,
160 {
161 // For safety let's (deep) copy A
163 // First pre-multiply by the scalars
164 mv *= beta;
165 C *= alpha;
166
167 // Find the number of vectors in each multivector
168 int mv_n_vector = mv.nvector();
169 int A_n_vector = A.nvector();
170 // Find the number of rows
171 int n_row_local = A.nrow_local();
172 // Now need to multiply by the matrix and add the result
173 // to the appropriate entry in the multivector
174 for (int i = 0; i < n_row_local; ++i)
175 {
176 for (int j = 0; j < mv_n_vector; j++)
177 {
178 double res = 0.0;
179 for (int k = 0; k < A_n_vector; k++)
180 {
181 res += C(k, i) * B(k, j);
182 }
183 mv(j, i) += res;
184 }
185 }
186 }
187
188 /// \brief Replace \c mv with \f$\alpha A + \beta B\f$.
189 static void MvAddMv(const double alpha,
191 const double beta,
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 {
210 mv *= alpha;
211 }
212
213 /*! \brief Scale each element of the \c i-th vector in \c mv with \c
214 * alpha[i].
215 */
217 const std::vector<double>& alpha)
218 {
219 const unsigned n_vector = mv.nvector();
220 const unsigned n_row_local = mv.nrow_local();
221 for (unsigned v = 0; v < n_vector; v++)
222 {
223 for (unsigned i = 0; i < n_row_local; i++)
224 {
225 mv(v, i) *= alpha[v];
226 }
227 }
228 }
229
230 /*! \brief Compute a dense matrix \c B through the matrix-matrix multiply
231 * \f$ \alpha A^Hmv \f$.
232 */
233 static void MvTransMv(const double alpha,
235 const oomph::DoubleMultiVector& mv,
236 Teuchos::SerialDenseMatrix<int, double>& B)
237 {
238 const unsigned A_nvec = A.nvector();
239 const unsigned A_nrow_local = A.nrow_local();
240 const unsigned mv_nvec = mv.nvector();
241 // const unsigned mv_nrow_local = mv.nrow_local();
242
243 for (unsigned v1 = 0; v1 < A_nvec; v1++)
244 {
245 for (unsigned v2 = 0; v2 < mv_nvec; v2++)
246 {
247 B(v1, v2) = 0.0;
248 for (unsigned i = 0; i < A_nrow_local; i++)
249 {
250 B(v1, v2) += A(v1, i) * mv(v2, i);
251 }
252 B(v1, v2) *= alpha;
253 }
254 }
255
256 // This will need to be communicated if we're in parallel and distributed
257#ifdef OOMPH_HAS_MPI
258 if (A.distributed() &&
259 A.distribution_pt()->communicator_pt()->nproc() > 1)
260 {
261 const int n_total_val = A_nvec * mv_nvec;
262 // Pack entries into a vector
263 double b[n_total_val];
264 double b2[n_total_val];
265 unsigned count = 0;
266 for (unsigned v1 = 0; v1 < A_nvec; v1++)
267 {
268 for (unsigned v2 = 0; v2 < mv_nvec; v2++)
269 {
270 b[count] = B(v1, v2);
271 b2[count] = b[count];
272 ++count;
273 }
274 }
275
276
277 MPI_Allreduce(&b,
278 &b2,
279 n_total_val,
280 MPI_DOUBLE,
281 MPI_SUM,
282 A.distribution_pt()->communicator_pt()->mpi_comm());
283
284 count = 0;
285 for (unsigned v1 = 0; v1 < A_nvec; v1++)
286 {
287 for (unsigned v2 = 0; v2 < mv_nvec; v2++)
288 {
289 B(v1, v2) = b2[count];
290 ++count;
291 }
292 }
293 }
294#endif
295 }
296
297
298 /*! \brief Compute a vector \c b where the components are the individual
299 * dot-products of the \c i-th columns of \c A and \c mv, i.e.
300 * \f$ b[i] = * A[i]^Hmv[i] \f$.
301 */
302 static void MvDot(const oomph::DoubleMultiVector& mv,
304 std::vector<double>& b)
305 {
306 mv.dot(A, b);
307 }
308
309
310 /*! \brief Compute the 2-norm of each individual vector of \c mv.
311 Upon return, \c normvec[i] holds the value of \f$||mv_i||_2\f$, the \c
312 i-th column of \c mv.
313 */
314 static void MvNorm(const oomph::DoubleMultiVector& mv,
315 std::vector<double>& normvec)
316 {
317 mv.norm(normvec);
318 }
319
320 //@}
321
322 //! @name Initialization methods
323 //@{
324 /*! \brief Copy the vectors in \c A to a set of vectors in \c mv indicated
325 by the indices given in \c index.
326
327 The \c numvecs vectors in \c A are copied to a subset of vectors in \c mv
328 indicated by the indices given in \c index, i.e.<tt> mv[index[i]] =
329 A[i]</tt>.
330 */
331 static void SetBlock(const oomph::DoubleMultiVector& A,
332 const std::vector<int>& index,
334 {
335 // Check some stuff
336 const unsigned n_index = index.size();
337 if (n_index == 0)
338 {
339 return;
340 }
341 const unsigned n_row_local = mv.nrow_local();
342 for (unsigned v = 0; v < n_index; v++)
343 {
344 for (unsigned i = 0; i < n_row_local; i++)
345 {
346 mv(index[v], i) = A(v, i);
347 }
348 }
349 }
350
351 /// \brief Deep copy of A into specified columns of mv
352 ///
353 /// (Deeply) copy the first <tt>index.size()</tt> columns of \c A
354 /// into the columns of \c mv specified by the given index range.
355 ///
356 /// Postcondition: <tt>mv[i] = A[i - index.lbound()]</tt>
357 /// for all <tt>i</tt> in <tt>[index.lbound(), index.ubound()]</tt>
358 ///
359 /// \param A [in] Source multivector
360 /// \param index [in] Inclusive index range of columns of mv;
361 /// index set of the target
362 /// \param mv [out] Target multivector
363 static void SetBlock(const oomph::DoubleMultiVector& A,
364 const Teuchos::Range1D& index,
366 {
367 // Check some stuff
368 const unsigned n_index = index.size();
369 if (n_index == 0)
370 {
371 return;
372 }
373 const unsigned n_row_local = mv.nrow_local();
374 unsigned range_index = index.lbound();
375 for (unsigned v = 0; v < n_index; v++)
376 {
377 for (unsigned i = 0; i < n_row_local; i++)
378 {
379 mv(range_index, i) = A(v, i);
380 }
381 ++range_index;
382 }
383 }
384
385 /// \brief mv := A
386 ///
387 /// Assign (deep copy) A into mv.
388 static void Assign(const oomph::DoubleMultiVector& A,
390 {
391 mv = A;
392 }
393
394 /*! \brief Replace the vectors in \c mv with random vectors.
395 */
397 {
398 const unsigned n_vector = mv.nvector();
399 const unsigned n_row_local = mv.nrow_local();
400 for (unsigned v = 0; v < n_vector; v++)
401 {
402 for (unsigned i = 0; i < n_row_local; i++)
403 {
404 mv(v, i) = Teuchos::ScalarTraits<double>::random();
405 }
406 }
407 }
408
409 /*! \brief Replace each element of the vectors in \c mv with \c alpha.
410 */
411 static void MvInit(
413 const double alpha = Teuchos::ScalarTraits<double>::zero())
414 {
415 mv.initialise(alpha);
416 }
417
418 //@}
419
420 //! @name Print method
421 //@{
422
423 /*! \brief Print the \c mv multi-vector to the \c os output stream.
424 */
425 static void MvPrint(const oomph::DoubleMultiVector& mv, std::ostream& os)
426 {
427 mv.output(os);
428 }
429
430 //@}
431 };
432
433
434} // namespace Anasazi
435
436namespace oomph
437{
438 //===================================================================
439 /// Base class for Oomph-lib's Vector Operator classes that will be
440 /// used with the DoubleMultiVector
441 //==================================================================
443 {
444 public:
445 /// Empty constructor
447
448 /// virtual destructor
450
451 /// The apply interface
452 virtual void apply(const DoubleMultiVector& x,
453 DoubleMultiVector& y) const = 0;
454 };
455
456} // namespace oomph
457
458namespace Anasazi
459{
460 //======================================================================
461 /// Anasazi Traits class that specialises the apply operator based on
462 /// oomph-lib's DoubleVector class and double as the primitive type
463 //======================================================================
464 template<>
465 class OperatorTraits<double,
466 oomph::DoubleMultiVector,
468 {
469 public:
470 /// Matrix operator application method
474 {
475 Op.apply(x, y);
476 }
477
478 }; // End of the specialised traits class
479
480} // namespace Anasazi
481
482
483namespace oomph
484{
485 //====================================================================
486 /// Class for the shift invert operation
487 //====================================================================
489 {
490 private:
491 // Pointer to the problem
493
494 // Pointer to the linear solver used
496
497 // Storage for the shift
498 double Sigma;
499
500 // Storage for the matrices
502
503 public:
505 LinearSolver* const& linear_solver_pt,
506 const double& sigma = 0.0)
507 : Problem_pt(problem_pt), Linear_solver_pt(linear_solver_pt), Sigma(sigma)
508 {
509 // Before we get into the Arnoldi loop solve the shifted matrix problem
510 // Allocated Row compressed matrices for the mass matrix and shifted main
511 // matrix
512 M_pt = new CRDoubleMatrix(problem_pt->dof_distribution_pt());
513 AsigmaM_pt = new CRDoubleMatrix(problem_pt->dof_distribution_pt());
514
515 // Assemble the matrices
517
518 // Do not report the time taken
520 }
521
522
523 // Now specify how to apply the operator
525 {
526 const unsigned n_vec = x.nvector();
527 const unsigned n_row_local = x.nrow_local();
528 if (n_vec > 1)
529 {
531 }
532 // Solve the first vector
533
535 // Premultiply by mass matrix
536 M_pt->multiply(x.doublevector(0), X);
537 // For some reason I need to create a new vector and copy here
538 // This is odd and not expected, examine carefully
541 // Need to synchronise
542 //#ifdef OOMPH_HAS_MPI
543 // Problem_pt->synchronise_all_dofs();
544 //#endif
545 for (unsigned i = 0; i < n_row_local; i++)
546 {
547 y(0, i) = Y[i];
548 }
549
550 // Now loop over the others and resolve
551 for (unsigned v = 1; v < n_vec; ++v)
552 {
553 M_pt->multiply(x.doublevector(v), X);
555 //#ifdef OOMPH_HAS_MPI
556 // Problem_pt->synchronise_all_dofs();
557 //#endif
558 for (unsigned i = 0; i < n_row_local; i++)
559 {
560 y(v, i) = Y[i];
561 }
562 }
563 }
564 };
565
566
567 //=====================================================================
568 /// Class for the Anasazi eigensolver
569 //=====================================================================
570 class ANASAZI : public EigenSolver
571 {
572 private:
573 typedef double ST;
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;
580
581 /// Pointer to output manager
582 Anasazi::OutputManager<ST>* Output_manager_pt;
583
584 /// Pointer to a linear solver
586
587 /// Pointer to a default linear solver
589
590 /// Integer to set whether the real, imaginary or magnitude is
591 /// required
592 /// to be small or large.
594
595 /// Number of Arnoldi vectors to compute
597
598 /// Set the shifted value
599 double Sigma;
600
601 /// Boolean to set which part of the spectrum left (default) or right
602 /// of the shifted value.
603 bool Small;
604
605 /// Boolean to indicate whether or not to compute the eigenvectors
607
608
609 public:
610 /// Constructor
612 : Linear_solver_pt(0),
614 Spectrum(0),
615 NArnoldi(10),
616 Sigma(0.0),
618 {
619 Output_manager_pt = new Anasazi::BasicOutputManager<ST>();
620 // Set verbosity level
621 int verbosity =
622 Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
623 Output_manager_pt->setVerbosity(verbosity);
624
625 // print greeting
626 Output_manager_pt->stream(Anasazi::Warnings)
627 << Anasazi::Anasazi_Version() << std::endl
628 << std::endl;
629 }
630
631 /// Empty copy constructor
632 ANASAZI(const ANASAZI&) : Sigma(0.0) {}
633
634 /// Destructor, delete the linear solver
635 virtual ~ANASAZI() {}
636
637 /// Access function for the number of Arnoldi vectors
638 int& narnoldi()
639 {
640 return NArnoldi;
641 }
642
643 /// Access function for the number of Arnoldi vectors (const version)
644 const int& narnoldi() const
645 {
646 return NArnoldi;
647 }
648
649 /// Set to enable the computation of the eigenvectors (default)
651 {
653 }
654
655 /// Set to disable the computation of the eigenvectors
657 {
658 Compute_eigenvectors = false;
659 }
660
661 /// Solve the eigen problem
662 void solve_eigenproblem(Problem* const& problem_pt,
663 const int& n_eval,
664 Vector<std::complex<double>>& eigenvalue,
665 Vector<DoubleVector>& eigenvector)
666 {
667 // No access to sigma, so set from sigma real
669
670 // Initially be dumb here
671 Linear_solver_pt = problem_pt->linear_solver_pt();
672
673 // Let's make the initial vector
674 Teuchos::RCP<DoubleMultiVector> initial = Teuchos::rcp(
675 new DoubleMultiVector(1, problem_pt->dof_distribution_pt()));
676 Anasazi::MultiVecTraits<double, DoubleMultiVector>::MvRandom(*initial);
677
678 // Make the operator
679 Teuchos::RCP<DoubleMultiVectorOperator> Op_pt =
680 Teuchos::rcp(new ProblemBasedShiftInvertOperator(
681 problem_pt, this->linear_solver_pt(), Sigma));
682
683 // Create the basic problem
684 Teuchos::RCP<Anasazi::BasicEigenproblem<double,
687 anasazi_pt = Teuchos::rcp(
688 new Anasazi::BasicEigenproblem<double,
691 initial));
692
693 // Think I have it?
694 anasazi_pt->setHermitian(false);
695
696 // set the number of eigenvalues requested
697 anasazi_pt->setNEV(n_eval);
698
699 // Inform the eigenproblem that you are done passing it information
700 if (anasazi_pt->setProblem() != true)
701 {
702 oomph_info << "Anasazi::BasicEigenproblem::setProblem() had an error."
703 << std::endl
704 << "End Result: TEST FAILED" << std::endl;
705 }
706
707 // Create the solver manager
708 // No need to have ncv specificed, Triliinos has a sensible default
709 // int ncv = 10;
710 MT tol = 1.0e-10;
711 int verbosity =
712 Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
713
714 Teuchos::ParameterList MyPL;
715 MyPL.set("Which", "LM");
716 MyPL.set("Block Size", 1);
717 // MyPL.set( "Num Blocks", ncv );
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);
726
727 // Solve the problem to the specified tolerances or length
728 Anasazi::ReturnType ret = BKS.solve();
729 bool testFailed = false;
730 if (ret != Anasazi::Converged)
731 {
732 testFailed = true;
733 }
734
735 if (testFailed)
736 {
737 oomph_info << "Eigensolver not converged\n";
738 }
739
740 // Get the eigenvalues and eigenvectors from the eigenproblem
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;
745
746 eigenvalue.resize(evals.size());
747 eigenvector.resize(evals.size());
748
749 for (unsigned i = 0; i < evals.size(); i++)
750 {
751 // Undo shift and invert
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);
756
757 // Now set the eigenvectors, I hope
758 eigenvector[i].build(evecs->distribution_pt());
759 unsigned nrow_local = evecs->nrow_local();
760 // Would be faster with pointers, but I'll sort that out later!
761 for (unsigned n = 0; n < nrow_local; n++)
762 {
763 eigenvector[i][n] = (*evecs)(i, n);
764 }
765 }
766 }
767
768 /// Set the desired eigenvalues to be left of the shift
770 {
771 Small = true;
772 }
773
774 /// Set the desired eigenvalues to be right of the shift
776 {
777 Small = false;
778 }
779
780 /// Set the real part to be the quantity of interest (default)
782 {
783 Spectrum = 1;
784 }
785
786 /// Set the imaginary part fo the quantity of interest
788 {
789 Spectrum = -1;
790 }
791
792 /// Set the magnitude to be the quantity of interest
794 {
795 Spectrum = 0;
796 }
797
798 /// Return a pointer to the linear solver object
800 {
801 return Linear_solver_pt;
802 }
803
804 /// Return a pointer to the linear solver object (const version)
806 {
807 return Linear_solver_pt;
808 }
809 };
810
811} // namespace oomph
812
813
814#endif
cstr elem_len * i
Definition: cfortran.h:603
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)
SCT::magnitudeType MT
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.
Definition: matrices.h:888
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1782
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.
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....
Definition: double_vector.h:58
Base class for all EigenProblem solves. This simply defines standard interfaces so that different sol...
Definition: eigen_solver.h:61
double Sigma_real
Double value that represents the real part of the shift in shifted eigensolvers.
Definition: eigen_solver.h:65
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
Base class for all linear solvers. This merely defines standard interfaces for linear solvers,...
Definition: linear_solver.h:68
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)
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
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:8368
LinearAlgebraDistribution *const & dof_distribution_pt() const
Return the pointer to the dof distribution (read-only)
Definition: problem.h:1668
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1466
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...