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:
00090         Element () :
00091           rowIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ()),
00092           colIndex_ (Teuchos::OrdinalTraits<Ordinal>::invalid ()),
00093           value_ (Teuchos::ScalarTraits<Scalar>::zero ())
00094         {}
00095 
00097         Element (const Ordinal i, const Ordinal j, const Scalar& Aij) :
00098           rowIndex_ (i), colIndex_ (j), value_ (Aij) {}
00099 
00101         bool operator== (const Element& rhs) {
00102           return rowIndex_ == rhs.rowIndex_ && colIndex_ == rhs.colIndex_;
00103         }
00104 
00106         bool operator!= (const Element& rhs) {
00107           return ! (*this == rhs);
00108         }
00109 
00111         bool operator< (const Element& rhs) const {
00112           if (rowIndex_ < rhs.rowIndex_)
00113             return true;
00114           else if (rowIndex_ > rhs.rowIndex_)
00115             return false;
00116           else { // equal
00117             return colIndex_ < rhs.colIndex_;
00118           }
00119         }
00120 
00126         template<class BinaryFunction>
00127         void merge (const Element& rhs, const BinaryFunction& f) {
00128           TEUCHOS_TEST_FOR_EXCEPTION(
00129             rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
00130             std::invalid_argument,
00131             "Attempt to merge elements at different locations in the sparse "
00132             "matrix.  The current element is at (" << rowIndex() << ", "
00133             << colIndex() << ") and the element you asked me to merge with it "
00134             "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << ").  This "
00135             "probably indicates a bug in the sparse matrix reader.");
00136 
00137           value_ = f (rhs.value_, value_);
00138         }
00139 
00146         void merge (const Element& rhs, const bool replace=false) {
00147           TEUCHOS_TEST_FOR_EXCEPTION(
00148             rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex(),
00149             std::invalid_argument,
00150             "Attempt to merge elements at different locations in the sparse "
00151             "matrix.  The current element is at (" << rowIndex() << ", "
00152             << colIndex() << ") and the element you asked me to merge with it "
00153             "is at (" << rhs.rowIndex() << ", " << rhs.colIndex() << ").  This "
00154             "probably indicates a bug in the sparse matrix reader.");
00155 
00156           if (replace) {
00157             value_ = rhs.value_;
00158           }
00159           else {
00160             value_ += rhs.value_;
00161           }
00162         }
00163 
00165         Ordinal rowIndex() const { return rowIndex_; }
00166 
00168         Ordinal colIndex() const { return colIndex_; }
00169 
00171         Scalar value() const { return value_; }
00172 
00173       private:
00174         Ordinal rowIndex_, colIndex_;
00175         Scalar value_;
00176       };
00177 
00188       template<class Scalar, class Ordinal>
00189       std::ostream&
00190       operator<< (std::ostream& out, const Element<Scalar, Ordinal>& elt)
00191       {
00192         typedef ScalarTraits<Scalar> STS;
00193         // Non-Ordinal types are floating-point types.  In order not to
00194         // lose information when we print a floating-point type, we have
00195         // to set the number of digits to print.  C++ standard behavior
00196         // in the default locale seems to be to print only five decimal
00197         // digits after the decimal point; this does not suffice for
00198         // double precision.  We solve the problem of how many digits to
00199         // print more generally below.  It's a rough solution so please
00200         // feel free to audit and revise it.
00201         //
00202         // FIXME (mfh 01 Feb 2011)
00203         // This really calls for the following approach:
00204         //
00205         // Guy L. Steele and Jon L. White, "How to print floating-point
00206         // numbers accurately", 20 Years of the ACM/SIGPLAN Conference
00207         // on Programming Language Design and Implementation
00208         // (1979-1999): A Selection, 2003.
00209         if (! STS::isOrdinal) {
00210           // std::scientific, std::fixed, and default are the three
00211           // output states for floating-point numbers.  A reasonable
00212           // user-defined floating-point type should respect these
00213           // flags; hopefully it does.
00214           out << std::scientific;
00215 
00216           // Decimal output is standard for Matrix Market format.
00217           out << std::setbase (10);
00218 
00219           // Compute the number of decimal digits required for expressing
00220           // a Scalar, by comparing with IEEE 754 double precision (16
00221           // decimal digits, 53 binary digits).  This would be easier if
00222           // Teuchos exposed std::numeric_limits<T>::digits10, alas.
00223           const double numDigitsAsDouble =
00224             16 * ((double) STS::t() / (double) ScalarTraits<double>::t());
00225           // Adding 0.5 and truncating is a portable "floor".
00226           const int numDigits = static_cast<int> (numDigitsAsDouble + 0.5);
00227 
00228           // Precision to which a Scalar should be written.  Add one
00229           // for good measure, since 17 is necessary for IEEE 754
00230           // doubles.
00231           out << std::setprecision (numDigits + 1);
00232         }
00233         out << elt.rowIndex () << " " << elt.colIndex () << " ";
00234         if (STS::isComplex) {
00235           out << STS::real (elt.value ()) << " " << STS::imag (elt.value ());
00236         }
00237         else {
00238           out << elt.value ();
00239         }
00240         return out;
00241       }
00242 
00279       template<class Scalar, class Ordinal>
00280       class Adder {
00281       public:
00282         typedef Ordinal index_type;
00283         typedef Scalar value_type;
00284         typedef Element<Scalar, Ordinal> element_type;
00285         typedef typename std::vector<element_type>::size_type size_type;
00286 
00297         Adder () :
00298           expectedNumRows_(0),
00299           expectedNumCols_(0),
00300           expectedNumEntries_(0),
00301           seenNumRows_(0),
00302           seenNumCols_(0),
00303           seenNumEntries_(0),
00304           tolerant_ (true),
00305           debug_ (false)
00306         {}
00307 
00325         Adder (const Ordinal expectedNumRows,
00326                const Ordinal expectedNumCols,
00327                const Ordinal expectedNumEntries,
00328                const bool tolerant=false,
00329                const bool debug=false) :
00330           expectedNumRows_(expectedNumRows),
00331           expectedNumCols_(expectedNumCols),
00332           expectedNumEntries_(expectedNumEntries),
00333           seenNumRows_(0),
00334           seenNumCols_(0),
00335           seenNumEntries_(0),
00336           tolerant_ (tolerant),
00337           debug_ (debug)
00338         {}
00339 
00360         void
00361         operator() (const Ordinal i,
00362                     const Ordinal j,
00363                     const Scalar& Aij,
00364                     const bool countAgainstTotal=true)
00365         {
00366           if (! tolerant_) {
00367             const bool indexPairOutOfRange = i < 1 || j < 1 ||
00368               i > expectedNumRows_ || j > expectedNumCols_;
00369 
00370             TEUCHOS_TEST_FOR_EXCEPTION(indexPairOutOfRange,
00371               std::invalid_argument, "Matrix is " << expectedNumRows_ << " x "
00372               << expectedNumCols_ << ", so entry A(" << i << "," << j << ") = "
00373               << Aij << " is out of range.");
00374             if (countAgainstTotal) {
00375               TEUCHOS_TEST_FOR_EXCEPTION(seenNumEntries_ >= expectedNumEntries_,
00376                 std::invalid_argument, "Cannot add entry A(" << i << "," << j
00377                 << ") = " << Aij << " to matrix; already have expected number "
00378                 "of entries " << expectedNumEntries_ << ".");
00379             }
00380           }
00381           // i and j are 1-based indices, but we store them as 0-based.
00382           elts_.push_back (element_type (i-1, j-1, Aij));
00383 
00384           // Keep track of the rightmost column containing a matrix
00385           // entry, and the bottommost row containing a matrix entry.
00386           // This gives us a lower bound for the dimensions of the
00387           // matrix, and a check for the reported dimensions of the
00388           // matrix in the Matrix Market file.
00389           seenNumRows_ = std::max (seenNumRows_, i);
00390           seenNumCols_ = std::max (seenNumCols_, j);
00391           if (countAgainstTotal) {
00392             ++seenNumEntries_;
00393           }
00394         }
00395 
00405         void
00406         print (std::ostream& out, const bool doMerge, const bool replace=false)
00407         {
00408           if (doMerge) {
00409             merge (replace);
00410           } else {
00411             std::sort (elts_.begin(), elts_.end());
00412           }
00413           // Print out the results, delimited by newlines.
00414           typedef std::ostream_iterator<element_type> iter_type;
00415           std::copy (elts_.begin(), elts_.end(), iter_type (out, "\n"));
00416         }
00417 
00440         std::pair<size_type, size_type>
00441         merge (const bool replace=false)
00442         {
00443           typedef typename std::vector<element_type>::iterator iter_type;
00444 
00445           // Start with a sorted container.  Element objects sort in
00446           // lexicographic order of their (row, column) indices, for
00447           // easy conversion to CSR format.  If you expect that the
00448           // elements will usually be sorted in the desired order, you
00449           // can check first whether they are already sorted.  We have
00450           // no such expectation, so we don't even bother to spend the
00451           // extra O(# entries) operations to check.
00452           std::sort (elts_.begin(), elts_.end());
00453 
00454           // Walk through the array of elements in place, merging
00455           // duplicates and pushing unique elements up to the front of
00456           // the array.  We can't use std::unique for this because it
00457           // doesn't let us merge duplicate elements; it only removes
00458           // them from the sequence.
00459           size_type numUnique = 0;
00460           iter_type cur = elts_.begin();
00461           if (cur == elts_.end()) { // No elements to merge
00462             return std::make_pair (numUnique, size_type (0));
00463           }
00464           else {
00465             iter_type next = cur;
00466             ++next; // There is one unique element
00467             ++numUnique;
00468             while (next != elts_.end()) {
00469               if (*cur == *next) {
00470                 // Merge in the duplicated element *next
00471                 cur->merge (*next, replace);
00472               } else {
00473                 // *cur is already a unique element.  Move over one to
00474                 // *make space for the new unique element.
00475                 ++cur;
00476                 *cur = *next; // Add the new unique element
00477                 ++numUnique;
00478               }
00479               // Look at the "next" not-yet-considered element
00480               ++next;
00481             }
00482             // Remember how many elements we removed before resizing.
00483             const size_type numRemoved = elts_.size() - numUnique;
00484             elts_.resize (numUnique);
00485             return std::make_pair (numUnique, numRemoved);
00486           }
00487         }
00488 
00531         void
00532         mergeAndConvertToCSR (size_type& numUniqueElts,
00533                               size_type& numRemovedElts,
00534                               Teuchos::ArrayRCP<Ordinal>& rowptr,
00535                               Teuchos::ArrayRCP<Ordinal>& colind,
00536                               Teuchos::ArrayRCP<Scalar>& values,
00537                               const bool replace=false)
00538         {
00539           using Teuchos::arcp;
00540           using Teuchos::ArrayRCP;
00541 
00542           std::pair<size_type, size_type> mergeResult = merge (replace);
00543 
00544           // At this point, elts_ is already in CSR order.
00545           // Now we can allocate and fill the ind and val arrays.
00546           ArrayRCP<Ordinal> ind = arcp<Ordinal> (elts_.size ());
00547           ArrayRCP<Scalar> val = arcp<Scalar> (elts_.size ());
00548 
00549           // Number of rows in the matrix.
00550           const Ordinal nrows = tolerant_ ? seenNumRows_ : expectedNumRows_;
00551           ArrayRCP<Ordinal> ptr = arcp<Ordinal> (nrows + 1);
00552 
00553           // Copy over the elements, and fill in the ptr array with
00554           // offsets.  Note that merge() sorted the entries by row
00555           // index, so we can assume the row indices are increasing in
00556           // the list of entries.
00557           Ordinal curRow = 0;
00558           Ordinal curInd = 0;
00559           typedef typename std::vector<element_type>::const_iterator iter_type;
00560           ptr[0] = 0; // ptr always has at least one entry.
00561           for (iter_type it = elts_.begin(); it != elts_.end(); ++it) {
00562             const Ordinal i = it->rowIndex ();
00563             const Ordinal j = it->colIndex ();
00564             const Scalar Aij = it->value ();
00565 
00566             TEUCHOS_TEST_FOR_EXCEPTION(i < curRow, std::logic_error, "The "
00567               "current matrix entry's row index " << i << " is less then what "
00568               "should be the current row index lower bound " << curRow << ".");
00569             for (Ordinal k = curRow+1; k <= i; ++k) {
00570               ptr[k] = curInd;
00571             }
00572             curRow = i;
00573 
00574             TEUCHOS_TEST_FOR_EXCEPTION(
00575               static_cast<size_t> (curInd) >= elts_.size (),
00576               std::logic_error, "The current index " << curInd << " into ind "
00577               "and val is >= the number of matrix entries " << elts_.size ()
00578               << ".");
00579             ind[curInd] = j;
00580             val[curInd] = Aij;
00581             ++curInd;
00582           }
00583           for (Ordinal k = curRow+1; k <= nrows; ++k) {
00584             ptr[k] = curInd;
00585           }
00586 
00587           // Assign to outputs here, to ensure the strong exception
00588           // guarantee (assuming that ArrayRCP's operator= doesn't
00589           // throw).
00590           rowptr = ptr;
00591           colind = ind;
00592           values = val;
00593           numUniqueElts = mergeResult.first;
00594           numRemovedElts = mergeResult.second;
00595         }
00596 
00598         const std::vector<element_type>& getEntries() const {
00599           return elts_;
00600         }
00601 
00603         void clear() {
00604           seenNumRows_ = 0;
00605           seenNumCols_ = 0;
00606           seenNumEntries_ = 0;
00607           elts_.resize (0);
00608         }
00609 
00613         const Ordinal numRows() const { return seenNumRows_; }
00614 
00618         const Ordinal numCols() const { return seenNumCols_; }
00619 
00620       private:
00621         Ordinal expectedNumRows_, expectedNumCols_, expectedNumEntries_;
00622         Ordinal seenNumRows_, seenNumCols_, seenNumEntries_;
00623         bool tolerant_;
00624         bool debug_;
00625 
00627         std::vector<element_type> elts_;
00628       };
00629     } // namespace Raw
00630   } // namespace MatrixMarket
00631 } // namespace Teuchos
00632 
00633 #endif // #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines