Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_MatrixMarket_Raw_Graph_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_Graph_Adder_hpp
43 #define __Teuchos_MatrixMarket_Raw_Graph_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 {
85  template<class Ordinal>
86  class GraphElement {
87  public:
90  rowIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ()),
91  colIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ())
92  {}
93 
95  GraphElement (const Ordinal i, const Ordinal j) :
96  rowIndex_ (i), colIndex_ (j) {}
97 
99  bool operator== (const GraphElement& rhs) {
100  return rowIndex_ == rhs.rowIndex_ && colIndex_ == rhs.colIndex_;
101  }
102 
104  bool operator!= (const GraphElement& rhs) {
105  return ! (*this == rhs);
106  }
107 
109  bool operator< (const GraphElement& rhs) const {
110  if (rowIndex_ < rhs.rowIndex_)
111  return true;
112  else if (rowIndex_ > rhs.rowIndex_)
113  return false;
114  else { // equal
115  return colIndex_ < rhs.colIndex_;
116  }
117  }
118 
120  Ordinal rowIndex() const { return rowIndex_; }
121 
123  Ordinal colIndex() const { return colIndex_; }
124 
125  private:
126  Ordinal rowIndex_, colIndex_;
127  };
128 
133  template<class Ordinal>
134  std::ostream&
135  operator<< (std::ostream& out, const GraphElement<Ordinal>& elt)
136  {
137  out << elt.rowIndex () << " " << elt.colIndex ();
138  return out;
139  }
140 
160  template<class Ordinal>
161  class GraphAdder {
162  public:
163  typedef Ordinal index_type;
165  typedef typename std::vector<element_type>::size_type size_type;
166 
178  expectedNumRows_(0),
179  expectedNumCols_(0),
180  expectedNumEntries_(0),
181  seenNumRows_(0),
182  seenNumCols_(0),
183  seenNumEntries_(0),
184  tolerant_ (true),
185  debug_ (false)
186  {}
187 
205  GraphAdder (const Ordinal expectedNumRows,
206  const Ordinal expectedNumCols,
207  const Ordinal expectedNumEntries,
208  const bool tolerant=false,
209  const bool debug=false) :
210  expectedNumRows_(expectedNumRows),
211  expectedNumCols_(expectedNumCols),
212  expectedNumEntries_(expectedNumEntries),
213  seenNumRows_(0),
214  seenNumCols_(0),
215  seenNumEntries_(0),
216  tolerant_ (tolerant),
217  debug_ (debug)
218  {}
219 
240  void
241  operator() (const Ordinal i,
242  const Ordinal j,
243  const bool countAgainstTotal=true)
244  {
245  if (! tolerant_) {
246  const bool indexPairOutOfRange = i < 1 || j < 1 ||
247  i > expectedNumRows_ || j > expectedNumCols_;
248 
249  TEUCHOS_TEST_FOR_EXCEPTION(indexPairOutOfRange,
250  std::invalid_argument, "Graph is " << expectedNumRows_ << " x "
251  << expectedNumCols_ << ", so entry A(" << i << "," << j
252  << ") is out of range.");
253  if (countAgainstTotal) {
254  TEUCHOS_TEST_FOR_EXCEPTION(seenNumEntries_ >= expectedNumEntries_,
255  std::invalid_argument, "Cannot add entry A(" << i << "," << j
256  << ") to graph; already have expected number "
257  "of entries " << expectedNumEntries_ << ".");
258  }
259  }
260  // i and j are 1-based indices, but we store them as 0-based.
261  elts_.push_back (element_type (i-1, j-1));
262 
263  // Keep track of the rightmost column containing a matrix
264  // entry, and the bottommost row containing a matrix entry.
265  // This gives us a lower bound for the dimensions of the
266  // graph, and a check for the reported dimensions of the
267  // graph in the Matrix Market file.
268  seenNumRows_ = std::max (seenNumRows_, i);
269  seenNumCols_ = std::max (seenNumCols_, j);
270  if (countAgainstTotal) {
271  ++seenNumEntries_;
272  }
273  }
274 
291  void
292  print (std::ostream& out, const bool doMerge, const bool replace=false)
293  {
294  if (doMerge) {
296  (replace, std::logic_error, "replace = true not implemented!");
297  //merge (replace);
298  merge ();
299  } else {
300  std::sort (elts_.begin(), elts_.end());
301  }
302  // Print out the results, delimited by newlines.
303  typedef std::ostream_iterator<element_type> iter_type;
304  std::copy (elts_.begin(), elts_.end(), iter_type (out, "\n"));
305  }
306 
320  std::pair<size_type, size_type>
321  merge ()
322  {
323  typedef typename std::vector<element_type>::iterator iter_type;
324 
325  // Start with a sorted container. GraphElement objects sort in
326  // lexicographic order of their (row, column) indices, for
327  // easy conversion to CSR format. If you expect that the
328  // elements will usually be sorted in the desired order, you
329  // can check first whether they are already sorted. We have
330  // no such expectation, so we don't even bother to spend the
331  // extra O(# entries) operations to check.
332  std::sort (elts_.begin(), elts_.end());
333 
334  // Remove duplicate elements from the sequence
335  iter_type it;
336  it = std::unique(elts_.begin(), elts_.end());
337  size_type numUnique = std::distance(elts_.begin(),it);
338  const size_type numRemoved = elts_.size() - numUnique;
339  elts_.resize( std::distance(elts_.begin(),it) );
340  elts_.resize (numUnique);
341  return std::make_pair (numUnique, numRemoved);
342  }
343 
375  void
376  mergeAndConvertToCSR (size_type& numUniqueElts,
377  size_type& numRemovedElts,
380  {
381  using Teuchos::arcp;
382  using Teuchos::ArrayRCP;
383 
384  std::pair<size_type, size_type> mergeResult = merge();
385 
386  // At this point, elts_ is already in CSR order.
387  // Now we can allocate and fill the ind array.
388  ArrayRCP<Ordinal> ind = arcp<Ordinal> (elts_.size ());
389 
390  // Number of rows in the graph.
391  const Ordinal nrows = tolerant_ ? seenNumRows_ : expectedNumRows_;
392  ArrayRCP<Ordinal> ptr = arcp<Ordinal> (nrows + 1);
393 
394  // Copy over the elements, and fill in the ptr array with
395  // offsets. Note that merge() sorted the entries by row
396  // index, so we can assume the row indices are increasing in
397  // the list of entries.
398  Ordinal curRow = 0;
399  Ordinal curInd = 0;
400  typedef typename std::vector<element_type>::const_iterator iter_type;
401  ptr[0] = 0; // ptr always has at least one entry.
402  for (iter_type it = elts_.begin(); it != elts_.end(); ++it) {
403  const Ordinal i = it->rowIndex ();
404  const Ordinal j = it->colIndex ();
405 
406  TEUCHOS_TEST_FOR_EXCEPTION(i < curRow, std::logic_error, "The "
407  "current graph entry's row index " << i << " is less then what "
408  "should be the current row index lower bound " << curRow << ".");
409  for (Ordinal k = curRow+1; k <= i; ++k) {
410  ptr[k] = curInd;
411  }
412  curRow = i;
413 
415  static_cast<size_t> (curInd) >= elts_.size (),
416  std::logic_error, "The current index " << curInd << " into ind "
417  "is >= the number of matrix entries " << elts_.size ()
418  << ".");
419  ind[curInd] = j;
420  ++curInd;
421  }
422  for (Ordinal k = curRow+1; k <= nrows; ++k) {
423  ptr[k] = curInd;
424  }
425 
426  // Assign to outputs here, to ensure the strong exception
427  // guarantee (assuming that ArrayRCP's operator= doesn't
428  // throw).
429  rowptr = ptr;
430  colind = ind;
431  numUniqueElts = mergeResult.first;
432  numRemovedElts = mergeResult.second;
433  }
434 
436  const std::vector<element_type>& getEntries() const {
437  return elts_;
438  }
439 
441  void clear() {
442  seenNumRows_ = 0;
443  seenNumCols_ = 0;
444  seenNumEntries_ = 0;
445  elts_.resize (0);
446  }
447 
451  const Ordinal numRows() const { return seenNumRows_; }
452 
456  const Ordinal numCols() const { return seenNumCols_; }
457 
458  private:
459  Ordinal expectedNumRows_, expectedNumCols_, expectedNumEntries_;
460  Ordinal seenNumRows_, seenNumCols_, seenNumEntries_;
461  bool tolerant_;
462  bool debug_;
463 
465  std::vector<element_type> elts_;
466  };
467  } // namespace Raw
468  } // namespace MatrixMarket
469 } // namespace Teuchos
470 
471 #endif // #ifndef __Teuchos_MatrixMarket_Raw_Graph_Adder_hpp
bool operator==(const GraphElement &rhs)
Compare row and column indices.
To be used with Checker for &quot;raw&quot; sparse matrix input.
bool operator<(const GraphElement &rhs) const
Lexicographic order first by row index, then by column index.
void clear()
Clear all the added graph entries and reset metadata.
bool operator!=(const GraphElement &rhs)
Compare row and column indices.
const std::vector< element_type > & getEntries() const
A temporary const view of the entries of the graph.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
const Ordinal numRows() const
Computed number of rows.
const Ordinal numCols() const
Computed number of columns.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
GraphElement()
Default constructor: an invalid entry of the graph.
Templated Parameter List class.
Ordinal colIndex() const
Column index (zero-based) of this GraphElement.
GraphElement(const Ordinal i, const Ordinal j)
Create a sparse graph entry at (i,j).
void operator()(const Ordinal i, const Ordinal j, const bool countAgainstTotal=true)
Add an entry to the sparse graph.
This structure defines some basic traits for the ordinal field type.
void print(std::ostream &out, const bool doMerge, const bool replace=false)
Print the sparse graph data.
GraphAdder(const Ordinal expectedNumRows, const Ordinal expectedNumCols, const Ordinal expectedNumEntries, const bool tolerant=false, const bool debug=false)
Standard constructor.
Ordinal rowIndex() const
Row index (zero-based) of this GraphElement.
void mergeAndConvertToCSR(size_type &numUniqueElts, size_type &numRemovedElts, Teuchos::ArrayRCP< Ordinal > &rowptr, Teuchos::ArrayRCP< Ordinal > &colind)
Merge duplicate elements and convert to zero-based CSR.
std::pair< size_type, size_type > merge()
Merge duplicate elements.
Reference-counted smart pointer for managing arrays.