Ifpack2_IlukGraph.hpp

Go to the documentation of this file.

/*@HEADER
// ***********************************************************************
//
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2009) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@HEADER
 */


#ifndef IFPACK2_ILUK_GRAPH_HPP
#define IFPACK2_ILUK_GRAPH_HPP

#include <algorithm>
#include <vector>

#include "KokkosSparse_spiluk.hpp"

#include <Ifpack2_ConfigDefs.hpp>
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_CommHelpers.hpp>
#include <Tpetra_CrsGraph.hpp>
#include <Tpetra_Details_WrappedDualView.hpp>
#include <Tpetra_Import.hpp>
#include <Ifpack2_CreateOverlapGraph.hpp>
#include <Ifpack2_Parameters.hpp>

namespace Ifpack2 {

template<class GraphType, class KKHandleType>
class IlukGraph : public Teuchos::Describable {
public:
  typedef typename GraphType::local_ordinal_type local_ordinal_type;
  typedef typename GraphType::global_ordinal_type global_ordinal_type;
  typedef typename GraphType::node_type node_type;

  typedef Tpetra::RowGraph<local_ordinal_type,
                           global_ordinal_type,
                           node_type> row_graph_type;
  typedef Tpetra::CrsGraph<local_ordinal_type,
                           global_ordinal_type,
                           node_type> crs_graph_type;



  typedef typename crs_graph_type::nonconst_global_inds_host_view_type nonconst_global_inds_host_view_type;
  typedef typename crs_graph_type::nonconst_local_inds_host_view_type nonconst_local_inds_host_view_type;
  typedef typename crs_graph_type::global_inds_host_view_type global_inds_host_view_type;
  typedef typename crs_graph_type::local_inds_host_view_type  local_inds_host_view_type;

  IlukGraph (const Teuchos::RCP<const GraphType>& G,
             const int levelFill,
             const int levelOverlap,
             const double overalloc = 2.);

  virtual ~IlukGraph ();

  void setParameters (const Teuchos::ParameterList& parameterlist);

  void initialize();
  void initialize(const Teuchos::RCP<KKHandleType>& KernelHandle);

  int getLevelFill () const { return LevelFill_; }

  int getLevelOverlap () const { return LevelOverlap_; }

  Teuchos::RCP<const GraphType> getA_Graph () const {
    return Graph_;
  }

  Teuchos::RCP<crs_graph_type> getL_Graph () const {
    return L_Graph_;
  }

  Teuchos::RCP<crs_graph_type> getU_Graph () const {
    return U_Graph_;
  }

  Teuchos::RCP<const crs_graph_type> getOverlapGraph () const {
    return OverlapGraph_;
  }

  size_t getNumGlobalDiagonals() const { return NumGlobalDiagonals_; }

private:
  typedef typename GraphType::map_type map_type;  

  IlukGraph (const IlukGraph<GraphType, KKHandleType>&);

  IlukGraph& operator= (const IlukGraph<GraphType, KKHandleType>&);

  void constructOverlapGraph();

  Teuchos::RCP<const GraphType> Graph_;
  Teuchos::RCP<const crs_graph_type> OverlapGraph_;
  int LevelFill_;
  int LevelOverlap_;
  const double Overalloc_;
  Teuchos::RCP<crs_graph_type> L_Graph_;
  Teuchos::RCP<crs_graph_type> U_Graph_;
  size_t NumMyDiagonals_;
  size_t NumGlobalDiagonals_;
};


template<class GraphType, class KKHandleType>
IlukGraph<GraphType, KKHandleType>::
IlukGraph (const Teuchos::RCP<const GraphType>& G,
           const int levelFill,
           const int levelOverlap,
           const double overalloc)
  : Graph_ (G),
    LevelFill_ (levelFill),
    LevelOverlap_ (levelOverlap),
    Overalloc_ (overalloc),
    NumMyDiagonals_ (0),
    NumGlobalDiagonals_ (0)
{
  TEUCHOS_TEST_FOR_EXCEPTION(Overalloc_ <= 1., std::runtime_error,
    "Ifpack2::IlukGraph: FATAL: overalloc must be greater than 1.")
}


template<class GraphType, class KKHandleType>
IlukGraph<GraphType, KKHandleType>::~IlukGraph()
{}


template<class GraphType, class KKHandleType>
void IlukGraph<GraphType, KKHandleType>::
setParameters (const Teuchos::ParameterList& parameterlist)
{
  getParameter (parameterlist, "iluk level-of-fill", LevelFill_);
  getParameter (parameterlist, "iluk level-of-overlap", LevelOverlap_);
}


template<class GraphType, class KKHandleType>
void IlukGraph<GraphType, KKHandleType>::constructOverlapGraph () {
  // FIXME (mfh 22 Dec 2013) This won't do if we want
  // RILUK::initialize() to do the right thing (that is,
  // unconditionally recompute the "symbolic factorization").
  if (OverlapGraph_ == Teuchos::null) {
    OverlapGraph_ = createOverlapGraph<GraphType> (Graph_, LevelOverlap_);
  }
}


template<class GraphType, class KKHandleType>
void IlukGraph<GraphType, KKHandleType>::initialize()
{
  using Teuchos::Array;
  using Teuchos::ArrayView;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::REDUCE_SUM;
  using Teuchos::reduceAll;

  size_t NumIn, NumL, NumU;
  bool DiagFound;

  constructOverlapGraph();

  // Get Maximum Row length
  const int MaxNumIndices = OverlapGraph_->getLocalMaxNumRowEntries ();

  // FIXME (mfh 23 Dec 2013) Use size_t or whatever
  // getLocalNumElements() returns, instead of ptrdiff_t.
  const int NumMyRows = OverlapGraph_->getRowMap ()->getLocalNumElements ();

  using device_type = typename node_type::device_type;
  using execution_space = typename device_type::execution_space;
  using dual_view_type = Kokkos::DualView<size_t*,device_type>;
  dual_view_type numEntPerRow_dv("numEntPerRow",NumMyRows);
  Tpetra::Details::WrappedDualView<dual_view_type> numEntPerRow(numEntPerRow_dv);

  const auto overalloc = Overalloc_;
  const auto levelfill = LevelFill_;
  {
    // Scoping for the  localOverlapGraph access
    auto numEntPerRow_d = numEntPerRow.getDeviceView(Tpetra::Access::OverwriteAll);
    auto localOverlapGraph = OverlapGraph_->getLocalGraphDevice();
    Kokkos::parallel_for("CountOverlapGraphRowEntries",
                         Kokkos::RangePolicy<execution_space>(0, NumMyRows),
                         KOKKOS_LAMBDA(const int i)
                         {
                           // Heuristic to get the maximum number of entries per row.
                           int RowMaxNumIndices = localOverlapGraph.rowConst(i).length;
                           numEntPerRow_d(i) = (levelfill == 0) ? RowMaxNumIndices  // No additional storage needed
                             : Kokkos::ceil(static_cast<double>(RowMaxNumIndices)
                                    * Kokkos::pow(overalloc, levelfill));
                         });
   
  };

  bool insertError;  // No error found yet while inserting entries
  do {
    insertError = false;
    Teuchos::ArrayView<const size_t> a_numEntPerRow(numEntPerRow.getHostView(Tpetra::Access::ReadOnly).data(),NumMyRows);
    L_Graph_ = rcp (new crs_graph_type (OverlapGraph_->getRowMap (),
                                        OverlapGraph_->getRowMap (),
                                        a_numEntPerRow));
    U_Graph_ = rcp (new crs_graph_type (OverlapGraph_->getRowMap (),
                                        OverlapGraph_->getRowMap (),
                                        a_numEntPerRow));

    Array<local_ordinal_type> L (MaxNumIndices);
    Array<local_ordinal_type> U (MaxNumIndices);

    // First we copy the user's graph into L and U, regardless of fill level

    NumMyDiagonals_ = 0;

    for (int i = 0; i< NumMyRows; ++i) {
      local_inds_host_view_type my_indices;
      OverlapGraph_->getLocalRowView (i, my_indices);

      // Split into L and U (we don't assume that indices are ordered).

      NumL = 0;
      NumU = 0;
      DiagFound = false;
      NumIn = my_indices.size();

      for (size_t j = 0; j < NumIn; ++j) {
        const local_ordinal_type k = my_indices[j];

        if (k<NumMyRows) { // Ignore column elements that are not in the square matrix

          if (k==i) {
            DiagFound = true;
          }
          else if (k < i) {
            L[NumL] = k;
            NumL++;
          }
          else {
            U[NumU] = k;
            NumU++;
          }
        }
      }

      // Check in things for this row of L and U

      if (DiagFound) {
        ++NumMyDiagonals_;
      }
      if (NumL) {
        L_Graph_->insertLocalIndices (i, NumL, L.data());
      }
      if (NumU) {
        U_Graph_->insertLocalIndices (i, NumU, U.data());
      }
    }

    if (LevelFill_ > 0) {
      // Complete Fill steps
      RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
      RCP<const map_type> L_RangeMap = Graph_->getRangeMap ();
      RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
      RCP<const map_type> U_RangeMap = OverlapGraph_->getRowMap ();
      RCP<Teuchos::ParameterList> params = Teuchos::parameterList ();
      params->set ("Optimize Storage",false);
      L_Graph_->fillComplete (L_DomainMap, L_RangeMap, params);
      U_Graph_->fillComplete (U_DomainMap, U_RangeMap, params);
      L_Graph_->resumeFill ();
      U_Graph_->resumeFill ();

      // At this point L_Graph and U_Graph are filled with the pattern of input graph,
      // sorted and have redundant indices (if any) removed.  Indices are zero based.
      // LevelFill is greater than zero, so continue...

      int MaxRC = NumMyRows;
      std::vector<std::vector<int> > Levels(MaxRC);
      std::vector<int> LinkList(MaxRC);
      std::vector<int> CurrentLevel(MaxRC);
      Array<local_ordinal_type> CurrentRow (MaxRC + 1);
      std::vector<int> LevelsRowU(MaxRC);

      try {
        for (int i = 0; i < NumMyRows; ++i) {
          int First, Next;

          // copy column indices of row into workspace and sort them

          size_t LenL = L_Graph_->getNumEntriesInLocalRow(i);
          size_t LenU = U_Graph_->getNumEntriesInLocalRow(i);
          size_t Len = LenL + LenU + 1;
          CurrentRow.resize(Len);
          nonconst_local_inds_host_view_type CurrentRow_view(CurrentRow.data(),CurrentRow.size());
          L_Graph_->getLocalRowCopy(i, CurrentRow_view, LenL);  // Get L Indices
          CurrentRow[LenL] = i;                              // Put in Diagonal
          if (LenU > 0) {
            ArrayView<local_ordinal_type> URowView = CurrentRow.view (LenL+1,LenU);            
            nonconst_local_inds_host_view_type URowView_v(URowView.data(),URowView.size());

            // Get U Indices
            U_Graph_->getLocalRowCopy (i, URowView_v, LenU);
          }

          // Construct linked list for current row

          for (size_t j=0; j<Len-1; j++) {
            LinkList[CurrentRow[j]] = CurrentRow[j+1];
            CurrentLevel[CurrentRow[j]] = 0;
          }

          LinkList[CurrentRow[Len-1]] = NumMyRows;
          CurrentLevel[CurrentRow[Len-1]] = 0;

          // Merge List with rows in U

          First = CurrentRow[0];
          Next = First;
          while (Next < i) {
            int PrevInList = Next;
            int NextInList = LinkList[Next];
            int RowU = Next;
            // Get Indices for this row of U
            local_inds_host_view_type IndicesU;
            U_Graph_->getLocalRowView (RowU, IndicesU);
            // FIXME (mfh 23 Dec 2013) size() returns ptrdiff_t, not int.
            int LengthRowU = IndicesU.size ();

            int ii;

            // Scan RowU

            for (ii = 0; ii < LengthRowU; /*nop*/) {
              int CurInList = IndicesU[ii];
              if (CurInList < NextInList) {
                // new fill-in
                int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
                if (NewLevel <= LevelFill_) {
                  LinkList[PrevInList]  = CurInList;
                  LinkList[CurInList] = NextInList;
                  PrevInList = CurInList;
                  CurrentLevel[CurInList] = NewLevel;
                }
                ii++;
              }
              else if (CurInList == NextInList) {
                PrevInList = NextInList;
                NextInList = LinkList[PrevInList];
                int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
                CurrentLevel[CurInList] = std::min (CurrentLevel[CurInList],
                                                    NewLevel);
                ii++;
              }
              else { // (CurInList > NextInList)
                PrevInList = NextInList;
                NextInList = LinkList[PrevInList];
              }
            }
            Next = LinkList[Next];
          }

          // Put pattern into L and U
          CurrentRow.resize(0);

          Next = First;

          // Lower
          while (Next < i) {
            CurrentRow.push_back(Next);
            Next = LinkList[Next];
          }

          // FIXME (mfh 23 Dec 2013) It's not clear to me that
          // removeLocalIndices works like people expect it to work.  In
          // particular, it does not actually change the column Map.
          L_Graph_->removeLocalIndices (i); // Delete current set of Indices
          if (CurrentRow.size() > 0) {
            L_Graph_->insertLocalIndices (i, CurrentRow.size(),CurrentRow.data());
          }

          // Diagonal

          TEUCHOS_TEST_FOR_EXCEPTION(
            Next != i, std::runtime_error,
            "Ifpack2::IlukGraph::initialize: FATAL: U has zero diagonal")

          LevelsRowU[0] = CurrentLevel[Next];
          Next = LinkList[Next];

          // Upper
          CurrentRow.resize(0);
          LenU = 0;

          while (Next < NumMyRows) {
            LevelsRowU[LenU+1] = CurrentLevel[Next];
            CurrentRow.push_back (Next);
            ++LenU;
            Next = LinkList[Next];
          }

          // FIXME (mfh 23 Dec 2013) It's not clear to me that
          // removeLocalIndices works like people expect it to work.  In
          // particular, it does not actually change the column Map.

          U_Graph_->removeLocalIndices (i); // Delete current set of Indices
          if (LenU > 0) {
            U_Graph_->insertLocalIndices (i, CurrentRow.size(),CurrentRow.data());
          }

          // Allocate and fill Level info for this row
          Levels[i] = std::vector<int> (LenU+1);
          for (size_t jj=0; jj<LenU+1; jj++) {
            Levels[i][jj] = LevelsRowU[jj];
          }
        }
      }
      catch (std::runtime_error &e) {
        insertError = true;
        auto numEntPerRow_d = numEntPerRow.getDeviceView(Tpetra::Access::OverwriteAll);
        Kokkos::parallel_for("CountOverlapGraphRowEntries",
          Kokkos::RangePolicy<execution_space>(0, NumMyRows),
          KOKKOS_LAMBDA(const int i)
          {
            const auto numRowEnt = numEntPerRow_d(i);
            numEntPerRow_d(i) = ceil(static_cast<double>((numRowEnt != 0 ? numRowEnt : 1)) * overalloc);
          });
      }
      const int localInsertError = insertError ? 1 : 0;
      int globalInsertError = 0;
      reduceAll (* (OverlapGraph_->getRowMap ()->getComm ()), REDUCE_SUM, 1,
                 &localInsertError, &globalInsertError);
      insertError = globalInsertError > 0;
    }
  } while (insertError);  // do until all insertions complete successfully

  // Complete Fill steps
  RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
  RCP<const map_type> L_RangeMap  = Graph_->getRangeMap ();
  RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
  RCP<const map_type> U_RangeMap  = OverlapGraph_->getRowMap ();
  L_Graph_->fillComplete (L_DomainMap, L_RangeMap);//DoOptimizeStorage is default here...
  U_Graph_->fillComplete (U_DomainMap, U_RangeMap);//DoOptimizeStorage is default here...

  reduceAll<int, size_t> (* (L_DomainMap->getComm ()), REDUCE_SUM, 1,
                          &NumMyDiagonals_, &NumGlobalDiagonals_);
}


template<class GraphType, class KKHandleType>
void IlukGraph<GraphType, KKHandleType>::initialize(const Teuchos::RCP<KKHandleType>& KernelHandle)
{
  using Teuchos::Array;
  using Teuchos::ArrayView;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::REDUCE_SUM;
  using Teuchos::reduceAll;

  typedef typename crs_graph_type::local_graph_device_type local_graph_device_type;
  typedef typename local_graph_device_type::size_type      size_type;
  typedef typename local_graph_device_type::data_type      data_type;
  typedef typename local_graph_device_type::array_layout   array_layout;
  typedef typename local_graph_device_type::device_type    device_type;

  typedef typename Kokkos::View<size_type*, array_layout, device_type> lno_row_view_t;
  typedef typename Kokkos::View<data_type*, array_layout, device_type> lno_nonzero_view_t;
  
  constructOverlapGraph();

  // FIXME (mfh 23 Dec 2013) Use size_t or whatever
  // getLocalNumElements() returns, instead of ptrdiff_t.
  const int NumMyRows = OverlapGraph_->getRowMap()->getLocalNumElements();
  auto localOverlapGraph = OverlapGraph_->getLocalGraphDevice();

  if (KernelHandle->get_spiluk_handle()->get_nrows() < static_cast<size_type>(NumMyRows)) {
    KernelHandle->get_spiluk_handle()->reset_handle(NumMyRows,
                                                    KernelHandle->get_spiluk_handle()->get_nnzL(),
                                                    KernelHandle->get_spiluk_handle()->get_nnzU());
  }

  lno_row_view_t     L_row_map("L_row_map", NumMyRows + 1);
  lno_nonzero_view_t L_entries("L_entries", KernelHandle->get_spiluk_handle()->get_nnzL());
  lno_row_view_t     U_row_map("U_row_map", NumMyRows + 1);
  lno_nonzero_view_t U_entries("U_entries", KernelHandle->get_spiluk_handle()->get_nnzU());

  bool symbolicError;
  do {
    symbolicError = false;
    try {
      KokkosSparse::Experimental::spiluk_symbolic( KernelHandle.getRawPtr(), LevelFill_,
                                                   localOverlapGraph.row_map, localOverlapGraph.entries,
                                                   L_row_map, L_entries, U_row_map, U_entries );
    }
    catch (std::runtime_error &e) {
      symbolicError = true;
      data_type nnzL = static_cast<data_type>(Overalloc_)*L_entries.extent(0);
      data_type nnzU = static_cast<data_type>(Overalloc_)*U_entries.extent(0);    
      KernelHandle->get_spiluk_handle()->reset_handle(NumMyRows, nnzL, nnzU);
      Kokkos::resize(L_entries, KernelHandle->get_spiluk_handle()->get_nnzL());
      Kokkos::resize(U_entries, KernelHandle->get_spiluk_handle()->get_nnzU());
    }
    const int localSymbolicError = symbolicError ? 1 : 0;
    int globalSymbolicError = 0;
    reduceAll (* (OverlapGraph_->getRowMap ()->getComm ()), REDUCE_SUM, 1,
                   &localSymbolicError, &globalSymbolicError);
    symbolicError = globalSymbolicError > 0;
  } while (symbolicError);

  Kokkos::resize(L_entries, KernelHandle->get_spiluk_handle()->get_nnzL());
  Kokkos::resize(U_entries, KernelHandle->get_spiluk_handle()->get_nnzU());
  
  RCP<Teuchos::ParameterList> params = Teuchos::parameterList ();
  params->set ("Optimize Storage",false);
  
  L_Graph_ = rcp (new crs_graph_type (OverlapGraph_->getRowMap (),
                                      OverlapGraph_->getRowMap (), 
                                      L_row_map, L_entries));
  U_Graph_ = rcp (new crs_graph_type (OverlapGraph_->getRowMap (),
                                      OverlapGraph_->getRowMap (), 
                                      U_row_map, U_entries));
                    
  RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
  RCP<const map_type> L_RangeMap  = Graph_->getRangeMap ();
  RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
  RCP<const map_type> U_RangeMap  = OverlapGraph_->getRowMap ();
  L_Graph_->fillComplete (L_DomainMap, L_RangeMap, params);
  U_Graph_->fillComplete (U_DomainMap, U_RangeMap, params);
}

} // namespace Ifpack2

#endif /* IFPACK2_ILUK_GRAPH_HPP */