Merge pull request #3249 from trilinos/Fix-3235

Tpetra,Belos: Fix #3235, & add unit test

packages/belos/tpetra/test/CMakeLists.txt

@@ -13,3 +13,9 @@ IF(${PACKAGE_NAME}_ENABLE_Experimental)
   ADD_SUBDIRECTORY(BlockGCRODR)
 ENDIF() # ${PACKAGE_NAME}_ENABLE_Experimental
 
+TRIBITS_ADD_EXECUTABLE_AND_TEST(
+  Issue_3235
+  SOURCES Issue_3235.cpp
+  COMM mpi
+  NUM_MPI_PROCS 2
+  )

packages/belos/tpetra/test/Issue_3235.cpp

@@ -0,0 +1,80 @@
+#include "BelosLinearProblem.hpp"
+#include "BelosTpetraAdapter.hpp"
+#include "BelosBlockGmresSolMgr.hpp"
+
+#include "Tpetra_Core.hpp"
+#include "Tpetra_CrsMatrix.hpp"
+#include "Tpetra_Map.hpp"
+#include "Tpetra_MultiVector.hpp"
+#include "Tpetra_Operator.hpp"
+
+int main(int argc, char *argv[]) {
+  using Teuchos::Array;
+  using Teuchos::ParameterList;
+  using Teuchos::RCP;
+  using Teuchos::rcp;
+
+  using SC = double;
+  using crs_matrix_type = Tpetra::CrsMatrix<SC>;
+  using map_type = Tpetra::Map<>;
+  using OP = Tpetra::Operator<SC>;
+  using MV = Tpetra::MultiVector<SC>;
+  using LO = map_type::local_ordinal_type; // int LO;
+  using GO = map_type::global_ordinal_type; // int GO;
+  // typedef KokkosClassic::DefaultNode::DefaultNodeType NO;
+
+  Tpetra::ScopeGuard tpetraScope (&argc, &argv);
+  {
+    auto TeuchosComm = Tpetra::getDefaultComm ();
+
+    int NumMyElements = 0;
+    if (TeuchosComm->getRank()==0) {
+      NumMyElements = 10;
+    }
+    Array<GO> uniqueMapArray(NumMyElements);
+    for (LO i=0; i<uniqueMapArray.size(); i++) {
+      uniqueMapArray[i] = i;
+    }
+
+    RCP<const map_type> UniqueMap = rcp (new map_type (-1, uniqueMapArray, 0, TeuchosComm));
+
+    RCP<crs_matrix_type> K = rcp (new crs_matrix_type (UniqueMap,1));
+
+    for (LO i = 0; i < static_cast<LO>(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<MV> Solution = rcp(new MV(UniqueMap,1));
+    RCP<MV> RightHandSide = rcp(new MV(UniqueMap,1));
+
+    Solution->putScalar(0.0);
+    RightHandSide->putScalar(1.0);
+
+    using Belos::LinearProblem;
+    RCP<LinearProblem<SC,MV,OP> > belosLinearProblem =
+      rcp (new LinearProblem<SC,MV,OP> (K, Solution, RightHandSide));
+    belosLinearProblem->setProblem();
+
+    RCP<ParameterList> solverParameterList = rcp(new ParameterList());
+    solverParameterList->set("Convergence Tolerance",1.0e-12);
+    solverParameterList->set("Verbosity",47);
+    solverParameterList->set("Output Frequency",1);
+    solverParameterList->set("Output Style",1);
+
+    Belos::BlockGmresSolMgr<SC,MV,OP> solver(belosLinearProblem,
+                                             solverParameterList);
+    solver.solve();
+  }
+
+  return(EXIT_SUCCESS);
+}

packages/tpetra/core/src/Tpetra_Details_Blas.hpp

@@ -99,7 +99,10 @@ getStride2DView (const ViewType& A)
                  "IndexType must be a built-in integer type.");
   IndexType stride[8];
   A.stride (stride);
-  return (A.extent (1) > 1) ? stride[1] : A.extent (0);
+  // BLAS implementations do not like zero LDA, even if (e.g.,) the
+  // number of rows is actually zero.  See e.g., GitHub Issue #3235.
+  const auto LDA = (A.extent (1) > 1) ? stride[1] : A.extent (0);
+  return LDA == 0 ? IndexType (1) : LDA;
 }
 
 namespace Impl {