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           if (rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex()) {
00131             throw std::invalid_argument ("Can only merge elements at the same "
00132                                          "location in the sparse matrix");
00133           }
00134           else {
00135             value_ = f (rhs.value_, value_);
00136           }
00137         }
00138 
00144         void merge (const Element& rhs, const bool replace=false) {
00145           if (rowIndex() != rhs.rowIndex() || colIndex() != rhs.colIndex())
00146             throw std::logic_error("Can only merge elements at the same "
00147                                    "location in the sparse matrix");
00148           else if (replace)
00149             value_ = rhs.value_;
00150           else
00151             value_ += rhs.value_;
00152         }
00153 
00155         Ordinal rowIndex() const { return rowIndex_; }
00156 
00158         Ordinal colIndex() const { return colIndex_; }
00159 
00161         Scalar value() const { return value_; }
00162 
00163       private:
00164         Ordinal rowIndex_, colIndex_;
00165         Scalar value_;
00166       };
00167 
00178       template<class Scalar, class Ordinal>
00179       std::ostream&
00180       operator<< (std::ostream& out, const Element<Scalar, Ordinal>& elt)
00181       {
00182         typedef ScalarTraits<Scalar> STS;
00183         // Non-Ordinal types are floating-point types.  In order not to
00184         // lose information when we print a floating-point type, we have
00185         // to set the number of digits to print.  C++ standard behavior
00186         // in the default locale seems to be to print only five decimal
00187         // digits after the decimal point; this does not suffice for
00188         // double precision.  We solve the problem of how many digits to
00189         // print more generally below.  It's a rough solution so please
00190         // feel free to audit and revise it.
00191         //
00192         // FIXME (mfh 01 Feb 2011)
00193         // This really calls for the following approach:
00194         //
00195         // Guy L. Steele and Jon L. White, "How to print floating-point
00196         // numbers accurately", 20 Years of the ACM/SIGPLAN Conference
00197         // on Programming Language Design and Implementation
00198         // (1979-1999): A Selection, 2003.
00199         if (! STS::isOrdinal) {
00200           // std::scientific, std::fixed, and default are the three
00201           // output states for floating-point numbers.  A reasonable
00202           // user-defined floating-point type should respect these
00203           // flags; hopefully it does.
00204           out << std::scientific;
00205 
00206           // Decimal output is standard for Matrix Market format.
00207           out << std::setbase (10);
00208 
00209           // Compute the number of decimal digits required for expressing
00210           // a Scalar, by comparing with IEEE 754 double precision (16
00211           // decimal digits, 53 binary digits).  This would be easier if
00212           // Teuchos exposed std::numeric_limits<T>::digits10, alas.
00213           const double numDigitsAsDouble =
00214             16 * ((double) STS::t() / (double) ScalarTraits<double>::t());
00215           // Adding 0.5 and truncating is a portable "floor".
00216           const int numDigits = static_cast<int> (numDigitsAsDouble + 0.5);
00217 
00218           // Precision to which a Scalar should be written.  Add one
00219           // for good measure, since 17 is necessary for IEEE 754
00220           // doubles.
00221           out << std::setprecision (numDigits + 1);
00222         }
00223         out << elt.rowIndex () << " " << elt.colIndex () << " ";
00224         if (STS::isComplex) {
00225           out << STS::real (elt.value ()) << " " << STS::imag (elt.value ());
00226         }
00227         else {
00228           out << elt.value ();
00229         }
00230         return out;
00231       }
00232 
00245       template<class Scalar, class Ordinal>
00246       class Adder {
00247       public:
00248         typedef Ordinal index_type;
00249         typedef Scalar value_type;
00250         typedef Element<Scalar, Ordinal> element_type;
00251         typedef typename std::vector<element_type>::size_type size_type;
00252 
00263         Adder () :
00264           expectedNumRows_(0),
00265           expectedNumCols_(0),
00266           expectedNumEntries_(0),
00267           seenNumRows_(0),
00268           seenNumCols_(0),
00269           seenNumEntries_(0),
00270           tolerant_ (true),
00271           debug_ (false)
00272         {}
00273 
00291         Adder (const Ordinal expectedNumRows,
00292                const Ordinal expectedNumCols,
00293                const Ordinal expectedNumEntries,
00294                const bool tolerant=false,
00295                const bool debug=false) :
00296           expectedNumRows_(expectedNumRows),
00297           expectedNumCols_(expectedNumCols),
00298           expectedNumEntries_(expectedNumEntries),
00299           seenNumRows_(0),
00300           seenNumCols_(0),
00301           seenNumEntries_(0),
00302           tolerant_ (tolerant),
00303           debug_ (debug)
00304         {}
00305 
00318         void
00319         operator() (const Ordinal i, const Ordinal j, const Scalar& Aij)
00320         {
00321           if (! tolerant_) {
00322             const bool indexPairOutOfRange = i < 1 || j < 1 ||
00323               i > expectedNumRows_ || j > expectedNumCols_;
00324 
00325             TEUCHOS_TEST_FOR_EXCEPTION(indexPairOutOfRange,
00326               std::invalid_argument, "Matrix is " << expectedNumRows_ << " x "
00327               << expectedNumCols_ << ", so entry A(" << i << "," << j << ") = "
00328               << Aij << " is out of range.");
00329             TEUCHOS_TEST_FOR_EXCEPTION(seenNumEntries_ >= expectedNumEntries_,
00330               std::invalid_argument, "Cannot add entry A(" << i << "," << j
00331               << ") = " << Aij << " to matrix; already have expected number "
00332               "of entries " << expectedNumEntries_ << ".");
00333           }
00334           // i and j are 1-based indices, but we store them as 0-based.
00335           elts_.push_back (element_type (i-1, j-1, Aij));
00336 
00337           // Keep track of the rightmost column containing a matrix
00338           // entry, and the bottommost row containing a matrix entry.
00339           // This gives us a lower bound for the dimensions of the
00340           // matrix, and a check for the reported dimensions of the
00341           // matrix in the Matrix Market file.
00342           seenNumRows_ = std::max (seenNumRows_, i);
00343           seenNumCols_ = std::max (seenNumCols_, j);
00344           seenNumEntries_++;
00345         }
00346 
00356         void
00357         print (std::ostream& out, const bool doMerge, const bool replace=false)
00358         {
00359           if (doMerge) {
00360             merge (replace);
00361           } else {
00362             std::sort (elts_.begin(), elts_.end());
00363           }
00364           // Print out the results, delimited by newlines.
00365           typedef std::ostream_iterator<element_type> iter_type;
00366           std::copy (elts_.begin(), elts_.end(), iter_type (out, "\n"));
00367         }
00368 
00391         std::pair<size_type, size_type>
00392         merge (const bool replace=false)
00393         {
00394           typedef typename std::vector<element_type>::iterator iter_type;
00395 
00396           // Start with a sorted container.  Element objects sort in
00397           // lexicographic order of their (row, column) indices, for
00398           // easy conversion to CSR format.  If you expect that the
00399           // elements will usually be sorted in the desired order, you
00400           // can check first whether they are already sorted.  We have
00401           // no such expectation, so we don't even bother to spend the
00402           // extra O(# entries) operations to check.
00403           std::sort (elts_.begin(), elts_.end());
00404 
00405           // Walk through the array of elements in place, merging
00406           // duplicates and pushing unique elements up to the front of
00407           // the array.  We can't use std::unique for this because it
00408           // doesn't let us merge duplicate elements; it only removes
00409           // them from the sequence.
00410           size_type numUnique = 0;
00411           iter_type cur = elts_.begin();
00412           if (cur == elts_.end()) { // No elements to merge
00413             return std::make_pair (numUnique, size_type (0));
00414           }
00415           else {
00416             iter_type next = cur;
00417             ++next; // There is one unique element
00418             ++numUnique;
00419             while (next != elts_.end()) {
00420               if (*cur == *next) {
00421                 // Merge in the duplicated element *next
00422                 cur->merge (*next, replace);
00423               } else {
00424                 // *cur is already a unique element.  Move over one to
00425                 // *make space for the new unique element.
00426                 ++cur;
00427                 *cur = *next; // Add the new unique element
00428                 ++numUnique;
00429               }
00430               // Look at the "next" not-yet-considered element
00431               ++next;
00432             }
00433             // Remember how many elements we removed before resizing.
00434             const size_type numRemoved = elts_.size() - numUnique;
00435             elts_.resize (numUnique);
00436             return std::make_pair (numUnique, numRemoved);
00437           }
00438         }
00439 
00482         void
00483         mergeAndConvertToCSR (size_type& numUniqueElts,
00484                               size_type& numRemovedElts,
00485                               Teuchos::ArrayRCP<Ordinal>& rowptr,
00486                               Teuchos::ArrayRCP<Ordinal>& colind,
00487                               Teuchos::ArrayRCP<Scalar>& values,
00488                               const bool replace=false)
00489         {
00490           using Teuchos::arcp;
00491           using Teuchos::ArrayRCP;
00492 
00493           std::pair<size_type, size_type> mergeResult = merge (replace);
00494 
00495           // At this point, elts_ is already in CSR order.
00496           // Now we can allocate and fill the ind and val arrays.
00497           ArrayRCP<Ordinal> ind = arcp<Ordinal> (elts_.size ());
00498           ArrayRCP<Scalar> val = arcp<Scalar> (elts_.size ());
00499 
00500           // Number of rows in the matrix.
00501           const Ordinal nrows = tolerant_ ? seenNumRows_ : expectedNumRows_;
00502           ArrayRCP<Ordinal> ptr = arcp<Ordinal> (nrows + 1);
00503 
00504           // Copy over the elements, and fill in the ptr array with
00505           // offsets.  Note that merge() sorted the entries by row
00506           // index, so we can assume the row indices are increasing in
00507           // the list of entries.
00508           Ordinal curRow = 0;
00509           Ordinal curInd = 0;
00510           typedef typename std::vector<element_type>::const_iterator iter_type;
00511           ptr[0] = 0; // ptr always has at least one entry.
00512           for (iter_type it = elts_.begin(); it != elts_.end(); ++it) {
00513             const Ordinal i = it->rowIndex ();
00514             const Ordinal j = it->colIndex ();
00515             const Scalar Aij = it->value ();
00516 
00517             TEUCHOS_TEST_FOR_EXCEPTION(i < curRow, std::logic_error, "The "
00518               "current matrix entry's row index " << i << " is less then what "
00519               "should be the current row index lower bound " << curRow << ".");
00520             for (Ordinal k = curRow+1; k <= i; ++k) {
00521               ptr[k] = curInd;
00522             }
00523             curRow = i;
00524 
00525             TEUCHOS_TEST_FOR_EXCEPTION(
00526               static_cast<size_t> (curInd) >= elts_.size (),
00527               std::logic_error, "The current index " << curInd << " into ind "
00528               "and val is >= the number of matrix entries " << elts_.size ()
00529               << ".");
00530             ind[curInd] = j;
00531             val[curInd] = Aij;
00532             ++curInd;
00533           }
00534           for (Ordinal k = curRow+1; k <= nrows; ++k) {
00535             ptr[k] = curInd;
00536           }
00537 
00538           // Assign to outputs here, to ensure the strong exception
00539           // guarantee (assuming that ArrayRCP's operator= doesn't
00540           // throw).
00541           rowptr = ptr;
00542           colind = ind;
00543           values = val;
00544           numUniqueElts = mergeResult.first;
00545           numRemovedElts = mergeResult.second;
00546         }
00547 
00549         const std::vector<element_type>& getEntries() const {
00550           return elts_;
00551         }
00552 
00554         void clear() {
00555           seenNumRows_ = 0;
00556           seenNumCols_ = 0;
00557           seenNumEntries_ = 0;
00558           elts_.resize (0);
00559         }
00560 
00562         const Ordinal numRows() const { return seenNumRows_; }
00563 
00565         const Ordinal numCols() const { return seenNumCols_; }
00566 
00567       private:
00568         Ordinal expectedNumRows_, expectedNumCols_, expectedNumEntries_;
00569         Ordinal seenNumRows_, seenNumCols_, seenNumEntries_;
00570         bool tolerant_;
00571         bool debug_;
00572         std::vector<element_type> elts_;
00573       };
00574     } // namespace Raw
00575   } // namespace MatrixMarket
00576 } // namespace Teuchos
00577 
00578 #endif // #ifndef __Teuchos_MatrixMarket_Raw_Adder_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines