Belos: BLAS error in GMRES - matrix with no entries on a rank
Created by: chochmuth
Hello,
I'm encountering a problem using Belos GMRES (Flexible GMRES too) with a matrix, that has an empty row, or several.
cblas_dgemm BLAS routine throws the following error:
lda must be >= MAX(K,1): lda=0 K=0ldb must be >= MAX(K,1): ldb=0 K=0BLAS error: Parameter number 9 passed to cblas_dgemm had an invalid value
With parameter number 9(lda) standing for the number of rows.
Is this a desired/necessary behaviour?
You can find an example for a matrix with entries on the first mpi-rank only below. It runs fine if executed on one mpi-rank but with more than one mpi-rank i get the above error.
#include <mpi.h>
#include <Teuchos_XMLParameterListCoreHelpers.hpp>
#include <Xpetra_Operator.hpp>
#include <Xpetra_Matrix_fwd.hpp>
#include <Xpetra_MatrixMatrix.hpp>
#include "Xpetra_MapFactory.hpp"
#include "Xpetra_MultiVectorFactory.hpp"
#include "Xpetra_MatrixFactory.hpp"
#include <BelosOperatorT.hpp>
#include <BelosXpetraAdapter.hpp>
#include <BelosSolverFactory.hpp>
#include <Xpetra_Export.hpp>
typedef unsigned UN;
typedef double SC;
typedef int LO;
typedef int GO;
typedef KokkosClassic::DefaultNode::DefaultNodeType EpetraNode;
typedef EpetraNode NO;
using namespace std;
using namespace Teuchos;
using namespace Xpetra;
using namespace Belos;
int main(int argc, char *argv[]) {
MPI_Init(&argc,&argv);
{
RCP<const Teuchos::Comm<int> > TeuchosComm = rcp(new MpiComm<int> (MPI_COMM_WORLD));
RCP<Map<LO,GO,NO> > UniqueMap;
RCP<Matrix<SC,LO,GO,NO> > K;
int NumMyElements = 0;
if (TeuchosComm->getRank()==0) {
NumMyElements = 10;
}
Array<GO> uniqueMapArray(NumMyElements);
for (LO i=0; i<uniqueMapArray.size(); i++) {
uniqueMapArray[i] = i;
}
UniqueMap = MapFactory<LO,GO,NO>::Build(UseTpetra,-1,uniqueMapArray(),0,TeuchosComm);
K = MatrixFactory<SC,LO,GO,NO>::Build(UniqueMap,1);
for (LO i=0; i<UniqueMap->getNodeNumElements(); i++) {
LO numEntries = 10-i;
Array<GO> indicesArray(numEntries);
Array<SC> valuesArray(numEntries);
for (LO k=0; k<numEntries; k++) {
indicesArray[k] = k+i;
valuesArray[k] = 1;
}
K->insertGlobalValues(UniqueMap->getGlobalElement(i),indicesArray(),valuesArray());
}
K->fillComplete();
//Solve with GMRES
RCP<MultiVector<SC,LO,GO,NO> > xSolution = MultiVectorFactory<SC,LO,GO,NO>::Build(UniqueMap,1);
RCP<MultiVector<SC,LO,GO,NO> > xRightHandSide = MultiVectorFactory<SC,LO,GO,NO>::Build(UniqueMap,1);
xSolution->putScalar(0.0);
xRightHandSide->putScalar(1.0);
RCP<OperatorT<MultiVector<SC,LO,GO,NO> > > OpK = rcp(new XpetraOp<SC, LO, GO, NO>(K));
RCP<LinearProblem<SC,MultiVector<SC,LO,GO,NO>,OperatorT<MultiVector<SC,LO,GO,NO> > > > belosLinearProblem(new LinearProblem<SC,MultiVector<SC,LO,GO,NO>,OperatorT<MultiVector<SC,LO,GO,NO> > >(OpK,xSolution,xRightHandSide));
SolverFactory<SC,MultiVector<SC,LO,GO,NO>,OperatorT<MultiVector<SC,LO,GO,NO> > > belosFactory;
RCP<ParameterList> solverParameterList = rcp(new ParameterList());
solverParameterList->set("Convergence Tolerance",1.e-12);
solverParameterList->set("Verbosity",47);
solverParameterList->set("Output Frequency",1);
solverParameterList->set("Output Style",1);
RCP<SolverManager<SC,MultiVector<SC,LO,GO,NO>,OperatorT<MultiVector<SC,LO,GO,NO> > > > belosSoverManager = belosFactory.create("GMRES",solverParameterList);
belosSoverManager->setProblem(belosLinearProblem);
belosLinearProblem->setProblem(xSolution,xRightHandSide);
belosSoverManager->solve();
}
MPI_Finalize();
return(EXIT_SUCCESS);
}
@trilinos/belos Belos