Skip to content

Commit

Permalink
Add Epetra noncontigous gid
Browse files Browse the repository at this point in the history 
Signed-off-by: maxfirmbach <[email protected]>
  • Loading branch information
@maxfirmbach
maxfirmbach committed Dec 28, 2024
1 parent ece6cfc commit 4b3742f
Show file tree
Hide file tree
Showing 3 changed files with 188 additions and 9 deletions.
36 changes: 33 additions & 3 deletions packages/amesos2/src/Amesos2_EpetraCrsMatrix_MatrixAdapter_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ namespace Amesos2 {
{}

Teuchos::RCP<const MatrixAdapter<Epetra_CrsMatrix> >
ConcreteMatrixAdapter<Epetra_CrsMatrix>::get_impl(const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > map, EDistribution /* distribution */) const
ConcreteMatrixAdapter<Epetra_CrsMatrix>::get_impl(const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > map, EDistribution distribution) const
{
using Teuchos::as;
using Teuchos::rcp;
Expand All @@ -36,13 +36,43 @@ namespace Amesos2 {
o_map = rcpFromRef(this->mat_->RowMap());
t_map = Util::tpetra_map_to_epetra_map<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(*map);

int maxRowNNZ = 0; // this->getMaxRowNNZ();
const int maxRowNNZ = 0;
RCP<Epetra_CrsMatrix> t_mat = rcp(new Epetra_CrsMatrix(Copy, *t_map, maxRowNNZ));

Epetra_Import importer(*t_map, *o_map);
t_mat->Import(*(this->mat_), importer, Insert);
t_mat->FillComplete();

t_mat->FillComplete(); // Must be in local form for later extraction of rows
// Case for non-contiguous GIDs
if ( distribution == CONTIGUOUS_AND_ROOTED ) {

auto myRank = map->getComm()->getRank();

const int global_num_contiguous_entries = t_mat->NumGlobalRows();
const int local_num_contiguous_entries = (myRank == 0) ? t_mat->NumGlobalRows() : 0;

RCP<const Epetra_Map> contiguousRowMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) );
RCP<const Epetra_Map> contiguousColMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) );
RCP<const Epetra_Map> contiguousDomainMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) );
RCP<const Epetra_Map> contiguousRangeMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) );

RCP<Epetra_CrsMatrix> contiguous_t_mat = rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy, *contiguousRowMap, *contiguousColMap, t_mat->MaxNumEntries()) );

// fill local sparse matrix on rank zero
if(myRank == 0) {
int num_entries;
int *indices;
double *values;
for (int row = 0; row < t_mat->NumMyRows(); row++) {
t_mat->ExtractMyRowView(row, num_entries, values, indices);
contiguous_t_mat->InsertMyValues(row, num_entries, values, indices);
}
}

contiguous_t_mat->FillComplete(*contiguousDomainMap, *contiguousRangeMap);

return rcp (new ConcreteMatrixAdapter<Epetra_CrsMatrix> (contiguous_t_mat));
}

return( rcp(new ConcreteMatrixAdapter<Epetra_CrsMatrix>(t_mat)) );
}
Expand Down
8 changes: 2 additions & 6 deletions packages/amesos2/src/Amesos2_EpetraRowMatrix_AbstractMatrixAdapter_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,12 +274,8 @@ namespace Amesos2 {
MatrixTraits<Epetra_RowMatrix>::node_t> >
AbstractConcreteMatrixAdapter<Epetra_RowMatrix, DerivedMat>::getMap_impl() const
{
// Should Map() be used (returns Epetra_BlockMap)
/*
printf("Amesos2_EpetraRowMatrix_AbstractMatrixAdapter: Epetra does not support a 'getMap()' method. Returning rcp(Teuchos::null). \
If your map contains non-contiguous GIDs please use Tpetra instead of Epetra. \n");
*/
return( Teuchos::null );
const Epetra_BlockMap map = this->mat_->Map();
return( Util::epetra_map_to_tpetra_map<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(map) );
}

template <class DerivedMat>
Expand Down
153 changes: 153 additions & 0 deletions packages/amesos2/test/solvers/KLU2_UnitTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,17 @@
#include "Amesos2.hpp"
#include "Amesos2_Meta.hpp"

#if defined(HAVE_AMESOS2_EPETRA)
#ifdef HAVE_MPI
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif
#include <Epetra_Map.h>
#include <Epetra_MultiVector.h>
#include <Epetra_CrsMatrix.h>
#endif

namespace {

using std::cout;
Expand Down Expand Up @@ -596,6 +607,143 @@ namespace {
} // end if numProcs = 2
}

#if defined(HAVE_AMESOS2_EPETRA)
//! @test Test for Epetra non-contiguous GIDs.
TEUCHOS_UNIT_TEST_TEMPLATE_3_DECL( KLU2, NonContigGIDEpetra, SCALAR, LO, GO )
{
typedef Epetra_CrsMatrix MAT;
typedef ScalarTraits<SCALAR> ST;
typedef Epetra_MultiVector MV;
typedef typename ST::magnitudeType Mag;
typedef Epetra_Map map_type;

using Tpetra::global_size_t;
using Teuchos::tuple;
using Teuchos::RCP;
using Teuchos::rcp;
using Scalar = SCALAR;

RCP<const Comm<int> > comm = Tpetra::getDefaultComm();

size_t myRank = comm->getRank();
const global_size_t numProcs = comm->getSize();

// Unit test created for 2 processes
if ( numProcs == 2 ) {

const global_size_t numVectors = 1;
const global_size_t nrows = 6;

const int numGlobalEntries = nrows;
const LO numLocalEntries = nrows / numProcs;

// Create non-contiguous Map
// This example: np 2 leads to GIDS: proc0 - 0,2,4 proc 1 - 1,3,5
Teuchos::Array<int> elementList(numLocalEntries);
for ( LO k = 0; k < numLocalEntries; ++k ) {
elementList[k] = myRank + k*numProcs + 4*myRank;
}

#ifdef HAVE_MPI
const Epetra_MpiComm comm(MPI_COMM_WORLD);
#else
const Epetra_SerialComm comm;
#endif

RCP< const map_type > map = rcp( new Epetra_Map( numGlobalEntries, numLocalEntries, elementList.data(), 0, comm ));
RCP<MAT> A = rcp( new MAT(Epetra_DataAccess::Copy, *map, numLocalEntries) );

/*
* We will solve a system with a known solution, for which we will be using
* the following matrix:
*
* GID 0 2 4 5 7 9
* [ 0 [ 7, 0, -3, 0, -1, 0 ]
* 2 [ 2, 8, 0, 0, 0, 0 ]
* 4 [ 0, 0, 1, 0, 0, 0 ]
* 5 [-3, 0, 0, 5, 0, 0 ]
* 7 [ 0, -1, 0, 0, 4, 0 ]
* 9 [ 0, 0, 0, -2, 0, 6 ] ]
*
*/

// Construct matrix
if( myRank == 0 ){
A->InsertGlobalValues(0,3,tuple<Scalar>(7,-3,-1).data(),tuple<GO>(0,4,7).data());
A->InsertGlobalValues(2,2,tuple<Scalar>(2,8).data(),tuple<GO>(0,2).data());
A->InsertGlobalValues(4,1,tuple<Scalar>(1).data(),tuple<GO>(4).data());
} else {
A->InsertGlobalValues(5,2,tuple<Scalar>(-3,5).data(),tuple<GO>(0,5).data());
A->InsertGlobalValues(7,2,tuple<Scalar>(-1,4).data(),tuple<GO>(2,7).data());
A->InsertGlobalValues(9,2,tuple<Scalar>(-2,6).data(),tuple<GO>(5,9).data());
}
A->FillComplete();

// Create X with random values
RCP<MV> X = rcp(new MV(*map, numVectors));
X->Random();

/* Create B, use same GIDs
*
* Use RHS:
*
* [[-7]
* [18]
* [ 3]
* [17]
* [18]
* [28]]
*/
RCP<MV> B = rcp(new MV(*map, numVectors));
GO rowids[nrows] = {0,2,4,5,7,9};
Scalar data[nrows] = {-7,18,3,17,18,28};
for( global_size_t i = 0; i < nrows; ++i ){
if( B->Map().MyGID(rowids[i]) ){
B->ReplaceGlobalValue(rowids[i],0,data[i]);
}
}

// Create solver interface with Amesos2 factory method
RCP<Amesos2::Solver<MAT,MV> > solver = Amesos2::create<MAT,MV>("KLU2", A, X, B);

// Create a Teuchos::ParameterList to hold solver parameters
Teuchos::ParameterList amesos2_params("Amesos2");
amesos2_params.sublist("KLU2").set("IsContiguous", false, "Are GIDs Contiguous");

solver->setParameters( Teuchos::rcpFromRef(amesos2_params) );

solver->symbolicFactorization().numericFactorization().solve();

/* Check the solution
*
* Should be:
*
* [[1]
* [2]
* [3]
* [4]
* [5]
* [6]]
*/
// Solution Vector for later comparison
RCP<MV> Xhat = rcp(new MV(*map, numVectors));
GO rowids_soln[nrows] = {0,2,4,5,7,9};
Scalar data_soln[nrows] = {1,2,3,4,5,6};
for( global_size_t i = 0; i < nrows; ++i ){
if( Xhat->Map().MyGID(rowids_soln[i]) ){
Xhat->ReplaceGlobalValue(rowids_soln[i],0,data_soln[i]);
}
}

// Check result of solve
Array<Mag> xhatnorms(numVectors), xnorms(numVectors);
Xhat->Norm2(xhatnorms().data());
X->Norm2(xnorms().data());
TEST_COMPARE_FLOATING_ARRAYS( xhatnorms, xnorms, 0.005 );
} // end if numProcs = 2
}
#endif

/*
* Unit Tests for Complex data types
*/
Expand Down Expand Up @@ -824,6 +972,11 @@ namespace {
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( KLU2, SolveTrans, SCALAR, LO, GO ) \
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( KLU2, NonContigGID, SCALAR, LO, GO )

#ifdef HAVE_AMESOS2_EPETRA
#define UNIT_TEST_GROUP_EPETRA( LO, GO, SCALAR) \
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( KLU2, NonContigGIDEpetra, SCALAR, LO, GO )
#endif

#define UNIT_TEST_GROUP_ORDINAL( ORDINAL ) \
UNIT_TEST_GROUP_ORDINAL_ORDINAL( ORDINAL, ORDINAL )

Expand Down

0 comments on commit 4b3742f

Please sign in to comment.