From fb8dc4ab715d1ddf01abcdb9d57d8a4a854eee82 Mon Sep 17 00:00:00 2001
From: Mark Hoemmen <mhoemme@sandia.gov>
Date: Tue, 7 Aug 2018 12:54:26 -0600
Subject: [PATCH] Tpetra,Belos: Fix #3235, & add unit test

@trilinos/tpetra @trilinos/belos

Tpetra::MultiVector::multiply calls Tpetra::Details::Blas::gemm.  The
latter was sending LDA=0 into DGEMM, in the case where some process
has a (multi)vector with zero rows.  BLAS implementations do not like
LDA=0, even in this case.  This commit fixes the issue by ensuring
that LDA is never zero.
---
 packages/belos/tpetra/test/CMakeLists.txt     |  6 ++
 packages/belos/tpetra/test/Issue_3235.cpp     | 80 +++++++++++++++++++
 .../tpetra/core/src/Tpetra_Details_Blas.hpp   |  5 +-
 3 files changed, 90 insertions(+), 1 deletion(-)
 create mode 100644 packages/belos/tpetra/test/Issue_3235.cpp

diff --git a/packages/belos/tpetra/test/CMakeLists.txt b/packages/belos/tpetra/test/CMakeLists.txt
index a3d5233c95ea..65296876cb82 100644
--- a/packages/belos/tpetra/test/CMakeLists.txt
+++ b/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
+  )
diff --git a/packages/belos/tpetra/test/Issue_3235.cpp b/packages/belos/tpetra/test/Issue_3235.cpp
new file mode 100644
index 000000000000..d6d9ed2ff319
--- /dev/null
+++ b/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);
+}
diff --git a/packages/tpetra/core/src/Tpetra_Details_Blas.hpp b/packages/tpetra/core/src/Tpetra_Details_Blas.hpp
index 2a377f343077..2e8a1156c998 100644
--- a/packages/tpetra/core/src/Tpetra_Details_Blas.hpp
+++ b/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 {