Amesos2: `Did not get the expected number of non-zero vals` for non-0-based maps
Created by: nschloe
@trilinos/amesos2 @jdbooth @srajama1 When running
#include <Teuchos_DefaultComm.hpp>
#include <Tpetra_CrsMatrix.hpp>
#include <Teuchos_VerboseObject.hpp>
#include <Amesos2.hpp>
using SC = double;
using LO = int;
using GO = int;
using MV = Tpetra::MultiVector<SC,LO,GO>;
using OP = Tpetra::CrsMatrix<SC,LO,GO>;
using NormType = MV::mag_type;
int main ( int argc, char *argv[] )
{
Teuchos::GlobalMPISession session(&argc, &argv, NULL);
auto comm = Teuchos::DefaultComm<int>::getComm();
// std::vector<int> ids = {0, 1}; const int base = 0;// no problems
// std::vector<int> ids = {1, 2}; const int base = 1; // no problems
std::vector<int> ids = {1, 2}; const int base = 0;
auto map = Teuchos::rcp(new Tpetra::Map<int>(
Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
ids,
base,
comm
));
// insert indices from all processes
//auto graph = Teuchos::rcp(new Tpetra::CrsGraph<int, int>(map, map, 0));
auto graph = Teuchos::rcp(new Tpetra::CrsGraph<int, int>(map, 0));
for (int k = 0; k < 2; k++) {
graph->insertGlobalIndices(ids[k], Teuchos::tuple(ids[k]));
}
graph->fillComplete();
auto x = Teuchos::rcp(new Tpetra::Vector<double,int>(map));
auto b = Teuchos::rcp(new Tpetra::Vector<double,int>(map));
auto A = Teuchos::rcp(new Tpetra::CrsMatrix<double,int>(graph));
// fill matrix, rhs
for (int k = 0; k < 2; k++) {
A->replaceGlobalValues(ids[k], Teuchos::tuple(ids[k]), Teuchos::tuple(1.0));
b->replaceGlobalValue(ids[k], 3.14);
}
// let's see what we did
auto out = Teuchos::VerboseObjectBase::getDefaultOStream();
out->setOutputToRootOnly(-1);
A->describe(*out, Teuchos::VERB_EXTREME);
b->describe(*out, Teuchos::VERB_EXTREME);
// solve
auto solver = Amesos2::create<OP,MV>("Superlu", A, x, b);
solver->symbolicFactorization().numericFactorization().solve();
// show the solution
std::cout << std::endl;
x->describe(*out, Teuchos::VERB_EXTREME);
return EXIT_SUCCESS;
}
(even on only one process), one gets
terminate called after throwing an instance of 'std::runtime_error'
what(): /build/trilinos-LuoN8k/trilinos-12.5~20160107030412/packages/amesos2/src/Amesos2_Superlu_def.hpp:716:
Throw number = 2
Throw test that evaluated to true: nnz_ret != as<int>(this->globalNumNonZeros_)
Did not get the expected number of non-zero vals
The reason is the global index offset 0
together with the IDs {1,2}
. If the offset is set to 1
or the IDs to {0,1}
, everything runs fine. It's probably getting tricky with non-contiguous maps.
Same as https://software.sandia.gov/bugzilla/show_bug.cgi?id=6340 (@aprokop).