From 9007b0d8b7432369fa736438e152d624d84bf317 Mon Sep 17 00:00:00 2001 From: maxfirmbach Date: Fri, 27 Dec 2024 11:50:43 +0100 Subject: [PATCH] Add Epetra noncontigous gid --- ...sos2_EpetraCrsMatrix_MatrixAdapter_def.hpp | 36 ++++- ...traRowMatrix_AbstractMatrixAdapter_def.hpp | 8 +- .../amesos2/test/solvers/KLU2_UnitTests.cpp | 149 ++++++++++++++++++ 3 files changed, 184 insertions(+), 9 deletions(-) diff --git a/packages/amesos2/src/Amesos2_EpetraCrsMatrix_MatrixAdapter_def.hpp b/packages/amesos2/src/Amesos2_EpetraCrsMatrix_MatrixAdapter_def.hpp index 7dc7cd4babe5..d256ab9864d5 100644 --- a/packages/amesos2/src/Amesos2_EpetraCrsMatrix_MatrixAdapter_def.hpp +++ b/packages/amesos2/src/Amesos2_EpetraCrsMatrix_MatrixAdapter_def.hpp @@ -26,7 +26,7 @@ namespace Amesos2 { {} Teuchos::RCP > - ConcreteMatrixAdapter::get_impl(const Teuchos::Ptr > map, EDistribution /* distribution */) const + ConcreteMatrixAdapter::get_impl(const Teuchos::Ptr > map, EDistribution distribution) const { using Teuchos::as; using Teuchos::rcp; @@ -36,13 +36,43 @@ namespace Amesos2 { o_map = rcpFromRef(this->mat_->RowMap()); t_map = Util::tpetra_map_to_epetra_map(*map); - int maxRowNNZ = 0; // this->getMaxRowNNZ(); + const int maxRowNNZ = 0; RCP 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 contiguousRowMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) ); + RCP contiguousColMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) ); + RCP contiguousDomainMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) ); + RCP contiguousRangeMap = rcp( new Epetra_Map(global_num_contiguous_entries, local_num_contiguous_entries, 0, (t_mat->Comm() ) ) ); + + RCP 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 (contiguous_t_mat)); + } return( rcp(new ConcreteMatrixAdapter(t_mat)) ); } diff --git a/packages/amesos2/src/Amesos2_EpetraRowMatrix_AbstractMatrixAdapter_def.hpp b/packages/amesos2/src/Amesos2_EpetraRowMatrix_AbstractMatrixAdapter_def.hpp index de9a9db8cde5..34a35f2eb4ed 100644 --- a/packages/amesos2/src/Amesos2_EpetraRowMatrix_AbstractMatrixAdapter_def.hpp +++ b/packages/amesos2/src/Amesos2_EpetraRowMatrix_AbstractMatrixAdapter_def.hpp @@ -274,12 +274,8 @@ namespace Amesos2 { MatrixTraits::node_t> > AbstractConcreteMatrixAdapter::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(map) ); } template diff --git a/packages/amesos2/test/solvers/KLU2_UnitTests.cpp b/packages/amesos2/test/solvers/KLU2_UnitTests.cpp index 700c3cb997a3..7025a9c2db40 100644 --- a/packages/amesos2/test/solvers/KLU2_UnitTests.cpp +++ b/packages/amesos2/test/solvers/KLU2_UnitTests.cpp @@ -28,6 +28,18 @@ #include "Amesos2.hpp" #include "Amesos2_Meta.hpp" +#if defined(HAVE_AMESOS2_EPETRA) +#ifdef HAVE_MPI +#include +#include +#else +#include +#endif +#include +#include +#include +#endif + namespace { using std::cout; @@ -596,6 +608,138 @@ 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 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 > 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 elementList(numLocalEntries); + for ( LO k = 0; k < numLocalEntries; ++k ) { + elementList[k] = myRank + k*numProcs + 4*myRank; + } + + RCP< const map_type > map = rcp( new Epetra_Map(numGlobalEntries, numLocalEntries, elementList.data(), 0, Epetra_MpiComm(MPI_COMM_WORLD) )); + RCP 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(7,-3,-1).data(),tuple(0,4,7).data()); + A->InsertGlobalValues(2,2,tuple(2,8).data(),tuple(0,2).data()); + A->InsertGlobalValues(4,1,tuple(1).data(),tuple(4).data()); + } else { + A->InsertGlobalValues(5,2,tuple(-3,5).data(),tuple(0,5).data()); + A->InsertGlobalValues(7,2,tuple(-1,4).data(),tuple(2,7).data()); + A->InsertGlobalValues(9,2,tuple(-2,6).data(),tuple(5,9).data()); + } + A->FillComplete(); + + // Create X with random values + RCP X = rcp(new MV(*map, numVectors)); + //X->setObjectLabel("X"); + X->Random(); + + /* Create B, use same GIDs + * + * Use RHS: + * + * [[-7] + * [18] + * [ 3] + * [17] + * [18] + * [28]] + */ + RCP 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 > solver = Amesos2::create("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 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 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 */ @@ -824,6 +968,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 )