Tpetra_Details_determineLocalTriangularStructure.hpp

Go to the documentation of this file.

/*
// @HEADER
// ***********************************************************************
//
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// 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 TPETRA_DETAILS_DETERMINELOCALTRIANGULARSTRUCTURE_HPP
#define TPETRA_DETAILS_DETERMINELOCALTRIANGULARSTRUCTURE_HPP


#include "Kokkos_Core.hpp"
#include "Tpetra_Details_OrdinalTraits.hpp"

namespace Tpetra {
namespace Details {

template<class LO>
struct LocalTriangularStructureResult {
  LO diagCount;
  LO maxNumRowEnt;
  bool couldBeLowerTriangular;
  bool couldBeUpperTriangular;
};

namespace Impl {
  template<class LocalGraphType, class LocalMapType>
  class DetermineLocalTriangularStructure {
  public:
    // Result can't be more than the number of local rows, so
    // local_ordinal_type is appropriate.
    using result_type =
      LocalTriangularStructureResult<typename LocalMapType::local_ordinal_type>;

    DetermineLocalTriangularStructure (const LocalGraphType& G,
                                       const LocalMapType& rowMap,
                                       const LocalMapType& colMap,
                                       const bool ignoreMapsForTriangularStructure) :
      G_ (G),
      rowMap_ (rowMap),
      colMap_ (colMap),
      ignoreMapsForTriangularStructure_ (ignoreMapsForTriangularStructure)
    {}

    KOKKOS_INLINE_FUNCTION void init (result_type& dst) const
    {
      dst.diagCount = 0;
      dst.maxNumRowEnt = 0;
      dst.couldBeLowerTriangular = true; // well, we don't know yet, do we?
      dst.couldBeUpperTriangular = true; // ditto
    }

    KOKKOS_INLINE_FUNCTION void
    join (volatile result_type& dst,
          const volatile result_type& src) const
    {
      dst.diagCount += src.diagCount;
      dst.maxNumRowEnt = (src.maxNumRowEnt > dst.maxNumRowEnt) ?
        src.maxNumRowEnt : dst.maxNumRowEnt;
      dst.couldBeLowerTriangular &= src.couldBeLowerTriangular;
      dst.couldBeUpperTriangular &= src.couldBeUpperTriangular;
    }

    KOKKOS_INLINE_FUNCTION void
    operator () (const typename LocalMapType::local_ordinal_type lclRow,
                 result_type& result) const
    {
      using LO = typename LocalMapType::local_ordinal_type;
      using GO = typename LocalMapType::global_ordinal_type;
      using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;

      auto G_row = G_.rowConst (lclRow);
      const LO numEnt = G_row.length;
      if (numEnt != 0) {
        result.maxNumRowEnt = (numEnt > result.maxNumRowEnt) ?
          numEnt : result.maxNumRowEnt;
        // Use global row and column indices to find the diagonal
        // entry.  Caller promises that local row index is in the row
        // Map on the calling process.
        const GO gblDiagCol = rowMap_.getGlobalElement (lclRow);
        const LO lclDiagCol = colMap_.getLocalElement (gblDiagCol);
        // If it's not in the column Map, then there's no diagonal entry.
        if (lclDiagCol != LOT::invalid ()) {
          // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
          // the sorted case, but note that it requires operator[].
          bool foundDiag = false; // don't count duplicates

          if (ignoreMapsForTriangularStructure_) {
            for (LO k = 0; k < numEnt && ! foundDiag; ++k) {
              const LO lclCol = G_row(k);
              if (lclCol == lclDiagCol) {
                foundDiag = true;
              }
            }
            // mfh 30 Apr 2018: See GitHub Issue #2658.  Per
            // current Tpetra::CrsGraph::computeLocalConstants
            // behavior, assume that local column indices are
            // sorted in each row.
            if (numEnt > LO (0)) {
              const LO smallestLclCol = G_row(0);
              const LO largestLclCol = G_row(numEnt-1); // could be same

              if (smallestLclCol < lclRow) {
                result.couldBeUpperTriangular = false;
              }
              if (lclRow < largestLclCol) {
                result.couldBeLowerTriangular = false;
              }
            }
          }
          else {
            for (LO k = 0; k < numEnt &&
                   ((! foundDiag) ||
                    result.couldBeLowerTriangular ||
                    result.couldBeUpperTriangular);
                 ++k) {
              const LO lclCol = G_row(k);
              if (lclCol == lclDiagCol) {
                foundDiag = true;
              }
              else {
                const GO gblCol = colMap_.getGlobalElement (lclCol);
                if (gblCol < gblDiagCol) {
                  result.couldBeUpperTriangular = false;
                }
                if (gblDiagCol < gblCol) {
                  result.couldBeLowerTriangular = false;
                }
              }
            } // for each entry in lclRow
          } // if-else ignoreMapsForTriangularStructure

          if (foundDiag) {
            ++(result.diagCount);
          }
        }
      }
    }

  private:
    LocalGraphType G_;
    LocalMapType rowMap_;
    LocalMapType colMap_;
    bool ignoreMapsForTriangularStructure_;
  };

} // namespace Impl

template<class LocalGraphType, class LocalMapType>
LocalTriangularStructureResult<typename LocalMapType::local_ordinal_type>
determineLocalTriangularStructure (const LocalGraphType& G,
                                   const LocalMapType& rowMap,
                                   const LocalMapType& colMap,
                                   const bool ignoreMapsForTriangularStructure)
{
  using LO = typename LocalMapType::local_ordinal_type;
  using execution_space = typename LocalGraphType::device_type::execution_space;
  using range_type = Kokkos::RangePolicy<execution_space, LO>;
  using functor_type =
    Impl::DetermineLocalTriangularStructure<LocalGraphType, LocalMapType>;

  LocalTriangularStructureResult<LO> result {0, 0, true, true};
  Kokkos::parallel_reduce ("Tpetra::Details::determineLocalTriangularStructure",
                           range_type (0, G.numRows ()),
                           functor_type (G, rowMap, colMap,
                                         ignoreMapsForTriangularStructure),
                           result);
  return result;
}

} // namespace Details
} // namespace Tpetra

#endif // TPETRA_DETAILS_DETERMINELOCALTRIANGULARSTRUCTURE_HPP