Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
46#include "Teuchos_ArrayRCP.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 ()),
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 () :
304 seenNumRows_(0),
305 seenNumCols_(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) :
336 seenNumRows_(0),
337 seenNumCols_(0),
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
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) {
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.
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 {
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.
487 elts_.resize (numUnique);
488 return std::make_pair (numUnique, numRemoved);
489 }
490 }
491
534 void
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.
551
552 // Number of rows in the matrix.
553 const Ordinal nrows = tolerant_ ? seenNumRows_ : expectedNumRows_;
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;
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:
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.
int size(const Comm< Ordinal > &comm)
Get the number of processes in the communicator.
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.
std::vector< element_type >::size_type size_type
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.
std::vector< element_type > elts_
The actual matrix entries, stored as an array of structs.
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.
Concrete serial communicator subclass.
#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.
std::ostream & operator<<(std::ostream &out, const Element< Scalar, Ordinal > &elt)
Print out an Element to the given output stream.
This structure defines some basic traits for the ordinal field type.
This structure defines some basic traits for a scalar field type.