Tpetra_CrsSingletonFilter_LinearProblem_decl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef TPETRA_CRSSINGLETONFILTER_LINEARPROBLEM_DECL_HPP
#define TPETRA_CRSSINGLETONFILTER_LINEARPROBLEM_DECL_HPP
 
 
#include "Tpetra_LinearProblem.hpp"
#include "Tpetra_Transform.hpp"
 
#include "Teuchos_DataAccess.hpp"
 
#include "Kokkos_Core.hpp"
#include "Kokkos_Sort.hpp"
 
#include "Tpetra_Vector.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_RowMatrix.hpp"
 
namespace Tpetra {
 
 
template <class Scalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node>
class CrsSingletonFilter_LinearProblem : public SameTypeTransform<Tpetra::LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node> > {
 public:
 
 
  using scalar_type         = Scalar;
  using local_ordinal_type  = LocalOrdinal;
  using global_ordinal_type = GlobalOrdinal;
 
  using map_type            = Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>;
  using crs_matrix_type     = Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
  using row_matrix_type     = Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
  using multivector_type    = Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
  using vector_type         = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
  using vector_type_int     = Tpetra::Vector<int, LocalOrdinal, GlobalOrdinal, Node>;
  using vector_type_LO      = Tpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>;
  using linear_problem_type = Tpetra::LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
  using import_type         = Tpetra::Import<LocalOrdinal, GlobalOrdinal, Node>;
  using export_type         = Tpetra::Export<LocalOrdinal, GlobalOrdinal, Node>;
 
  using OriginalType      = typename Transform<linear_problem_type, linear_problem_type>::OriginalType;
  using OriginalConstType = typename Transform<linear_problem_type, linear_problem_type>::OriginalConstType;
  using NewType           = typename Transform<linear_problem_type, linear_problem_type>::NewType;
  using NewConstType      = typename Transform<linear_problem_type, linear_problem_type>::NewConstType;
 
  using nonconst_local_inds_host_view_type  = typename row_matrix_type::nonconst_local_inds_host_view_type;
  using nonconst_global_inds_host_view_type = typename row_matrix_type::nonconst_global_inds_host_view_type;
  using nonconst_values_host_view_type      = typename row_matrix_type::nonconst_values_host_view_type;
 
  using local_map_type = typename map_type::local_map_type;
 
  using local_matrix_type         = typename crs_matrix_type::local_matrix_device_type;
  using local_ptr_view_type       = typename crs_matrix_type::row_ptrs_device_view_type;
  using local_ind_view_type       = typename crs_matrix_type::local_inds_device_view_type;
  using local_val_view_type       = typename crs_matrix_type::local_matrix_device_type::values_type;
  using const_local_val_view_type = typename local_val_view_type::const_type;
 
  using impl_scalar_type             = typename multivector_type::impl_scalar_type;
  using local_multivector_type       = typename multivector_type::dual_view_type::t_dev;
  using const_local_multivector_type = typename local_multivector_type::const_type;
 
  using local_vector_int_type = typename vector_type_int::dual_view_type::t_dev;
 
  using device_type             = typename Node::device_type;
  using vector_view_type_int    = Kokkos::View<local_ordinal_type*, device_type>;
  using vector_view_type_scalar = Kokkos::View<impl_scalar_type*, device_type>;
 
  using execution_space = typename vector_type_int::execution_space;
  using range_policy    = Kokkos::RangePolicy<execution_space>;
 
 
 
  CrsSingletonFilter_LinearProblem(bool run_on_host = false, bool verbose = false);
 
  virtual ~CrsSingletonFilter_LinearProblem() = default;
 
  NewType operator()(const OriginalType& originalLinearProblem);
 
  void analyze(const OriginalType& originalLinearProblem);
 
  NewType construct();
 
  void fwd();
 
  void rvs();
 
 
  void Analyze(const Teuchos::RCP<row_matrix_type>& FullMatrix);
 
  bool SingletonsDetected() const {
    if (!AnalysisDone_)
      return (false);
    else
      return (NumSingletons() > 0);
  }
 
 
 
  void ConstructReducedProblem(const Teuchos::RCP<linear_problem_type>& Problem);
 
  void UpdateReducedProblem(const Teuchos::RCP<linear_problem_type>& Problem);
 
 
 
  void ComputeFullSolution();
 
 
 
  int NumSingletonRows() const { return (globalNumSingletonRows_); }
 
  int NumSingletonCols() const { return (globalNumSingletonCols_); }
 
  int NumSingletons() const { return (NumSingletonCols() + NumSingletonRows()); }
 
  double RatioOfDimensions() const { return (RatioOfDimensions_); }
 
  double RatioOfNonzeros() const { return (RatioOfNonzeros_); }
 
 
 
  Teuchos::RCP<linear_problem_type> FullProblem() const { return (FullProblem_); }
 
  Teuchos::RCP<linear_problem_type> ReducedProblem() const { return (ReducedProblem_); }
 
  Teuchos::RCP<row_matrix_type> FullMatrix() const { return (FullMatrix_); }
 
  Teuchos::RCP<crs_matrix_type> ReducedMatrix() const { return ReducedMatrix_; }
 
  Teuchos::RCP<const map_type> ReducedMatrixRowMap() const { return ReducedMatrixRowMap_; }
 
  Teuchos::RCP<const map_type> ReducedMatrixColMap() const { return ReducedMatrixColMap_; }
 
  Teuchos::RCP<const map_type> ReducedMatrixDomainMap() const { return ReducedMatrixDomainMap_; }
 
  Teuchos::RCP<const map_type> ReducedMatrixRangeMap() const { return ReducedMatrixRangeMap_; }
 
 protected:
  //  This RCP will be null if full matrix is not a CrsMatrix.
  Teuchos::RCP<crs_matrix_type> FullCrsMatrix() const { return (FullCrsMatrix_); }
 
  Teuchos::RCP<const map_type> FullMatrixRowMap() const { return (FullMatrix()->getRowMap()); }
  Teuchos::RCP<const map_type> FullMatrixColMap() const { return (FullMatrix()->getColMap()); }
  Teuchos::RCP<const map_type> FullMatrixDomainMap() const { return (FullMatrix()->getDomainMap()); }
  Teuchos::RCP<const map_type> FullMatrixRangeMap() const { return (FullMatrix()->getRangeMap()); }
 
  // int ComputeEliminateMaps();
  int Setup(const Teuchos::RCP<linear_problem_type>& Problem);
  void InitFullMatrixAccess();
  // void GetRow(int localRow, int & NumIndices, int * & Indices);
  void GetRow(local_ordinal_type localRow, size_t& NumIndices,
              Teuchos::Array<local_ordinal_type>& localIndices);
  void GetRow(LocalOrdinal Row, size_t& NumIndices, Teuchos::ArrayView<const Scalar>& Values,
              Teuchos::ArrayView<const LocalOrdinal>& Indices);
 
  void GetRowGCIDs(LocalOrdinal Row, size_t& NumIndices, Teuchos::ArrayView<const Scalar>& Values,
                   Teuchos::Array<GlobalOrdinal>& GlobalIndices);
 
  void CreatePostSolveArrays(vector_type_LO localRowIDofSingletonCol,
                             vector_type_LO ColProfiles,
                             vector_type_LO NewColProfiles,
                             vector_type_LO ColHasRowWithSingleton);
 
  void ConstructRedistributeExporter(Teuchos::RCP<const map_type> SourceMap, Teuchos::RCP<const map_type> TargetMap,
                                     Teuchos::RCP<export_type>& RedistributeExporter,
                                     Teuchos::RCP<const map_type>& RedistributeMap);
 
  Teuchos::RCP<const map_type> GenerateReducedMap(const Teuchos::RCP<const map_type>& originalMap,
                                                  const Teuchos::RCP<vector_type_int>& mapColors,
                                                  int color = 0, bool locally_sort_gids = true);
 
  Teuchos::RCP<linear_problem_type> FullProblem_;
  Teuchos::RCP<linear_problem_type> ReducedProblem_;
  Teuchos::RCP<row_matrix_type> FullMatrix_;
  Teuchos::RCP<crs_matrix_type> FullCrsMatrix_;
  Teuchos::RCP<crs_matrix_type> ReducedMatrix_;
  Teuchos::RCP<multivector_type> ReducedRHS_;
  Teuchos::RCP<multivector_type> ReducedLHS_;
 
  Teuchos::RCP<const map_type> ReducedMatrixRowMap_;
  Teuchos::RCP<const map_type> ReducedMatrixColMap_;
  Teuchos::RCP<const map_type> ReducedMatrixDomainMap_;
  Teuchos::RCP<const map_type> ReducedMatrixRangeMap_;
  Teuchos::RCP<const map_type> OrigReducedMatrixDomainMap_;
  Teuchos::RCP<import_type> Full2ReducedRHSImporter_;
  Teuchos::RCP<import_type> Full2ReducedLHSImporter_;
  Teuchos::RCP<export_type> RedistributeDomainExporter_;
 
  vector_view_type_int ColSingletonRowLIDs_;
  vector_view_type_int ColSingletonColLIDs_;
  vector_view_type_int ColSingletonPivotLIDs_;
  vector_view_type_scalar ColSingletonPivots_;
 
  local_ordinal_type localNumSingletonRows_;  
  local_ordinal_type localNumSingletonCols_;  
  local_ordinal_type globalNumSingletonRows_;
  local_ordinal_type globalNumSingletonCols_;
  double RatioOfDimensions_;
  double RatioOfNonzeros_;
 
  bool HaveReducedProblem_;
  // bool UserDefinedEliminateMaps_;
  bool AnalysisDone_;
  bool SymmetricElimination_;
 
  Teuchos::RCP<multivector_type> tempExportX_;
  Teuchos::RCP<multivector_type> tempX_;
  Teuchos::RCP<multivector_type> tempB_;
  // Teuchos::RCP<multivector_type> RedistributeReducedLHS_;
 
  // Maximum number of entries in any row of the matrix, on this process.
  local_ordinal_type localMaxNumRowEntries_;
 
  Teuchos::RCP<vector_type_int> RowMapColors_;
  Teuchos::RCP<vector_type_int> ColMapColors_;
  bool FullMatrixIsCrsMatrix_;
 
  bool run_on_host_;
  bool verbose_;
 
 private:
  CrsSingletonFilter_LinearProblem(const Teuchos::RCP<CrsSingletonFilter_LinearProblem>& /* Problem */) {}
};
 
}  // namespace Tpetra
 
#endif  //  TPETRA_CRSSINGLETONFILTER_LINEARPROBLEM_DECL_HPP