Teuchos - Trilinos Tools Package Version of the Day
Teuchos_MatrixMarket_Raw_Adder.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
00043 #define __Teuchos_MatrixMarket_Raw_Adder_hpp
00044 
00045 #include "Teuchos_ConfigDefs.hpp"
00046 #include "Teuchos_ArrayRCP.hpp"
00047 #include "Teuchos_CommHelpers.hpp"
00048 #include "Teuchos_ParameterList.hpp"
00049 #include "Teuchos_MatrixMarket_Banner.hpp"
00050 #include "Teuchos_MatrixMarket_CoordDataReader.hpp"
00051 
00052 #include <algorithm>
00053 #include <fstream>
00054 #include <iostream>
00055 #include <iterator>
00056 #include <vector>
00057 #include <stdexcept>
00058 
00059 
00060 namespace Teuchos {
00061   namespace MatrixMarket {
00062     namespace Raw {
00086       template<class Scalar, class Ordinal>
00087       class Element {
00088       public:
00096         Element () : rowIndex_ (-1), colIndex_ (-1), value_ (0) {}
00097 
00099         Element (const Ordinal i, const Ordinal j, const Scalar& Aij) :
00100           rowIndex_ (i), colIndex_ (j), value_ (Aij) {}
00101 
00103         bool operator== (const Element& rhs) {
00104           return rowIndex_ == rhs.rowIndex_ && colIndex_ == rhs.colIndex_;
00105         }
00106 
00108         bool operator!= (const Element& rhs) {
00109           return ! (*this == rhs);
00110         }
00111 
00113         bool operator< (const Element& rhs) const {
00114           if (rowIndex_ < rhs.rowIndex_)
00115             return true;
00116           else if (rowIndex_ > rhs.rowIndex_)
00117             return false;
00118           else { // equal
00119             return colIndex_ < rhs.colIndex_;
00120           }
00121         }
00122 
00128         template<class BinaryFunction>
00129         void merge (const Element& rhs, const BinaryFunction& f) {
00130           TEUCHOS_TEST_FOR_EXCEPTION(
00131             rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
00132             std::invalid_argument,
00133             "Attempt to merge elements at different locations in the sparse "
00134             "matrix.  The current element is at (" << rowIndex() << ", "
00135             << colIndex() << ") and the element you asked me to merge with it "
00136             "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << ").  This "
00137             "probably indicates a bug in the sparse matrix reader.");
00138 
00139           value_ = f (rhs.value_, value_);
00140         }
00141 
00148         void merge (const Element& rhs, const bool replace=false) {
00149           TEUCHOS_TEST_FOR_EXCEPTION(
00150             rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
00151             std::invalid_argument,
00152             "Attempt to merge elements at different locations in the sparse "
00153             "matrix.  The current element is at (" << rowIndex() << ", "
00154             << colIndex() << ") and the element you asked me to merge with it "
00155             "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << ").  This "
00156             "probably indicates a bug in the sparse matrix reader.");
00157 
00158           if (replace) {
00159             value_ = rhs.value_;
00160           }
00161           else {
00162             value_ += rhs.value_;
00163           }
00164         }
00165 
00167         Ordinal rowIndex() const { return rowIndex_; }
00168 
00170         Ordinal colIndex() const { return colIndex_; }
00171 
00173         Scalar value() const { return value_; }
00174 
00175       private:
00176         Ordinal rowIndex_, colIndex_;
00177         Scalar value_;
00178       };
00179 
00190       template<class Scalar, class Ordinal>
00191       std::ostream&
00192       operator<< (std::ostream& out, const Element<Scalar, Ordinal>& elt)
00193       {
00194         typedef ScalarTraits<Scalar> STS;
00195         // Non-Ordinal types are floating-point types.  In order not to
00196         // lose information when we print a floating-point type, we have
00197         // to set the number of digits to print.  C++ standard behavior
00198         // in the default locale seems to be to print only five decimal
00199         // digits after the decimal point; this does not suffice for
00200         // double precision.  We solve the problem of how many digits to
00201         // print more generally below.  It's a rough solution so please
00202         // feel free to audit and revise it.
00203         //
00204         // FIXME (mfh 01 Feb 2011)
00205         // This really calls for the following approach:
00206         //
00207         // Guy L. Steele and Jon L. White, "How to print floating-point
00208         // numbers accurately", 20 Years of the ACM/SIGPLAN Conference
00209         // on Programming Language Design and Implementation
00210         // (1979-1999): A Selection, 2003.
00211         if (! STS::isOrdinal) {
00212           // std::scientific, std::fixed, and default are the three
00213           // output states for floating-point numbers.  A reasonable
00214           // user-defined floating-point type should respect these
00215           // flags; hopefully it does.
00216           out << std::scientific;
00217 
00218           // Decimal output is standard for Matrix Market format.
00219           out << std::setbase (10);
00220 
00221           // Compute the number of decimal digits required for expressing
00222           // a Scalar, by comparing with IEEE 754 double precision (16
00223           // decimal digits, 53 binary digits).  This would be easier if
00224           // Teuchos exposed std::numeric_limits<T>::digits10, alas.
00225           const double numDigitsAsDouble =
00226             16 * ((double) STS::t() / (double) ScalarTraits<double>::t());
00227           // Adding 0.5 and truncating is a portable "floor".
00228           const int numDigits = static_cast<int> (numDigitsAsDouble + 0.5);
00229 
00230           // Precision to which a Scalar should be written.  Add one
00231           // for good measure, since 17 is necessary for IEEE 754
00232           // doubles.
00233           out << std::setprecision (numDigits + 1);
00234         }
00235         out << elt.rowIndex () << " " << elt.colIndex () << " ";
00236         if (STS::isComplex) {
00237           out << STS::real (elt.value ()) << " " << STS::imag (elt.value ());
00238         }
00239         else {
00240           out << elt.value ();
00241         }
00242         return out;
00243       }
00244 
00281       template<class Scalar, class Ordinal>
00282       class Adder {
00283       public:
00284         typedef Ordinal index_type;
00285         typedef Scalar value_type;
00286         typedef Element<Scalar, Ordinal> element_type;
00287         typedef typename std::vector<element_type>::size_type size_type;
00288 
00299         Adder () :
00300           expectedNumRows_(0),
00301           expectedNumCols_(0),
00302           expectedNumEntries_(0),
00303           seenNumRows_(0),
00304           seenNumCols_(0),
00305           seenNumEntries_(0),
00306           tolerant_ (true),
00307           debug_ (false)
00308         {}
00309 
00327         Adder (const Ordinal expectedNumRows,
00328                const Ordinal expectedNumCols,
00329                const Ordinal expectedNumEntries,
00330                const bool tolerant=false,
00331                const bool debug=false) :
00332           expectedNumRows_(expectedNumRows),
00333           expectedNumCols_(expectedNumCols),
00334           expectedNumEntries_(expectedNumEntries),
00335           seenNumRows_(0),
00336           seenNumCols_(0),
00337           seenNumEntries_(0),
00338           tolerant_ (tolerant),
00339           debug_ (debug)
00340         {}
00341 
00362         void
00363         operator() (const Ordinal i,
00364                     const Ordinal j,
00365                     const Scalar& Aij,
00366                     const bool countAgainstTotal=true)
00367         {
00368           if (! tolerant_) {
00369             const bool indexPairOutOfRange = i < 1 || j < 1 ||
00370               i > expectedNumRows_ || j > expectedNumCols_;
00371 
00372             TEUCHOS_TEST_FOR_EXCEPTION(indexPairOutOfRange,
00373               std::invalid_argument, "Matrix is " << expectedNumRows_ << " x "
00374               << expectedNumCols_ << ", so entry A(" << i << "," << j << ") = "
00375               << Aij << " is out of range.");
00376             if (countAgainstTotal) {
00377               TEUCHOS_TEST_FOR_EXCEPTION(seenNumEntries_ >= expectedNumEntries_,
00378                 std::invalid_argument, "Cannot add entry A(" << i << "," << j
00379                 << ") = " << Aij << " to matrix; already have expected number "
00380                 "of entries " << expectedNumEntries_ << ".");
00381             }
00382           }
00383           // i and j are 1-based indices, but we store them as 0-based.
00384           elts_.push_back (element_type (i-1, j-1, Aij));
00385 
00386           // Keep track of the rightmost column containing a matrix
00387           // entry, and the bottommost row containing a matrix entry.
00388           // This gives us a lower bound for the dimensions of the
00389           // matrix, and a check for the reported dimensions of the
00390           // matrix in the Matrix Market file.
00391           seenNumRows_ = std::max (seenNumRows_, i);
00392           seenNumCols_ = std::max (seenNumCols_, j);
00393           if (countAgainstTotal) {
00394             ++seenNumEntries_;
00395           }
00396         }
00397 
00407         void
00408         print (std::ostream& out, const bool doMerge, const bool replace=false)
00409         {
00410           if (doMerge) {
00411             merge (replace);
00412           } else {
00413             std::sort (elts_.begin(), elts_.end());
00414           }
00415           // Print out the results, delimited by newlines.
00416           typedef std::ostream_iterator<element_type> iter_type;
00417           std::copy (elts_.begin(), elts_.end(), iter_type (out, "\n"));
00418         }
00419 
00442         std::pair<size_type, size_type>
00443         merge (const bool replace=false)
00444         {
00445           typedef typename std::vector<element_type>::iterator iter_type;
00446 
00447           // Start with a sorted container.  Element objects sort in
00448           // lexicographic order of their (row, column) indices, for
00449           // easy conversion to CSR format.  If you expect that the
00450           // elements will usually be sorted in the desired order, you
00451           // can check first whether they are already sorted.  We have
00452           // no such expectation, so we don't even bother to spend the
00453           // extra O(# entries) operations to check.
00454           std::sort (elts_.begin(), elts_.end());
00455 
00456           // Walk through the array of elements in place, merging
00457           // duplicates and pushing unique elements up to the front of
00458           // the array.  We can't use std::unique for this because it
00459           // doesn't let us merge duplicate elements; it only removes
00460           // them from the sequence.
00461           size_type numUnique = 0;
00462           iter_type cur = elts_.begin();
00463           if (cur == elts_.end()) { // No elements to merge
00464             return std::make_pair (numUnique, size_type (0));
00465           }
00466           else {
00467             iter_type next = cur;
00468             ++next; // There is one unique element
00469             ++numUnique;
00470             while (next != elts_.end()) {
00471               if (*cur == *next) {
00472                 // Merge in the duplicated element *next
00473                 cur->merge (*next, replace);
00474               } else {
00475                 // *cur is already a unique element.  Move over one to
00476                 // *make space for the new unique element.
00477                 ++cur;
00478                 *cur = *next; // Add the new unique element
00479                 ++numUnique;
00480               }
00481               // Look at the "next" not-yet-considered element
00482               ++next;
00483             }
00484             // Remember how many elements we removed before resizing.
00485             const size_type numRemoved = elts_.size() - numUnique;
00486             elts_.resize (numUnique);
00487             return std::make_pair (numUnique, numRemoved);
00488           }
00489         }
00490 
00533         void
00534         mergeAndConvertToCSR (size_type& numUniqueElts,
00535                               size_type& numRemovedElts,
00536                               Teuchos::ArrayRCP<Ordinal>& rowptr,
00537                               Teuchos::ArrayRCP<Ordinal>& colind,
00538                               Teuchos::ArrayRCP<Scalar>& values,
00539                               const bool replace=false)
00540         {
00541           using Teuchos::arcp;
00542           using Teuchos::ArrayRCP;
00543 
00544           std::pair<size_type, size_type> mergeResult = merge (replace);
00545 
00546           // At this point, elts_ is already in CSR order.
00547           // Now we can allocate and fill the ind and val arrays.
00548           ArrayRCP<Ordinal> ind = arcp<Ordinal> (elts_.size ());
00549           ArrayRCP<Scalar> val = arcp<Scalar> (elts_.size ());
00550 
00551           // Number of rows in the matrix.
00552           const Ordinal nrows = tolerant_ ? seenNumRows_ : expectedNumRows_;
00553           ArrayRCP<Ordinal> ptr = arcp<Ordinal> (nrows + 1);
00554 
00555           // Copy over the elements, and fill in the ptr array with
00556           // offsets.  Note that merge() sorted the entries by row
00557           // index, so we can assume the row indices are increasing in
00558           // the list of entries.
00559           Ordinal curRow = 0;
00560           Ordinal curInd = 0;
00561           typedef typename std::vector<element_type>::const_iterator iter_type;
00562           ptr[0] = 0; // ptr always has at least one entry.
00563           for (iter_type it = elts_.begin(); it != elts_.end(); ++it) {
00564             const Ordinal i = it->rowIndex ();
00565             const Ordinal j = it->colIndex ();
00566             const Scalar Aij = it->value ();
00567 
00568             TEUCHOS_TEST_FOR_EXCEPTION(i < curRow, std::logic_error, "The "
00569               "current matrix entry's row index " << i << " is less then what "
00570               "should be the current row index lower bound " << curRow << ".");
00571             for (Ordinal k = curRow+1; k <= i; ++k) {
00572               ptr[k] = curInd;
00573             }
00574             curRow = i;
00575 
00576             TEUCHOS_TEST_FOR_EXCEPTION(
00577               static_cast<size_t> (curInd) >= elts_.size (),
00578               std::logic_error, "The current index " << curInd << " into ind "
00579               "and val is >= the number of matrix entries " << elts_.size ()
00580               << ".");
00581             ind[curInd] = j;
00582             val[curInd] = Aij;
00583             ++curInd;
00584           }
00585           for (Ordinal k = curRow+1; k <= nrows; ++k) {
00586             ptr[k] = curInd;
00587           }
00588 
00589           // Assign to outputs here, to ensure the strong exception
00590           // guarantee (assuming that ArrayRCP's operator= doesn't
00591           // throw).
00592           rowptr = ptr;
00593           colind = ind;
00594           values = val;
00595           numUniqueElts = mergeResult.first;
00596           numRemovedElts = mergeResult.second;
00597         }
00598 
00600         const std::vector<element_type>& getEntries() const {
00601           return elts_;
00602         }
00603 
00605         void clear() {
00606           seenNumRows_ = 0;
00607           seenNumCols_ = 0;
00608           seenNumEntries_ = 0;
00609           elts_.resize (0);
00610         }
00611 
00615         const Ordinal numRows() const { return seenNumRows_; }
00616 
00620         const Ordinal numCols() const { return seenNumCols_; }
00621 
00622       private:
00623         Ordinal expectedNumRows_, expectedNumCols_, expectedNumEntries_;
00624         Ordinal seenNumRows_, seenNumCols_, seenNumEntries_;
00625         bool tolerant_;
00626         bool debug_;
00627 
00629         std::vector<element_type> elts_;
00630       };
00631     } // namespace Raw
00632   } // namespace MatrixMarket
00633 } // namespace Teuchos
00634 
00635 #endif // #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines