46#ifndef ANASAZI_EPETRA_ADAPTER_HPP
47#define ANASAZI_EPETRA_ADAPTER_HPP
50#include "Anasaziepetra_DLLExportMacro.h"
56#include "Teuchos_Assert.hpp"
57#include "Teuchos_SerialDenseMatrix.hpp"
58#include "Teuchos_FancyOStream.hpp"
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"
67#if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
68# include <Tpetra_ConfigDefs.hpp>
69# if defined(HAVE_TPETRA_EPETRA)
70# include <Epetra_TsqrAdaptor.hpp>
153 EpetraMultiVec(
const Epetra_BlockMap& Map_in,
double * array,
const int numvecs,
const int stride=0);
162 EpetraMultiVec(Epetra_DataAccess CV,
const Epetra_MultiVector& P_vec,
const std::vector<int>& index);
217 if ( Map().GlobalIndicesLongLong() )
218 return static_cast<ptrdiff_t
>( GlobalLength64() );
220 return static_cast<ptrdiff_t
>( GlobalLength() );
234 const Teuchos::SerialDenseMatrix<int,double>& B,
244 void MvTransMv (
double alpha,
const MultiVec<double>& A, Teuchos::SerialDenseMatrix<int,double>& B
245#ifdef HAVE_ANASAZI_EXPERIMENTAL
253#ifdef HAVE_ANASAZI_EXPERIMENTAL
262 "Anasazi::EpetraMultiVec::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
267 void MvScale (
const std::vector<double>& alpha );
276 void MvNorm ( std::vector<double> & normvec )
const {
277 if (((
int)normvec.size() >= GetNumberVecs()) ) {
279 "Anasazi::EpetraMultiVec::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
297 "Anasazi::EpetraMultiVec::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
304 "Anasazi::EpetraMultiVec::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
323 void MvPrint( std::ostream& os )
const { os << *
this << std::endl; };
348 EpetraOp(
const Teuchos::RCP<Epetra_Operator> &Op );
368#pragma warning(disable:4251)
370 Teuchos::RCP<Epetra_Operator> Epetra_Op;
396 class ANASAZIEPETRA_LIB_DLL_EXPORT
EpetraGenOp :
public virtual Operator<double>,
public virtual Epetra_Operator {
402 EpetraGenOp(
const Teuchos::RCP<Epetra_Operator> &AOp,
403 const Teuchos::RCP<Epetra_Operator> &MOp,
404 bool isAInverse =
true );
417 int Apply(
const Epetra_MultiVector &X, Epetra_MultiVector &Y)
const;
422 int ApplyInverse(
const Epetra_MultiVector &X, Epetra_MultiVector &Y)
const;
425 const char*
Label()
const {
return "Epetra_Operator applying A^{-1}M"; };
440 const Epetra_Comm&
Comm()
const {
return Epetra_AOp->Comm(); };
455#pragma warning(disable:4251)
457 Teuchos::RCP<Epetra_Operator> Epetra_AOp;
458 Teuchos::RCP<Epetra_Operator> Epetra_MOp;
482 class ANASAZIEPETRA_LIB_DLL_EXPORT
EpetraSymOp :
public virtual Operator<double>,
public virtual Epetra_Operator {
487 EpetraSymOp(
const Teuchos::RCP<Epetra_Operator> &Op,
bool isTrans =
false );
500 int Apply(
const Epetra_MultiVector &X, Epetra_MultiVector &Y)
const;
506 int ApplyInverse(
const Epetra_MultiVector &X, Epetra_MultiVector &Y)
const;
509 const char*
Label()
const {
return "Epetra_Operator applying A^TA or AA^T"; };
524 const Epetra_Comm&
Comm()
const {
return Epetra_Op->Comm(); };
537#pragma warning(disable:4251)
539 Teuchos::RCP<Epetra_Operator> Epetra_Op;
571 EpetraSymMVOp(
const Teuchos::RCP<const Epetra_MultiVector> &MV,
572 bool isTrans =
false );
588#pragma warning(disable:4251)
590 Teuchos::RCP<const Epetra_MultiVector> Epetra_MV;
591 Teuchos::RCP<const Epetra_Map> MV_localmap;
592 Teuchos::RCP<const Epetra_BlockMap> MV_blockmap;
622 const Teuchos::RCP<Epetra_Operator> &OP );
637#pragma warning(disable:4251)
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;
671 const Teuchos::RCP<Epetra_Operator> &OP );
686#pragma warning(disable:4251)
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;
727 static Teuchos::RCP<Epetra_MultiVector>
728 Clone (
const Epetra_MultiVector& mv,
const int outNumVecs)
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.");
738 return Teuchos::rcp (
new Epetra_MultiVector (mv.Map(), outNumVecs));
745 static Teuchos::RCP<Epetra_MultiVector>
748 return Teuchos::rcp (
new Epetra_MultiVector (mv));
756 static Teuchos::RCP<Epetra_MultiVector>
757 CloneCopy (
const Epetra_MultiVector& mv,
const std::vector<int>& index)
760 const int outNumVecs = index.size();
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"
767 if (outNumVecs > inNumVecs)
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());
783 const int minIndex = *std::min_element (index.begin(), index.end());
784 const int maxIndex = *std::max_element (index.begin(), index.end());
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());
797 if (maxIndex >= inNumVecs)
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());
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()));
817 static Teuchos::RCP<Epetra_MultiVector>
818 CloneCopy (
const Epetra_MultiVector& mv,
const Teuchos::Range1D& index)
821 const int outNumVecs = index.size();
822 const bool validRange = outNumVecs > 0 && index.lbound() >= 0 &&
823 index.ubound() < inNumVecs;
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.");
838 return Teuchos::rcp (
new Epetra_MultiVector (Epetra_DataAccess::Copy, mv, index.lbound(), index.size()));
846 static Teuchos::RCP<Epetra_MultiVector>
850 const int outNumVecs = index.size();
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)
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());
873 const int minIndex = *std::min_element (index.begin(), index.end());
874 const int maxIndex = *std::max_element (index.begin(), index.end());
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());
887 if (maxIndex >= inNumVecs)
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());
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()));
907 static Teuchos::RCP<Epetra_MultiVector>
910 const bool validRange = index.size() > 0 &&
911 index.lbound() >= 0 &&
912 index.ubound() < mv.NumVectors();
915 std::ostringstream os;
916 os <<
"Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
917 "NonConst(mv,index=[" << index.lbound() <<
", " << index.ubound()
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.");
929 return Teuchos::rcp (
new Epetra_MultiVector (Epetra_DataAccess::View, mv, index.lbound(), index.size()));
937 static Teuchos::RCP<const Epetra_MultiVector>
938 CloneView (
const Epetra_MultiVector& mv,
const std::vector<int>& index)
941 const int outNumVecs = index.size();
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)
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());
964 const int minIndex = *std::min_element (index.begin(), index.end());
965 const int maxIndex = *std::max_element (index.begin(), index.end());
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());
978 if (maxIndex >= inNumVecs)
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());
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()));
998 static Teuchos::RCP<Epetra_MultiVector>
999 CloneView (
const Epetra_MultiVector& mv,
const Teuchos::Range1D& index)
1001 const bool validRange = index.size() > 0 &&
1002 index.lbound() >= 0 &&
1003 index.ubound() < mv.NumVectors();
1006 std::ostringstream os;
1007 os <<
"Anasazi::MultiVecTraits<double,Epetra_MultiVector>::CloneView"
1008 "(mv,index=[" << index.lbound() <<
", " << index.ubound()
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.");
1020 return Teuchos::rcp (
new Epetra_MultiVector(Epetra_DataAccess::View, mv, index.lbound(), index.size()));
1031 if (mv.Map().GlobalIndicesLongLong())
1032 return static_cast<ptrdiff_t
>( mv.GlobalLength64() );
1034 return static_cast<ptrdiff_t
>( mv.GlobalLength() );
1039 {
return mv.NumVectors(); }
1041 static bool HasConstantStride(
const Epetra_MultiVector& mv )
1042 {
return mv.ConstantStride(); }
1051 const Teuchos::SerialDenseMatrix<int,double>& B,
1052 double beta, Epetra_MultiVector& mv )
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());
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.");
1063 static void MvAddMv(
double alpha,
const Epetra_MultiVector& A,
double beta,
const Epetra_MultiVector& B, Epetra_MultiVector& mv )
1097 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,0.0) returned a nonzero value.");
1100 else if (alpha == 0.0) {
1108 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(beta,B,0.0) returned a nonzero value.");
1114 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvAddMv call to Epetra_MultiVector::Update(alpha,A,beta,B,0.0) returned a nonzero value.");
1120 static void MvTransMv(
double alpha,
const Epetra_MultiVector& A,
const Epetra_MultiVector& mv, Teuchos::SerialDenseMatrix<int,double>& B
1121#ifdef HAVE_ANASAZI_EXPERIMENTAL
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());
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.");
1135 static void MvDot(
const Epetra_MultiVector& A,
const Epetra_MultiVector& B, std::vector<double> &b
1136#ifdef HAVE_ANASAZI_EXPERIMENTAL
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.");
1148 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvDot(A,B,b) call to Epetra_MultiVector::Dot() returned a nonzero value.");
1158 static void MvNorm(
const Epetra_MultiVector& mv, std::vector<double> &normvec )
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.");
1165 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvNorm call to Epetra_MultiVector::Norm2() returned a nonzero value.");
1176 const std::vector<int>& index,
1177 Epetra_MultiVector& mv)
1180 const int outNumVecs = index.size();
1188 if (inNumVecs != outNumVecs)
1190 std::ostringstream os;
1191 os <<
"Belos::MultiVecTraits<double,Epetra_MultiVector>::"
1192 "SetBlock(A, mv, index = {";
1195 for (
int k = 0; k < outNumVecs - 1; ++k)
1196 os << index[k] <<
", ";
1197 os << index[outNumVecs-1];
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());
1209 Teuchos::RCP<const Epetra_MultiVector> A_view;
1210 if (outNumVecs == inNumVecs)
1211 A_view = Teuchos::rcpFromRef (A);
1213 A_view =
CloneView (A, Teuchos::Range1D(0, outNumVecs - 1));
1226 SetBlock (
const Epetra_MultiVector& A,
1227 const Teuchos::Range1D& index,
1228 Epetra_MultiVector& mv)
1230 const int numColsA = A.NumVectors();
1231 const int numColsMv = mv.NumVectors();
1233 const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
1235 const bool validSource = index.size() <= numColsA;
1237 if (! validIndex || ! validSource)
1239 std::ostringstream os;
1240 os <<
"Anasazi::MultiVecTraits<double, Epetra_MultiVector>::SetBlock"
1241 "(A, index=[" << index.lbound() <<
", " << index.ubound() <<
"], "
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!");
1260 Teuchos::RCP<Epetra_MultiVector> mv_view;
1261 if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
1262 mv_view = Teuchos::rcpFromRef (mv);
1269 Teuchos::RCP<const Epetra_MultiVector> A_view;
1270 if (index.size() == numColsA)
1271 A_view = Teuchos::rcpFromRef (A);
1273 A_view =
CloneView (A, Teuchos::Range1D(0, index.size()-1));
1286 Assign (
const Epetra_MultiVector& A,
1287 Epetra_MultiVector& mv)
1291 if (numColsA > numColsMv)
1293 std::ostringstream os;
1294 os <<
"Anasazi::MultiVecTraits<double, Epetra_MultiVector>::Assign"
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!");
1303 Teuchos::RCP<Epetra_MultiVector> mv_view;
1304 if (numColsMv == numColsA)
1305 mv_view = Teuchos::rcpFromRef (mv);
1307 mv_view =
CloneView (mv, Teuchos::Range1D(0, numColsA - 1));
1321 static void MvScale ( Epetra_MultiVector& mv,
double alpha )
1324 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale(mv,double alpha) returned a nonzero value.");
1329 static void MvScale ( Epetra_MultiVector& mv,
const std::vector<double>& alpha )
1332 int numvecs = mv.NumVectors();
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.")
1337 for (
int i=0; i<numvecs; i++) {
1339 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvScale call to Epetra_MultiVector::Scale() returned a nonzero value.");
1348 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvRandom call to Epetra_MultiVector::Random() returned a nonzero value.");
1353 static void MvInit( Epetra_MultiVector& mv,
double alpha = Teuchos::ScalarTraits<double>::zero() )
1356 "Anasazi::MultiVecTraits<double, Epetra_MultiVector>::MvInit call to Epetra_MultiVector::PutScalar() returned a nonzero value.");
1366 static void MvPrint(
const Epetra_MultiVector& mv, std::ostream& os )
1367 { os << mv << std::endl; }
1371#if defined(HAVE_ANASAZI_TPETRA) && defined(HAVE_ANASAZI_TSQR)
1372# if defined(HAVE_TPETRA_EPETRA)
1378 typedef Epetra::TsqrAdaptor tsqr_adaptor_type;
1408 static void Apply (
const Epetra_Operator& Op,
1409 const Epetra_MultiVector& x,
1410 Epetra_MultiVector& y )
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.");
1416 int ret = Op.Apply(x,y);
1418 "Anasazi::OperatorTraits<double,Epetra_Multivector,Epetra_Operator>::Apply(): Error in Epetra_Operator::Apply(). Code " << ret);
1424 struct OutputStreamTraits<Epetra_Operator>
1426 static Teuchos::RCP<Teuchos::FancyOStream>
1427 getOutputStream (
const Epetra_Operator& op,
int rootRank = 0)
1429 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
1430 const Epetra_Comm & comm = op.Comm();
1433 int myRank = comm.MyPID();
1434 int numProcs = comm.NumProc();
1437 comm.MinAll( &myRank, &rootRank, 1 );
1443 fos->setProcRankAndSize (myRank, numProcs);
1444 fos->setOutputToRootOnly (rootRank);
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.
~EpetraSymMVOp()
Destructor.
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...
~EpetraW2SymMVOp()
Destructor.
Adapter class for creating a weighted operator from an Epetra_MultiVector and Epetra_Operator.
~EpetraWSymMVOp()
Destructor.
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.