Teuchos - Trilinos Tools Package  Version of the Day
Teuchos_MatrixMarket_Raw_Adder.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 
42 #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
43 #define __Teuchos_MatrixMarket_Raw_Adder_hpp
44 
45 #include "Teuchos_ConfigDefs.hpp"
46 #include "Teuchos_ArrayRCP.hpp"
47 #include "Teuchos_CommHelpers.hpp"
49 #include "Teuchos_MatrixMarket_Banner.hpp"
50 #include "Teuchos_MatrixMarket_CoordDataReader.hpp"
51 
52 #include <algorithm>
53 #include <fstream>
54 #include <iostream>
55 #include <iterator>
56 #include <vector>
57 #include <stdexcept>
58 
59 
60 namespace Teuchos {
61  namespace MatrixMarket {
62  namespace Raw {
86  template<class Scalar, class Ordinal>
87  class Element {
88  public:
90  Element () :
91  rowIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ()),
92  colIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ()),
93  value_ (Teuchos::ScalarTraits<Scalar>::zero ())
94  {}
95 
97  Element (const Ordinal i, const Ordinal j, const Scalar& Aij) :
98  rowIndex_ (i), colIndex_ (j), value_ (Aij) {}
99 
101  bool operator== (const Element& rhs) {
102  return rowIndex_ == rhs.rowIndex_ && colIndex_ == rhs.colIndex_;
103  }
104 
106  bool operator!= (const Element& rhs) {
107  return ! (*this == rhs);
108  }
109 
111  bool operator< (const Element& rhs) const {
112  if (rowIndex_ < rhs.rowIndex_)
113  return true;
114  else if (rowIndex_ > rhs.rowIndex_)
115  return false;
116  else { // equal
117  return colIndex_ < rhs.colIndex_;
118  }
119  }
120 
126  template<class BinaryFunction>
127  void merge (const Element& rhs, const BinaryFunction& f) {
129  rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
130  std::invalid_argument,
131  "Attempt to merge elements at different locations in the sparse "
132  "matrix. The current element is at (" << rowIndex() << ", "
133  << colIndex() << ") and the element you asked me to merge with it "
134  "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << "). This "
135  "probably indicates a bug in the sparse matrix reader.");
136 
137  value_ = f (rhs.value_, value_);
138  }
139 
146  void merge (const Element& rhs, const bool replace=false) {
148  rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
149  std::invalid_argument,
150  "Attempt to merge elements at different locations in the sparse "
151  "matrix. The current element is at (" << rowIndex() << ", "
152  << colIndex() << ") and the element you asked me to merge with it "
153  "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << "). This "
154  "probably indicates a bug in the sparse matrix reader.");
155 
156  if (replace) {
157  value_ = rhs.value_;
158  }
159  else {
160  value_ += rhs.value_;
161  }
162  }
163 
165  Ordinal rowIndex() const { return rowIndex_; }
166 
168  Ordinal colIndex() const { return colIndex_; }
169 
171  Scalar value() const { return value_; }
172 
173  private:
174  Ordinal rowIndex_, colIndex_;
175  Scalar value_;
176  };
177 
188  template<class Scalar, class Ordinal>
189  std::ostream&
190  operator<< (std::ostream& out, const Element<Scalar, Ordinal>& elt)
191  {
192  typedef ScalarTraits<Scalar> STS;
193  std::ios::fmtflags f( out.flags() );
194  // Non-Ordinal types are floating-point types. In order not to
195  // lose information when we print a floating-point type, we have
196  // to set the number of digits to print. C++ standard behavior
197  // in the default locale seems to be to print only five decimal
198  // digits after the decimal point; this does not suffice for
199  // double precision. We solve the problem of how many digits to
200  // print more generally below. It's a rough solution so please
201  // feel free to audit and revise it.
202  //
203  // FIXME (mfh 01 Feb 2011)
204  // This really calls for the following approach:
205  //
206  // Guy L. Steele and Jon L. White, "How to print floating-point
207  // numbers accurately", 20 Years of the ACM/SIGPLAN Conference
208  // on Programming Language Design and Implementation
209  // (1979-1999): A Selection, 2003.
210  if (! STS::isOrdinal) {
211  // std::scientific, std::fixed, and default are the three
212  // output states for floating-point numbers. A reasonable
213  // user-defined floating-point type should respect these
214  // flags; hopefully it does.
215  out << std::scientific;
216 
217  // Decimal output is standard for Matrix Market format.
218  out << std::setbase (10);
219 
220  // Compute the number of decimal digits required for expressing
221  // a Scalar, by comparing with IEEE 754 double precision (16
222  // decimal digits, 53 binary digits). This would be easier if
223  // Teuchos exposed std::numeric_limits<T>::digits10, alas.
224  const double numDigitsAsDouble =
225  16 * ((double) STS::t() / (double) ScalarTraits<double>::t());
226  // Adding 0.5 and truncating is a portable "floor".
227  const int numDigits = static_cast<int> (numDigitsAsDouble + 0.5);
228 
229  // Precision to which a Scalar should be written. Add one
230  // for good measure, since 17 is necessary for IEEE 754
231  // doubles.
232  out << std::setprecision (numDigits + 1);
233  }
234  out << elt.rowIndex () << " " << elt.colIndex () << " ";
235  if (STS::isComplex) {
236  out << STS::real (elt.value ()) << " " << STS::imag (elt.value ());
237  }
238  else {
239  out << elt.value ();
240  }
241  // Restore flags
242  out.flags( f );
243  return out;
244  }
245 
282  template<class Scalar, class Ordinal>
283  class Adder {
284  public:
285  typedef Ordinal index_type;
286  typedef Scalar value_type;
288  typedef typename std::vector<element_type>::size_type size_type;
289 
300  Adder () :
301  expectedNumRows_(0),
302  expectedNumCols_(0),
303  expectedNumEntries_(0),
304  seenNumRows_(0),
305  seenNumCols_(0),
306  seenNumEntries_(0),
307  tolerant_ (true),
308  debug_ (false)
309  {}
310 
328  Adder (const Ordinal expectedNumRows,
329  const Ordinal expectedNumCols,
330  const Ordinal expectedNumEntries,
331  const bool tolerant=false,
332  const bool debug=false) :
333  expectedNumRows_(expectedNumRows),
334  expectedNumCols_(expectedNumCols),
335  expectedNumEntries_(expectedNumEntries),
336  seenNumRows_(0),
337  seenNumCols_(0),
338  seenNumEntries_(0),
339  tolerant_ (tolerant),
340  debug_ (debug)
341  {}
342 
363  void
364  operator() (const Ordinal i,
365  const Ordinal j,
366  const Scalar& Aij,
367  const bool countAgainstTotal=true)
368  {
369  if (! tolerant_) {
370  const bool indexPairOutOfRange = i < 1 || j < 1 ||
371  i > expectedNumRows_ || j > expectedNumCols_;
372 
373  TEUCHOS_TEST_FOR_EXCEPTION(indexPairOutOfRange,
374  std::invalid_argument, "Matrix is " << expectedNumRows_ << " x "
375  << expectedNumCols_ << ", so entry A(" << i << "," << j << ") = "
376  << Aij << " is out of range.");
377  if (countAgainstTotal) {
378  TEUCHOS_TEST_FOR_EXCEPTION(seenNumEntries_ >= expectedNumEntries_,
379  std::invalid_argument, "Cannot add entry A(" << i << "," << j
380  << ") = " << Aij << " to matrix; already have expected number "
381  "of entries " << expectedNumEntries_ << ".");
382  }
383  }
384  // i and j are 1-based indices, but we store them as 0-based.
385  elts_.push_back (element_type (i-1, j-1, Aij));
386 
387  // Keep track of the rightmost column containing a matrix
388  // entry, and the bottommost row containing a matrix entry.
389  // This gives us a lower bound for the dimensions of the
390  // matrix, and a check for the reported dimensions of the
391  // matrix in the Matrix Market file.
392  seenNumRows_ = std::max (seenNumRows_, i);
393  seenNumCols_ = std::max (seenNumCols_, j);
394  if (countAgainstTotal) {
395  ++seenNumEntries_;
396  }
397  }
398 
408  void
409  print (std::ostream& out, const bool doMerge, const bool replace=false)
410  {
411  if (doMerge) {
412  merge (replace);
413  } else {
414  std::sort (elts_.begin(), elts_.end());
415  }
416  // Print out the results, delimited by newlines.
417  typedef std::ostream_iterator<element_type> iter_type;
418  std::copy (elts_.begin(), elts_.end(), iter_type (out, "\n"));
419  }
420 
443  std::pair<size_type, size_type>
444  merge (const bool replace=false)
445  {
446  typedef typename std::vector<element_type>::iterator iter_type;
447 
448  // Start with a sorted container. Element objects sort in
449  // lexicographic order of their (row, column) indices, for
450  // easy conversion to CSR format. If you expect that the
451  // elements will usually be sorted in the desired order, you
452  // can check first whether they are already sorted. We have
453  // no such expectation, so we don't even bother to spend the
454  // extra O(# entries) operations to check.
455  std::sort (elts_.begin(), elts_.end());
456 
457  // Walk through the array of elements in place, merging
458  // duplicates and pushing unique elements up to the front of
459  // the array. We can't use std::unique for this because it
460  // doesn't let us merge duplicate elements; it only removes
461  // them from the sequence.
462  size_type numUnique = 0;
463  iter_type cur = elts_.begin();
464  if (cur == elts_.end()) { // No elements to merge
465  return std::make_pair (numUnique, size_type (0));
466  }
467  else {
468  iter_type next = cur;
469  ++next; // There is one unique element
470  ++numUnique;
471  while (next != elts_.end()) {
472  if (*cur == *next) {
473  // Merge in the duplicated element *next
474  cur->merge (*next, replace);
475  } else {
476  // *cur is already a unique element. Move over one to
477  // *make space for the new unique element.
478  ++cur;
479  *cur = *next; // Add the new unique element
480  ++numUnique;
481  }
482  // Look at the "next" not-yet-considered element
483  ++next;
484  }
485  // Remember how many elements we removed before resizing.
486  const size_type numRemoved = elts_.size() - numUnique;
487  elts_.resize (numUnique);
488  return std::make_pair (numUnique, numRemoved);
489  }
490  }
491 
534  void
535  mergeAndConvertToCSR (size_type& numUniqueElts,
536  size_type& numRemovedElts,
540  const bool replace=false)
541  {
542  using Teuchos::arcp;
543  using Teuchos::ArrayRCP;
544 
545  std::pair<size_type, size_type> mergeResult = merge (replace);
546 
547  // At this point, elts_ is already in CSR order.
548  // Now we can allocate and fill the ind and val arrays.
549  ArrayRCP<Ordinal> ind = arcp<Ordinal> (elts_.size ());
550  ArrayRCP<Scalar> val = arcp<Scalar> (elts_.size ());
551 
552  // Number of rows in the matrix.
553  const Ordinal nrows = tolerant_ ? seenNumRows_ : expectedNumRows_;
554  ArrayRCP<Ordinal> ptr = arcp<Ordinal> (nrows + 1);
555 
556  // Copy over the elements, and fill in the ptr array with
557  // offsets. Note that merge() sorted the entries by row
558  // index, so we can assume the row indices are increasing in
559  // the list of entries.
560  Ordinal curRow = 0;
561  Ordinal curInd = 0;
562  typedef typename std::vector<element_type>::const_iterator iter_type;
563  ptr[0] = 0; // ptr always has at least one entry.
564  for (iter_type it = elts_.begin(); it != elts_.end(); ++it) {
565  const Ordinal i = it->rowIndex ();
566  const Ordinal j = it->colIndex ();
567  const Scalar Aij = it->value ();
568 
569  TEUCHOS_TEST_FOR_EXCEPTION(i < curRow, std::logic_error, "The "
570  "current matrix entry's row index " << i << " is less then what "
571  "should be the current row index lower bound " << curRow << ".");
572  for (Ordinal k = curRow+1; k <= i; ++k) {
573  ptr[k] = curInd;
574  }
575  curRow = i;
576 
578  static_cast<size_t> (curInd) >= elts_.size (),
579  std::logic_error, "The current index " << curInd << " into ind "
580  "and val is >= the number of matrix entries " << elts_.size ()
581  << ".");
582  ind[curInd] = j;
583  val[curInd] = Aij;
584  ++curInd;
585  }
586  for (Ordinal k = curRow+1; k <= nrows; ++k) {
587  ptr[k] = curInd;
588  }
589 
590  // Assign to outputs here, to ensure the strong exception
591  // guarantee (assuming that ArrayRCP's operator= doesn't
592  // throw).
593  rowptr = ptr;
594  colind = ind;
595  values = val;
596  numUniqueElts = mergeResult.first;
597  numRemovedElts = mergeResult.second;
598  }
599 
601  const std::vector<element_type>& getEntries() const {
602  return elts_;
603  }
604 
606  void clear() {
607  seenNumRows_ = 0;
608  seenNumCols_ = 0;
609  seenNumEntries_ = 0;
610  elts_.resize (0);
611  }
612 
616  const Ordinal numRows() const { return seenNumRows_; }
617 
621  const Ordinal numCols() const { return seenNumCols_; }
622 
623  private:
624  Ordinal expectedNumRows_, expectedNumCols_, expectedNumEntries_;
625  Ordinal seenNumRows_, seenNumCols_, seenNumEntries_;
626  bool tolerant_;
627  bool debug_;
628 
630  std::vector<element_type> elts_;
631  };
632  } // namespace Raw
633  } // namespace MatrixMarket
634 } // namespace Teuchos
635 
636 #endif // #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
std::pair< size_type, size_type > merge(const bool replace=false)
Merge duplicate elements.
const std::vector< element_type > & getEntries() const
A temporary const view of the entries of the matrix.
bool operator<(const Element &rhs) const
Lexicographic order first by row index, then by column index.
Element(const Ordinal i, const Ordinal j, const Scalar &Aij)
Create a sparse matrix entry at (i,j) with value Aij.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
"Raw" input of sparse matrices from Matrix Market files.
Scalar value() const
Value (A(rowIndex(), colIndex()) of this Element.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
void clear()
Clear all the added matrix entries and reset metadata.
void mergeAndConvertToCSR(size_type &numUniqueElts, size_type &numRemovedElts, Teuchos::ArrayRCP< Ordinal > &rowptr, Teuchos::ArrayRCP< Ordinal > &colind, Teuchos::ArrayRCP< Scalar > &values, const bool replace=false)
Merge duplicate elements and convert to zero-based CSR.
Element()
Default constructor: an invalid entry of the matrix.
This structure defines some basic traits for a scalar field type.
Templated Parameter List class.
Ordinal rowIndex() const
Row index (zero-based) of this Element.
bool operator==(const Element &rhs)
Ignore the matrix value for comparisons.
This structure defines some basic traits for the ordinal field type.
const Ordinal numRows() const
Computed number of rows.
void print(std::ostream &out, const bool doMerge, const bool replace=false)
Print the sparse matrix data.
Stores one entry of a sparse matrix.
const Ordinal numCols() const
Computed number of columns.
void merge(const Element &rhs, const BinaryFunction &f)
Merge rhs into this Element, using custom binary function.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
bool operator!=(const Element &rhs)
Ignore the matrix value for comparisons.
void merge(const Element &rhs, const bool replace=false)
Merge rhs into this Element, either by addition or replacement.
Ordinal colIndex() const
Column index (zero-based) of this Element.
Adder(const Ordinal expectedNumRows, const Ordinal expectedNumCols, const Ordinal expectedNumEntries, const bool tolerant=false, const bool debug=false)
Standard constructor.
Matrix Market file utilities.
Reference-counted smart pointer for managing arrays.
To be used with Checker for "raw" sparse matrix input.