Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_MatrixMarket_Raw_Adder.hpp
Go to the documentation of this file.
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"
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:
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;
288  typedef typename std::vector<element_type>::size_type size_type;
289 
300  Adder () :
301  expectedNumRows_(0),
302  expectedNumCols_(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 ||
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) {
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
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 
626  const Ordinal numEntries() const { return seenNumEntries_; }
627 
628 
629  private:
632  bool tolerant_;
633  bool debug_;
634 
636  std::vector<element_type> elts_;
637  };
638  } // namespace Raw
639  } // namespace MatrixMarket
640 } // namespace Teuchos
641 
642 #endif // #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
std::pair< size_type, size_type > merge(const bool replace=false)
Merge duplicate elements.
std::vector< element_type > elts_
The actual matrix entries, stored as an array of structs.
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.
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.
void operator()(const Ordinal i, const Ordinal j, const Scalar &Aij, const bool countAgainstTotal=true)
Add an entry to the 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.
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.
const Ordinal numEntries() const
Computed number of columns.
Adder(const Ordinal expectedNumRows, const Ordinal expectedNumCols, const Ordinal expectedNumEntries, const bool tolerant=false, const bool debug=false)
Standard constructor.
std::vector< element_type >::size_type size_type
Reference-counted smart pointer for managing arrays.
To be used with Checker for &quot;raw&quot; sparse matrix input.