GitLab

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