Teuchos - Trilinos Tools Package Version of the Day
Loading...
Searching...
No Matches
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
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
60namespace Teuchos {
61 namespace MatrixMarket {
62 namespace Raw {
86 template<class Scalar, class Ordinal>
87 class Element {
88 public:
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
626 const Ordinal numEntries() const { return seenNumEntries_; }
627
628
629 private:
630 Ordinal expectedNumRows_, expectedNumCols_, expectedNumEntries_;
631 Ordinal seenNumRows_, seenNumCols_, seenNumEntries_;
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
Teuchos header file which uses auto-configuration information to include necessary C++ headers.
Templated Parameter List class.
Reference-counted smart pointer for managing arrays.
To be used with Checker for "raw" sparse matrix input.
void operator()(const Ordinal i, const Ordinal j, const Scalar &Aij, const bool countAgainstTotal=true)
Add an entry to the sparse matrix.
void print(std::ostream &out, const bool doMerge, const bool replace=false)
Print the sparse matrix data.
const std::vector< element_type > & getEntries() const
A temporary const view of the entries of the matrix.
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::pair< size_type, size_type > merge(const bool replace=false)
Merge duplicate elements.
const Ordinal numCols() const
Computed number of columns.
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.
void clear()
Clear all the added matrix entries and reset metadata.
const Ordinal numRows() const
Computed number of rows.
Stores one entry of a sparse matrix.
Scalar value() const
Value (A(rowIndex(), colIndex()) of this Element.
bool operator<(const Element &rhs) const
Lexicographic order first by row index, then by column index.
void merge(const Element &rhs, const BinaryFunction &f)
Merge rhs into this Element, using custom binary function.
Ordinal colIndex() const
Column index (zero-based) of this Element.
Ordinal rowIndex() const
Row index (zero-based) of this Element.
void merge(const Element &rhs, const bool replace=false)
Merge rhs into this Element, either by addition or replacement.
Element()
Default constructor: an invalid entry of the matrix.
bool operator==(const Element &rhs)
Ignore the matrix value for comparisons.
bool operator!=(const Element &rhs)
Ignore the matrix value for comparisons.
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.
Matrix Market file utilities.
"Raw" input of sparse matrices from Matrix Market files.
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos,...
This structure defines some basic traits for the ordinal field type.
This structure defines some basic traits for a scalar field type.