Merge Pull Request #10829 from trilinos/Trilinos/csiefer-39fa682

Automatically Merged using Trilinos Pull Request AutoTester
PR Title: MueLu: Enabling BlockCrs through the hierarchy
PR Author: csiefer2

packages/ifpack2/src/Ifpack2_BlockRelaxation_def.hpp

@@ -1053,6 +1053,11 @@ description () const
   } else {
     out << "INVALID";
   }
+
+  // BlockCrs if we have that
+  if(hasBlockCrsMatrix_)
+    out<<", BlockCrs";
+
   // Print the approximate # rows per part
   int approx_rows_per_part = A_->getLocalNumRows()/Partitioner_->numLocalParts();
   out <<", blocksize: "<<approx_rows_per_part;

packages/muelu/adapters/tpetra/MueLu_CreateTpetraPreconditioner.hpp

@@ -171,6 +171,25 @@ namespace MueLu {
     MueLu::ReuseXpetraPreconditioner<SC,LO,GO,NO>(A, H);
   }
 
+  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
+  void ReuseTpetraPreconditioner(const Teuchos::RCP<Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& inA,
+                                 MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op) {
+    typedef Scalar          SC;
+    typedef LocalOrdinal    LO;
+    typedef GlobalOrdinal   GO;
+    typedef Node            NO;
+
+    typedef Xpetra::Matrix<SC,LO,GO,NO>     Matrix;
+    typedef MueLu ::Hierarchy<SC,LO,GO,NO>  Hierarchy;
+
+    RCP<Hierarchy> H = Op.GetHierarchy();
+    RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > temp = rcp(new Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO>(inA));
+    TEUCHOS_TEST_FOR_EXCEPTION(temp==Teuchos::null, Exceptions::RuntimeError, "ReuseTpetraPreconditioner: cast from Tpetra::BlockCrsMatrix to Xpetra::TpetraBlockCrsMatrix failed.");
+    RCP<Matrix> A = rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(temp));
+
+    MueLu::ReuseXpetraPreconditioner<SC,LO,GO,NO>(A, H);
+  }
+
 
 
 } //namespace

packages/muelu/src/Graph/MatrixTransformation/MueLu_AmalgamationFactory_def.hpp

@@ -75,11 +75,30 @@ namespace MueLu {
 
     RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
 
+    /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
+       fullblocksize is the number of storage blocks that must kept together during the amalgamation process.
+
+       Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
+
+       numPDEs = fullblocksize * storageblocksize.
+       
+       If numPDEs==1
+         Matrix is point storage (classical CRS storage).  storageblocksize=1 and fullblocksize=1
+         No other values makes sense.
+
+       If numPDEs>1
+         If matrix uses point storage, then storageblocksize=1  and fullblockssize=numPDEs.
+         If matrix uses block storage, with block size of n, then storageblocksize=n, and fullblocksize=numPDEs/n.  
+         Thus far, only storageblocksize=numPDEs and fullblocksize=1 has been tested.
+    */
+
+
     LO fullblocksize    = 1;   // block dim for fixed size blocks
     GO offset           = 0;   // global offset of dof gids
     LO blockid          = -1;  // block id in strided map
     LO nStridedOffset   = 0;   // DOF offset for strided block id "blockid" (default = 0)
     LO stridedblocksize = fullblocksize; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
+    LO storageblocksize = A->GetStorageBlockSize();
     // GO indexBase        = A->getRowMap()->getIndexBase();  // index base for maps (unused)
 
     // 1) check for blocking/striding information
@@ -101,13 +120,20 @@ namespace MueLu {
       } else {
         stridedblocksize = fullblocksize;
       }
+      // Correct for the storageblocksize
+      // NOTE:  Before this point fullblocksize is actually numPDEs
+      TEUCHOS_TEST_FOR_EXCEPTION(fullblocksize % storageblocksize != 0,Exceptions::RuntimeError,"AmalgamationFactory: fullblocksize needs to be a multiple of A->GetStorageBlockSize()");
+      fullblocksize /= storageblocksize;
+      stridedblocksize /= storageblocksize;
+
       oldView = A->SwitchToView(oldView);
       GetOStream(Runtime1) << "AmalagamationFactory::Build():" << " found fullblocksize=" << fullblocksize << " and stridedblocksize=" << stridedblocksize << " from strided maps. offset=" << offset << std::endl;
 
     } else {
       GetOStream(Warnings0) << "AmalagamationFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
     }
 
+
     // build node row map (uniqueMap) and node column map (nonUniqueMap)
     // the arrays rowTranslation and colTranslation contain the local node id
     // given a local dof id. They are only necessary for the CoalesceDropFactory if
@@ -166,7 +192,7 @@ namespace MueLu {
     container               filter;
 
     GO offset = 0;
-    LO blkSize = A.GetFixedBlockSize();
+    LO blkSize = A.GetFixedBlockSize() / A.GetStorageBlockSize();
     if (A.IsView("stridedMaps") == true) {
       Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
       Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);

packages/muelu/src/Graph/MatrixTransformation/MueLu_AmalgamationFactory_kokkos_def.hpp

@@ -80,11 +80,29 @@ namespace MueLu {
 
     RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
 
+   /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
+       fullblocksize is the number of storage blocks that must kept together during the amalgamation process.
+
+       Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
+
+       numPDEs = fullblocksize * storageblocksize.
+       
+       If numPDEs==1
+         Matrix is point storage (classical CRS storage).  storageblocksize=1 and fullblocksize=1
+         No other values makes sense.
+
+       If numPDEs>1
+         If matrix uses point storage, then storageblocksize=1  and fullblockssize=numPDEs.
+         If matrix uses block storage, with block size of n, then storageblocksize=n, and fullblocksize=numPDEs/n.  
+         Thus far, only storageblocksize=numPDEs and fullblocksize=1 has been tested.
+    */
+
     LO fullblocksize    = 1;   // block dim for fixed size blocks
     GO offset           = 0;   // global offset of dof gids
     LO blockid          = -1;  // block id in strided map
     LO nStridedOffset   = 0;   // DOF offset for strided block id "blockid" (default = 0)
     LO stridedblocksize = fullblocksize; // size of strided block id "blockid" (default = fullblocksize, only if blockid!=-1 stridedblocksize <= fullblocksize)
+    LO storageblocksize = A->GetStorageBlockSize();
     // GO indexBase        = A->getRowMap()->getIndexBase();  // index base for maps (unused)
 
     // 1) check for blocking/striding information
@@ -106,6 +124,12 @@ namespace MueLu {
       } else {
         stridedblocksize = fullblocksize;
       }
+      // Correct for the storageblocksize
+      // NOTE:  Before this point fullblocksize is actually numPDEs
+      TEUCHOS_TEST_FOR_EXCEPTION(fullblocksize % storageblocksize != 0,Exceptions::RuntimeError,"AmalgamationFactory::Build(): fullblocksize needs to be a multiple of A->GetStorageBlockSize()");
+      fullblocksize /= storageblocksize;
+      stridedblocksize /= storageblocksize;
+
       oldView = A->SwitchToView(oldView);
       GetOStream(Runtime1) << "AmalagamationFactory::Build():" << " found fullblocksize=" << fullblocksize << " and stridedblocksize=" << stridedblocksize << " from strided maps. offset=" << offset << std::endl;
 
@@ -172,7 +196,7 @@ namespace MueLu {
     container               filter; // TODO:  replace std::set with an object having faster lookup/insert, hashtable for instance
 
     GO offset = 0;
-    LO blkSize = A.GetFixedBlockSize();
+    LO blkSize = A.GetFixedBlockSize() / A.GetStorageBlockSize();
     if (A.IsView("stridedMaps") == true) {
       Teuchos::RCP<const Map> myMap = A.getRowMap("stridedMaps");
       Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);

packages/muelu/src/Graph/MatrixTransformation/MueLu_AmalgamationInfo_def.hpp

@@ -132,7 +132,7 @@ namespace MueLu {
 
   template <class LocalOrdinal, class GlobalOrdinal, class Node>
   void AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node>::UnamalgamateAggregatesLO(const Aggregates& aggregates,
-        Teuchos::ArrayRCP<LO>& aggStart, Teuchos::ArrayRCP<LO>& aggToRowMap) const {
+                      Teuchos::ArrayRCP<LO>& aggStart, Teuchos::ArrayRCP<LO>& aggToRowMap) const {
 
     int myPid = aggregates.GetMap()->getComm()->getRank();
     Teuchos::ArrayView<const GO> nodeGlobalElts = aggregates.GetMap()->getLocalElementList();

packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp

@@ -370,7 +370,7 @@ namespace MueLu {
           distanceLaplacianAlgo = scaled_cut_symmetric;
         else
           TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut), not \"" << distanceLaplacianAlgoStr << "\"");
-        GetOStream(Runtime0) << "algorithm = \"" << algo << "\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize() << std::endl;
+        GetOStream(Runtime0) << "algorithm = \"" << algo << "\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr << "\": threshold = " << threshold << ", blocksize = " << A->GetFixedBlockSize()<< std::endl;
       } else if (algo == "classical") {
         if (classicalAlgoStr == "default")
           classicalAlgo = defaultAlgo;
@@ -396,6 +396,27 @@ namespace MueLu {
       GO numDropped = 0, numTotal = 0;
       std::string graphType = "unamalgamated"; //for description purposes only
 
+
+      /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
+       BlockSize is the number of storage blocks that must kept together during the amalgamation process.
+
+       Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
+
+       numPDEs = BlockSize * storageblocksize.
+       
+       If numPDEs==1
+         Matrix is point storage (classical CRS storage).  storageblocksize=1 and BlockSize=1
+         No other values makes sense.
+
+       If numPDEs>1
+         If matrix uses point storage, then storageblocksize=1  and BlockSize=numPDEs.
+         If matrix uses block storage, with block size of n, then storageblocksize=n, and BlockSize=numPDEs/n.  
+         Thus far, only storageblocksize=numPDEs and BlockSize=1 has been tested.
+      */      
+      TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0,Exceptions::RuntimeError,"A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
+      const LO BlockSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
+
+
       /************************** RS or SA-style Classical Dropping (and variants) **************************/
       if (algo == "classical") {
         if (predrop_ == null) {
@@ -417,7 +438,7 @@ namespace MueLu {
         // At this points we either have
         //     (predrop_ != null)
         // Therefore, it is sufficient to check only threshold
-        if (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && !useSignedClassicalRS && !useSignedClassicalSA && A->hasCrsGraph()) {
+        if ( BlockSize==1 && threshold == STS::zero() && !useSignedClassicalRS && !useSignedClassicalSA && A->hasCrsGraph()) {
           // Case 1:  scalar problem, no dropping => just use matrix graph
           RCP<GraphBase> graph = rcp(new Graph(A->getCrsGraph(), "graph of A"));
           // Detect and record rows that correspond to Dirichlet boundary conditions
@@ -442,10 +463,10 @@ namespace MueLu {
           Set(currentLevel, "DofsPerNode", 1);
           Set(currentLevel, "Graph", graph);
 
-        } else if ( (A->GetFixedBlockSize() == 1 && threshold != STS::zero()) ||
-                    (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && !A->hasCrsGraph()) ||
-                    (A->GetFixedBlockSize() == 1 && useSignedClassicalRS) ||
-                    (A->GetFixedBlockSize() == 1 && useSignedClassicalSA) ) {	  
+        } else if ( (BlockSize == 1 && threshold != STS::zero()) ||
+                    (BlockSize == 1 && threshold == STS::zero() && !A->hasCrsGraph()) ||
+                    (BlockSize == 1 && useSignedClassicalRS) ||
+                    (BlockSize == 1 && useSignedClassicalSA) ) {	  
           // Case 2:  scalar problem with dropping => record the column indices of undropped entries, but still use original
           //                                          graph's map information, e.g., whether index is local
           // OR a matrix without a CrsGraph
@@ -721,7 +742,7 @@ namespace MueLu {
            }
 #endif
           }//end generateColoringGraph
-        } else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
+        } else if (BlockSize > 1 && threshold == STS::zero()) {
           // Case 3:  Multiple DOF/node problem without dropping
           const RCP<const Map> rowMap = A->getRowMap();
           const RCP<const Map> colMap = A->getColMap();
@@ -853,7 +874,7 @@ namespace MueLu {
           Set(currentLevel, "Graph",       graph);
           Set(currentLevel, "DofsPerNode", blkSize); // full block size
 
-        } else if (A->GetFixedBlockSize() > 1 && threshold != STS::zero()) {
+        } else if (BlockSize > 1 && threshold != STS::zero()) {
           // Case 4:  Multiple DOF/node problem with dropping
           const RCP<const Map> rowMap = A->getRowMap();
           const RCP<const Map> colMap = A->getColMap();

packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_kokkos_def.hpp

@@ -506,7 +506,27 @@ namespace MueLu {
     const MT zero = Teuchos::ScalarTraits<MT>::zero();
 
     auto A         = Get< RCP<Matrix> >(currentLevel, "A");
-    LO   blkSize   = A->GetFixedBlockSize();
+
+
+    /* NOTE: storageblocksize (from GetStorageBlockSize()) is the size of a block in the chosen storage scheme.
+       blkSize is the number of storage blocks that must kept together during the amalgamation process.
+
+       Both of these quantities may be different than numPDEs (from GetFixedBlockSize()), but the following must always hold:
+
+       numPDEs = blkSize * storageblocksize.
+       
+       If numPDEs==1
+         Matrix is point storage (classical CRS storage).  storageblocksize=1 and  blkSize=1
+         No other values makes sense.
+
+       If numPDEs>1
+         If matrix uses point storage, then storageblocksize=1  and blkSize=numPDEs.
+         If matrix uses block storage, with block size of n, then storageblocksize=n, and blkSize=numPDEs/n.  
+         Thus far, only storageblocksize=numPDEs and blkSize=1 has been tested.
+      */      
+
+    TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0,Exceptions::RuntimeError,"A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
+    LO   blkSize   = A->GetFixedBlockSize() / A->GetStorageBlockSize();
 
     auto amalInfo = Get< RCP<AmalgamationInfo_kokkos> >(currentLevel, "UnAmalgamationInfo");
 
@@ -542,7 +562,7 @@ namespace MueLu {
       boundaryNodes = Utilities_kokkos::DetectDirichletRows(*A, dirichletThreshold);
 
       // Trivial LWGraph construction
-      graph = rcp(new LWGraph_kokkos(A->getLocalMatrixDevice().graph, A->getRowMap(), A->getColMap(), "graph of A"));
+      graph = rcp(new LWGraph_kokkos(A->getCrsGraph()->getLocalGraphDevice(), A->getRowMap(), A->getColMap(), "graph of A"));
       graph->getLocalLWGraph().SetBoundaryNodeMap(boundaryNodes);
 
       numTotal = A->getLocalNumEntries();

packages/muelu/src/Misc/MueLu_RAPFactory_def.hpp

@@ -56,6 +56,7 @@
 #include <Xpetra_TripleMatrixMultiply.hpp>
 #include <Xpetra_Vector.hpp>
 #include <Xpetra_VectorFactory.hpp>
+#include <Xpetra_IO.hpp>
 
 #include "MueLu_RAPFactory_decl.hpp"
 
@@ -281,8 +282,7 @@ namespace MueLu {
           Xpetra::TripleMatrixMultiply<SC,LO,GO,NO>::
             MultiplyRAP(*P, doTranspose, *A, !doTranspose, *P, !doTranspose, *Ac, doFillComplete,
                         doOptimizeStorage, labelstr+std::string("MueLu::R*A*P-implicit-")+levelstr.str(),
-                        RAPparams);
-
+                        RAPparams);         
         } else {
           RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
           Ac = MatrixFactory::Build(R->getRowMap(), Teuchos::as<LO>(0));

packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp

@@ -1245,9 +1245,10 @@ namespace MueLu {
           break;
         }
 
-        Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
+        LO storageblocksize=Am->GetStorageBlockSize();
+        Xpetra::global_size_t nnz = Am->getGlobalNumEntries()*storageblocksize*storageblocksize;
         nnzPerLevel     .push_back(nnz);
-        rowsPerLevel    .push_back(Am->getGlobalNumRows());
+        rowsPerLevel    .push_back(Am->getGlobalNumRows()*storageblocksize);
         numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
       }
 
@@ -1434,8 +1435,9 @@ namespace MueLu {
     }
 
     GetOStream(Runtime1) << "Replacing coordinate map" << std::endl;
+    TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0, Exceptions::RuntimeError, "Hierarchy::ReplaceCoordinateMap: Storage block size does not evenly divide fixed block size");
 
-    size_t blkSize = A->GetFixedBlockSize();
+    size_t blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
 
     RCP<const Map> nodeMap = A->getRowMap();
     if (blkSize > 1) {

packages/muelu/src/Smoothers/MueLu_Ifpack2Smoother_def.hpp

@@ -67,6 +67,7 @@
 #include <Xpetra_CrsMatrixWrap.hpp>
 #include <Xpetra_TpetraBlockCrsMatrix.hpp>
 #include <Xpetra_Matrix.hpp>
+#include <Xpetra_MatrixMatrix.hpp>
 #include <Xpetra_MultiVectorFactory.hpp>
 #include <Xpetra_TpetraMultiVector.hpp>
 
@@ -234,16 +235,21 @@ namespace MueLu {
         if(Acrs.is_null())
           throw std::runtime_error("Ifpack2Smoother: Cannot extract CrsMatrix from matrix A.");
         RCP<TpetraCrsMatrix> At = rcp_dynamic_cast<TpetraCrsMatrix>(Acrs);
-        if(At.is_null())
-          throw std::runtime_error("Ifpack2Smoother: Cannot extract TpetraCrsMatrix from matrix A.");
-
-        RCP<Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> > blockCrs = Tpetra::convertToBlockCrsMatrix(*At->getTpetra_CrsMatrix(),blocksize);
-        RCP<CrsMatrix> blockCrs_as_crs  = rcp(new TpetraBlockCrsMatrix(blockCrs));
-        RCP<CrsMatrixWrap> blockWrap = rcp(new CrsMatrixWrap(blockCrs_as_crs));
-        A_ = blockWrap;
-        this->GetOStream(Statistics0) << "Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
-
-        paramList.remove("smoother: use blockcrsmatrix storage");
+        if(At.is_null()) {
+          if(!Xpetra::Helpers<Scalar,LO,GO,Node>::isTpetraBlockCrs(matA))
+            throw std::runtime_error("Ifpack2Smoother: Cannot extract CrsMatrix or BlockCrsMatrix from matrix A.");          
+          this->GetOStream(Statistics0) << "Ifpack2Smoother: Using (native) BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
+          paramList.remove("smoother: use blockcrsmatrix storage");
+        }
+        else {
+          RCP<Tpetra::BlockCrsMatrix<Scalar, LO, GO, Node> > blockCrs = Tpetra::convertToBlockCrsMatrix(*At->getTpetra_CrsMatrix(),blocksize);
+          RCP<CrsMatrix> blockCrs_as_crs  = rcp(new TpetraBlockCrsMatrix(blockCrs));
+          RCP<CrsMatrixWrap> blockWrap = rcp(new CrsMatrixWrap(blockCrs_as_crs));
+          A_ = blockWrap;
+          this->GetOStream(Statistics0) << "Ifpack2Smoother: Using BlockCrsMatrix storage with blocksize "<<blocksize<<std::endl;
+
+          paramList.remove("smoother: use blockcrsmatrix storage");
+        }
       }
     }