Anasazi Version of the Day
Loading...
Searching...
No Matches
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 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
46#ifndef ANASAZI_EPETRA_ADAPTER_HPP
47#define ANASAZI_EPETRA_ADAPTER_HPP
48
49#include "AnasaziConfigDefs.hpp"
50#include "Anasaziepetra_DLLExportMacro.h"
51#include "AnasaziTypes.hpp"
52#include "AnasaziMultiVec.hpp"
53#include "AnasaziOperator.hpp"
55
56#include "Teuchos_Assert.hpp"
57#include "Teuchos_SerialDenseMatrix.hpp"
58#include "Teuchos_FancyOStream.hpp"
59
60#include "Epetra_MultiVector.h"
61#include "Epetra_Vector.h"
62#include "Epetra_Operator.h"
63#include "Epetra_Map.h"
64#include "Epetra_LocalMap.h"
65#include "Epetra_Comm.h"
66
67#if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
68# include <Tpetra_ConfigDefs.hpp> // HAVE_TPETRA_EPETRA
69# if defined(HAVE_TPETRA_EPETRA)
70# include <Epetra_TsqrAdaptor.hpp>
71# endif // defined(HAVE_TPETRA_EPETRA)
72#endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
73
74namespace Anasazi {
75
77
78
82 class EpetraMultiVecFailure : public AnasaziError {public:
83 EpetraMultiVecFailure(const std::string& what_arg) : AnasaziError(what_arg)
84 {}};
85
89 class EpetraOpFailure : public AnasaziError {public:
90 EpetraOpFailure(const std::string& what_arg) : AnasaziError(what_arg)
91 {}};
92
94
96
97
102
103 public:
106
108 virtual Epetra_MultiVector* GetEpetraMultiVec() { return 0; }
109
111 virtual const Epetra_MultiVector* GetEpetraMultiVec() const { return 0; }
112 };
113
115
117 //
118 //--------template class AnasaziEpetraMultiVec-----------------
119 //
121
128 class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraMultiVec : public MultiVec<double>, public Epetra_MultiVector, public EpetraMultiVecAccessor {
129 public:
131
132
134
139 EpetraMultiVec(const Epetra_BlockMap& Map_in, const int numvecs);
140
142 EpetraMultiVec(const Epetra_MultiVector & P_vec);
143
145
153 EpetraMultiVec(const Epetra_BlockMap& Map_in, double * array, const int numvecs, const int stride=0);
154
156
162 EpetraMultiVec(Epetra_DataAccess CV, const Epetra_MultiVector& P_vec, const std::vector<int>& index);
163
165 virtual ~EpetraMultiVec() {};
166
168
170
171
176 MultiVec<double> * Clone ( const int numvecs ) const;
177
183 MultiVec<double> * CloneCopy () const;
184
192 MultiVec<double> * CloneCopy ( const std::vector<int>& index ) const;
193
201 MultiVec<double> * CloneViewNonConst ( const std::vector<int>& index );
202
210 const MultiVec<double> * CloneView ( const std::vector<int>& index ) const;
211
213
215 ptrdiff_t GetGlobalLength () const
216 {
217 if ( Map().GlobalIndicesLongLong() )
218 return static_cast<ptrdiff_t>( GlobalLength64() );
219 else
220 return static_cast<ptrdiff_t>( GlobalLength() );
221 }
222
225 int GetNumberVecs () const { return NumVectors(); }
226
228
230
231
233 void MvTimesMatAddMv ( double alpha, const MultiVec<double>& A,
234 const Teuchos::SerialDenseMatrix<int,double>& B,
235 double beta );
236
239 void MvAddMv ( double alpha, const MultiVec<double>& A,
240 double beta, const MultiVec<double>& B);
241
244 void MvTransMv ( double alpha, const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B
245#ifdef HAVE_ANASAZI_EXPERIMENTAL
246 , ConjType conj = Anasazi::CONJ
247#endif
248 ) const;
249
252 void MvDot ( const MultiVec<double>& A, std::vector<double> &b
253#ifdef HAVE_ANASAZI_EXPERIMENTAL
254 , ConjType conj = Anasazi::CONJ
255#endif
256 ) const;
257
260 void MvScale ( double alpha ) {
261 TEUCHOS_TEST_FOR_EXCEPTION( this->Scale( alpha )!=0, EpetraMultiVecFailure,
262 "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
263 }
264
267 void MvScale ( const std::vector<double>& alpha );
268
270
272
276 void MvNorm ( std::vector<double> & normvec ) const {
277 if (((int)normvec.size() >= GetNumberVecs()) ) {
278 TEUCHOS_TEST_FOR_EXCEPTION( this->Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
279 "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
280 }
281 };
283
285
286
291 void SetBlock ( const MultiVec<double>& A, const std::vector<int>& index );
292
295 void MvRandom() {
296 TEUCHOS_TEST_FOR_EXCEPTION( this->Random()!=0, EpetraMultiVecFailure,
297 "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
298 };
299
302 void MvInit ( double alpha ) {
303 TEUCHOS_TEST_FOR_EXCEPTION( this->PutScalar( alpha )!=0, EpetraMultiVecFailure,
304 "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
305 };
306
308
309
311 Epetra_MultiVector* GetEpetraMultiVec() { return this; };
312
314 const Epetra_MultiVector* GetEpetraMultiVec() const { return this; };
315
317
319
321
323 void MvPrint( std::ostream& os ) const { os << *this << std::endl; };
325
326 private:
327 };
328 //-------------------------------------------------------------
329
331 //
332 //--------template class AnasaziEpetraOp---------------------
333 //
335
342 class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraOp : public virtual Operator<double> {
343 public:
345
346
348 EpetraOp(const Teuchos::RCP<Epetra_Operator> &Op );
349
351 ~EpetraOp();
353
355
356
360 void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
362
363 private:
364//use pragmas to disable some false-positive warnings for windows
365// sharedlibs export
366#ifdef _MSC_VER
367#pragma warning(push)
368#pragma warning(disable:4251)
369#endif
370 Teuchos::RCP<Epetra_Operator> Epetra_Op;
371#ifdef _MSC_VER
372#pragma warning(pop)
373#endif
374 };
375 //-------------------------------------------------------------
376
378 //
379 //--------template class AnasaziEpetraGenOp--------------------
380 //
382
396 class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraGenOp : public virtual Operator<double>, public virtual Epetra_Operator {
397 public:
399
402 EpetraGenOp(const Teuchos::RCP<Epetra_Operator> &AOp,
403 const Teuchos::RCP<Epetra_Operator> &MOp,
404 bool isAInverse = true );
405
407 ~EpetraGenOp();
408
410
412 void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
413
415
417 int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
418
420
422 int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
423
425 const char* Label() const { return "Epetra_Operator applying A^{-1}M"; };
426
428 bool UseTranspose() const { return (false); };
429
431 int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
432
434 bool HasNormInf() const { return (false); };
435
437 double NormInf() const { return (-1.0); };
438
440 const Epetra_Comm& Comm() const { return Epetra_AOp->Comm(); };
441
443 const Epetra_Map& OperatorDomainMap() const { return Epetra_AOp->OperatorDomainMap(); };
444
446 const Epetra_Map& OperatorRangeMap() const { return Epetra_AOp->OperatorRangeMap(); };
447
448 private:
449 bool isAInverse;
450
451//use pragmas to disable some false-positive warnings for windows
452// sharedlibs export
453#ifdef _MSC_VER
454#pragma warning(push)
455#pragma warning(disable:4251)
456#endif
457 Teuchos::RCP<Epetra_Operator> Epetra_AOp;
458 Teuchos::RCP<Epetra_Operator> Epetra_MOp;
459#ifdef _MSC_VER
460#pragma warning(pop)
461#endif
462 };
463
465 //
466 //--------template class AnasaziEpetraSymOp--------------------
467 //
469
482 class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymOp : public virtual Operator<double>, public virtual Epetra_Operator {
483 public:
485
487 EpetraSymOp(const Teuchos::RCP<Epetra_Operator> &Op, bool isTrans = false );
488
490 ~EpetraSymOp();
491
493
495 void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
496
498
500 int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
501
503
506 int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const;
507
509 const char* Label() const { return "Epetra_Operator applying A^TA or AA^T"; };
510
512 bool UseTranspose() const { return (false); };
513
515 int SetUseTranspose(bool /*UseTranspose_in*/) { return 0; };
516
518 bool HasNormInf() const { return (false); };
519
521 double NormInf() const { return (-1.0); };
522
524 const Epetra_Comm& Comm() const { return Epetra_Op->Comm(); };
525
527 const Epetra_Map& OperatorDomainMap() const { return Epetra_Op->OperatorDomainMap(); };
528
530 const Epetra_Map& OperatorRangeMap() const { return Epetra_Op->OperatorRangeMap(); };
531
532 private:
533
534//use pragmas to disable false-positive warnings in generating windows sharedlib exports
535#ifdef _MSC_VER
536#pragma warning(push)
537#pragma warning(disable:4251)
538#endif
539 Teuchos::RCP<Epetra_Operator> Epetra_Op;
540#ifdef _MSC_VER
541#pragma warning(pop)
542#endif
543
544 bool isTrans_;
545 };
546
547
549 //
550 //--------template class AnasaziEpetraSymMVOp---------------------
551 //
553
566 class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraSymMVOp : public virtual Operator<double> {
567 public:
569
571 EpetraSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
572 bool isTrans = false );
573
576
578
580 void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
581
582 private:
583
584//use pragmas to disable some false-positive warnings for windows
585// sharedlibs export
586#ifdef _MSC_VER
587#pragma warning(push)
588#pragma warning(disable:4251)
589#endif
590 Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
591 Teuchos::RCP<const Epetra_Map> MV_localmap;
592 Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
593#ifdef _MSC_VER
594#pragma warning(pop)
595#endif
596
597 bool isTrans_;
598 };
599
601 //
602 //--------template class AnasaziEpetraWSymMVOp---------------------
603 //
605
618 class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraWSymMVOp : public virtual Operator<double> {
619 public:
621 EpetraWSymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
622 const Teuchos::RCP<Epetra_Operator> &OP );
623
626
628
630 void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
631
632 private:
633//use pragmas to disable some false-positive warnings for windows
634// sharedlibs export
635#ifdef _MSC_VER
636#pragma warning(push)
637#pragma warning(disable:4251)
638#endif
639 Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
640 Teuchos::RCP<Epetra_Operator> Epetra_OP;
641 Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
642 Teuchos::RCP<const Epetra_Map> MV_localmap;
643 Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
644#ifdef _MSC_VER
645#pragma warning(pop)
646#endif
647 };
648
650 //
651 //--------template class AnasaziEpetraW2SymMVOp---------------------
652 //
654
667 class ANASAZIEPETRA_LIB_DLL_EXPORT EpetraW2SymMVOp : public virtual Operator<double> {
668 public:
670 EpetraW2SymMVOp(const Teuchos::RCP<const Epetra_MultiVector> &MV,
671 const Teuchos::RCP<Epetra_Operator> &OP );
672
675
677
679 void Apply ( const MultiVec<double>& X, MultiVec<double>& Y ) const;
680
681 private:
682//use pragmas to disable some false-positive warnings for windows
683// sharedlibs export
684#ifdef _MSC_VER
685#pragma warning(push)
686#pragma warning(disable:4251)
687#endif
688 Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
689 Teuchos::RCP<Epetra_Operator> Epetra_OP;
690 Teuchos::RCP<Epetra_MultiVector> Epetra_WMV;
691 Teuchos::RCP<const Epetra_Map> MV_localmap;
692 Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
693#ifdef _MSC_VER
694#pragma warning(pop)
695#endif
696 };
697
698
700 //
701 // Implementation of the Anasazi::MultiVecTraits for Epetra::MultiVector.
702 //
704
715 template<>
716 class MultiVecTraits<double, Epetra_MultiVector>
717 {
718 public:
719
721
722
727 static Teuchos::RCP<Epetra_MultiVector>
728 Clone (const Epetra_MultiVector& mv, const int outNumVecs)
729 {
730 TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs <= 0, std::invalid_argument,
731 "Belos::MultiVecTraits<double, Epetra_MultiVector>::"
732 "Clone(mv, outNumVecs = " << outNumVecs << "): "
733 "outNumVecs must be positive.");
734 // FIXME (mfh 13 Jan 2011) Anasazi currently lets Epetra fill in
735 // the entries of the returned multivector with zeros, but Belos
736 // does not. We retain this different behavior for now, but the
737 // two versions will need to be reconciled.
738 return Teuchos::rcp (new Epetra_MultiVector (mv.Map(), outNumVecs));
739 }
740
745 static Teuchos::RCP<Epetra_MultiVector>
746 CloneCopy (const Epetra_MultiVector& mv)
747 {
748 return Teuchos::rcp (new Epetra_MultiVector (mv));
749 }
750
756 static Teuchos::RCP<Epetra_MultiVector>
757 CloneCopy (const Epetra_MultiVector& mv, const std::vector<int>& index)
758 {
759 const int inNumVecs = GetNumberVecs (mv);
760 const int outNumVecs = index.size();
761
762 // Simple, inexpensive tests of the index vector.
763 TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
764 "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
765 "CloneCopy(mv, index = {}): At least one vector must be"
766 " cloned from mv.");
767 if (outNumVecs > inNumVecs)
768 {
769 std::ostringstream os;
770 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
771 "CloneCopy(mv, index = {";
772 for (int k = 0; k < outNumVecs - 1; ++k)
773 os << index[k] << ", ";
774 os << index[outNumVecs-1] << "}): There are " << outNumVecs
775 << " indices to copy, but only " << inNumVecs << " columns of mv.";
776 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
777 }
778#ifdef TEUCHOS_DEBUG
779 // In debug mode, we perform more expensive tests of the index
780 // vector, to ensure all the elements are in range.
781 // Dereferencing the iterator is valid because index has length
782 // > 0.
783 const int minIndex = *std::min_element (index.begin(), index.end());
784 const int maxIndex = *std::max_element (index.begin(), index.end());
785
786 if (minIndex < 0)
787 {
788 std::ostringstream os;
789 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
790 "CloneCopy(mv, index = {";
791 for (int k = 0; k < outNumVecs - 1; ++k)
792 os << index[k] << ", ";
793 os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
794 "the smallest index " << minIndex << " is negative.";
795 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
796 }
797 if (maxIndex >= inNumVecs)
798 {
799 std::ostringstream os;
800 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
801 "CloneCopy(mv, index = {";
802 for (int k = 0; k < outNumVecs - 1; ++k)
803 os << index[k] << ", ";
804 os << index[outNumVecs-1] << "}): Indices must be strictly less than "
805 "the number of vectors " << inNumVecs << " in mv; the largest index "
806 << maxIndex << " is out of bounds.";
807 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
808 }
809#endif // TEUCHOS_DEBUG
810 // Cast to nonconst, because Epetra_MultiVector's constructor
811 // wants a nonconst int array argument. It doesn't actually
812 // change the entries of the array.
813 std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
814 return Teuchos::rcp (new Epetra_MultiVector (Epetra_DataAccess::Copy, mv, &tmpind[0], index.size()));
815 }
816
817 static Teuchos::RCP<Epetra_MultiVector>
818 CloneCopy (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
819 {
820 const int inNumVecs = GetNumberVecs (mv);
821 const int outNumVecs = index.size();
822 const bool validRange = outNumVecs > 0 && index.lbound() >= 0 &&
823 index.ubound() < inNumVecs;
824 if (! validRange)
825 {
826 std::ostringstream os;
827 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::Clone(mv,"
828 "index=[" << index.lbound() << ", " << index.ubound() << "]): ";
829 TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
830 os.str() << "Column index range must be nonempty.");
831 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
832 os.str() << "Column index range must be nonnegative.");
833 TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= inNumVecs, std::invalid_argument,
834 os.str() << "Column index range must not exceed "
835 "number of vectors " << inNumVecs << " in the "
836 "input multivector.");
837 }
838 return Teuchos::rcp (new Epetra_MultiVector (Epetra_DataAccess::Copy, mv, index.lbound(), index.size()));
839 }
840
846 static Teuchos::RCP<Epetra_MultiVector>
847 CloneViewNonConst (Epetra_MultiVector& mv, const std::vector<int>& index)
848 {
849 const int inNumVecs = GetNumberVecs (mv);
850 const int outNumVecs = index.size();
851
852 // Simple, inexpensive tests of the index vector.
853 TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
854 "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
855 "CloneViewNonConst(mv, index = {}): The output view "
856 "must have at least one column.");
857 if (outNumVecs > inNumVecs)
858 {
859 std::ostringstream os;
860 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
861 "CloneViewNonConst(mv, index = {";
862 for (int k = 0; k < outNumVecs - 1; ++k)
863 os << index[k] << ", ";
864 os << index[outNumVecs-1] << "}): There are " << outNumVecs
865 << " indices to view, but only " << inNumVecs << " columns of mv.";
866 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
867 }
868#ifdef TEUCHOS_DEBUG
869 // In debug mode, we perform more expensive tests of the index
870 // vector, to ensure all the elements are in range.
871 // Dereferencing the iterator is valid because index has length
872 // > 0.
873 const int minIndex = *std::min_element (index.begin(), index.end());
874 const int maxIndex = *std::max_element (index.begin(), index.end());
875
876 if (minIndex < 0)
877 {
878 std::ostringstream os;
879 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
880 "CloneViewNonConst(mv, index = {";
881 for (int k = 0; k < outNumVecs - 1; ++k)
882 os << index[k] << ", ";
883 os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
884 "the smallest index " << minIndex << " is negative.";
885 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
886 }
887 if (maxIndex >= inNumVecs)
888 {
889 std::ostringstream os;
890 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::"
891 "CloneViewNonConst(mv, index = {";
892 for (int k = 0; k < outNumVecs - 1; ++k)
893 os << index[k] << ", ";
894 os << index[outNumVecs-1] << "}): Indices must be strictly less than "
895 "the number of vectors " << inNumVecs << " in mv; the largest index "
896 << maxIndex << " is out of bounds.";
897 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
898 }
899#endif // TEUCHOS_DEBUG
900 // Cast to nonconst, because Epetra_MultiVector's constructor
901 // wants a nonconst int array argument. It doesn't actually
902 // change the entries of the array.
903 std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
904 return Teuchos::rcp (new Epetra_MultiVector (Epetra_DataAccess::View, mv, &tmpind[0], index.size()));
905 }
906
907 static Teuchos::RCP<Epetra_MultiVector>
908 CloneViewNonConst (Epetra_MultiVector& mv, const Teuchos::Range1D& index)
909 {
910 const bool validRange = index.size() > 0 &&
911 index.lbound() >= 0 &&
912 index.ubound() < mv.NumVectors();
913 if (! validRange)
914 {
915 std::ostringstream os;
916 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
917 "NonConst(mv,index=[" << index.lbound() << ", " << index.ubound()
918 << "]): ";
919 TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
920 os.str() << "Column index range must be nonempty.");
921 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
922 os.str() << "Column index range must be nonnegative.");
923 TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(),
924 std::invalid_argument,
925 os.str() << "Column index range must not exceed "
926 "number of vectors " << mv.NumVectors() << " in "
927 "the input multivector.");
928 }
929 return Teuchos::rcp (new Epetra_MultiVector (Epetra_DataAccess::View, mv, index.lbound(), index.size()));
930 }
931
937 static Teuchos::RCP<const Epetra_MultiVector>
938 CloneView (const Epetra_MultiVector& mv, const std::vector<int>& index)
939 {
940 const int inNumVecs = GetNumberVecs (mv);
941 const int outNumVecs = index.size();
942
943 // Simple, inexpensive tests of the index vector.
944 TEUCHOS_TEST_FOR_EXCEPTION(outNumVecs == 0, std::invalid_argument,
945 "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
946 "CloneView(mv, index = {}): The output view "
947 "must have at least one column.");
948 if (outNumVecs > inNumVecs)
949 {
950 std::ostringstream os;
951 os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
952 "CloneView(mv, index = {";
953 for (int k = 0; k < outNumVecs - 1; ++k)
954 os << index[k] << ", ";
955 os << index[outNumVecs-1] << "}): There are " << outNumVecs
956 << " indices to view, but only " << inNumVecs << " columns of mv.";
957 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
958 }
959#ifdef TEUCHOS_DEBUG
960 // In debug mode, we perform more expensive tests of the index
961 // vector, to ensure all the elements are in range.
962 // Dereferencing the iterator is valid because index has length
963 // > 0.
964 const int minIndex = *std::min_element (index.begin(), index.end());
965 const int maxIndex = *std::max_element (index.begin(), index.end());
966
967 if (minIndex < 0)
968 {
969 std::ostringstream os;
970 os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
971 "CloneView(mv, index = {";
972 for (int k = 0; k < outNumVecs - 1; ++k)
973 os << index[k] << ", ";
974 os << index[outNumVecs-1] << "}): Indices must be nonnegative, but "
975 "the smallest index " << minIndex << " is negative.";
976 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
977 }
978 if (maxIndex >= inNumVecs)
979 {
980 std::ostringstream os;
981 os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
982 "CloneView(mv, index = {";
983 for (int k = 0; k < outNumVecs - 1; ++k)
984 os << index[k] << ", ";
985 os << index[outNumVecs-1] << "}): Indices must be strictly less than "
986 "the number of vectors " << inNumVecs << " in mv; the largest index "
987 << maxIndex << " is out of bounds.";
988 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
989 }
990#endif // TEUCHOS_DEBUG
991 // Cast to nonconst, because Epetra_MultiVector's constructor
992 // wants a nonconst int array argument. It doesn't actually
993 // change the entries of the array.
994 std::vector<int>& tmpind = const_cast< std::vector<int>& > (index);
995 return Teuchos::rcp (new Epetra_MultiVector (Epetra_DataAccess::View, mv, &tmpind[0], index.size()));
996 }
997
998 static Teuchos::RCP<Epetra_MultiVector>
999 CloneView (const Epetra_MultiVector& mv, const Teuchos::Range1D& index)
1000 {
1001 const bool validRange = index.size() > 0 &&
1002 index.lbound() >= 0 &&
1003 index.ubound() < mv.NumVectors();
1004 if (! validRange)
1005 {
1006 std::ostringstream os;
1007 os << "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
1008 "(mv,index=[" << index.lbound() << ", " << index.ubound()
1009 << "]): ";
1010 TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
1011 os.str() << "Column index range must be nonempty.");
1012 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
1013 os.str() << "Column index range must be nonnegative.");
1014 TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= mv.NumVectors(),
1015 std::invalid_argument,
1016 os.str() << "Column index range must not exceed "
1017 "number of vectors " << mv.NumVectors() << " in "
1018 "the input multivector.");
1019 }
1020 return Teuchos::rcp (new Epetra_MultiVector(Epetra_DataAccess::View, mv, index.lbound(), index.size()));
1021 }
1022
1024
1026
1027
1029 static ptrdiff_t GetGlobalLength( const Epetra_MultiVector& mv )
1030 {
1031 if (mv.Map().GlobalIndicesLongLong())
1032 return static_cast<ptrdiff_t>( mv.GlobalLength64() );
1033 else
1034 return static_cast<ptrdiff_t>( mv.GlobalLength() );
1035 }
1036
1038 static int GetNumberVecs( const Epetra_MultiVector& mv )
1039 { return mv.NumVectors(); }
1040
1041 static bool HasConstantStride( const Epetra_MultiVector& mv )
1042 { return mv.ConstantStride(); }
1044
1046
1047
1050 static void MvTimesMatAddMv( double alpha, const Epetra_MultiVector& A,
1051 const Teuchos::SerialDenseMatrix<int,double>& B,
1052 double beta, Epetra_MultiVector& mv )
1053 {
1054 Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
1055 Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
1056
1057 TEUCHOS_TEST_FOR_EXCEPTION( mv.Multiply( 'N', 'N', alpha, A, B_Pvec, beta )!=0, EpetraMultiVecFailure,
1058 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTimesMatAddMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
1059 }
1060
1063 static void MvAddMv( double alpha, const Epetra_MultiVector& A, double beta, const Epetra_MultiVector& B, Epetra_MultiVector& mv )
1064 {
1065 // epetra mv.Update(alpha,A,beta,B,gamma) will check
1066 // alpha == 0.0
1067 // and
1068 // beta == 0.0
1069 // and based on this will call
1070 // mv.Update(beta,B,gamma)
1071 // or
1072 // mv.Update(alpha,A,gamma)
1073 //
1074 // mv.Update(alpha,A,gamma)
1075 // will then check for one of
1076 // gamma == 0
1077 // or
1078 // gamma == 1
1079 // or
1080 // alpha == 1
1081 // in that order. however, it will not look for the combination
1082 // alpha == 1 and gamma = 0
1083 // which is a common use case when we wish to assign
1084 // mv = A (in which case alpha == 1, beta == gamma == 0)
1085 // or
1086 // mv = B (in which case beta == 1, alpha == gamma == 0)
1087 //
1088 // therefore, we will check for these use cases ourselves
1089 if (beta == 0.0) {
1090 if (alpha == 1.0) {
1091 // assign
1092 mv = A;
1093 }
1094 else {
1095 // single update
1096 TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, 0.0 )!=0, EpetraMultiVecFailure,
1097 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
1098 }
1099 }
1100 else if (alpha == 0.0) {
1101 if (beta == 1.0) {
1102 // assign
1103 mv = B;
1104 }
1105 else {
1106 // single update
1107 TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( beta, B, 0.0 )!=0, EpetraMultiVecFailure,
1108 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
1109 }
1110 }
1111 else {
1112 // double update
1113 TEUCHOS_TEST_FOR_EXCEPTION( mv.Update( alpha, A, beta, B, 0.0 )!=0, EpetraMultiVecFailure,
1114 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
1115 }
1116 }
1117
1120 static void MvTransMv( double alpha, const Epetra_MultiVector& A, const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
1121#ifdef HAVE_ANASAZI_EXPERIMENTAL
1122 , ConjType conj = Anasazi::CONJ
1123#endif
1124 )
1125 {
1126 Epetra_LocalMap LocalMap(B.numRows(), 0, mv.Map().Comm());
1127 Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
1128
1129 TEUCHOS_TEST_FOR_EXCEPTION( B_Pvec.Multiply( 'T', 'N', alpha, A, mv, 0.0 )!=0, EpetraMultiVecFailure,
1130 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvTransMv call to Epetra_MultiVector::Multiply() returned a nonzero value.");
1131 }
1132
1135 static void MvDot( const Epetra_MultiVector& A, const Epetra_MultiVector& B, std::vector<double> &b
1136#ifdef HAVE_ANASAZI_EXPERIMENTAL
1137 , ConjType conj = Anasazi::CONJ
1138#endif
1139 )
1140 {
1141#ifdef TEUCHOS_DEBUG
1142 TEUCHOS_TEST_FOR_EXCEPTION(A.NumVectors() != B.NumVectors(),std::invalid_argument,
1143 "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): A and B must have the same number of vectors.");
1144 TEUCHOS_TEST_FOR_EXCEPTION(b.size() != (unsigned int)A.NumVectors(),std::invalid_argument,
1145 "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvDot(A,B,b): b must have room for all dot products.");
1146#endif
1147 TEUCHOS_TEST_FOR_EXCEPTION( A.Dot( B, &b[0] )!=0, EpetraMultiVecFailure,
1148 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");
1149 }
1150
1152
1154
1158 static void MvNorm( const Epetra_MultiVector& mv, std::vector<double> &normvec )
1159 {
1160#ifdef TEUCHOS_DEBUG
1161 TEUCHOS_TEST_FOR_EXCEPTION((unsigned int)mv.NumVectors() != normvec.size(),std::invalid_argument,
1162 "Anasazi::MultiVecTraits<double,Epetra_MultiVector>::MvNorm(mv,normvec): normvec must be the same size of mv.");
1163#endif
1164 TEUCHOS_TEST_FOR_EXCEPTION( mv.Norm2(&normvec[0])!=0, EpetraMultiVecFailure,
1165 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
1166 }
1167
1169
1171
1172
1174 static void
1175 SetBlock (const Epetra_MultiVector& A,
1176 const std::vector<int>& index,
1177 Epetra_MultiVector& mv)
1178 {
1179 const int inNumVecs = GetNumberVecs (A);
1180 const int outNumVecs = index.size();
1181
1182 // FIXME (mfh 13 Jan 2011) Belos allows A to have more columns
1183 // than index.size(), in which case we just take the first
1184 // index.size() columns of A. Anasazi requires that A have the
1185 // same number of columns as index.size(). Changing Anasazi's
1186 // behavior should not break existing Anasazi solvers, but the
1187 // tests need to be done.
1188 if (inNumVecs != outNumVecs)
1189 {
1190 std::ostringstream os;
1191 os << "Belos::MultiVecTraits<double,Epetra_MultiVector>::"
1192 "SetBlock(A, mv, index = {";
1193 if (outNumVecs > 0)
1194 {
1195 for (int k = 0; k < outNumVecs - 1; ++k)
1196 os << index[k] << ", ";
1197 os << index[outNumVecs-1];
1198 }
1199 os << "}): A has only " << inNumVecs << " columns, but there are "
1200 << outNumVecs << " indices in the index vector.";
1201 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
1202 }
1203 // Make a view of the columns of mv indicated by the index std::vector.
1204 Teuchos::RCP<Epetra_MultiVector> mv_view = CloneViewNonConst (mv, index);
1205
1206 // View of columns [0, outNumVecs-1] of the source multivector A.
1207 // If A has fewer columns than mv_view, then create a view of
1208 // the first outNumVecs columns of A.
1209 Teuchos::RCP<const Epetra_MultiVector> A_view;
1210 if (outNumVecs == inNumVecs)
1211 A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
1212 else
1213 A_view = CloneView (A, Teuchos::Range1D(0, outNumVecs - 1));
1214
1215 // Assignment calls Epetra_MultiVector::Assign(), which deeply
1216 // copies the data directly, ignoring the underlying
1217 // Epetra_Map(s). If A and mv don't have the same data
1218 // distribution (Epetra_Map), this may result in incorrect or
1219 // undefined behavior. Epetra_MultiVector::Update() also
1220 // ignores the Epetra_Maps, so we might as well just use the
1221 // (perhaps slightly cheaper) Assign() method via operator=().
1222 *mv_view = *A_view;
1223 }
1224
1225 static void
1226 SetBlock (const Epetra_MultiVector& A,
1227 const Teuchos::Range1D& index,
1228 Epetra_MultiVector& mv)
1229 {
1230 const int numColsA = A.NumVectors();
1231 const int numColsMv = mv.NumVectors();
1232 // 'index' indexes into mv; it's the index set of the target.
1233 const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
1234 // We can't take more columns out of A than A has.
1235 const bool validSource = index.size() <= numColsA;
1236
1237 if (! validIndex || ! validSource)
1238 {
1239 std::ostringstream os;
1240 os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::SetBlock"
1241 "(A, index=[" << index.lbound() << ", " << index.ubound() << "], "
1242 "mv): ";
1243 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
1244 os.str() << "Range lower bound must be nonnegative.");
1245 TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
1246 os.str() << "Range upper bound must be less than "
1247 "the number of columns " << numColsA << " in the "
1248 "'mv' output argument.");
1249 TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
1250 os.str() << "Range must have no more elements than"
1251 " the number of columns " << numColsA << " in the "
1252 "'A' input argument.");
1253 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
1254 }
1255
1256 // View of columns [index.lbound(), index.ubound()] of the
1257 // target multivector mv. We avoid view creation overhead by
1258 // only creating a view if the index range is different than [0,
1259 // (# columns in mv) - 1].
1260 Teuchos::RCP<Epetra_MultiVector> mv_view;
1261 if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
1262 mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
1263 else
1264 mv_view = CloneViewNonConst (mv, index);
1265
1266 // View of columns [0, index.size()-1] of the source multivector
1267 // A. If A has fewer columns than mv_view, then create a view
1268 // of the first index.size() columns of A.
1269 Teuchos::RCP<const Epetra_MultiVector> A_view;
1270 if (index.size() == numColsA)
1271 A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
1272 else
1273 A_view = CloneView (A, Teuchos::Range1D(0, index.size()-1));
1274
1275 // Assignment calls Epetra_MultiVector::Assign(), which deeply
1276 // copies the data directly, ignoring the underlying
1277 // Epetra_Map(s). If A and mv don't have the same data
1278 // distribution (Epetra_Map), this may result in incorrect or
1279 // undefined behavior. Epetra_MultiVector::Update() also
1280 // ignores the Epetra_Maps, so we might as well just use the
1281 // (perhaps slightly cheaper) Assign() method via operator=().
1282 *mv_view = *A_view;
1283 }
1284
1285 static void
1286 Assign (const Epetra_MultiVector& A,
1287 Epetra_MultiVector& mv)
1288 {
1289 const int numColsA = GetNumberVecs (A);
1290 const int numColsMv = GetNumberVecs (mv);
1291 if (numColsA > numColsMv)
1292 {
1293 std::ostringstream os;
1294 os << "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::Assign"
1295 "(A, mv): ";
1296 TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
1297 os.str() << "Input multivector 'A' has "
1298 << numColsA << " columns, but output multivector "
1299 "'mv' has only " << numColsMv << " columns.");
1300 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
1301 }
1302 // View of the first [0, numColsA-1] columns of mv.
1303 Teuchos::RCP<Epetra_MultiVector> mv_view;
1304 if (numColsMv == numColsA)
1305 mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
1306 else // numColsMv > numColsA
1307 mv_view = CloneView (mv, Teuchos::Range1D(0, numColsA - 1));
1308
1309 // Assignment calls Epetra_MultiVector::Assign(), which deeply
1310 // copies the data directly, ignoring the underlying
1311 // Epetra_Map(s). If A and mv don't have the same data
1312 // distribution (Epetra_Map), this may result in incorrect or
1313 // undefined behavior. Epetra_MultiVector::Update() also
1314 // ignores the Epetra_Maps, so we might as well just use the
1315 // (perhaps slightly cheaper) Assign() method via operator=().
1316 *mv_view = A;
1317 }
1318
1321 static void MvScale ( Epetra_MultiVector& mv, double alpha )
1322 {
1323 TEUCHOS_TEST_FOR_EXCEPTION( mv.Scale( alpha )!=0, EpetraMultiVecFailure,
1324 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value.");
1325 }
1326
1329 static void MvScale ( Epetra_MultiVector& mv, const std::vector<double>& alpha )
1330 {
1331 // Check to make sure the vector is as long as the multivector has columns.
1332 int numvecs = mv.NumVectors();
1333#ifdef TEUCHOS_DEBUG
1334 TEUCHOS_TEST_FOR_EXCEPTION( alpha.size() != (unsigned int)numvecs, std::invalid_argument,
1335 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale(mv,vector alpha): size of alpha inconsistent with number of vectors in mv.")
1336#endif
1337 for (int i=0; i<numvecs; i++) {
1338 TEUCHOS_TEST_FOR_EXCEPTION( mv(i)->Scale(alpha[i])!=0, EpetraMultiVecFailure,
1339 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
1340 }
1341 }
1342
1345 static void MvRandom( Epetra_MultiVector& mv )
1346 {
1347 TEUCHOS_TEST_FOR_EXCEPTION( mv.Random()!=0, EpetraMultiVecFailure,
1348 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
1349 }
1350
1353 static void MvInit( Epetra_MultiVector& mv, double alpha = Teuchos::ScalarTraits<double>::zero() )
1354 {
1355 TEUCHOS_TEST_FOR_EXCEPTION( mv.PutScalar(alpha)!=0, EpetraMultiVecFailure,
1356 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
1357 }
1358
1360
1362
1363
1366 static void MvPrint( const Epetra_MultiVector& mv, std::ostream& os )
1367 { os << mv << std::endl; }
1368
1370
1371#if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
1372# if defined(HAVE_TPETRA_EPETRA)
1378 typedef Epetra::TsqrAdaptor tsqr_adaptor_type;
1379# endif // defined(HAVE_TPETRA_EPETRA)
1380#endif // defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
1381 };
1382
1384 //
1385 // Implementation of the Anasazi::OperatorTraits for Epetra::Operator.
1386 //
1388
1400 template <>
1401 class OperatorTraits < double, Epetra_MultiVector, Epetra_Operator >
1402 {
1403 public:
1404
1408 static void Apply ( const Epetra_Operator& Op,
1409 const Epetra_MultiVector& x,
1410 Epetra_MultiVector& y )
1411 {
1412#ifdef TEUCHOS_DEBUG
1413 TEUCHOS_TEST_FOR_EXCEPTION(x.NumVectors() != y.NumVectors(),std::invalid_argument,
1414 "Anasazi::OperatorTraits<double,Epetra_MultiVector,Epetra_Operator>::Apply(Op,x,y): x and y must have the same number of columns.");
1415#endif
1416 int ret = Op.Apply(x,y);
1417 TEUCHOS_TEST_FOR_EXCEPTION(ret != 0, OperatorError,
1418 "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
1419 }
1420
1421 };
1422
1423 template<>
1424 struct OutputStreamTraits<Epetra_Operator>
1425 {
1426 static Teuchos::RCP<Teuchos::FancyOStream>
1427 getOutputStream (const Epetra_Operator& op, int rootRank = 0)
1428 {
1429 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
1430 const Epetra_Comm & comm = op.Comm();
1431
1432 // Select minimum MPI rank as the root rank for printing if provided rank is less than 0.
1433 int myRank = comm.MyPID();
1434 int numProcs = comm.NumProc();
1435 if (rootRank < 0)
1436 {
1437 comm.MinAll( &myRank, &rootRank, 1 );
1438 }
1439
1440 // This is irreversible, but that's only a problem if the input std::ostream
1441 // is actually a Teuchos::FancyOStream on which this method has been
1442 // called before, with a different root rank.
1443 fos->setProcRankAndSize (myRank, numProcs);
1444 fos->setOutputToRootOnly (rootRank);
1445 return fos;
1446 }
1447 };
1448
1449} // end of Anasazi namespace
1450
1451#endif
1452// end of file ANASAZI_EPETRA_ADAPTER_HPP
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Interface for multivectors used by Anasazi' linear solvers.
Templated virtual class for creating operators that can interface with the Anasazi::OperatorTraits cl...
Abstract class definition for Anasazi output stream.
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
Adapter class for creating an operators often used in solving generalized eigenproblems.
bool UseTranspose() const
Returns the current UseTranspose setting [always false for this operator].
bool HasNormInf() const
Returns true if this object can provide an approximate inf-norm [always false for this operator].
int SetUseTranspose(bool)
If set true, the transpose of this operator will be applied [not functional for this operator].
const Epetra_Comm & Comm() const
Returns the Epetra_Comm communicator associated with this operator.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
const char * Label() const
Returns a character string describing the operator.
double NormInf() const
Returns the infinity norm of the global matrix [not functional for this operator].
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
EpetraMultiVecAccessor is an interfaceto allow any Anasazi::MultiVec implementation that is based on ...
virtual ~EpetraMultiVecAccessor()
Destructor.
virtual const Epetra_MultiVector * GetEpetraMultiVec() const
Return the pointer to the Epetra_MultiVector object.
virtual Epetra_MultiVector * GetEpetraMultiVec()
Return the pointer to the Epetra_MultiVector object.
EpetraMultiVecFailure is thrown when a return value from an Epetra call on an Epetra_MultiVector is n...
Basic adapter class for Anasazi::MultiVec that uses Epetra_MultiVector.
virtual ~EpetraMultiVec()
Destructor.
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
void MvRandom()
Fill the vectors in *this with random numbers.
void MvScale(double alpha)
Scale each element of the vectors in *this with alpha.
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 ...
void MvPrint(std::ostream &os) const
Print *this EpetraMultiVec.
const Epetra_MultiVector * GetEpetraMultiVec() const
Return the pointer to the Epetra_MultiVector object.
Epetra_MultiVector * GetEpetraMultiVec()
Return the pointer to the Epetra_MultiVector object.
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
void MvInit(double alpha)
Replace each element of the vectors in *this with alpha.
EpetraOpFailure is thrown when a return value from an Epetra call on an Epetra_Operator is non-zero.
Basic adapter class for Anasazi::Operator that uses Epetra_Operator.
Adapter class for creating a symmetric operator from an Epetra_MultiVector.
Adapter class for creating a symmetric operator from an Epetra_Operator.
const Epetra_Comm & Comm() const
Returns the Epetra_Comm communicator associated with this operator.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
double NormInf() const
Returns the infinity norm of the global matrix [not functional for this operator].
bool UseTranspose() const
Returns the current UseTranspose setting [always false for this operator].
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
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].
int SetUseTranspose(bool)
If set true, the transpose of this operator will be applied [not functional for this operator].
Adapter class for creating a weighted symmetric operator from an Epetra_MultiVector and Epetra_Operat...
Adapter class for creating a weighted operator from an Epetra_MultiVector and Epetra_Operator.
static void MvInit(Epetra_MultiVector &mv, double alpha=Teuchos::ScalarTraits< double >::zero())
Replace each element of the vectors in mv with alpha.
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 void MvAddMv(double alpha, const Epetra_MultiVector &A, double beta, const Epetra_MultiVector &B, Epetra_MultiVector &mv)
Replace mv with .
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 .
static ptrdiff_t GetGlobalLength(const Epetra_MultiVector &mv)
Obtain the vector length of mv.
static void MvTimesMatAddMv(double alpha, const Epetra_MultiVector &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta, Epetra_MultiVector &mv)
Update mv with .
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).
static Teuchos::RCP< Epetra_MultiVector > Clone(const Epetra_MultiVector &mv, const int outNumVecs)
Creates a new empty Epetra_MultiVector containing numVecs columns.
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).
static void MvPrint(const Epetra_MultiVector &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
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).
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...
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 ...
static void MvScale(Epetra_MultiVector &mv, const std::vector< double > &alpha)
Scale each element of the i-th vector in mv with alpha[i].
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.
static void MvRandom(Epetra_MultiVector &mv)
Replace the vectors in mv with random vectors.
static void MvScale(Epetra_MultiVector &mv, double alpha)
Scale each element of the vectors in mv with alpha.
static int GetNumberVecs(const Epetra_MultiVector &mv)
Obtain the number of vectors in mv.
Traits class which defines basic operations on multivectors.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void Assign(const MV &A, MV &mv)
mv := A
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index.
Interface for multivectors used by Anasazi's linear solvers.
Exceptions thrown to signal error in operator application.
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 ...
Virtual base class which defines basic traits for the operator type.
Anasazi's templated virtual class for constructing an operator that can interface with the OperatorTr...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ConjType
Enumerated types used to specify conjugation arguments.