Anasazi  Version of the Day
AnasaziEpetraAdapter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
33 #ifndef ANASAZI_EPETRA_ADAPTER_HPP
34 #define ANASAZI_EPETRA_ADAPTER_HPP
35 
36 #include "AnasaziConfigDefs.hpp"
37 #include "Anasaziepetra_DLLExportMacro.h"
38 #include "AnasaziTypes.hpp"
39 #include "AnasaziMultiVec.hpp"
40 #include "AnasaziOperator.hpp"
41 
42 #include "Teuchos_Assert.hpp"
43 #include "Teuchos_SerialDenseMatrix.hpp"
44 #include "Epetra_MultiVector.h"
45 #include "Epetra_Vector.h"
46 #include "Epetra_Operator.h"
47 #include "Epetra_Map.h"
48 #include "Epetra_LocalMap.h"
49 
50 #if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
51 # include <Tpetra_ConfigDefs.hpp> // HAVE_TPETRA_EPETRA
52 # if defined(HAVE_TPETRA_EPETRA)
53 # include <Epetra_TsqrAdaptor.hpp>
54 # endif // defined(HAVE_TPETRA_EPETRA)
55 #endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
56 
57 namespace Anasazi {
58 
60 
61 
65  class EpetraMultiVecFailure : public AnasaziError {public:
66  EpetraMultiVecFailure(const std::string& what_arg) : AnasaziError(what_arg)
67  {}};
68 
72  class EpetraOpFailure : public AnasaziError {public:
73  EpetraOpFailure(const std::string& what_arg) : AnasaziError(what_arg)
74  {}};
75 
77 
79 
80 
85 
86  public:
89 
91  virtual Epetra_MultiVector* GetEpetraMultiVec() { return 0; }
92 
94  virtual const Epetra_MultiVector* GetEpetraMultiVec() const { return 0; }
95  };
96 
98 
100  //
101  //--------template class AnasaziEpetraMultiVec-----------------
102  //
104 
111  class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector, public EpetraMultiVecAccessor {
112  public:
114 
115 
117 
122  EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs);
123 
125  EpetraMultiVec(const Epetra_MultiVector & P_vec);
126 
128 
136  EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride=0);
137 
139 
145  EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
146 
148  virtual ~EpetraMultiVec() {};
149 
151 
153 
154 
159  MultiVec<double> * Clone ( const int numvecs ) const;
160 
166  MultiVec<double> * CloneCopy () const;
167 
175  MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
176 
184  MultiVec<double> * CloneViewNonConst ( const std::vector<int>& index );
185 
193  const MultiVec<double> * CloneView ( const std::vector<int>& index ) const;
194 
196 
198  ptrdiff_t GetGlobalLength () const
199  {
200  if ( Map().GlobalIndicesLongLong() )
201  return static_cast<ptrdiff_t>( GlobalLength64() );
202  else
203  return static_cast<ptrdiff_t>( GlobalLength() );
204  }
205 
208  int GetNumberVecs () const { return NumVectors(); }
209 
211 
213 
214 
216  void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A,
217  const Teuchos::SerialDenseMatrix<int,double>& B,
218  double beta );
219 
222  void MvAddMv ( double alpha, const MultiVec<double>& A,
223  double beta, const MultiVec<double>& B);
224 
227  void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B
228 #ifdef HAVE_ANASAZI_EXPERIMENTAL
229  , ConjType conj = Anasazi::CONJ
230 #endif
231  ) const;
232 
235  void MvDot ( const MultiVec<double>& A, std::vector<double> &b
236 #ifdef HAVE_ANASAZI_EXPERIMENTAL
237  , ConjType conj = Anasazi::CONJ
238 #endif
239  ) const;
240 
243  void MvScale ( double alpha ) {
244  TEUCHOS_TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
245  "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
246  }
247 
250  void MvScale ( const std::vector<double>& alpha );
251 
253 
255 
259  void MvNorm ( std::vector<double> & normvec ) const {
260  if (((int)normvec.size() >= GetNumberVecs()) ) {
261  TEUCHOS_TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
262  "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
263  }
264  };
266 
268 
269 
274  void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
275 
278  void MvRandom() {
279  TEUCHOS_TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
280  "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
281  };
282 
285  void MvInit ( double alpha ) {
286  TEUCHOS_TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
287  "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
288  };
289 
291 
292 
294  Epetra_MultiVector* GetEpetraMultiVec() { return this; };
295 
297  const Epetra_MultiVector* GetEpetraMultiVec() const { return this; };
298 
300 
302 
304 
306  void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
308 
309  private:
310  };
311  //-------------------------------------------------------------
312 
314  //
315  //--------template class AnasaziEpetraOp---------------------
316  //
318 
325  class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraOp : public virtual Operator<double> {
326  public:
328 
329 
331  EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
332 
334  ~EpetraOp();
336 
338 
339 
343  void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
345 
346  private:
347 //use pragmas to disable some false-positive warnings for windows
348 // sharedlibs export
349 #ifdef _MSC_VER
350 #pragma warning(push)
351 #pragma warning(disable:4251)
352 #endif
353  Teuchos::RCP<Epetra_Operator> Epetra_Op;
354 #ifdef _MSC_VER
355 #pragma warning(pop)
356 #endif
357  };
358  //-------------------------------------------------------------
359 
361  //
362  //--------template class AnasaziEpetraGenOp--------------------
363  //
365 
379  class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
380  public:
382 
385  EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp,
386  const Teuchos::RCP<Epetra_Operator> &MOp,
387  bool isAInverse = true );
388 
390  ~EpetraGenOp();
391 
393 
395  void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
396 
398 
400  int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
401 
403 
405  int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
406 
408  const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
409 
411  bool UseTranspose() const { return (false); };
412 
414  int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
415 
417  bool HasNormInf() const { return (false); };
418 
420  double NormInf() const { return (-1.0); };
421 
423  const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
424 
426  const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
427 
429  const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
430 
431  private:
432  bool isAInverse;
433 
434 //use pragmas to disable some false-positive warnings for windows
435 // sharedlibs export
436 #ifdef _MSC_VER
437 #pragma warning(push)
438 #pragma warning(disable:4251)
439 #endif
440  Teuchos::RCP<Epetra_Operator> Epetra_AOp;
441  Teuchos::RCP<Epetra_Operator> Epetra_MOp;
442 #ifdef _MSC_VER
443 #pragma warning(pop)
444 #endif
445  };
446 
448  //
449  //--------template class AnasaziEpetraSymOp--------------------
450  //
452 
465  class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
466  public:
468 
470  EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
471 
473  ~EpetraSymOp();
474 
476 
478  void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
479 
481 
483  int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
484 
486 
489  int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
490 
492  const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
493 
495  bool UseTranspose() const { return (false); };
496 
498  int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
499 
501  bool HasNormInf() const { return (false); };
502 
504  double NormInf() const { return (-1.0); };
505 
507  const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
508 
510  const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
511 
513  const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
514 
515  private:
516 
517 //use pragmas to disable false-positive warnings in generating windows sharedlib exports
518 #ifdef _MSC_VER
519 #pragma warning(push)
520 #pragma warning(disable:4251)
521 #endif
522  Teuchos::RCP<Epetra_Operator> Epetra_Op;
523 #ifdef _MSC_VER
524 #pragma warning(pop)
525 #endif
526 
527  bool isTrans_;
528  };
529 
530 
532  //
533  //--------template class AnasaziEpetraSymMVOp---------------------
534  //
536 
549  class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymMVOp : public virtual Operator<double> {
550  public:
552 
554  EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
555  bool isTrans = false );
556 
559 
561 
563  void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
564 
565  private:
566 
567 //use pragmas to disable some false-positive warnings for windows
568 // sharedlibs export
569 #ifdef _MSC_VER
570 #pragma warning(push)
571 #pragma warning(disable:4251)
572 #endif
573  Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
574  Teuchos::RCP<const Epetra_Map> MV_localmap;
575  Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
576 #ifdef _MSC_VER
577 #pragma warning(pop)
578 #endif
579 
580  bool isTrans_;
581  };
582 
584  //
585  //--------template class AnasaziEpetraWSymMVOp---------------------
586  //
588 
601  class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraWSymMVOp : public virtual Operator<double> {
602  public:
604  EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
605  const Teuchos::RCP<Epetra_Operator> &OP );
606 
609 
611 
613  void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
614 
615  private:
616 //use pragmas to disable some false-positive warnings for windows
617 // sharedlibs export
618 #ifdef _MSC_VER
619 #pragma warning(push)
620 #pragma warning(disable:4251)
621 #endif
622  Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
623  Teuchos::RCP<Epetra_Operator> Epetra_OP;
624  Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
625  Teuchos::RCP<const Epetra_Map> MV_localmap;
626  Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
627 #ifdef _MSC_VER
628 #pragma warning(pop)
629 #endif
630  };
631 
633  //
634  //--------template class AnasaziEpetraW2SymMVOp---------------------
635  //
637 
650  class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraW2SymMVOp : public virtual Operator<double> {
651  public:
653  EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
654  const Teuchos::RCP<Epetra_Operator> &OP );
655 
658 
660 
662  void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
663 
664  private:
665 //use pragmas to disable some false-positive warnings for windows
666 // sharedlibs export
667 #ifdef _MSC_VER
668 #pragma warning(push)
669 #pragma warning(disable:4251)
670 #endif
671  Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
672  Teuchos::RCP<Epetra_Operator> Epetra_OP;
673  Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
674  Teuchos::RCP<const Epetra_Map> MV_localmap;
675  Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
676 #ifdef _MSC_VER
677 #pragma warning(pop)
678 #endif
679  };
680 
681 
683  //
684  // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
685  //
687 
698  template<>
699  class MultiVecTraits<double, Epetra_MultiVector>
700  {
701  public:
702 
704 
705 
710  static Teuchos::RCP<Epetra_MultiVector>
711  Clone (const Epetra_MultiVector& mv, const int outNumVecs)
712  {
713  TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs <= 0, std::invalid_argument,
714  "Belos::MultiVecTraits<double, Epetra_MultiVector>::"
715  "Clone(mv, outNumVecs = " << outNumVecs << "): "
716  "outNumVecs must be positive.");
717  // FIXME (mfh 13 Jan 2011) Anasazi currently lets Epetra fill in
718  // the entries of the returned multivector with zeros, but Belos
719  // does not. We retain this different behavior for now, but the
720  // two versions will need to be reconciled.
721  return Teuchos::rcp (new Epetra_MultiVector (mv.Map(), outNumVecs));
722  }
723 
728  static Teuchos::RCP<Epetra_MultiVector>
729  CloneCopy (const Epetra_MultiVector& mv)
730  {
731  return Teuchos::rcp (new Epetra_MultiVector (mv));
732  }
733 
739  static Teuchos::RCP<Epetra_MultiVector>
740  CloneCopy (const Epetra_MultiVector& mv, const std::vector<int>& index)
741  {
742  const int inNumVecs = GetNumberVecs (mv);
743  const int outNumVecs = index.size();
744 
745  // Simple, inexpensive tests of the index vector.
746  TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
747  "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
748  "CloneCopy(mv, index = {}): At least one vector must be"
749  " cloned from mv.");
750  if (outNumVecs > inNumVecs)
751  {
752  std::ostringstream os;
753  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
754  "CloneCopy(mv, index = {";
755  for (int k = 0; k < outNumVecs - 1; ++k)
756  os << index[k] << ", ";
757  os << index[outNumVecs-1] << "}): There are " << outNumVecs
758  << " indices to copy, but only " << inNumVecs << " columns of mv.";
759  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
760  }
761 #ifdef TEUCHOS_DEBUG
762  // In debug mode, we perform more expensive tests of the index
763  // vector, to ensure all the elements are in range.
764  // Dereferencing the iterator is valid because index has length
765  // > 0.
766  const int minIndex = *std::min_element (index.begin(), index.end());
767  const int maxIndex = *std::max_element (index.begin(), index.end());
768 
769  if (minIndex < 0)
770  {
771  std::ostringstream os;
772  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
773  "CloneCopy(mv, index = {";
774  for (int k = 0; k < outNumVecs - 1; ++k)
775  os << index[k] << ", ";
776  os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
777  "the smallest index " << minIndex << " is negative.";
778  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
779  }
780  if (maxIndex >= inNumVecs)
781  {
782  std::ostringstream os;
783  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
784  "CloneCopy(mv, index = {";
785  for (int k = 0; k < outNumVecs - 1; ++k)
786  os << index[k] << ", ";
787  os << index[outNumVecs-1] << "}): Indices must be strictly less than "
788  "the number of vectors " << inNumVecs << " in mv; the largest index "
789  << maxIndex << " is out of bounds.";
790  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
791  }
792 #endif // TEUCHOS_DEBUG
793  // Cast to nonconst, because Epetra_MultiVector's constructor
794  // wants a nonconst int array argument. It doesn't actually
795  // change the entries of the array.
796  std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
797  return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, &tmpind[0], index.size()));
798  // return Teuchos::rcp (new Epetra_MultiVector (::Copy, mv, &tmpind[0], index.size()));
799  }
800 
801  static Teuchos::RCP<Epetra_MultiVector>
802  CloneCopy (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
803  {
804  const int inNumVecs = GetNumberVecs (mv);
805  const int outNumVecs = index.size();
806  const bool validRange = outNumVecs > 0 && index.lbound() >= 0 &&
807  index.ubound() < inNumVecs;
808  if (! validRange)
809  {
810  std::ostringstream os;
811  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,"
812  "index=[" << index.lbound() << ", " << index.ubound() << "]): ";
813  TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
814  os.str() << "Column index range must be nonempty.");
815  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
816  os.str() << "Column index range must be nonnegative.");
817  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= inNumVecs, std::invalid_argument,
818  os.str() << "Column index range must not exceed "
819  "number of vectors " << inNumVecs << " in the "
820  "input multivector.");
821  }
822  return Teuchos::rcp (new Epetra_MultiVector (Copy, mv, index.lbound(), index.size()));
823  }
824 
830  static Teuchos::RCP<Epetra_MultiVector>
831  CloneViewNonConst (Epetra_MultiVector& mv, const std::vector<int>& index)
832  {
833  const int inNumVecs = GetNumberVecs (mv);
834  const int outNumVecs = index.size();
835 
836  // Simple, inexpensive tests of the index vector.
837  TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
838  "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
839  "CloneViewNonConst(mv, index = {}): The output view "
840  "must have at least one column.");
841  if (outNumVecs > inNumVecs)
842  {
843  std::ostringstream os;
844  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
845  "CloneViewNonConst(mv, index = {";
846  for (int k = 0; k < outNumVecs - 1; ++k)
847  os << index[k] << ", ";
848  os << index[outNumVecs-1] << "}): There are " << outNumVecs
849  << " indices to view, but only " << inNumVecs << " columns of mv.";
850  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
851  }
852 #ifdef TEUCHOS_DEBUG
853  // In debug mode, we perform more expensive tests of the index
854  // vector, to ensure all the elements are in range.
855  // Dereferencing the iterator is valid because index has length
856  // > 0.
857  const int minIndex = *std::min_element (index.begin(), index.end());
858  const int maxIndex = *std::max_element (index.begin(), index.end());
859 
860  if (minIndex < 0)
861  {
862  std::ostringstream os;
863  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
864  "CloneViewNonConst(mv, index = {";
865  for (int k = 0; k < outNumVecs - 1; ++k)
866  os << index[k] << ", ";
867  os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
868  "the smallest index " << minIndex << " is negative.";
869  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
870  }
871  if (maxIndex >= inNumVecs)
872  {
873  std::ostringstream os;
874  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
875  "CloneViewNonConst(mv, index = {";
876  for (int k = 0; k < outNumVecs - 1; ++k)
877  os << index[k] << ", ";
878  os << index[outNumVecs-1] << "}): Indices must be strictly less than "
879  "the number of vectors " << inNumVecs << " in mv; the largest index "
880  << maxIndex << " is out of bounds.";
881  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
882  }
883 #endif // TEUCHOS_DEBUG
884  // Cast to nonconst, because Epetra_MultiVector's constructor
885  // wants a nonconst int array argument. It doesn't actually
886  // change the entries of the array.
887  std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
888  return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
889  }
890 
891  static Teuchos::RCP<Epetra_MultiVector>
892  CloneViewNonConst (Epetra_MultiVector& mv, const Teuchos::Range1D& index)
893  {
894  const bool validRange = index.size() > 0 &&
895  index.lbound() >= 0 &&
896  index.ubound() < mv.NumVectors();
897  if (! validRange)
898  {
899  std::ostringstream os;
900  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
901  "NonConst(mv,index=[" << index.lbound() << ", " << index.ubound()
902  << "]): ";
903  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
904  os.str() << "Column index range must be nonempty.");
905  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
906  os.str() << "Column index range must be nonnegative.");
907  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(),
908  std::invalid_argument,
909  os.str() << "Column index range must not exceed "
910  "number of vectors " << mv.NumVectors() << " in "
911  "the input multivector.");
912  }
913  return Teuchos::rcp (new Epetra_MultiVector (View, mv, index.lbound(), index.size()));
914  }
915 
921  static Teuchos::RCP<const Epetra_MultiVector>
922  CloneView (const Epetra_MultiVector& mv, const std::vector<int>& index)
923  {
924  const int inNumVecs = GetNumberVecs (mv);
925  const int outNumVecs = index.size();
926 
927  // Simple, inexpensive tests of the index vector.
928  TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
929  "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
930  "CloneView(mv, index = {}): The output view "
931  "must have at least one column.");
932  if (outNumVecs > inNumVecs)
933  {
934  std::ostringstream os;
935  os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
936  "CloneView(mv, index = {";
937  for (int k = 0; k < outNumVecs - 1; ++k)
938  os << index[k] << ", ";
939  os << index[outNumVecs-1] << "}): There are " << outNumVecs
940  << " indices to view, but only " << inNumVecs << " columns of mv.";
941  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
942  }
943 #ifdef TEUCHOS_DEBUG
944  // In debug mode, we perform more expensive tests of the index
945  // vector, to ensure all the elements are in range.
946  // Dereferencing the iterator is valid because index has length
947  // > 0.
948  const int minIndex = *std::min_element (index.begin(), index.end());
949  const int maxIndex = *std::max_element (index.begin(), index.end());
950 
951  if (minIndex < 0)
952  {
953  std::ostringstream os;
954  os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
955  "CloneView(mv, index = {";
956  for (int k = 0; k < outNumVecs - 1; ++k)
957  os << index[k] << ", ";
958  os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
959  "the smallest index " << minIndex << " is negative.";
960  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
961  }
962  if (maxIndex >= inNumVecs)
963  {
964  std::ostringstream os;
965  os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
966  "CloneView(mv, index = {";
967  for (int k = 0; k < outNumVecs - 1; ++k)
968  os << index[k] << ", ";
969  os << index[outNumVecs-1] << "}): Indices must be strictly less than "
970  "the number of vectors " << inNumVecs << " in mv; the largest index "
971  << maxIndex << " is out of bounds.";
972  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
973  }
974 #endif // TEUCHOS_DEBUG
975  // Cast to nonconst, because Epetra_MultiVector's constructor
976  // wants a nonconst int array argument. It doesn't actually
977  // change the entries of the array.
978  std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
979  return Teuchos::rcp (new Epetra_MultiVector (View, mv, &tmpind[0], index.size()));
980  }
981 
982  static Teuchos::RCP<Epetra_MultiVector>
983  CloneView (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
984  {
985  const bool validRange = index.size() > 0 &&
986  index.lbound() >= 0 &&
987  index.ubound() < mv.NumVectors();
988  if (! validRange)
989  {
990  std::ostringstream os;
991  os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
992  "(mv,index=[" << index.lbound() << ", " << index.ubound()
993  << "]): ";
994  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
995  os.str() << "Column index range must be nonempty.");
996  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
997  os.str() << "Column index range must be nonnegative.");
998  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(),
999  std::invalid_argument,
1000  os.str() << "Column index range must not exceed "
1001  "number of vectors " << mv.NumVectors() << " in "
1002  "the input multivector.");
1003  }
1004  return Teuchos::rcp (new Epetra_MultiVector(View, mv, index.lbound(), index.size()));
1005  }
1006 
1008 
1010 
1011 
1013  static ptrdiff_t GetGlobalLength( const Epetra_MultiVector& mv )
1014  {
1015  if (mv.Map().GlobalIndicesLongLong())
1016  return static_cast<ptrdiff_t>( mv.GlobalLength64() );
1017  else
1018  return static_cast<ptrdiff_t>( mv.GlobalLength() );
1019  }
1020 
1022  static int GetNumberVecs( const Epetra_MultiVector& mv )
1023  { return mv.NumVectors(); }
1024 
1025  static bool HasConstantStride( const Epetra_MultiVector& mv )
1026  { return mv.ConstantStride(); }
1028 
1030 
1031 
1034  static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A,
1035  const Teuchos::SerialDenseMatrix<int,double>& B,
1036  double beta, Epetra_MultiVector& mv )
1037  {
1038  Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
1039  Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
1040 
1041  TEUCHOS_TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
1042  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTimesMatAddMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
1043  }
1044 
1047  static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
1048  {
1049  // epetra mv.Update(alpha,A,beta,B,gamma) will check
1050  // alpha == 0.0
1051  // and
1052  // beta == 0.0
1053  // and based on this will call
1054  // mv.Update(beta,B,gamma)
1055  // or
1056  // mv.Update(alpha,A,gamma)
1057  //
1058  // mv.Update(alpha,A,gamma)
1059  // will then check for one of
1060  // gamma == 0
1061  // or
1062  // gamma == 1
1063  // or
1064  // alpha == 1
1065  // in that order. however, it will not look for the combination
1066  // alpha == 1 and gamma = 0
1067  // which is a common use case when we wish to assign
1068  // mv = A (in which case alpha == 1, beta == gamma == 0)
1069  // or
1070  // mv = B (in which case beta == 1, alpha == gamma == 0)
1071  //
1072  // therefore, we will check for these use cases ourselves
1073  if (beta == 0.0) {
1074  if (alpha == 1.0) {
1075  // assign
1076  mv = A;
1077  }
1078  else {
1079  // single update
1080  TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
1081  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
1082  }
1083  }
1084  else if (alpha == 0.0) {
1085  if (beta == 1.0) {
1086  // assign
1087  mv = B;
1088  }
1089  else {
1090  // single update
1091  TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
1092  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
1093  }
1094  }
1095  else {
1096  // double update
1097  TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
1098  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
1099  }
1100  }
1101 
1104  static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
1105 #ifdef HAVE_ANASAZI_EXPERIMENTAL
1106  , ConjType conj = Anasazi::CONJ
1107 #endif
1108  )
1109  {
1110  Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
1111  Epetra_MultiVector B_Pvec(::View, LocalMap, B.values(), B.stride(), B.numCols());
1112 
1113  TEUCHOS_TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
1114  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
1115  }
1116 
1119  static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
1120 #ifdef HAVE_ANASAZI_EXPERIMENTAL
1121  , ConjType conj = Anasazi::CONJ
1122 #endif
1123  )
1124  {
1125 #ifdef TEUCHOS_DEBUG
1126  TEUCHOS_TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
1127  "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
1128  TEUCHOS_TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
1129  "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
1130 #endif
1131  TEUCHOS_TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
1132  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");
1133  }
1134 
1136 
1138 
1142  static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
1143  {
1144 #ifdef TEUCHOS_DEBUG
1145  TEUCHOS_TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
1146  "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
1147 #endif
1148  TEUCHOS_TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
1149  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
1150  }
1151 
1153 
1155 
1156 
1158  static void
1159  SetBlock (const Epetra_MultiVector& A,
1160  const std::vector<int>& index,
1161  Epetra_MultiVector& mv)
1162  {
1163  const int inNumVecs = GetNumberVecs (A);
1164  const int outNumVecs = index.size();
1165 
1166  // FIXME (mfh 13 Jan 2011) Belos allows A to have more columns
1167  // than index.size(), in which case we just take the first
1168  // index.size() columns of A. Anasazi requires that A have the
1169  // same number of columns as index.size(). Changing Anasazi's
1170  // behavior should not break existing Anasazi solvers, but the
1171  // tests need to be done.
1172  if (inNumVecs != outNumVecs)
1173  {
1174  std::ostringstream os;
1175  os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
1176  "SetBlock(A, mv, index = {";
1177  if (outNumVecs > 0)
1178  {
1179  for (int k = 0; k < outNumVecs - 1; ++k)
1180  os << index[k] << ", ";
1181  os << index[outNumVecs-1];
1182  }
1183  os << "}): A has only " << inNumVecs << " columns, but there are "
1184  << outNumVecs << " indices in the index vector.";
1185  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
1186  }
1187  // Make a view of the columns of mv indicated by the index std::vector.
1188  Teuchos::RCP<Epetra_MultiVector> mv_view = CloneViewNonConst (mv, index);
1189 
1190  // View of columns [0, outNumVecs-1] of the source multivector A.
1191  // If A has fewer columns than mv_view, then create a view of
1192  // the first outNumVecs columns of A.
1193  Teuchos::RCP<const Epetra_MultiVector> A_view;
1194  if (outNumVecs == inNumVecs)
1195  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
1196  else
1197  A_view = CloneView (A, Teuchos::Range1D(0, outNumVecs - 1));
1198 
1199  // Assignment calls Epetra_MultiVector::Assign(), which deeply
1200  // copies the data directly, ignoring the underlying
1201  // Epetra_Map(s). If A and mv don't have the same data
1202  // distribution (Epetra_Map), this may result in incorrect or
1203  // undefined behavior. Epetra_MultiVector::Update() also
1204  // ignores the Epetra_Maps, so we might as well just use the
1205  // (perhaps slightly cheaper) Assign() method via operator=().
1206  *mv_view = *A_view;
1207  }
1208 
1209  static void
1210  SetBlock (const Epetra_MultiVector& A,
1211  const Teuchos::Range1D& index,
1212  Epetra_MultiVector& mv)
1213  {
1214  const int numColsA = A.NumVectors();
1215  const int numColsMv = mv.NumVectors();
1216  // 'index' indexes into mv; it's the index set of the target.
1217  const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
1218  // We can't take more columns out of A than A has.
1219  const bool validSource = index.size() <= numColsA;
1220 
1221  if (! validIndex || ! validSource)
1222  {
1223  std::ostringstream os;
1224  os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::SetBlock"
1225  "(A, index=[" << index.lbound() << ", " << index.ubound() << "], "
1226  "mv): ";
1227  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
1228  os.str() << "Range lower bound must be nonnegative.");
1229  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
1230  os.str() << "Range upper bound must be less than "
1231  "the number of columns " << numColsA << " in the "
1232  "'mv' output argument.");
1233  TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
1234  os.str() << "Range must have no more elements than"
1235  " the number of columns " << numColsA << " in the "
1236  "'A' input argument.");
1237  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
1238  }
1239 
1240  // View of columns [index.lbound(), index.ubound()] of the
1241  // target multivector mv. We avoid view creation overhead by
1242  // only creating a view if the index range is different than [0,
1243  // (# columns in mv) - 1].
1244  Teuchos::RCP<Epetra_MultiVector> mv_view;
1245  if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
1246  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
1247  else
1248  mv_view = CloneViewNonConst (mv, index);
1249 
1250  // View of columns [0, index.size()-1] of the source multivector
1251  // A. If A has fewer columns than mv_view, then create a view
1252  // of the first index.size() columns of A.
1253  Teuchos::RCP<const Epetra_MultiVector> A_view;
1254  if (index.size() == numColsA)
1255  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
1256  else
1257  A_view = CloneView (A, Teuchos::Range1D(0, index.size()-1));
1258 
1259  // Assignment calls Epetra_MultiVector::Assign(), which deeply
1260  // copies the data directly, ignoring the underlying
1261  // Epetra_Map(s). If A and mv don't have the same data
1262  // distribution (Epetra_Map), this may result in incorrect or
1263  // undefined behavior. Epetra_MultiVector::Update() also
1264  // ignores the Epetra_Maps, so we might as well just use the
1265  // (perhaps slightly cheaper) Assign() method via operator=().
1266  *mv_view = *A_view;
1267  }
1268 
1269  static void
1270  Assign (const Epetra_MultiVector& A,
1271  Epetra_MultiVector& mv)
1272  {
1273  const int numColsA = GetNumberVecs (A);
1274  const int numColsMv = GetNumberVecs (mv);
1275  if (numColsA > numColsMv)
1276  {
1277  std::ostringstream os;
1278  os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::Assign"
1279  "(A, mv): ";
1280  TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
1281  os.str() << "Input multivector 'A' has "
1282  << numColsA << " columns, but output multivector "
1283  "'mv' has only " << numColsMv << " columns.");
1284  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
1285  }
1286  // View of the first [0, numColsA-1] columns of mv.
1287  Teuchos::RCP<Epetra_MultiVector> mv_view;
1288  if (numColsMv == numColsA)
1289  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
1290  else // numColsMv > numColsA
1291  mv_view = CloneView (mv, Teuchos::Range1D(0, numColsA - 1));
1292 
1293  // Assignment calls Epetra_MultiVector::Assign(), which deeply
1294  // copies the data directly, ignoring the underlying
1295  // Epetra_Map(s). If A and mv don't have the same data
1296  // distribution (Epetra_Map), this may result in incorrect or
1297  // undefined behavior. Epetra_MultiVector::Update() also
1298  // ignores the Epetra_Maps, so we might as well just use the
1299  // (perhaps slightly cheaper) Assign() method via operator=().
1300  *mv_view = A;
1301  }
1302 
1305  static void MvScale ( Epetra_MultiVector& mv, double alpha )
1306  {
1307  TEUCHOS_TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
1308  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value.");
1309  }
1310 
1313  static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
1314  {
1315  // Check to make sure the vector is as long as the multivector has columns.
1316  int numvecs = mv.NumVectors();
1317 #ifdef TEUCHOS_DEBUG
1318  TEUCHOS_TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
1319  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
1320 #endif
1321  for (int i=0; i<numvecs; i++) {
1322  TEUCHOS_TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i])!=0, EpetraMultiVecFailure,
1323  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
1324  }
1325  }
1326 
1329  static void MvRandom( Epetra_MultiVector& mv )
1330  {
1331  TEUCHOS_TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
1332  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
1333  }
1334 
1337  static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
1338  {
1339  TEUCHOS_TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
1340  "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
1341  }
1342 
1344 
1346 
1347 
1350  static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
1351  { os << mv << std::endl; }
1352 
1354 
1355 #if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
1356 # if defined(HAVE_TPETRA_EPETRA)
1357  typedef Epetra::TsqrAdaptor tsqr_adaptor_type;
1363 # endif // defined(HAVE_TPETRA_EPETRA)
1364 #endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
1365  };
1366 
1368  //
1369  // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
1370  //
1372 
1384  template <>
1385  class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
1386  {
1387  public:
1388 
1392  static void Apply ( const Epetra_Operator& Op,
1393  const Epetra_MultiVector& x,
1394  Epetra_MultiVector& y )
1395  {
1396 #ifdef TEUCHOS_DEBUG
1397  TEUCHOS_TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
1398  "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
1399 #endif
1400  int ret = Op.Apply(x,y);
1401  TEUCHOS_TEST_FOR_EXCEPTION(ret != 0, OperatorError,
1402  "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
1403  }
1404 
1405  };
1406 
1407 } // end of Anasazi namespace
1408 
1409 #endif
1410 // end of file ANASAZI_EPETRA_ADAPTER_HPP
void MvRandom()
Fill the vectors in *this with random numbers.
Adapter class for creating an operators often used in solving generalized eigenproblems.
EpetraMultiVecAccessor is an interfaceto allow any Anasazi::MultiVec implementation that is based on ...
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
static void MvAddMv(double alpha, const Epetra_MultiVector &A, double beta, const Epetra_MultiVector &B, Epetra_MultiVector &mv)
Replace mv with .
void MvInit(double alpha)
Replace each element of the vectors in *this with alpha.
static void SetBlock(const Epetra_MultiVector &A, const std::vector< int > &index, Epetra_MultiVector &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
virtual ~EpetraMultiVecAccessor()
Destructor.
Adapter class for creating a weighted symmetric operator from an Epetra_MultiVector and Epetra_Operat...
static void MvInit(Epetra_MultiVector &mv, double alpha=Teuchos::ScalarTraits< double >::zero())
Replace each element of the vectors in mv with alpha.
Virtual base class which defines basic traits for the operator type.
const char * Label() const
Returns a character string describing the operator.
void MvNorm(std::vector< double > &normvec) const
Compute the 2-norm of each individual vector of *this. Upon return, normvec[i] holds the 2-norm of th...
virtual Epetra_MultiVector * GetEpetraMultiVec()
Return the pointer to the Epetra_MultiVector object.
static void MvPrint(const Epetra_MultiVector &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
static void MvRandom(Epetra_MultiVector &mv)
Replace the vectors in mv with random vectors.
static void MvTransMv(double alpha, const Epetra_MultiVector &A, const Epetra_MultiVector &mv, Teuchos::SerialDenseMatrix< int, double > &B)
Compute a dense matrix B through the matrix-matrix multiply .
An exception class parent to all Anasazi exceptions.
Interface for multivectors used by Anasazi&#39; linear solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Basic adapter class for Anasazi::Operator that uses Epetra_Operator.
static void MvDot(const Epetra_MultiVector &A, const Epetra_MultiVector &B, std::vector< double > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
const char * Label() const
Returns a character string describing the operator.
bool HasNormInf() const
Returns true if this object can provide an approximate inf-norm [always false for this operator]...
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
ConjType
Enumerated types used to specify conjugation arguments.
static Teuchos::RCP< Epetra_MultiVector > CloneCopy(const Epetra_MultiVector &mv, const std::vector< int > &index)
Creates a new Epetra_MultiVector and copies the selected contents of mv into the new vector (deep cop...
Epetra_MultiVector * GetEpetraMultiVec()
Return the pointer to the Epetra_MultiVector object.
bool HasNormInf() const
Returns true if this object can provide an approximate inf-norm [always false for this operator]...
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
static void MvScale(Epetra_MultiVector &mv, const std::vector< double > &alpha)
Scale each element of the i-th vector in mv with alpha[i].
const Epetra_Comm & Comm() const
Returns the Epetra_Comm communicator associated with this operator.
static void MvTimesMatAddMv(double alpha, const Epetra_MultiVector &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta, Epetra_MultiVector &mv)
Update mv with .
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< const Epetra_MultiVector > CloneView(const Epetra_MultiVector &mv, const std::vector< int > &index)
Creates a new const Epetra_MultiVector that shares the selected contents of mv (shallow copy)...
Adapter class for creating a weighted operator from an Epetra_MultiVector and Epetra_Operator.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
bool UseTranspose() const
Returns the current UseTranspose setting [always false for this operator].
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
const Epetra_MultiVector * GetEpetraMultiVec() const
Return the pointer to the Epetra_MultiVector object.
static Teuchos::RCP< Epetra_MultiVector > Clone(const Epetra_MultiVector &mv, const int outNumVecs)
Creates a new empty Epetra_MultiVector containing numVecs columns.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
double NormInf() const
Returns the infinity norm of the global matrix [not functional for this operator].
static ptrdiff_t GetGlobalLength(const Epetra_MultiVector &mv)
Obtain the vector length of mv.
int SetUseTranspose(bool)
If set true, the transpose of this operator will be applied [not functional for this operator]...
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Adapter class for creating a symmetric operator from an Epetra_Operator.
void MvScale(double alpha)
Scale each element of the vectors in *this with alpha.
static void MvScale(Epetra_MultiVector &mv, double alpha)
Scale each element of the vectors in mv with alpha.
static void Apply(const Epetra_Operator &Op, const Epetra_MultiVector &x, Epetra_MultiVector &y)
This method takes the Epetra_MultiVector x and applies the Epetra_Operator Op to it resulting in the ...
const Epetra_Comm & Comm() const
Returns the Epetra_Comm communicator associated with this operator.
Templated virtual class for creating operators that can interface with the Anasazi::OperatorTraits cl...
static Teuchos::RCP< Epetra_MultiVector > CloneViewNonConst(Epetra_MultiVector &mv, const std::vector< int > &index)
Creates a new Epetra_MultiVector that shares the selected contents of mv (shallow copy)...
EpetraMultiVecFailure is thrown when a return value from an Epetra call on an Epetra_MultiVector is n...
Types and exceptions used within Anasazi solvers and interfaces.
bool UseTranspose() const
Returns the current UseTranspose setting [always false for this operator].
virtual const Epetra_MultiVector * GetEpetraMultiVec() const
Return the pointer to the Epetra_MultiVector object.
Interface for multivectors used by Anasazi&#39;s linear solvers.
Adapter class for creating a symmetric operator from an Epetra_MultiVector.
static void MvNorm(const Epetra_MultiVector &mv, std::vector< double > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
static int GetNumberVecs(const Epetra_MultiVector &mv)
Obtain the number of vectors in mv.
static Teuchos::RCP< Epetra_MultiVector > CloneCopy(const Epetra_MultiVector &mv)
Creates a new Epetra_MultiVector and copies contents of mv into the new vector (deep copy)...
Anasazi&#39;s templated virtual class for constructing an operator that can interface with the OperatorTr...
virtual ~EpetraMultiVec()
Destructor.
int SetUseTranspose(bool)
If set true, the transpose of this operator will be applied [not functional for this operator]...
double NormInf() const
Returns the infinity norm of the global matrix [not functional for this operator].
EpetraOpFailure is thrown when a return value from an Epetra call on an Epetra_Operator is non-zero...
void MvPrint(std::ostream &os) const
Print *this EpetraMultiVec.
Exceptions thrown to signal error in operator application.
Basic adapter class for Anasazi::MultiVec that uses Epetra_MultiVector.