Teuchos_MatrixMarket_Raw_Adder.hpp

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
#define __Teuchos_MatrixMarket_Raw_Adder_hpp

#include "Teuchos_ConfigDefs.hpp"
#include "Teuchos_ArrayRCP.hpp"
#include "Teuchos_CommHelpers.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_MatrixMarket_Banner.hpp"
#include "Teuchos_MatrixMarket_CoordDataReader.hpp"

#include <algorithm>
#include <fstream>
#include <iostream>
#include <iterator>
#include <vector>
#include <stdexcept>


namespace Teuchos {
  namespace MatrixMarket {
    namespace Raw {
      template<class Scalar, class Ordinal>
      class Element {
      public:
        Element () :
          rowIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ()),
          colIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ()),
          value_ (Teuchos::ScalarTraits<Scalar>::zero ())
        {}

        Element (const Ordinal i, const Ordinal j, const Scalar& Aij) :
          rowIndex_ (i), colIndex_ (j), value_ (Aij) {}

        bool operator== (const Element& rhs) const {
          return rowIndex_ == rhs.rowIndex_ && colIndex_ == rhs.colIndex_;
        }

        bool operator!= (const Element& rhs) const {
          return ! (*this == rhs);
        }

        bool operator< (const Element& rhs) const {
          if (rowIndex_ < rhs.rowIndex_)
            return true;
          else if (rowIndex_ > rhs.rowIndex_)
            return false;
          else { // equal
            return colIndex_ < rhs.colIndex_;
          }
        }

        template<class BinaryFunction>
        void merge (const Element& rhs, const BinaryFunction& f) {
          TEUCHOS_TEST_FOR_EXCEPTION(
            rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
            std::invalid_argument,
            "Attempt to merge elements at different locations in the sparse "
            "matrix.  The current element is at (" << rowIndex() << ", "
            << colIndex() << ") and the element you asked me to merge with it "
            "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << ").  This "
            "probably indicates a bug in the sparse matrix reader.");

          value_ = f (rhs.value_, value_);
        }

        void merge (const Element& rhs, const bool replace=false) {
          TEUCHOS_TEST_FOR_EXCEPTION(
            rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
            std::invalid_argument,
            "Attempt to merge elements at different locations in the sparse "
            "matrix.  The current element is at (" << rowIndex() << ", "
            << colIndex() << ") and the element you asked me to merge with it "
            "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << ").  This "
            "probably indicates a bug in the sparse matrix reader.");

          if (replace) {
            value_ = rhs.value_;
          }
          else {
            value_ += rhs.value_;
          }
        }

        Ordinal rowIndex() const { return rowIndex_; }

        Ordinal colIndex() const { return colIndex_; }

        Scalar value() const { return value_; }

      private:
        Ordinal rowIndex_, colIndex_;
        Scalar value_;
      };

      template<class Scalar, class Ordinal>
      std::ostream&
      operator<< (std::ostream& out, const Element<Scalar, Ordinal>& elt)
      {
        typedef ScalarTraits<Scalar> STS;
        std::ios::fmtflags f( out.flags() );
        // Non-Ordinal types are floating-point types.  In order not to
        // lose information when we print a floating-point type, we have
        // to set the number of digits to print.  C++ standard behavior
        // in the default locale seems to be to print only five decimal
        // digits after the decimal point; this does not suffice for
        // double precision.  We solve the problem of how many digits to
        // print more generally below.  It's a rough solution so please
        // feel free to audit and revise it.
        //
        // FIXME (mfh 01 Feb 2011)
        // This really calls for the following approach:
        //
        // Guy L. Steele and Jon L. White, "How to print floating-point
        // numbers accurately", 20 Years of the ACM/SIGPLAN Conference
        // on Programming Language Design and Implementation
        // (1979-1999): A Selection, 2003.
        if (! STS::isOrdinal) {
          // std::scientific, std::fixed, and default are the three
          // output states for floating-point numbers.  A reasonable
          // user-defined floating-point type should respect these
          // flags; hopefully it does.
          out << std::scientific;

          // Decimal output is standard for Matrix Market format.
          out << std::setbase (10);

          // Compute the number of decimal digits required for expressing
          // a Scalar, by comparing with IEEE 754 double precision (16
          // decimal digits, 53 binary digits).  This would be easier if
          // Teuchos exposed std::numeric_limits<T>::digits10, alas.
          const double numDigitsAsDouble =
            16 * ((double) STS::t() / (double) ScalarTraits<double>::t());
          // Adding 0.5 and truncating is a portable "floor".
          const int numDigits = static_cast<int> (numDigitsAsDouble + 0.5);

          // Precision to which a Scalar should be written.  Add one
          // for good measure, since 17 is necessary for IEEE 754
          // doubles.
          out << std::setprecision (numDigits + 1);
        }
        out << elt.rowIndex () << " " << elt.colIndex () << " ";
        if (STS::isComplex) {
          out << STS::real (elt.value ()) << " " << STS::imag (elt.value ());
        }
        else {
          out << elt.value ();
        }
        // Restore flags
        out.flags( f );
        return out;
      }

      template<class Scalar, class Ordinal>
      class Adder {
      public:
        typedef Ordinal index_type;
        typedef Scalar value_type;
        typedef Element<Scalar, Ordinal> element_type;
        typedef typename std::vector<element_type>::size_type size_type;

        Adder () :
          expectedNumRows_(0),
          expectedNumCols_(0),
          expectedNumEntries_(0),
          seenNumRows_(0),
          seenNumCols_(0),
          seenNumEntries_(0),
          tolerant_ (true),
          debug_ (false)
        {}

        Adder (const Ordinal expectedNumRows,
               const Ordinal expectedNumCols,
               const Ordinal expectedNumEntries,
               const bool tolerant=false,
               const bool debug=false) :
          expectedNumRows_(expectedNumRows),
          expectedNumCols_(expectedNumCols),
          expectedNumEntries_(expectedNumEntries),
          seenNumRows_(0),
          seenNumCols_(0),
          seenNumEntries_(0),
          tolerant_ (tolerant),
          debug_ (debug)
        {}

        void
        operator() (const Ordinal i,
                    const Ordinal j,
                    const Scalar& Aij,
                    const bool countAgainstTotal=true)
        {
          if (! tolerant_) {
            const bool indexPairOutOfRange = i < 1 || j < 1 ||
              i > expectedNumRows_ || j > expectedNumCols_;

            TEUCHOS_TEST_FOR_EXCEPTION(indexPairOutOfRange,
              std::invalid_argument, "Matrix is " << expectedNumRows_ << " x "
              << expectedNumCols_ << ", so entry A(" << i << "," << j << ") = "
              << Aij << " is out of range.");
            if (countAgainstTotal) {
              TEUCHOS_TEST_FOR_EXCEPTION(seenNumEntries_ >= expectedNumEntries_,
                std::invalid_argument, "Cannot add entry A(" << i << "," << j
                << ") = " << Aij << " to matrix; already have expected number "
                "of entries " << expectedNumEntries_ << ".");
            }
          }
          // i and j are 1-based indices, but we store them as 0-based.
          elts_.push_back (element_type (i-1, j-1, Aij));

          // Keep track of the rightmost column containing a matrix
          // entry, and the bottommost row containing a matrix entry.
          // This gives us a lower bound for the dimensions of the
          // matrix, and a check for the reported dimensions of the
          // matrix in the Matrix Market file.
          seenNumRows_ = std::max (seenNumRows_, i);
          seenNumCols_ = std::max (seenNumCols_, j);
          if (countAgainstTotal) {
            ++seenNumEntries_;
          }
        }

        void
        print (std::ostream& out, const bool doMerge, const bool replace=false)
        {
          if (doMerge) {
            merge (replace);
          } else {
            std::sort (elts_.begin(), elts_.end());
          }
          // Print out the results, delimited by newlines.
          typedef std::ostream_iterator<element_type> iter_type;
          std::copy (elts_.begin(), elts_.end(), iter_type (out, "\n"));
        }

        std::pair<size_type, size_type>
        merge (const bool replace=false)
        {
          typedef typename std::vector<element_type>::iterator iter_type;

          // Start with a sorted container.  Element objects sort in
          // lexicographic order of their (row, column) indices, for
          // easy conversion to CSR format.  If you expect that the
          // elements will usually be sorted in the desired order, you
          // can check first whether they are already sorted.  We have
          // no such expectation, so we don't even bother to spend the
          // extra O(# entries) operations to check.
          std::sort (elts_.begin(), elts_.end());

          // Walk through the array of elements in place, merging
          // duplicates and pushing unique elements up to the front of
          // the array.  We can't use std::unique for this because it
          // doesn't let us merge duplicate elements; it only removes
          // them from the sequence.
          size_type numUnique = 0;
          iter_type cur = elts_.begin();
          if (cur == elts_.end()) { // No elements to merge
            return std::make_pair (numUnique, size_type (0));
          }
          else {
            iter_type next = cur;
            ++next; // There is one unique element
            ++numUnique;
            while (next != elts_.end()) {
              if (*cur == *next) {
                // Merge in the duplicated element *next
                cur->merge (*next, replace);
              } else {
                // *cur is already a unique element.  Move over one to
                // *make space for the new unique element.
                ++cur;
                *cur = *next; // Add the new unique element
                ++numUnique;
              }
              // Look at the "next" not-yet-considered element
              ++next;
            }
            // Remember how many elements we removed before resizing.
            const size_type numRemoved = elts_.size() - numUnique;
            elts_.resize (numUnique);
            return std::make_pair (numUnique, numRemoved);
          }
        }

        void
        mergeAndConvertToCSR (size_type& numUniqueElts,
                              size_type& numRemovedElts,
                              Teuchos::ArrayRCP<Ordinal>& rowptr,
                              Teuchos::ArrayRCP<Ordinal>& colind,
                              Teuchos::ArrayRCP<Scalar>& values,
                              const bool replace=false)
        {
          using Teuchos::arcp;
          using Teuchos::ArrayRCP;

          std::pair<size_type, size_type> mergeResult = merge (replace);

          // At this point, elts_ is already in CSR order.
          // Now we can allocate and fill the ind and val arrays.
          ArrayRCP<Ordinal> ind = arcp<Ordinal> (elts_.size ());
          ArrayRCP<Scalar> val = arcp<Scalar> (elts_.size ());

          // Number of rows in the matrix.
          const Ordinal nrows = tolerant_ ? seenNumRows_ : expectedNumRows_;
          ArrayRCP<Ordinal> ptr = arcp<Ordinal> (nrows + 1);

          // Copy over the elements, and fill in the ptr array with
          // offsets.  Note that merge() sorted the entries by row
          // index, so we can assume the row indices are increasing in
          // the list of entries.
          Ordinal curRow = 0;
          Ordinal curInd = 0;
          typedef typename std::vector<element_type>::const_iterator iter_type;
          ptr[0] = 0; // ptr always has at least one entry.
          for (iter_type it = elts_.begin(); it != elts_.end(); ++it) {
            const Ordinal i = it->rowIndex ();
            const Ordinal j = it->colIndex ();
            const Scalar Aij = it->value ();

            TEUCHOS_TEST_FOR_EXCEPTION(i < curRow, std::logic_error, "The "
              "current matrix entry's row index " << i << " is less then what "
              "should be the current row index lower bound " << curRow << ".");
            for (Ordinal k = curRow+1; k <= i; ++k) {
              ptr[k] = curInd;
            }
            curRow = i;

            TEUCHOS_TEST_FOR_EXCEPTION(
              static_cast<size_t> (curInd) >= elts_.size (),
              std::logic_error, "The current index " << curInd << " into ind "
              "and val is >= the number of matrix entries " << elts_.size ()
              << ".");
            ind[curInd] = j;
            val[curInd] = Aij;
            ++curInd;
          }
          for (Ordinal k = curRow+1; k <= nrows; ++k) {
            ptr[k] = curInd;
          }

          // Assign to outputs here, to ensure the strong exception
          // guarantee (assuming that ArrayRCP's operator= doesn't
          // throw).
          rowptr = ptr;
          colind = ind;
          values = val;
          numUniqueElts = mergeResult.first;
          numRemovedElts = mergeResult.second;
        }

        const std::vector<element_type>& getEntries() const {
          return elts_;
        }

        void clear() {
          seenNumRows_ = 0;
          seenNumCols_ = 0;
          seenNumEntries_ = 0;
          elts_.resize (0);
        }

        const Ordinal numRows() const { return seenNumRows_; }

        const Ordinal numCols() const { return seenNumCols_; }

        const Ordinal numEntries() const { return seenNumEntries_; }


      private:
        Ordinal expectedNumRows_, expectedNumCols_, expectedNumEntries_;
        Ordinal seenNumRows_, seenNumCols_, seenNumEntries_;
        bool tolerant_;
        bool debug_;

        std::vector<element_type> elts_;
      };
    } // namespace Raw
  } // namespace MatrixMarket
} // namespace Teuchos

#endif // #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp