Skip to content

Commit

Permalink
#67: First test build but raised errors to fix
Browse files Browse the repository at this point in the history
  • Loading branch information
antoine.meyer1 committed Aug 17, 2023
1 parent 2522bc0 commit 43b27c5
Showing 1 changed file with 72 additions and 86 deletions.
158 changes: 72 additions & 86 deletions packages/teko/tests/unit_tests/tStridedTpetraOperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,6 @@

// Thyra includes
#include "Thyra_VectorStdOps.hpp"
//#include "Thyra_EpetraThyraWrappers.hpp" // TO REMOVE
//#include "Thyra_EpetraLinearOp.hpp" // TO REMOVE
//#include "Thyra_EpetraOperatorWrapper.hpp" // TO REMOVE
#include "Thyra_DefaultBlockedLinearOp.hpp"
#include "Thyra_ProductVectorBase.hpp"
#include "Thyra_SpmdVectorSpaceBase.hpp"
Expand All @@ -72,45 +69,41 @@
// TriUtils includes
//#include "Trilinos_Util_CrsMatrixGallery.h" // TO REMOVE EPETRA

//#include "Teko_StridedEpetraOperator.hpp" // TO REMOVE
// Teko includes
#include "Teko_StridedTpetraOperator.hpp"
#include "Teko_Utilities.hpp"

#define SS_ECHO(ops) { std::stringstream ss; ss << ops; ECHO(ss.str()); };

// Teuchos using
using Teuchos::Comm;
using Teuchos::null;
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::rcp_dynamic_cast;

// Thyra using
using Thyra::VectorBase;
using Thyra::LinearOpBase;
using Thyra::createMember;
using Thyra::LinearOpTester;

// Tpetra using
using ST = typename Tpetra::Vector<double>::scalar_type;
using LO = typename Tpetra::Vector<>::local_ordinal_type;
using GO = typename Tpetra::Vector<>::global_ordinal_type;
using NT = typename Tpetra::Vector<>::node_type;
using tcrsmatrix_t = Tpetra::CrsMatrix<ST,LO,GO,NT>;
using tmap_t = Tpetra::Map<LO,GO,NT>;
using tmultivector_t = Tpetra::MultiVector<ST,LO,GO,NT>;

double tolerance = 1e-14;
ST tolerance = 1e-14;

const RCP<const Comm<int>> GetComm()
{
return Tpetra::getDefaultComm();
}

// TO REMOVE v
TEUCHOS_UNIT_TEST(t_AM_BLALBA, AM_BLALBLA)
{
int max = 20;
int min = 2;
TEST_ASSERT(max >= min);
}
// TO REMOVE ^

TEUCHOS_UNIT_TEST(tStridedTpetraOperator, test_numvars_constr)
{
// communicator
Expand All @@ -123,83 +116,76 @@ TEUCHOS_UNIT_TEST(tStridedTpetraOperator, test_numvars_constr)

// create a big matrix to play with
RCP<tmap_t> FGallery = rcp(new tmap_t(nx * ny, 0, comm));
RCP<tcrsmatrix_t> A = rcp(new tcrsmatrix_t(FGallery, 1), false);
double beforeNorm = A->getFrobeniusNorm(); // AM: before was Epetra_CrsMatrix->NormOne()
RCP<tcrsmatrix_t> A = rcp(new tcrsmatrix_t(FGallery, 1));
ST beforeNorm = A->getFrobeniusNorm();



/*
// create a strided tpetra operator
int vars = 3;
int width = 3;
Epetra_MultiVector x(A->OperatorDomainMap(),width);
Epetra_MultiVector ys(A->OperatorRangeMap(),width);
Epetra_MultiVector y(A->OperatorRangeMap(),width);
Teko::Epetra::StridedEpetraOperator shell(vars,A);
// test the operator against a lot of random vectors
int numtests = 50;
double max = 0.0;
double min = 1.0;
for(int i=0;i<numtests;i++) {
std::vector<double> norm(width);
std::vector<double> rel(width);
x.Random();
shell.Apply(x,y);
A->Apply(x,ys);
Epetra_MultiVector e(y);
e.Update(-1.0,ys,1.0);
e.Norm2(&norm[0]);
// compute relative error
ys.Norm2(&rel[0]);
for(int j=0;j<width;j++) {
max = max>norm[j]/rel[j] ? max : norm[j]/rel[j];
min = min<norm[j]/rel[j] ? min : norm[j]/rel[j];
}
tmultivector_t x(A->getDomainMap(), width);
tmultivector_t ys(A->getRangeMap(), width);
tmultivector_t y(A->getRangeMap(), width);
Teko::TpetraHelpers::StridedTpetraOperator shell(vars, A); // AM: TO FIX

// test the operator against a lot of random vectors
int numtests = 50;
ST max = 0.0;
ST min = 1.0;
for(int i=0; i<numtests; i++) {
std::vector<ST> norm(width);
std::vector<ST> rel(width);
x.randomize(); // AM: TO FIX

shell.apply(x,y); // AM: TO FIX
A->apply(x,ys);

tmultivector_t e(y, Teuchos::Copy);
e.update(-1.0, ys, 1.0);
e.norm2(Teuchos::ArrayView<ST>(norm));

// compute relative error
ys.norm2(Teuchos::ArrayView<ST>(rel));
for(int j=0;j<width;j++) {
max = max>norm[j]/rel[j] ? max : norm[j]/rel[j];
min = min<norm[j]/rel[j] ? min : norm[j]/rel[j];
}
}
TEST_ASSERT(max>=min);
TEST_ASSERT(max<=tolerance);
TEST_ASSERT(max >= min);
TEST_ASSERT(max <= tolerance);

int * indexOffset,* indicies;
double * values;
A->ExtractCrsDataPointers(indexOffset,indicies,values);
for(int i=0;i<A->NumMyNonzeros();i++)
values[i] *= 2.0; // square everything!
double afterNorm = A->NormOne();
TEST_ASSERT(beforeNorm!=afterNorm);
shell.RebuildOps();
// test the operator against a lot of random vectors
numtests = 50;
max = 0.0;
min = 1.0;
for(int i=0;i<numtests;i++) {
std::vector<double> norm(width);
std::vector<double> rel(width);
x.Random();
shell.Apply(x,y);
A->Apply(x,ys);
Epetra_MultiVector e(y);
e.Update(-1.0,ys,1.0);
e.Norm2(&norm[0]);
// compute relative error
ys.Norm2(&rel[0]);
for(int j=0;j<width;j++) {
max = max>norm[j]/rel[j] ? max : norm[j]/rel[j];
min = min<norm[j]/rel[j] ? min : norm[j]/rel[j];
}
}
TEST_ASSERT(max>=min);
TEST_ASSERT(max<=tolerance);
*/
// double everything
A->scale(2.0);

ST afterNorm = A->getFrobeniusNorm();
TEST_ASSERT(beforeNorm != afterNorm);

shell.RebuildOps();

// test the operator against a lot of random vectors
numtests = 50;
max = 0.0;
min = 1.0;
for(int i=0; i<numtests; i++) {
std::vector<ST> norm(width);
std::vector<ST> rel(width);
x.randomize(); // AM: PROBABLY TO FIX

shell.apply(x,y); // AM: PROBABLY TO FIX
A->apply(x,ys);

tmultivector_t e(y, Teuchos::Copy);
e.update(-1.0, ys, 1.0);
e.norm2(Teuchos::ArrayView<ST>(norm));

// compute relative error
ys.norm2(Teuchos::ArrayView<ST>(rel));
for(int j=0; j<width; j++) {
max = max>norm[j]/rel[j] ? max : norm[j]/rel[j];
min = min<norm[j]/rel[j] ? min : norm[j]/rel[j];
}
}
TEST_ASSERT(max >= min);
TEST_ASSERT(max <= tolerance);
}

/*
Expand Down

0 comments on commit 43b27c5

Please sign in to comment.