Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
Teuchos_MatrixMarket_Raw_Graph_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_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"
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:
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),
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 ||
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) {
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
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:
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 "raw" sparse matrix input.
void clear()
Clear all the added graph entries and reset metadata.
bool operator!=(const GraphElement &rhs)
Compare row and column indices.
const Ordinal numCols() const
Computed number of columns.
#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.
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
const std::vector< element_type > & getEntries() const
A temporary const view of the entries of the graph.
GraphElement()
Default constructor: an invalid entry of the graph.
bool operator<(const GraphElement &rhs) const
Lexicographic order first by row index, then by column index.
const Ordinal numRows() const
Computed number of rows.
Templated Parameter List class.
GraphElement(const Ordinal i, const Ordinal j)
Create a sparse graph entry at (i,j).
Ordinal colIndex() const
Column index (zero-based) of this GraphElement.
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.
std::vector< element_type > elts_
The actual matrix entries, stored as an array of structs.
GraphAdder(const Ordinal expectedNumRows, const Ordinal expectedNumCols, const Ordinal expectedNumEntries, const bool tolerant=false, const bool debug=false)
Standard constructor.
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.
Ordinal rowIndex() const
Row index (zero-based) of this GraphElement.
Matrix Market file utilities.
Reference-counted smart pointer for managing arrays.