MatrixMarket_Tpetra.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
// @HEADER

#ifndef __MatrixMarket_Tpetra_hpp
#define __MatrixMarket_Tpetra_hpp

#include "Tpetra_CrsMatrix.hpp"
#include "Teuchos_MatrixMarket_Raw_Adder.hpp"
#include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp"
#include "Teuchos_MatrixMarket_assignScalar.hpp"
#include "Teuchos_MatrixMarket_Banner.hpp"
#include "Teuchos_MatrixMarket_CoordDataReader.hpp"
#include "Teuchos_MatrixMarket_SetScientific.hpp"

#include <algorithm>
#include <fstream>
#include <iostream>
#include <iterator>
#include <vector>
#include <stdexcept>

namespace Tpetra {
  namespace MatrixMarket {

    template<class SparseMatrixType>
    class Reader {
    public:
      typedef SparseMatrixType sparse_matrix_type;
      typedef RCP<sparse_matrix_type> sparse_matrix_ptr;

      typedef typename SparseMatrixType::scalar_type scalar_type;
      typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
      typedef typename SparseMatrixType::global_ordinal_type
        global_ordinal_type;

      typedef typename SparseMatrixType::node_type node_type;

      typedef MultiVector<scalar_type,
                          local_ordinal_type,
                          global_ordinal_type,
                          node_type> multivector_type;

      typedef RCP<node_type> node_ptr;
      typedef Comm<int> comm_type;
      typedef RCP<const comm_type> comm_ptr;
      typedef Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
      typedef RCP<const map_type> map_ptr;

    private:

      typedef typename ArrayRCP<global_ordinal_type>::size_type size_type;

      static RCP<const map_type>
      makeRangeMap (const RCP<const comm_type>& pComm,
                    const RCP<node_type>& pNode,
                    const global_ordinal_type numRows)
      {
        // A conventional, uniformly partitioned, contiguous map.
        return rcp (new map_type (static_cast<global_size_t> (numRows),
                                  static_cast<global_ordinal_type> (0),
                                  pComm, GloballyDistributed, pNode));
      }

      static RCP<const map_type>
      makeRowMap (const RCP<const map_type>& pRowMap,
                  const RCP<const comm_type>& pComm,
                  const RCP<node_type>& pNode,
                  const global_ordinal_type numRows)
      {
        // If the caller didn't provide a map, return a conventional,
        // uniformly partitioned, contiguous map.
        if (pRowMap.is_null()) {
          return rcp (new map_type (static_cast<global_size_t> (numRows),
                                    static_cast<global_ordinal_type> (0),
                                    pComm, GloballyDistributed, pNode));
        } else {
          TEUCHOS_TEST_FOR_EXCEPTION(! pRowMap->isDistributed() && pComm->getSize() > 1,
                             std::invalid_argument,
                             "The specified row map is not distributed, but "
                             "the given communicator includes more than one "
                             "rank (in fact, there are " << pComm->getSize()
                             << " ranks).");
          TEUCHOS_TEST_FOR_EXCEPTION(pRowMap->getComm() != pComm,
                             std::invalid_argument,
                             "The specified row map's communicator (pRowMap->"
                             "getComm()) is different than the given separately "
                             "supplied communicator pComm.");
          TEUCHOS_TEST_FOR_EXCEPTION(pRowMap->getNode() != pNode,
                             std::invalid_argument,
                             "The specified row map's node (pRowMap->getNode()) "
                             "is different than the given separately supplied "
                             "node pNode.");
          return pRowMap;
        }
      }

      static map_ptr
      makeDomainMap (const map_ptr& pRangeMap,
                     const global_ordinal_type numRows,
                     const global_ordinal_type numCols)
      {
        // Abbreviations so that the map creation call isn't too long.
        typedef local_ordinal_type LO;
        typedef global_ordinal_type GO;
        typedef node_type Node;

        if (numRows == numCols) {
          return pRangeMap;
        } else {
          comm_ptr pComm = pRangeMap->getComm();
          node_ptr pNode = pRangeMap->getNode();
          return createUniformContigMapWithNode<LO,GO,Node> (numCols,
                                                             pComm,
                                                             pNode);
        }
      }

      static void
      distribute (ArrayRCP<size_t>& myNumEntriesPerRow,
                  ArrayRCP<size_t>& myRowPtr,
                  ArrayRCP<global_ordinal_type>& myColInd,
                  ArrayRCP<scalar_type>& myValues,
                  const RCP<const map_type>& pRowMap,
                  ArrayRCP<size_t>& numEntriesPerRow,
                  ArrayRCP<size_t>& rowPtr,
                  ArrayRCP<global_ordinal_type>& colInd,
                  ArrayRCP<scalar_type>& values,
                  const bool debug=false)
      {
         using Teuchos::CommRequest;
         using Teuchos::receive;
         using Teuchos::send;
         using std::cerr;
         using std::endl;

         const bool extraDebug = false;
         comm_ptr pComm = pRowMap->getComm();
         const int numProcs = Teuchos::size (*pComm);
         const int myRank = Teuchos::rank (*pComm);
         const int rootRank = 0;

         // Type abbreviations to make the code more concise.
         typedef global_ordinal_type GO;
         typedef local_ordinal_type LO;

         // List of the global indices of my rows.  They may or may
         // not be contiguous, and the row map need not be one-to-one.
         ArrayView<const GO> myRows = pRowMap->getNodeElementList();
         const size_type myNumRows = myRows.size();
         TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
                            pRowMap->getNodeNumElements(),
                            std::logic_error,
                            "pRowMap->getNodeElementList().size() = "
                            << myNumRows
                            << " != pRowMap->getNodeNumElements() = "
                            << pRowMap->getNodeNumElements() << ".  "
                            "Please report this bug to the Tpetra developers.");
         TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
                            std::logic_error,
                            "On Proc 0: numEntriesPerRow.size() = "
                            << numEntriesPerRow.size()
                            << " != pRowMap->getNodeElementList().size() = "
                            << myNumRows << ".  Please report this bug to the "
                            "Tpetra developers.");

         // Space for my proc's number of entries per row.
         // Will be filled in below.
         myNumEntriesPerRow = arcp<size_t> (myNumRows);

         // Teuchos::receive() returns an int; here is space for it.
         int recvResult = 0;

         if (myRank != rootRank)
           {
             // Tell the root how many rows we have.  If we're sending
             // none, then we don't have anything else to send, nor
             // does the root have to receive anything else.
             send (*pComm, myNumRows, rootRank);
             if (myNumRows != 0)
               {
                 // Now send my rows' global indices.  Hopefully the
                 // cast to int doesn't overflow.  This is unlikely,
                 // since it should fit in a LO, even
                 // though it is a GO.
                 send (*pComm, static_cast<int> (myNumRows),
                       myRows.getRawPtr(), rootRank);

                 // I (this proc) don't care if my global row indices
                 // are contiguous, though the root proc does (since
                 // otherwise it needs to pack noncontiguous data into
                 // contiguous storage before sending).  That's why we
                 // don't check for contiguousness here.

                 // Ask the root processor for my part of the array of the
                 // number of entries per row.
                 recvResult = receive (*pComm, rootRank,
                                       static_cast<int> (myNumRows),
                                       myNumEntriesPerRow.getRawPtr());

                 // Use the resulting array to figure out how many column
                 // indices and values for which I should ask from the root
                 // processor.
                 const local_ordinal_type myNumEntries =
                   std::accumulate (myNumEntriesPerRow.begin(),
                                    myNumEntriesPerRow.end(),
                                    0);

                 // Make space for my entries of the sparse matrix.
                 // Note that they don't have to be sorted by row
                 // index.  Iterating through all my rows requires
                 // computing a running sum over myNumEntriesPerRow.
                 myColInd = arcp<GO> (myNumEntries);
                 myValues = arcp<scalar_type> (myNumEntries);
                 if (myNumEntries > 0)
                   { // Ask for that many column indices and values,
                     // if there are any.
                     recvResult = receive (*pComm, rootRank,
                                           static_cast<int> (myNumEntries),
                                           myColInd.getRawPtr());
                     recvResult = receive (*pComm, rootRank,
                                           static_cast<int> (myNumEntries),
                                           myValues.getRawPtr());
                   }
               } // If I own at least one row
           } // If I am not the root processor
         else // I _am_ the root processor
           {
             if (debug)
               cerr << "-- Proc 0: Copying my data from global arrays" << endl;

             // Proc 0 (the root processor) still needs to (allocate,
             // if not done already) and fill its part of the matrix
             // (my*).
             for (size_type k = 0; k < myNumRows; ++k)
               {
                 const GO myCurRow = myRows[k];
                 const local_ordinal_type numEntriesInThisRow = numEntriesPerRow[myCurRow];
                 //myNumEntriesPerRow[k] = numEntriesPerRow[myCurRow];
                 myNumEntriesPerRow[k] = numEntriesInThisRow;
               }
             if (extraDebug && debug)
               {
                 cerr << "Proc " << Teuchos::rank (*(pRowMap->getComm()))
                      << ": myNumEntriesPerRow[0.." << (myNumRows-1) << "] = [";
                 for (size_type k = 0; k < myNumRows; ++k)
                   {
                     cerr << myNumEntriesPerRow[k];
                     if (k < myNumRows-1)
                       cerr << " ";
                   }
                 cerr << "]" << endl;
               }
             // The total number of matrix entries that my proc owns.
             const local_ordinal_type myNumEntries =
               std::accumulate (myNumEntriesPerRow.begin(),
                                myNumEntriesPerRow.end(),
                                0);
             if (debug)
               cerr << "-- Proc 0: I own " << myNumRows << " rows and "
                    << myNumEntries << " entries" << endl;

             myColInd = arcp<GO> (myNumEntries);
             myValues = arcp<scalar_type> (myNumEntries);

             // Copy Proc 0's part of the matrix into the my* arrays.
             // It's important that myCurPos be updated _before_ k,
             // otherwise myCurPos will get the wrong number of entries
             // per row (it should be for the row in the just-completed
             // iteration, not for the next iteration's row).
             local_ordinal_type myCurPos = 0;
             for (size_type k = 0; k < myNumRows;
                  myCurPos += myNumEntriesPerRow[k], ++k)
               {
                 const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
                 const GO myRow = myRows[k];
                 const size_t curPos = rowPtr[myRow];
                 if (extraDebug && debug)
                   {
                     cerr << "k = " << k << ", myRow = " << myRow << ": colInd("
                          << curPos << "," << curNumEntries << "), myColInd("
                          << myCurPos << "," << curNumEntries << ")" << endl;
                   }
                 // Only copy if there are entries to copy, in order
                 // not to construct empty ranges for the ArrayRCP
                 // views.
                 if (curNumEntries > 0)
                   {
                     ArrayView<GO> colIndView = colInd(curPos, curNumEntries);
                     ArrayView<GO> myColIndView =
                       myColInd(myCurPos, curNumEntries);
                     std::copy (colIndView.begin(), colIndView.end(),
                                myColIndView.begin());

                     ArrayView<scalar_type> valuesView =
                       values(curPos, curNumEntries);
                     ArrayView<scalar_type> myValuesView =
                       myValues(myCurPos, curNumEntries);
                     std::copy (valuesView.begin(), valuesView.end(),
                                myValuesView.begin());
                   }
               }

             // Proc 0 processes each other proc p in turn.
             for (int p = 1; p < numProcs; ++p)
               {
                 if (debug)
                   cerr << "-- Proc 0: Processing proc " << p << endl;

                 size_type theirNumRows = 0;
                 // Ask Proc p how many rows it has.  If it doesn't
                 // have any, we can move on to the next proc.  This
                 // has to be a standard receive so that we can avoid
                 // the degenerate case of sending zero data.
                 recvResult = receive (*pComm, p, &theirNumRows);
                 if (debug)
                   cerr << "-- Proc 0: Proc " << p << " owns "
                        << theirNumRows << " rows" << endl;
                 if (theirNumRows != 0)
                   { // Ask Proc p which rows it owns.  The resulting
                     // global row indices are not guaranteed to be
                     // contiguous or sorted.  Global row indices are
                     // themselves indices into the numEntriesPerRow
                     // array.
                     ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
                     recvResult = receive (*pComm, p,
                                           static_cast<int> (theirNumRows),
                                           theirRows.getRawPtr());
                     // Extra test to make sure that the rows we
                     // received are all sensible.  This is a good
                     // idea since we are going to use the global row
                     // indices we've received to index into the
                     // numEntriesPerRow array.  Better to catch any
                     // bugs here and print a sensible error message,
                     // rather than segfault and print a cryptic error
                     // message.
                     {
                       const global_size_t numRows =
                         pRowMap->getGlobalNumElements();
                       bool theirRowsValid = true;
                       for (size_type k = 0; k < theirNumRows; ++k)
                         {
                           // global_ordinal_t is generally a signed type.
                           if (theirRows[k] < 0)
                             theirRowsValid = false;
                           // Signed to unsigned cast should not
                           // overflow.  We cast in order to avoid
                           // compiler warnings about comparing signed
                           // and unsigned types, in case
                           // global_size_t is unsigned.
                           else if (static_cast<global_size_t> (theirRows[k]) >= numRows)
                             theirRowsValid = false;
                         }
                       if (! theirRowsValid)
                         {
                           std::ostringstream os;
                           std::copy (theirRows.begin(), theirRows.end(),
                                      std::ostream_iterator<GO>(os, " "));
                           TEUCHOS_TEST_FOR_EXCEPTION(! theirRowsValid,
                                              std::logic_error,
                                              "Proc " << p << " has at least "
                                              "one invalid row index.  Here are"
                                              " all of them: [ " << os.str()
                                              << "]");
                         }
                     }

                     // Perhaps we could save a little work if we
                     // check whether Proc p's row indices are
                     // contiguous.  That would make lookups in the
                     // global data arrays faster.  For now, we just
                     // implement the general case and don't
                     // prematurely optimize.  (Remember that you're
                     // making Proc 0 read the whole file, so you've
                     // already lost scalability.)

                     // Compute the number of entries in each of Proc
                     // p's rows.  (Proc p will compute its row
                     // pointer array on its own, after it gets the
                     // data from Proc 0.)
                     ArrayRCP<size_t> theirNumEntriesPerRow;
                     theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
                     for (size_type k = 0; k < theirNumRows; ++k)
                       theirNumEntriesPerRow[k] =
                         numEntriesPerRow[theirRows[k]];

                     // Tell Proc p the number of entries in each of
                     // its rows.  Hopefully the cast to int doesn't
                     // overflow.  This is unlikely, since it should
                     // fit in a LO, even though it is
                     // a GO.
                     send (*pComm, static_cast<int> (theirNumRows),
                           theirNumEntriesPerRow.getRawPtr(), p);

                     // Figure out how many entries Proc p owns.
                     const local_ordinal_type theirNumEntries =
                       std::accumulate (theirNumEntriesPerRow.begin(),
                                        theirNumEntriesPerRow.end(),
                                        0);

                     if (debug)
                       cerr << "-- Proc 0: Proc " << p << " owns "
                            << theirNumEntries << " entries" << endl;

                     // If there are no entries to send, then we're
                     // done with Proc p.
                     if (theirNumEntries == 0)
                       continue;

                     // Construct (views of) proc p's column indices
                     // and values.  Later, we might like to optimize
                     // for the (common) contiguous case, for which we
                     // don't need to copy data into separate "their*"
                     // arrays (we can just use contiguous views of
                     // the global arrays).
                     ArrayRCP<GO> theirColInd = arcp<GO> (theirNumEntries);
                     ArrayRCP<scalar_type> theirValues =
                       arcp<scalar_type> (theirNumEntries);
                     // Copy Proc p's part of the matrix into the
                     // their* arrays.  It's important that
                     // theirCurPos be updated _before_ k, otherwise
                     // theirCurPos will get the wrong number of
                     // entries per row (it should be for the row in
                     // the just-completed iteration, not for the next
                     // iteration's row).
                     local_ordinal_type theirCurPos = 0;
                     for (size_type k = 0; k < theirNumRows;
                          theirCurPos += theirNumEntriesPerRow[k], k++)
                       {
                         const local_ordinal_type curNumEntries = theirNumEntriesPerRow[k];
                         const GO theirRow = theirRows[k];
                         const local_ordinal_type curPos = rowPtr[theirRow];

                         // Only copy if there are entries to copy, in
                         // order not to construct empty ranges for
                         // the ArrayRCP views.
                         if (curNumEntries > 0)
                           {
                             ArrayView<GO> colIndView =
                               colInd(curPos, curNumEntries);
                             ArrayView<GO> theirColIndView =
                               theirColInd(theirCurPos, curNumEntries);
                             std::copy (colIndView.begin(), colIndView.end(),
                                        theirColIndView.begin());

                             ArrayView<scalar_type> valuesView =
                               values(curPos, curNumEntries);
                             ArrayView<scalar_type> theirValuesView =
                               theirValues(theirCurPos, curNumEntries);
                             std::copy (valuesView.begin(), valuesView.end(),
                                        theirValuesView.begin());
                           }
                       }
                     // Send Proc p its column indices and values.
                     // Hopefully the cast to int doesn't overflow.
                     // This is unlikely, since it should fit in a LO,
                     // even though it is a GO.
                     send (*pComm, static_cast<int> (theirNumEntries),
                           theirColInd.getRawPtr(), p);
                     send (*pComm, static_cast<int> (theirNumEntries),
                           theirValues.getRawPtr(), p);

                     if (debug)
                       cerr << "-- Proc 0: Finished with proc " << p << endl;
                   } // If proc p owns at least one row
               } // For each proc p not the root proc 0
           } // If I'm (not) the root proc 0

         // Invalidate the input data to save space, since we don't
         // need it anymore.
         numEntriesPerRow = null;
         rowPtr = null;
         colInd = null;
         values = null;

         if (debug && myRank == 0)
           cerr << "-- Proc 0: About to fill in myRowPtr" << endl;

         // Allocate and fill in myRowPtr (the row pointer array for
         // my rank's rows).  We delay this until the end because we
         // don't need it to compute anything else in distribute().
         // Each proc can do this work for itself, since it only needs
         // myNumEntriesPerRow to do so.
         myRowPtr = arcp<size_t> (myNumRows+1);
         myRowPtr[0] = 0;
         for (size_type k = 1; k < myNumRows+1; ++k) {
           myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
         }
         if (extraDebug && debug)
           {
             cerr << "Proc " << Teuchos::rank (*(pRowMap->getComm()))
                  << ": myRowPtr[0.." << myNumRows << "] = [";
             for (size_type k = 0; k < myNumRows+1; ++k)
               {
                 cerr << myRowPtr[k];
                 if (k < myNumRows)
                   cerr << " ";
               }
             cerr << "]" << endl << endl;
           }

         if (debug && myRank == 0)
           cerr << "-- Proc 0: Done with distribute" << endl;
      }

      static sparse_matrix_ptr
      makeMatrix (ArrayRCP<size_t>& myNumEntriesPerRow,
                  ArrayRCP<size_t>& myRowPtr,
                  ArrayRCP<global_ordinal_type>& myColInd,
                  ArrayRCP<scalar_type>& myValues,
                  const map_ptr& pRowMap,
                  const map_ptr& pRangeMap,
                  const map_ptr& pDomainMap,
                  const bool callFillComplete = true)
      {
        using std::cerr;
        using std::endl;
        // Typedef to make certain type declarations shorter.
        typedef global_ordinal_type GO;

        const bool extraDebug = false;
        const bool debug = false;

        // The row pointer array always has at least one entry, even
        // if the matrix has zero rows.  myNumEntriesPerRow, myColInd,
        // and myValues would all be empty arrays in that degenerate
        // case, but the row and domain maps would still be nonnull
        // (though they would be trivial maps).
        TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
                           "makeMatrix: myRowPtr array is null.  "
                           "Please report this bug to the Tpetra developers.");
        TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
                           "makeMatrix: domain map is null.  "
                           "Please report this bug to the Tpetra developers.");
        TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
                           "makeMatrix: range map is null.  "
                           "Please report this bug to the Tpetra developers.");
        TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
                           "makeMatrix: row map is null.  "
                           "Please report this bug to the Tpetra developers.");

        // Handy for debugging output; not needed otherwise.
        const int myRank = Teuchos::rank (*(pRangeMap->getComm()));

        if (extraDebug && debug)
          {
            cerr << "Proc " << myRank << ":" << endl
                 << "-- myRowPtr = [ ";
            std::copy (myRowPtr.begin(), myRowPtr.end(),
                       std::ostream_iterator<size_type>(cerr, " "));
            cerr << "]" << endl << "-- myColInd = [ ";
            std::copy (myColInd.begin(), myColInd.end(),
                       std::ostream_iterator<size_type>(cerr, " "));
            cerr << "]" << endl << endl;
          }

        // Go through all of my columns, and see if any are not in the
        // domain map.  This is possible if numProcs > 1, otherwise
        // not.
        if (extraDebug && debug)
          {
            size_type numRemote = 0;
            std::vector<GO> remoteGIDs;

            typedef typename ArrayRCP<GO>::const_iterator iter_type;
            for (iter_type it = myColInd.begin(); it != myColInd.end(); ++it)
              {
                if (! pDomainMap->isNodeGlobalElement (*it))
                  {
                    numRemote++;
                    remoteGIDs.push_back (*it);
                  }
              }
            if (numRemote > 0)
              {
                cerr << "Proc " << myRank << ": " << numRemote
                     << " remote GIDs = [ " << endl;
                std::copy (remoteGIDs.begin(), remoteGIDs.end(),
                           std::ostream_iterator<GO>(cerr, " "));
                cerr << "]" << endl;
              }
          }

        // Construct the CrsMatrix, using the row map, with the
        // constructor specifying the number of nonzeros for each row.
        // Create with DynamicProfile, so that the fillComplete() can
        // do first-touch reallocation (a NUMA (Non-Uniform Memory
        // Access) optimization on multicore CPUs).
        sparse_matrix_ptr A =
          rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
                                       DynamicProfile));
        TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), std::logic_error,
                           "makeMatrix: Initial allocation of CrsMatrix failed"
                           ".  Please report this bug to the Tpetra developers"
                           ".");
        // List of the global indices of my rows.
        // They may or may not be contiguous.
        ArrayView<const GO> myRows = pRowMap->getNodeElementList();
        const size_type myNumRows = myRows.size();

        // Add this processor's matrix entries to the CrsMatrix.
        for (size_type k = 0; k < myNumRows; ++k)
          {
            const size_type myCurPos = myRowPtr[k];
            const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];

            if (extraDebug && debug)
              {
                cerr << "Proc " << myRank << ": k = " << k
                     << ", myCurPos = " << myCurPos
                     << ", curNumEntries = " << curNumEntries
                     << endl;
              }
            // Avoid constructing empty views of ArrayRCP objects.
            if (curNumEntries > 0)
              A->insertGlobalValues (myRows[k],
                                     myColInd(myCurPos, curNumEntries),
                                     myValues(myCurPos, curNumEntries));
          }
        // We've entered in all our matrix entries, so we can delete
        // the original data.  This will save memory when we call
        // fillComplete(), so that we never keep more than two copies
        // of the matrix's data in memory at once.
        myNumEntriesPerRow = null;
        myRowPtr = null;
        myColInd = null;
        myValues = null;

        if (callFillComplete)
          A->fillComplete (pDomainMap, pRangeMap);
        return A;
      }

    private:

      static RCP<const Teuchos::MatrixMarket::Banner>
      readBanner (std::istream& in,
                  size_t& lineNumber,
                  const bool tolerant=false,
                  const bool debug=false)
      {
        using Teuchos::MatrixMarket::Banner;
        using std::cerr;
        using std::endl;
        typedef Teuchos::ScalarTraits<scalar_type> STS;

        RCP<Banner> pBanner; // On output, if successful: the read-in Banner.
        std::string line; // If read from stream successful: the Banner line

        // Try to read a line from the input stream.
        const bool readFailed = ! getline(in, line);
        TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
          "Failed to get Matrix Market banner line from input.");

        // We read a line from the input stream.
        lineNumber++;

        // Assume that the line we found is the Banner line.
        try {
          pBanner = rcp (new Banner (line, tolerant));
        } catch (std::exception& e) {
          TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
            "Matrix Market banner line contains syntax error(s): "
            << e.what());
        }
        TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() != "matrix",
          std::invalid_argument, "The Matrix Market file does not contain "
          "matrix data.  Its Banner (first) line says that its object type is \""
          << pBanner->matrixType() << "\", rather than the required \"matrix\".");

        // Validate the data type of the matrix, with respect to the
        // Scalar type of the CrsMatrix entries.
        TEUCHOS_TEST_FOR_EXCEPTION(
          ! STS::isComplex && pBanner->dataType() == "complex",
          std::invalid_argument,
          "The Matrix Market file contains complex-valued data, but you are "
          "trying to read it into a matrix containing entries of the real-"
          "valued Scalar type \""
          << Teuchos::TypeNameTraits<scalar_type>::name() << "\".");
        TEUCHOS_TEST_FOR_EXCEPTION(
          pBanner->dataType() != "real" &&
          pBanner->dataType() != "complex" &&
          pBanner->dataType() != "integer",
          std::invalid_argument,
          "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
          "Matrix Market file may not contain a \"pattern\" matrix.  A "
          "pattern matrix is really just a graph with no weights.  It "
          "should be stored in a CrsGraph, not a CrsMatrix.");

        return pBanner;
      }

      static Tuple<global_ordinal_type, 3>
      readCoordDims (std::istream& in,
                     size_t& lineNumber,
                     const RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
                     const comm_ptr& pComm,
                     const bool tolerant = false,
                     const bool debug = false)
      {
        using Teuchos::MatrixMarket::readCoordinateDimensions;

        // Packed coordinate matrix dimensions (numRows, numCols,
        // numNonzeros); computed on Rank 0 and broadcasted to all MPI
        // ranks.
        Tuple<global_ordinal_type, 3> dims;

        // Read in the coordinate matrix dimensions from the input
        // stream.  "success" tells us whether reading in the
        // coordinate matrix dimensions succeeded ("Guilty unless
        // proven innocent").
        bool success = false;
        if (pComm->getRank() == 0) {
          TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "coordinate",
            std::invalid_argument, "The Tpetra::CrsMatrix Matrix Market reader "
            "only accepts \"coordinate\" (sparse) matrix data.");
          // Unpacked coordinate matrix dimensions
          global_ordinal_type numRows, numCols, numNonzeros;
          // Only MPI Rank 0 reads from the input stream
          success = readCoordinateDimensions (in, numRows, numCols,
                                              numNonzeros, lineNumber,
                                              tolerant);
          // Pack up the data into a Tuple so we can send them with
          // one broadcast instead of three.
          dims[0] = numRows;
          dims[1] = numCols;
          dims[2] = numNonzeros;
        }
        // Only Rank 0 did the reading, so it decides success.
        //
        // FIXME (mfh 02 Feb 2011) Teuchos::broadcast doesn't know how
        // to send bools.  For now, we convert to/from int instead,
        // using the usual "true is 1, false is 0" encoding.
        {
          int the_success = success ? 1 : 0; // only matters on MPI Rank 0
          Teuchos::broadcast (*pComm, 0, &the_success);
          success = (the_success == 1);
        }
        if (success) {
          // Broadcast (numRows, numCols, numNonzeros) from Rank 0
          // to all the other MPI ranks.
          Teuchos::broadcast (*pComm, 0, dims);
        }
        else {
          // Perhaps in tolerant mode, we could set all the
          // dimensions to zero for now, and deduce correct
          // dimensions by reading all of the file's entries and
          // computing the max(row index) and max(column index).
          // However, for now we just error out in that case.
          TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
            "Error reading Matrix Market sparse matrix: failed to read "
            "coordinate matrix dimensions.");
        }
        return dims;
      }

      typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;

      static RCP<adder_type>
      makeAdder (const RCP<const Comm<int> >& pComm,
                 RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
                 const Tuple<global_ordinal_type, 3>& dims,
                 const bool tolerant=false,
                 const bool debug=false)
      {
        if (pComm->getRank() == 0) {
          typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> raw_adder_type;
          RCP<raw_adder_type> pRaw =
            rcp (new raw_adder_type (dims[0], dims[1], dims[2], tolerant, debug));
          return rcp (new adder_type (pRaw, pBanner->symmType()));
        }
        else {
          return null;
        }
      }

    public:

      static sparse_matrix_ptr
      readSparseFile (const std::string& filename,
                      const RCP<const Comm<int> >& pComm,
                      const RCP<node_type>& pNode,
                      const bool callFillComplete=true,
                      const bool tolerant=false,
                      const bool debug=false)
      {
        const int myRank = Teuchos::rank (*pComm);
        std::ifstream in;

        // Only open the file on Rank 0.
        if (myRank == 0)
          in.open (filename.c_str());
        return readSparse (in, pComm, pNode, callFillComplete, tolerant, debug);
        // We can rely on the destructor of the input stream to close
        // the file on scope exit, even if readSparse() throws an
        // exception.
      }

      static sparse_matrix_ptr
      readSparse (std::istream& in,
                  const RCP<const Comm<int> >& pComm,
                  const RCP<node_type>& pNode,
                  const bool callFillComplete=true,
                  const bool tolerant=false,
                  const bool debug=false)
      {
        using Teuchos::MatrixMarket::Banner;
        using Teuchos::broadcast;
        using Teuchos::ptr;
        using Teuchos::reduceAll;
        using std::cerr;
        using std::endl;
        typedef Teuchos::ScalarTraits<scalar_type> STS;

        const bool extraDebug = false;
        const int myRank = pComm->getRank ();
        const int rootRank = 0;

        // Current line number in the input stream.  Various calls
        // will modify this depending on the number of lines that are
        // read from the input stream.  Only Rank 0 modifies this.
        size_t lineNumber = 1;

        if (debug && myRank == rootRank) {
          cerr << "Matrix Market reader: readSparse:" << endl
               << "-- Reading banner line" << endl;
        }

        // The "Banner" tells you whether the input stream represents
        // a sparse matrix, the symmetry type of the matrix, and the
        // type of the data it contains.
        //
        // pBanner will only be nonnull on MPI Rank 0.  It will be
        // null on all other MPI processes.
        RCP<const Banner> pBanner;
        {
          // We read and validate the Banner on Proc 0, but broadcast
          // the validation result to all processes.
          // Teuchos::broadcast doesn't currently work with bool, so
          // we use int (true -> 1, false -> 0).
          int bannerIsCorrect = 1;
          std::ostringstream errMsg;

          if (myRank == rootRank) {
            // Read the Banner line from the input stream.
            try {
              pBanner = readBanner (in, lineNumber, tolerant, debug);
            }
            catch (std::exception& e) {
              errMsg << "Attempt to read the Matrix Market file's Banner line "
                "threw an exception: " << e.what();
              bannerIsCorrect = 0;
            }

            if (bannerIsCorrect) {
              // Validate the Banner for the case of a sparse matrix.
              // We validate on Proc 0, since it reads the Banner.

              // In intolerant mode, the matrix type must be "coordinate".
              if (! tolerant && pBanner->matrixType() != "coordinate") {
                bannerIsCorrect = 0;
                errMsg << "The Matrix Market input file must contain a "
                  "\"coordinate\"-format sparse matrix in order to create a "
                  "Tpetra::CrsMatrix object from it, but the file's matrix "
                  "type is \"" << pBanner->matrixType() << "\" instead.";
              }
              // In tolerant mode, we allow the matrix type to be
              // anything other than "array" (which would mean that
              // the file contains a dense matrix).
              if (tolerant && pBanner->matrixType() == "array") {
                bannerIsCorrect = 0;
                errMsg << "Matrix Market file must contain a \"coordinate\"-"
                  "format sparse matrix in order to create a Tpetra::CrsMatrix "
                  "object from it, but the file's matrix type is \"array\" "
                  "instead.  That probably means the file contains dense matrix "
                  "data.";
              }
            }
          } // Proc 0: Done reading the Banner, hopefully successfully.

          // Broadcast from Proc 0 whether the Banner was read correctly.
          broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));

          // If the Banner is invalid, all processes throw an
          // exception.  Only Proc 0 gets the exception message, but
          // that's OK, since the main point is to "stop the world"
          // (rather than throw an exception on one process and leave
          // the others hanging).
          TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
            std::invalid_argument, errMsg.str ());
        } // Done reading the Banner line and broadcasting success.
        if (debug && myRank == rootRank) {
          cerr << "-- Reading dimensions line" << endl;
        }

        // Read the matrix dimensions from the Matrix Market metadata.
        // dims = (numRows, numCols, numEntries).  Proc 0 does the
        // reading, but it broadcasts the results to all MPI
        // processes.  Thus, readCoordDims() is a collective
        // operation.  It does a collective check for correctness too.
        Tuple<global_ordinal_type, 3> dims =
          readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);

        if (debug && myRank == rootRank) {
          cerr << "-- Making Adder for collecting matrix data" << endl;
        }

        // "Adder" object for collecting all the sparse matrix entries
        // from the input stream.  This is only nonnull on Proc 0.
        RCP<adder_type> pAdder =
          makeAdder (pComm, pBanner, dims, tolerant, debug);

        if (debug && myRank == rootRank) {
          cerr << "-- Reading matrix data" << endl;
        }
        //
        // Read the matrix entries from the input stream on Proc 0.
        //
        {
          // We use readSuccess to broadcast the results of the read
          // (succeeded or not) to all MPI processes.  Since
          // Teuchos::broadcast doesn't currently know how to send
          // bools, we convert to int (true -> 1, false -> 0).
          int readSuccess = 1;
          std::ostringstream errMsg; // Exception message (only valid on Proc 0)
          if (myRank == rootRank) {
            try {
              // Reader for "coordinate" format sparse matrix data.
              typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
                global_ordinal_type, scalar_type, STS::isComplex> reader_type;
              reader_type reader (pAdder);

              // Read the sparse matrix entries.
              std::pair<bool, std::vector<size_t> > results =
                reader.read (in, lineNumber, tolerant, debug);
              readSuccess = results.first ? 1 : 0;
            }
            catch (std::exception& e) {
              readSuccess = 0;
              errMsg << e.what();
            }
          }
          broadcast (*pComm, rootRank, ptr (&readSuccess));

          // It would be nice to add a "verbose" flag, so that in
          // tolerant mode, we could log any bad line number(s) on
          // Proc 0.  For now, we just throw if the read fails to
          // succeed.
          //
          // Question: If we're in tolerant mode, and if the read did
          // not succeed, should we attempt to call fillComplete()?
          TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
            "Failed to read the Matrix Market sparse matrix file: "
            << errMsg.str());
        } // Done reading the matrix entries (stored on Proc 0 for now)

        if (debug && myRank == rootRank) {
          cerr << "-- Successfully read the Matrix Market data" << endl;
        }

        // In tolerant mode, we need to rebroadcast the matrix
        // dimensions, since they may be different after reading the
        // actual matrix data.  We only need to broadcast the number
        // of rows and columns.  Only Rank 0 needs to know the actual
        // global number of entries, since (a) we need to merge
        // duplicates on Rank 0 first anyway, and (b) when we
        // distribute the entries, each rank other than Rank 0 will
        // only need to know how many entries it owns, not the total
        // number of entries.
        if (tolerant) {
          if (debug && myRank == rootRank) {
            cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
                 << endl
                 << "----- Dimensions before: "
                 << dims[0] << " x " << dims[1]
                 << endl;
          }
          // Packed coordinate matrix dimensions (numRows, numCols).
          Tuple<global_ordinal_type, 2> updatedDims;
          if (myRank == rootRank) {
            // If one or more bottom rows of the matrix contain no
            // entries, then the Adder will report that the number
            // of rows is less than that specified in the
            // metadata.  We allow this case, and favor the
            // metadata so that the zero row(s) will be included.
            updatedDims[0] =
              std::max (dims[0], pAdder->getAdder()->numRows());
            updatedDims[1] = pAdder->getAdder()->numCols();
          }
          broadcast (*pComm, rootRank, updatedDims);
          dims[0] = updatedDims[0];
          dims[1] = updatedDims[1];
          if (debug && myRank == rootRank) {
            cerr << "----- Dimensions after: " << dims[0] << " x "
                 << dims[1] << endl;
          }
        }
        else {
          // In strict mode, we require that the matrix's metadata and
          // its actual data agree, at least somewhat.  In particular,
          // the number of rows must agree, since otherwise we cannot
          // distribute the matrix correctly.

          // Teuchos::broadcast() doesn't know how to broadcast bools,
          // so we use an int with the standard 1 == true, 0 == false
          // encoding.
          int dimsMatch = 1;
          if (myRank == rootRank) {
            // If one or more bottom rows of the matrix contain no
            // entries, then the Adder will report that the number of
            // rows is less than that specified in the metadata.  We
            // allow this case, and favor the metadata, but do not
            // allow the Adder to think there are more rows in the
            // matrix than the metadata says.
            if (dims[0] < pAdder->getAdder ()->numRows ()) {
              dimsMatch = 0;
            }
          }
          broadcast (*pComm, 0, ptr (&dimsMatch));
          if (dimsMatch == 0) {
            // We're in an error state anyway, so we might as well
            // work a little harder to print an informative error
            // message.
            //
            // Broadcast the Adder's idea of the matrix dimensions
            // from Proc 0 to all processes.
            Tuple<global_ordinal_type, 2> addersDims;
            if (myRank == rootRank) {
              addersDims[0] = pAdder->getAdder()->numRows();
              addersDims[1] = pAdder->getAdder()->numCols();
            }
            broadcast (*pComm, 0, addersDims);
            TEUCHOS_TEST_FOR_EXCEPTION(
              dimsMatch == 0, std::runtime_error,
              "The matrix metadata says that the matrix is " << dims[0] << " x "
              << dims[1] << ", but the actual data says that the matrix is "
              << addersDims[0] << " x " << addersDims[1] << ".  That means the "
              "data includes more rows than reported in the metadata.  This "
              "is not allowed when parsing in strict mode.  Parse the matrix in "
              "tolerant mode to ignore the metadata when it disagrees with the "
              "data.");
          }
        } // Matrix dimensions (# rows, # cols, # entries) agree.

        if (debug && myRank == rootRank) {
          cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
        }

        // Now that we've read in all the matrix entries from the
        // input stream into the adder on Proc 0, post-process them
        // into CSR format (still on Proc 0).  This will facilitate
        // distributing them to all the processors.
        //
        // These arrays represent the global matrix data as a CSR
        // matrix (with numEntriesPerRow as redundant but convenient
        // metadata, since it's computable from rowPtr and vice
        // versa).  They are valid only on Proc 0.
        ArrayRCP<size_t> numEntriesPerRow;
        ArrayRCP<size_t> rowPtr;
        ArrayRCP<global_ordinal_type> colInd;
        ArrayRCP<scalar_type> values;

        // Proc 0 first merges duplicate entries, and then converts
        // the coordinate-format matrix data to CSR.
        {
          int mergeAndConvertSucceeded = 1;
          std::ostringstream errMsg;

          if (myRank == rootRank) {
            try {
              typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
                global_ordinal_type> element_type;
              typedef typename std::vector<element_type>::const_iterator iter_type;

              // Number of rows in the matrix.  If we are in tolerant
              // mode, we've already synchronized dims with the actual
              // matrix data.  If in strict mode, we should use dims
              // (as read from the file's metadata) rather than the
              // matrix data to determine the dimensions.  (The matrix
              // data will claim fewer rows than the metadata, if one
              // or more rows have no entries stored in the file.)
              const size_type numRows = dims[0];

              // Additively merge duplicate matrix entries.
              pAdder->getAdder()->merge ();

              // Get a temporary const view of the merged matrix entries.
              const std::vector<element_type>& entries =
                pAdder->getAdder()->getEntries();

              // Number of matrix entries (after merging).
              const size_t numEntries = (size_t)entries.size();

              if (debug) {
                cerr << "----- Proc 0: Matrix has numRows=" << numRows
                     << " rows and numEntries=" << numEntries
                     << " entries." << endl;
              }

              // Make space for the CSR matrix data.  Converting to
              // CSR is easier if we fill numEntriesPerRow with zeros
              // at first.
              numEntriesPerRow = arcp<size_t> (numRows);
              std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
              rowPtr = arcp<size_t> (numRows+1);
              std::fill (rowPtr.begin(), rowPtr.end(), 0);
              colInd = arcp<global_ordinal_type> (numEntries);
              values = arcp<scalar_type> (numEntries);

              // Convert from array-of-structs coordinate format to CSR
              // (compressed sparse row) format.
              global_ordinal_type prvRow = 0;
              size_t curPos = 0;
              rowPtr[0] = 0;
              for (curPos = 0; curPos < numEntries; ++curPos) {
                const element_type& curEntry = entries[curPos];
                const global_ordinal_type curRow = curEntry.rowIndex();
                TEUCHOS_TEST_FOR_EXCEPTION(
                  curRow < prvRow, std::logic_error,
                  "Row indices are out of order, even though they are supposed "
                  "to be sorted.  curRow = " << curRow << ", prvRow = "
                  << prvRow << ", at curPos = " << curPos << ".  Please report "
                  "this bug to the Tpetra developers.");
                if (curRow > prvRow) {
                  for (size_type r = prvRow+1; r <= curRow; ++r) {
                    rowPtr[r] = curPos;
                  }
                  prvRow = curRow;
                }
                numEntriesPerRow[curRow]++;
                colInd[curPos] = curEntry.colIndex();
                values[curPos] = curEntry.value();
              }
              // rowPtr has one more entry than numEntriesPerRow.  The
              // last entry of rowPtr is the number of entries in
              // colInd and values.
              rowPtr[numRows] = numEntries;
            } // Finished conversion to CSR format
            catch (std::exception& e) {
              mergeAndConvertSucceeded = 0;
              errMsg << "Failed to merge sparse matrix entries and convert to "
                "CSR format: " << e.what();
            }

            if (debug && mergeAndConvertSucceeded) {
              // Number of rows in the matrix.
              const size_type numRows = dims[0];
              const size_type maxToDisplay = 100;

              cerr << "----- Proc 0: numEntriesPerRow[0.."
                   << (numEntriesPerRow.size()-1) << "] ";
              if (numRows > maxToDisplay) {
                cerr << "(only showing first and last few entries) ";
              }
              cerr << "= [";
              if (numRows > 0) {
                if (numRows > maxToDisplay) {
                  for (size_type k = 0; k < 2; ++k) {
                    cerr << numEntriesPerRow[k] << " ";
                  }
                  cerr << "... ";
                  for (size_type k = numRows-2; k < numRows-1; ++k) {
                    cerr << numEntriesPerRow[k] << " ";
                  }
                }
                else {
                  for (size_type k = 0; k < numRows-1; ++k) {
                    cerr << numEntriesPerRow[k] << " ";
                  }
                }
                cerr << numEntriesPerRow[numRows-1];
              } // numRows > 0
              cerr << "]" << endl;

              cerr << "----- Proc 0: rowPtr ";
              if (numRows > maxToDisplay) {
                cerr << "(only showing first and last few entries) ";
              }
              cerr << "= [";
              if (numRows > maxToDisplay) {
                for (size_type k = 0; k < 2; ++k) {
                  cerr << rowPtr[k] << " ";
                }
                cerr << "... ";
                for (size_type k = numRows-2; k < numRows; ++k) {
                  cerr << rowPtr[k] << " ";
                }
              }
              else {
                for (size_type k = 0; k < numRows; ++k) {
                  cerr << rowPtr[k] << " ";
                }
              }
              cerr << rowPtr[numRows] << "]" << endl;
            }
          } // if myRank == rootRank
        } // Done converting sparse matrix data to CSR format

        // Now we're done with the Adder, so we can release the
        // reference ("free" it) to save space.  This only actually
        // does anything on Rank 0, since pAdder is null on all the
        // other MPI processes.
        pAdder = null;

        if (debug && myRank == rootRank) {
          cerr << "-- Making range, domain, and row maps" << endl;
        }

        // Make the maps that describe the matrix's range and domain,
        // and the distribution of its rows.  Creating a Map is a
        // collective operation, so we don't have to do a broadcast of
        // a success Boolean.
        map_ptr pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
        map_ptr pDomainMap = makeDomainMap (pRangeMap, dims[0], dims[1]);
        map_ptr pRowMap = makeRowMap (null, pComm, pNode, dims[0]);

        if (debug && myRank == rootRank) {
          cerr << "-- Distributing the matrix data" << endl;
        }

        // Distribute the matrix data.  Each processor has to add the
        // rows that it owns.  If you try to make Proc 0 call
        // insertGlobalValues() for _all_ the rows, not just those it
        // owns, then fillComplete() will compute the number of
        // columns incorrectly.  That's why Proc 0 has to distribute
        // the matrix data and why we make all the processors (not
        // just Proc 0) call insertGlobalValues() on their own data.
        //
        // These arrays represent each processor's part of the matrix
        // data, in "CSR" format (sort of, since the row indices might
        // not be contiguous).
        ArrayRCP<size_t> myNumEntriesPerRow;
        ArrayRCP<size_t> myRowPtr;
        ArrayRCP<global_ordinal_type> myColInd;
        ArrayRCP<scalar_type> myValues;
        // Distribute the matrix data.  This is a collective operation.
        distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
                    numEntriesPerRow, rowPtr, colInd, values, debug);

        if (debug && myRank == rootRank) {
          cerr << "-- Inserting matrix entries on each processor";
          if (callFillComplete) {
            cerr << " and calling fillComplete()";
          }
          cerr << endl;
        }
        // Each processor inserts its part of the matrix data, and
        // then they all call fillComplete().  This method invalidates
        // the my* distributed matrix data before calling
        // fillComplete(), in order to save space.  In general, we
        // never store more than two copies of the matrix's entries in
        // memory at once, which is no worse than what Tpetra
        // promises.
        sparse_matrix_ptr pMatrix =
          makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
                      pRowMap, pRangeMap, pDomainMap, callFillComplete);
        // Only use a reduce-all in debug mode to check if pMatrix is
        // null.  Otherwise, just throw an exception.  We never expect
        // a null pointer here, so we can save a communication.
        if (debug) {
          int localIsNull = pMatrix.is_null () ? 1 : 0;
          int globalIsNull = 0;
          reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
          TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
            "Reader::makeMatrix() returned a null pointer on at least one "
            "process.  Please report this bug to the Tpetra developers.");
        }
        else {
          TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
            "Reader::makeMatrix() returned a null pointer.  "
            "Please report this bug to the Tpetra developers.");
        }

        // We can't get the dimensions of the matrix until after
        // fillComplete() is called.  Thus, we can't do the sanity
        // check (dimensions read from the Matrix Market data,
        // vs. dimensions reported by the CrsMatrix) unless the user
        // asked makeMatrix() to call fillComplete().
        //
        // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
        // what one might think it does, so you have to ask the range
        // resp. domain map for the number of rows resp. columns.
        if (callFillComplete) {
          const int numProcs = pComm->getSize ();

          if (extraDebug && debug) {
            const global_size_t globalNumRows =
              pRangeMap->getGlobalNumElements();
            const global_size_t globalNumCols =
              pDomainMap->getGlobalNumElements();
            if (myRank == rootRank) {
              cerr << "-- Matrix is "
                   << globalNumRows << " x " << globalNumCols
                   << " with " << pMatrix->getGlobalNumEntries()
                   << " entries, and index base "
                   << pMatrix->getIndexBase() << "." << endl;
            }
            pComm->barrier ();
            for (int p = 0; p < numProcs; ++p) {
              if (myRank == p) {
                cerr << "-- Proc " << p << " owns "
                     << pMatrix->getNodeNumCols() << " columns, and "
                     << pMatrix->getNodeNumEntries() << " entries." << endl;
              }
              pComm->barrier ();
            }
          } // if (extraDebug && debug)
        } // if (callFillComplete)

        if (debug && myRank == rootRank) {
          cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
               << endl;
        }
        return pMatrix;
      }

      static RCP<multivector_type>
      readDenseFile (const std::string& filename,
                     const RCP<const comm_type>& pComm,
                     const RCP<node_type>& pNode,
                     RCP<const map_type>& pMap,
                     const bool tolerant=false,
                     const bool debug=false)
      {
        const int myRank = Teuchos::rank (*pComm);
        std::ifstream in;

        // Only open the file on Rank 0.
        if (myRank == 0)
          in.open (filename.c_str());
        return readDense (in, pComm, pNode, pMap, tolerant, debug);
        // We can rely on the destructor of the input stream to close
        // the file on scope exit, even if readSparse() throws an
        // exception.
      }

      static RCP<multivector_type>
      readDense (std::istream& in,
                 const RCP<const comm_type>& pComm,
                 const RCP<node_type>& pNode,
                 RCP<const map_type>& pMap,
                 const bool tolerant=false,
                 const bool debug=false)
      {
        using Teuchos::MatrixMarket::Banner;
        using Teuchos::MatrixMarket::checkCommentLine;
        using std::cerr;
        using std::endl;

        // Abbreviations for typedefs, to make the code more concise.
        typedef scalar_type S;
        typedef local_ordinal_type LO;
        typedef global_ordinal_type GO;
        typedef node_type Node;
        typedef Teuchos::ScalarTraits<S> STS;
        typedef typename STS::magnitudeType M;
        typedef Teuchos::ScalarTraits<M> STM;

        // Rank 0 is the only (MPI) process allowed to read from the
        // input stream.
        const int myRank = pComm->getRank ();

        if (debug && myRank == 0) {
          cerr << "Matrix Market reader: readDense:" << endl;
        }

        // If pMap is nonnull, check the precondition that its
        // communicator resp. node equal pComm resp. pNode.  Checking
        // now avoids doing a lot of file reading before we detect the
        // violated precondition.
        TEUCHOS_TEST_FOR_EXCEPTION(! pMap.is_null() &&
                                   (pMap->getComm() != pComm ||
                                    pMap->getNode() != pNode),
          std::invalid_argument, "If you supply a nonnull Map, the Map's "
          "communicator and node must equal the supplied communicator resp. "
          "node.");

        // Rank 0 will read in the matrix dimensions from the file,
        // and broadcast them to all ranks in the given communicator.
        // There are only 2 dimensions in the matrix, but we use the
        // third element of the Tuple to encode the banner's reported
        // data type: "real" == 0, "complex" == 1, and "integer" == 0
        // (same as "real").  We don't allow pattern matrices (i.e.,
        // graphs) since they only make sense for sparse data.
        Tuple<GO, 3> dims;
        dims[0] = 0;
        dims[1] = 0;

        // Current line number in the input stream.  Only valid on
        // Rank 0.  Various calls will modify this depending on the
        // number of lines that are read from the input stream.
        size_t lineNumber = 1;

        // Only Rank 0 gets to read matrix data from the input stream.
        if (myRank == 0) {
          if (debug && myRank == 0) {
            cerr << "-- Reading banner line (dense)" << endl;
          }

          // The "Banner" tells you whether the input stream
          // represents a dense matrix, the symmetry type of the
          // matrix, and the type of the data it contains.
          RCP<const Banner> pBanner =
            readBanner (in, lineNumber, tolerant, debug);
          TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "array",
            std::invalid_argument, "The Matrix Market file does not contain "
            "dense matrix data.  Its banner (first) line says that its matrix "
            "type is \"" << pBanner->matrixType() << "\", rather than the "
            "required \"array\".");
          TEUCHOS_TEST_FOR_EXCEPTION(pBanner->dataType() == "pattern",
            std::invalid_argument, "The Matrix Market file's banner (first) "
            "line claims that the matrix's data type is \"pattern\".  This does "
            "not make sense for a dense matrix, yet the file reports the matrix "
            "as dense.  The only valid data types for a dense matrix are "
            "\"real\", \"complex\", and \"integer\".");
          // Encode the data type reported by the Banner as the third
          // element of the dimensions Tuple.
          dims[2] = encodeDataType (pBanner->dataType());

          if (debug && myRank == 0) {
            cerr << "-- Reading dimensions line (dense)" << endl;
          }
          // Keep reading lines from the input stream until we find a
          // non-comment line, or until we run out of lines.  The
          // latter is an error, since every "array" format Matrix
          // Market file must have a dimensions line after the banner
          // (even if the matrix has zero rows or columns, or zero
          // entries).
          std::string line;
          bool commentLine = true;
          while (commentLine) {
            // Test whether it is even valid to read from the
            // input stream wrapping the line.
            TEUCHOS_TEST_FOR_EXCEPTION(in.eof() || in.fail(), std::runtime_error,
              "Unable to get array dimensions line (at all) from line "
              << lineNumber << " of input stream.  The input stream claims that "
              "it is " << (in.eof() ? "at end-of-file." : "in a failed state."));
            // Try to get the next line from the input stream.
            if (getline(in, line)) {
              ++lineNumber; // We did actually read a line.
            }
            // Is the current line a comment line?  Ignore start and
            // size; they are only useful for reading the actual
            // matrix entries.  (We could use them here as an
            // optimization, but we've chosen not to.)
            size_t start = 0, size = 0;
            commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
          }
          //
          // Read in <numRows> <numCols> from input line.
          //
          std::istringstream istr (line);

          // Test whether it is even valid to read from the input
          // stream wrapping the line.
          TEUCHOS_TEST_FOR_EXCEPTION(istr.eof() || istr.fail(), std::runtime_error,
                                     "Unable to read any data from line " << lineNumber
                                     << " of input; the line should contain the matrix "
                                     << "dimensions \"<numRows> <numCols>\".");
          // Read in <numRows>.
          {
            GO theNumRows = 0;
            istr >> theNumRows;
            TEUCHOS_TEST_FOR_EXCEPTION(istr.fail(), std::runtime_error,
                                       "Failed to get number of rows from line "
                                       << lineNumber << " of input; the line should "
                                       "contain the matrix dimensions \"<numRows> "
                                       "<numCols>\".");
            // Capture the validly read result before checking for eof.
            dims[0] = theNumRows;
          }
          // There should be one more thing to read.
          TEUCHOS_TEST_FOR_EXCEPTION(istr.eof(), std::runtime_error,
                                     "No more data after number of rows on line "
                                     << lineNumber << " of input; the line should "
                                     "contain the matrix dimensions \"<numRows> "
                                     "<numCols>\".");
          // Read in <numCols>
          {
            GO theNumCols = 0;
            istr >> theNumCols;
            TEUCHOS_TEST_FOR_EXCEPTION(istr.fail(), std::runtime_error,
                                       "Failed to get number of columns from line "
                                       << lineNumber << " of input; the line should "
                                       "contain the matrix dimensions \"<numRows> "
                                       "<numCols>\".");
            // Capture the validly read result.
            dims[1] = theNumCols;
          }
        } // if (myRank == 0)

        // Broadcast matrix dimensions and the encoded data type from
        // MPI Rank 0 to all the MPI processes.
        Teuchos::broadcast (*pComm, 0, dims);

        // Tpetra objects want the matrix dimensions in these types.
        const global_size_t numRows = static_cast<global_size_t> (dims[0]);
        const size_t numCols = static_cast<size_t> (dims[1]);

        // Make an "everything on Rank 0" map pRank0Map.  We'll use
        // this to construct a "distributed" vector X owned entirely
        // by Rank 0.  Rank 0 will then read all the matrix entries
        // and put them in X.
        RCP<const map_type> pRank0Map =
          createContigMapWithNode<LO, GO, Node> (numRows,
                                                 (myRank == 0 ? numRows : 0),
                                                 pComm, pNode);

        // Make a multivector X owned entirely by Rank 0.
        RCP<multivector_type> X =
          createMultiVector<S, LO, GO, Node> (pRank0Map, numCols);

        //
        // On Rank 0, read the Matrix Market data from the input
        // stream into the multivector X.
        //
        if (myRank == 0) {
          if (debug && myRank == 0) {
            cerr << "-- Reading matrix data (dense)" << endl;
          }

          // Make sure that we can get a 1-D view of X.
          TEUCHOS_TEST_FOR_EXCEPTION(! X->isConstantStride (), std::logic_error,
                                     "Can't get a 1-D view of the entries of the "
                                     "MultiVector X on Rank 0, because the stride "
                                     "between the columns of X is not constant.  "
                                     "This shouldn't happen because we just created "
                                     "X and haven't filled it in yet.  Please report "
                                     "this bug to the Tpetra developers.");

          // Get a writeable 1-D view of the entries of X.  Rank 0
          // owns all of them.  The view will expire at the end of
          // scope, so (if necessary) it will be written back to X
          // at this time.
          ArrayRCP<S> X_view = X->get1dViewNonConst ();
          TEUCHOS_TEST_FOR_EXCEPTION(
            static_cast<global_size_t> (X_view.size()) < numRows * numCols,
            std::logic_error,
            "The view of X has size " << X_view << " which is not enough to "
            "accommodate the expected number of entries numRows*numCols = "
            << numRows << "*" << numCols << " = " << numRows*numCols << ".  "
            "Please report this bug to the Tpetra developers.");
          const size_t stride = X->getStride ();

          // The third element of the dimensions Tuple encodes the data
          // type reported by the Banner: "real" == 0, "complex" == 1,
          // "integer" == 0 (same as "real"), "pattern" == 2.  We do not
          // allow dense matrices to be pattern matrices, so dims[2] ==
          // 0 or 1.  We've already checked for this above.
          const bool isComplex = (dims[2] == 1);
          typename ArrayRCP<S>::size_type count = 0, curRow = 0, curCol = 0;

          std::string line;
          while (getline (in, line)) {
            ++lineNumber;
            // Is the current line a comment line?  If it's not,
            // line.substr(start,size) contains the data.
            size_t start = 0, size = 0;
            const bool commentLine =
              checkCommentLine (line, start, size, lineNumber, tolerant);
            if (! commentLine) {
              // Make sure we have room in which to put the new matrix
              // entry.  We check this only after checking for a
              // comment line, because there may be one or more
              // comment lines at the end of the file.  In tolerant
              // mode, we simply ignore any extra data.
              if (count >= X_view.size()) {
                if (tolerant) {
                  break;
                }
                else {
                  TEUCHOS_TEST_FOR_EXCEPTION(
                     count >= X_view.size(),
                     std::runtime_error,
                     "The Matrix Market input stream has more data in it than "
                     "its metadata reported.  Current line number is "
                     << lineNumber << ".");
                }
              }
              std::istringstream istr (line.substr (start, size));
              // Does the line contain anything at all?  Can we
              // safely read from the input stream wrapping the
              // line?
              if (istr.eof() || istr.fail()) {
                // In tolerant mode, simply ignore the line.
                if (tolerant) {
                  break;
                }
                // We repeat the full test here so the exception
                // message is more informative.
                TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && (istr.eof() || istr.fail()),
                  std::runtime_error, "Line " << lineNumber << " of the Matrix "
                  "Market file is empty, or we cannot read from it for some "
                  "other reason.");
              }
              // Current matrix entry to read in.
              S val = STS::zero();
              // Real and imaginary parts of the current matrix entry.
              // The imaginary part is zero if the matrix is
              // real-valued.
              M real = STM::zero(), imag = STM::zero();

              // isComplex refers to the input stream's data, not to
              // the scalar type S.  It's OK to read real-valued data
              // into a matrix storing complex-valued data; in that
              // case, all entries' imaginary parts are zero.
              if (isComplex) {
                // STS::real() and STS::imag() return a copy of their
                // respective components, not a writeable reference.
                // Otherwise we could just assign to them using the
                // istream extraction operator (>>).  That's why we
                // have separate magnitude type "real" and "imag"
                // variables.

                // Attempt to read the real part of the current matrix
                // entry.
                istr >> real;
                if (istr.fail()) {
                  TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
                                             std::runtime_error,
                                             "Failed to get real part of a "
                                             "complex-valued matrix entry "
                                             "from line " << lineNumber
                                             << " of Matrix Market file.");
                  // In tolerant mode, just skip bad lines.
                  if (tolerant) {
                    break;
                  }
                }
                else if (istr.eof()) {
                  TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.eof(),
                                             std::runtime_error,
                                             "Missing imaginary part of a "
                                             "complex-valued matrix entry "
                                             "on line " << lineNumber
                                             << " of Matrix Market file.");
                  // In tolerant mode, let any missing imaginary part
                  // be 0.
                }
                else {
                  // Attempt to read the imaginary part of the current
                  // matrix entry.
                  istr >> imag;
                  TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
                                             std::runtime_error,
                                             "Failed to get imaginary part "
                                             "of a complex-valued matrix "
                                             "entry from line " << lineNumber
                                             << " of Matrix Market file.");
                  // In tolerant mode, let any missing
                  // or corrupted imaginary part be 0.
                }
              }
              else { // Matrix Market file contains real-valued data.
                // Attempt to read the current matrix entry.
                istr >> real;
                TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
                                           std::runtime_error,
                                           "Failed to get a real-valued matrix "
                                           "entry from line " << lineNumber
                                           << " of Matrix Market file.");
                // In tolerant mode, simply ignore the line if
                // we failed to read a matrix entry.
                if (istr.fail() && tolerant) {
                  break;
                }
              }

              // In tolerant mode, we simply let pass through whatever
              // data we got.
              TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
                                         std::runtime_error,
                                         "Failed to read matrix data from line "
                                         << lineNumber << " of the Matrix Market "
                                         "file.");

              // val = S(real, imag), without potential badness
              // if S is a real type.
              Teuchos::MatrixMarket::details::assignScalar<S> (val, real, imag);

              curRow = count % numRows;
              curCol = count / numRows;
              X_view[curRow + curCol*stride] = val;
              ++count;
            } // if not a comment line
          } // while there are still lines in the file, get the next one

          TEUCHOS_TEST_FOR_EXCEPTION(! tolerant &&
                                     static_cast<global_size_t> (count) < numRows * numCols,
                                     std::runtime_error,
                                     "The Matrix Market metadata reports that the "
                                     "dense matrix is " << numRows <<  " x "
                                     << numCols << ", and thus has "
                                     << numRows*numCols << " total entries, but we "
                                     "only managed to find " << count << " entr"
                                     << (count == 1 ? "y" : "ies")
                                     << " in the file.");
        } // if (myRank == 0)

        // If there's only one MPI process and the user didn't supply
        // a Map (i.e., pMap is null), we're done.  Set pMap to the
        // Map used to distribute X, and return X.
        if (pComm->getSize() == 1 && pMap.is_null()) {
          pMap = pRank0Map;
          if (debug && myRank == 0) {
            cerr << "-- Done reading multivector" << endl;
          }
          return X;
        }

        if (debug && myRank == 0) {
          cerr << "-- Creating distributed Map and target MultiVector" << endl;
        }

        // If pMap is null, make a distributed map.  We've already
        // checked the preconditions above, in the case that pMap was
        // not null on input to this method.
        if (pMap.is_null()) {
          pMap = createUniformContigMapWithNode<LO, GO, Node> (numRows, pComm, pNode);
        }
        // Make a multivector Y with the distributed map pMap.
        RCP<multivector_type> Y = createMultiVector<S, LO, GO, Node> (pMap, numCols);

        // Make an Export object that will export X to Y.  First
        // argument is the source map, second argument is the target
        // map.
        Export<LO, GO, Node> exporter (pRank0Map, pMap);

        if (debug && myRank == 0) {
          cerr << "-- Exporting from X (owned by Rank 0) to globally owned Y"
               << endl;
        }
        // Export X into Y.
        Y->doExport (*X, exporter, INSERT);

        if (debug && myRank == 0) {
          cerr << "-- Done reading multivector" << endl;
        }

        // Y is distributed over all process(es) in the communicator.
        return Y;
      }

    private:

      static int
      encodeDataType (const std::string& dataType)
      {
        if (dataType == "real" || dataType == "integer") {
          return 0;
        }
        else if (dataType == "complex") {
          return 1;
        }
        else if (dataType == "pattern") {
          return 2;
        }
        else {
          // We should never get here, since Banner validates the
          // reported data type and ensures it is one of the accepted
          // values.
          TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
            "Unrecognized Matrix Market data type \"" << dataType << "\".");
        }
      }
    };

    template<class SparseMatrixType>
    class Writer {
    public:
      typedef SparseMatrixType sparse_matrix_type;
      typedef RCP<sparse_matrix_type> sparse_matrix_ptr;

      typedef typename SparseMatrixType::scalar_type scalar_type;
      typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
      typedef typename SparseMatrixType::global_ordinal_type global_ordinal_type;
      typedef typename SparseMatrixType::node_type node_type;

      typedef MultiVector<scalar_type,
                          local_ordinal_type,
                          global_ordinal_type,
                          node_type> multivector_type;
      typedef Map<local_ordinal_type, global_ordinal_type, node_type> map_type;

      static void
      writeSparseFile (const std::string& filename,
                       const RCP<const sparse_matrix_type>& pMatrix,
                       const std::string& matrixName,
                       const std::string& matrixDescription,
                       const bool debug=false)
      {
        const int myRank = Teuchos::rank (*(pMatrix->getComm()));
        std::ofstream out;

        // Only open the file on Rank 0.
        if (myRank == 0) out.open (filename.c_str());
        writeSparse (out, pMatrix, matrixName, matrixDescription, debug);
        // We can rely on the destructor of the output stream to close
        // the file on scope exit, even if writeSparse() throws an
        // exception.
      }

      static void
      writeSparseFile (const std::string& filename,
                       const RCP<const sparse_matrix_type>& pMatrix,
                       const bool debug=false)
      {
        writeSparseFile (filename, pMatrix, "", "", debug);
      }

    public:

      static void
      writeSparse (std::ostream& out,
                   const RCP<const sparse_matrix_type>& pMatrix,
                   const std::string& matrixName,
                   const std::string& matrixDescription,
                   const bool debug=false)
      {
        using std::cerr;
        using std::endl;
        typedef typename Teuchos::ScalarTraits<scalar_type> STS;
        typedef typename STS::magnitudeType magnitude_type;
        typedef typename Teuchos::ScalarTraits<magnitude_type> STM;
        typedef typename ArrayView<scalar_type>::size_type size_type;

        // Convenient abbreviations.
        typedef local_ordinal_type LO;
        typedef global_ordinal_type GO;

        // Make the output stream write floating-point numbers in
        // scientific notation.  It will politely put the output
        // stream back to its state on input, when this scope
        // terminates.
        Teuchos::MatrixMarket::details::SetScientific<scalar_type> sci (out);

        RCP<const Comm<int> > pComm = pMatrix->getComm();
        const int myRank = Teuchos::rank (*pComm);
        // Whether to print debugging output to stderr.
        const bool debugPrint = debug && myRank == 0;

        const global_size_t numRows =
          pMatrix->getRangeMap()->getGlobalNumElements();
        const global_size_t numCols =
          pMatrix->getDomainMap()->getGlobalNumElements();
        if (debugPrint) {
          cerr << "writeSparse:" << endl
               << "-- Input sparse matrix is:"
               << "---- " << numRows << " x " << numCols << " with "
               << pMatrix->getGlobalNumEntries() << " entries;" << endl
               << "---- "
               << (pMatrix->isGloballyIndexed() ? "Globally" : "Locally")
               << " indexed." << endl;
        }
        RCP<const map_type> pRowMap = pMatrix->getRowMap();
        if (debugPrint) {
          RCP<const map_type> pColMap = pMatrix->getColMap();
          cerr << "---- Its row map has " << pRowMap->getGlobalNumElements()
               << " elements." << endl
               << "---- Its col map has " << pColMap->getGlobalNumElements()
               << " elements." << endl;
        }
        // Make the "gather" row map, where Proc 0 owns all rows and
        // the other procs own no rows.
        const size_t localNumRows = (myRank == 0) ? numRows : 0;
        RCP<node_type> pNode = pRowMap->getNode();
        RCP<const map_type> pGatherRowMap =
          createContigMapWithNode<LO, GO, node_type> (numRows, localNumRows,
                                                      pComm, pNode);
        // Since the matrix may in general be non-square, we need to
        // make a column map as well.  In this case, the column map
        // contains all the columns of the original matrix, because we
        // are gathering the whole matrix onto Proc 0.
        const size_t localNumCols = (myRank == 0) ? numCols : 0;
        RCP<const map_type> pGatherColMap =
          createContigMapWithNode<LO, GO, node_type> (numCols, localNumCols,
                                                      pComm, pNode);

        // Current map is the source map, gather map is the target map.
        typedef Import<LO, GO, node_type> import_type;
        import_type importer (pRowMap, pGatherRowMap);

        // Create a new CrsMatrix to hold the result of the import.
        // The constructor needs a column map as well as a row map,
        // for the case that the matrix is not square.
        RCP<sparse_matrix_type> newMatrix =
          rcp (new sparse_matrix_type (pGatherRowMap,
                                       pGatherColMap,
                                       size_t(0)));
        // Import the sparse matrix onto Proc 0.
        newMatrix->doImport (*pMatrix, importer, INSERT);

        // fillComplete() needs the domain and range maps for the case
        // that the matrix is not square.
        {
          RCP<const map_type> pGatherDomainMap =
            createContigMapWithNode<LO, GO, node_type> (numCols, localNumCols,
                                                        pComm, pNode);
          RCP<const map_type> pGatherRangeMap =
            createContigMapWithNode<LO, GO, node_type> (numRows, localNumRows,
                                                        pComm, pNode);
          newMatrix->fillComplete (pGatherDomainMap, pGatherRangeMap);
        }

        if (debugPrint) {
          cerr << "-- Output sparse matrix is:"
               << "---- " << newMatrix->getRangeMap()->getGlobalNumElements()
               << " x "
               << newMatrix->getDomainMap()->getGlobalNumElements()
               << " with "
               << newMatrix->getGlobalNumEntries() << " entries;" << endl
               << "---- "
               << (newMatrix->isGloballyIndexed() ? "Globally" : "Locally")
               << " indexed." << endl
               << "---- Its row map has "
               << newMatrix->getRowMap()->getGlobalNumElements()
               << " elements." << endl
               << "---- Its col map has "
               << newMatrix->getColMap()->getGlobalNumElements()
               << " elements (may be different than the input matrix's column"
               << " map, if some columns of the matrix contain no entries)."
               << endl;
        }

        //
        // Print the metadata and the matrix entries on Rank 0.
        //
        if (myRank == 0) {
          // Print the Matrix Market banner line.  CrsMatrix stores
          // data nonsymmetrically ("general").  This implies that
          // readSparse() on a symmetrically stored input file,
          // followed by writeSparse() on the resulting sparse matrix,
          // will result in an output file with a different banner
          // line than the original input file.
          out << "%%MatrixMarket matrix coordinate "
              << (STS::isComplex ? "complex" : "real")
              << " general" << endl;

          // Print comments (the matrix name and / or description).
          if (matrixName != "")
            printAsComment (out, matrixName);
          if (matrixDescription != "")
            printAsComment (out, matrixDescription);

          // Print the Matrix Market header (# rows, # columns, #
          // nonzeros).  Use the range resp. domain map for the
          // number of rows resp. columns, since Tpetra::CrsMatrix
          // uses the column map for the number of columns.  That
          // only corresponds to the "linear-algebraic" number of
          // columns when the column map is 1-1, which only happens
          // if the matrix is (block) diagonal.
          out << newMatrix->getRangeMap()->getGlobalNumElements()
              << " "
              << newMatrix->getDomainMap()->getGlobalNumElements()
              << " "
              << newMatrix->getGlobalNumEntries()
              << endl;

          // Index base (0-based or 1-based) for the row map.
          const GO rowIndexBase = pGatherRowMap->getIndexBase();
          // Index base (0-based or 1-based) for the column map.
          const GO colIndexBase = newMatrix->getColMap()->getIndexBase();

          //
          // Print the entries of the matrix.
          //
          // newMatrix can never be globally indexed, since we
          // called fillComplete() on it.  We include code for both
          // cases (globally or locally indexed) just in case that
          // ever changes.
          if (newMatrix->isGloballyIndexed()) {
            for (GO globalRowIndex = pGatherRowMap->getMinAllGlobalIndex();
                 globalRowIndex <= pGatherRowMap->getMaxAllGlobalIndex();
                 ++globalRowIndex)
              {
                ArrayView<const GO> ind;
                ArrayView<const scalar_type> val;

                newMatrix->getGlobalRowView (globalRowIndex, ind, val);
                typename ArrayView<const GO>::const_iterator indIter =
                  ind.begin();
                typename ArrayView<const scalar_type>::const_iterator valIter =
                  val.begin();
                for (; indIter != ind.end() && valIter != val.end();
                     ++indIter, ++valIter)
                  {
                    const GO globalColIndex = *indIter;

                    // Matrix Market files use 1-based indices.
                    out << (globalRowIndex + 1 - rowIndexBase) << " "
                        << (globalColIndex + 1 - colIndexBase) << " ";
                    if (STS::isComplex)
                      out << STS::real(*valIter) << " " << STS::imag(*valIter);
                    else
                      out << *valIter;
                    out << endl;
                  }
              }
          }
          else // newMatrix is locally indexed
            {
              typedef OrdinalTraits<GO> OTG;
              RCP<const map_type> pColMap = newMatrix->getColMap ();

              for (LO localRowIndex = pGatherRowMap->getMinLocalIndex();
                   localRowIndex <= pGatherRowMap->getMaxLocalIndex();
                   ++localRowIndex)
                {
                  // Convert from local to global row index.
                  const GO globalRowIndex =
                    pGatherRowMap->getGlobalElement (localRowIndex);
                  TEUCHOS_TEST_FOR_EXCEPTION(globalRowIndex == OTG::invalid(),
                                     std::logic_error,
                                     "Failed to convert \"local\" row index "
                                     << localRowIndex << " into a global row "
                                     "index.  Please report this bug to the "
                                     "Tpetra developers.");
                  ArrayView<const LO> ind;
                  ArrayView<const scalar_type> val;

                  newMatrix->getLocalRowView (localRowIndex, ind, val);
                  typename ArrayView<const LO>::const_iterator indIter =
                    ind.begin();
                  typename ArrayView<const scalar_type>::const_iterator
                    valIter = val.begin();
                  for (; indIter != ind.end() && valIter != val.end();
                       ++indIter, ++valIter)
                    {
                      // Convert from local to global index.
                      const GO globalColIndex =
                        pColMap->getGlobalElement (*indIter);
                      TEUCHOS_TEST_FOR_EXCEPTION(globalColIndex == OTG::invalid(),
                                         std::logic_error,
                                         "Failed to convert \"local\" column "
                                         "index " << *indIter << " into a "
                                         "global column index.  Please report "
                                         "this bug to the Tpetra developers.");
                      // Matrix Market files use 1-based indices.
                      out << (globalRowIndex + 1 - rowIndexBase) << " "
                          << (globalColIndex + 1 - colIndexBase) << " ";
                      if (STS::isComplex)
                        out << STS::real(*valIter) << " "
                            << STS::imag(*valIter);
                      else
                        out << *valIter;
                      out << endl;
                    }
                }
            }
        }
      }

      static void
      writeSparse (std::ostream& out,
                   const RCP<const sparse_matrix_type>& pMatrix,
                   const bool debug=false)
      {
        writeSparse (out, pMatrix, "", "", debug);
      }

      static void
      writeDenseFile (const std::string& filename,
                      const RCP<const multivector_type>& X,
                      const std::string& matrixName,
                      const std::string& matrixDescription)
      {
        const int myRank = Teuchos::rank (*(X->getMap()->getComm()));
        std::ofstream out;

        if (myRank == 0) // Only open the file on Rank 0.
          out.open (filename.c_str());

        writeDense (out, X, matrixName, matrixDescription);
        // We can rely on the destructor of the output stream to close
        // the file on scope exit, even if writeDense() throws an
        // exception.
      }

      static void
      writeDenseFile (const std::string& filename,
                      const RCP<const multivector_type>& X)
      {
        writeDenseFile (filename, X, "", "");
      }

      static void
      writeDense (std::ostream& out,
                  const RCP<const multivector_type>& X,
                  const std::string& matrixName,
                  const std::string& matrixDescription)
      {
        using Teuchos::rcp;
        using Teuchos::rcp_const_cast;
        using std::endl;

        typedef typename Teuchos::ScalarTraits<scalar_type> STS;
        typedef typename STS::magnitudeType magnitude_type;
        typedef typename Teuchos::ScalarTraits<magnitude_type> STM;
        typedef typename ArrayView<scalar_type>::size_type size_type;

        // Convenient abbreviations.
        typedef local_ordinal_type LO;
        typedef global_ordinal_type GO;
        typedef node_type NT;

        // If true, when making the "gather" Map, use the same index
        // base as X's Map.  Otherwise, just use the default index
        // base for the gather Map.
        const bool keepTheSameIndexBase = false;

        // Make the output stream write floating-point numbers in
        // scientific notation.  It will politely put the output
        // stream back to its state on input, when this scope
        // terminates.
        Teuchos::MatrixMarket::details::SetScientific<scalar_type> sci (out);

        RCP<const Comm<int> > comm = X->getMap()->getComm();
        const int myRank = comm->getRank ();
        RCP<const map_type> map = X->getMap();
        const global_size_t numRows = map->getGlobalNumElements();
        // Promote to global_size_t.
        const global_size_t numCols = X->getNumVectors();

        // Make the "gather" map, where Proc 0 owns all rows of X, and
        // the other procs own no rows.
        const size_t localNumRows = (myRank == 0) ? numRows : 0;
        RCP<NT> node = map->getNode();
        RCP<const map_type> gatherMap = keepTheSameIndexBase ?
          rcp_const_cast<const map_type> (rcp (new map_type (numRows, localNumRows, map->getIndexBase(), comm, node))) :
          createContigMapWithNode<LO, GO, NT> (numRows, localNumRows, comm, node);

        // Create an Import object to import X's data into a
        // multivector Y owned entirely by Proc 0.  In the Import
        // constructor, X's map is the source map, and the "gather
        // map" is the target map.
        Import<LO, GO, NT> importer (map, gatherMap);

        // Create a new multivector Y to hold the result of the import.
        RCP<multivector_type> Y =
          createMultiVector<scalar_type, LO, GO, NT> (gatherMap, numCols);

        // Import the multivector onto Proc 0.
        Y->doImport (*X, importer, INSERT);

        //
        // Print the matrix in Matrix Market format on Rank 0.
        //
        if (myRank == 0) {
          // Print the Matrix Market banner line.  MultiVector
          // stores data nonsymmetrically, hence "general".
          out << "%%MatrixMarket matrix array "
              << (STS::isComplex ? "complex" : "real")
              << " general" << endl;

          // Print comments (the matrix name and / or description).
          if (matrixName != "") {
            printAsComment (out, matrixName);
          }
          if (matrixDescription != "") {
            printAsComment (out, matrixDescription);
          }
          // Print the Matrix Market dimensions header for dense matrices.
          out << numRows << " " << numCols << endl;
          //
          // Get a read-only view of the entries of Y.
          // Rank 0 owns all of the entries of Y.
          //
          // Y must have constant stride, since multivectors have
          // constant stride if created with a map and number of
          // vectors.  This should be the case even if X does not
          // have constant stride.  However, we check Y, just to
          // make sure.
          TEUCHOS_TEST_FOR_EXCEPTION(! Y->isConstantStride (), std::logic_error,
            "The multivector Y imported onto Rank 0 does not have constant "
            "stride.  Please report this bug to the Tpetra developers.");
          ArrayRCP<const scalar_type> Y_view = Y->get1dView();
          //
          // Print the entries of the matrix, in column-major order.
          //
          const global_size_t stride = Y->getStride ();
          for (global_size_t j = 0; j < numCols; ++j) {
            for (global_size_t i = 0; i < numRows; ++i) {
              const scalar_type Y_ij = Y_view[i + j*stride];
              if (STS::isComplex) {
                out << STS::real(Y_ij) << " " << STS::imag(Y_ij) << endl;
              }
              else {
                out << Y_ij << endl;
              }
            }
          }
        } // if (myRank == 0)
      }

      static void
      writeDense (std::ostream& out,
                  const RCP<const multivector_type>& X)
      {
        writeDense (out, X, "", "");
      }

    private:
      static void
      printAsComment (std::ostream& out, const std::string& str)
      {
        using std::endl;
        std::istringstream istream (str);
        std::string line;

        while (getline (istream, line)) {
          if (! line.empty()) {
            // Note that getline() doesn't store '\n', so we have to
            // append the endline ourselves.
            if (line[0] == '%') { // Line starts with a comment character.
              out << line << endl;
            }
            else { // Line doesn't start with a comment character.
              out << "%% " << line << endl;
            }
          }
        }
      }
    }; // class Writer

  } // namespace MatrixMarket
} // namespace Tpetra

#endif // __MatrixMarket_Tpetra_hpp