Skip to content
GitLab
Projects Groups Topics Snippets
  • /
  • Help
    • Help
    • Support
    • Community forum
    • Submit feedback
    • Contribute to GitLab
  • Sign in
  • T Trilinos
  • Project information
    • Project information
    • Activity
    • Labels
    • Members
  • Repository
    • Repository
    • Files
    • Commits
    • Branches
    • Tags
    • Contributor statistics
    • Graph
    • Compare revisions
  • Issues 936
    • Issues 936
    • List
    • Boards
    • Service Desk
    • Milestones
  • Merge requests 22
    • Merge requests 22
  • CI/CD
    • CI/CD
    • Pipelines
    • Jobs
    • Artifacts
    • Schedules
  • Deployments
    • Deployments
    • Environments
    • Releases
  • Packages and registries
    • Packages and registries
    • Model experiments
  • Monitor
    • Monitor
    • Incidents
  • Analytics
    • Analytics
    • Value stream
    • CI/CD
    • Repository
  • Wiki
    • Wiki
  • Snippets
    • Snippets
  • Activity
  • Graph
  • Create a new issue
  • Jobs
  • Commits
  • Issue Boards
Collapse sidebar
  • James Willenbring
  • Trilinos
  • Issues
  • #3235

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

Assignee
Assign to
Time tracking