Tpetra Matrix/Vector Services Version of the Day
MatrixMarket_Tpetra.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef __MatrixMarket_Tpetra_hpp
00043 #define __MatrixMarket_Tpetra_hpp
00044 
00056 #include "Tpetra_config.h"
00057 #include "Tpetra_CrsMatrix.hpp"
00058 #include "Tpetra_Vector.hpp"
00059 #include "Tpetra_ComputeGatherMap.hpp"
00060 #include "Teuchos_MatrixMarket_Raw_Adder.hpp"
00061 #include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp"
00062 #include "Teuchos_MatrixMarket_assignScalar.hpp"
00063 #include "Teuchos_MatrixMarket_Banner.hpp"
00064 #include "Teuchos_MatrixMarket_CoordDataReader.hpp"
00065 #include "Teuchos_MatrixMarket_SetScientific.hpp"
00066 
00067 #include <algorithm>
00068 #include <fstream>
00069 #include <iostream>
00070 #include <iterator>
00071 #include <vector>
00072 #include <stdexcept>
00073 
00074 
00075 // Macro that marks a function as "possibly unused," in order to
00076 // suppress build warnings.
00077 #if ! defined(TRILINOS_UNUSED_FUNCTION)
00078 #  if defined(__GNUC__) || defined(__INTEL_COMPILER)
00079 #    define TRILINOS_UNUSED_FUNCTION __attribute__((__unused__))
00080 #  elif defined(__clang__)
00081 #    if __has_attribute(unused)
00082 #      define TRILINOS_UNUSED_FUNCTION __attribute__((__unused__))
00083 #    else
00084 #      define TRILINOS_UNUSED_FUNCTION
00085 #    endif // Clang has 'unused' attribute
00086 #  elif defined(__IBMCPP__)
00087 // IBM's C++ compiler for Blue Gene/Q (V12.1) implements 'used' but not 'unused'.
00088 //
00089 // http://pic.dhe.ibm.com/infocenter/compbg/v121v141/index.jsp
00090 #    define TRILINOS_UNUSED_FUNCTION
00091 #  else // some other compiler
00092 #    define TRILINOS_UNUSED_FUNCTION
00093 #  endif
00094 #endif // ! defined(TRILINOS_UNUSED_FUNCTION)
00095 
00096 
00097 namespace Tpetra {
00126   namespace MatrixMarket {
00127 
00128     namespace {
00129 #ifdef HAVE_MPI
00130       // We're communicating integers of type IntType.  Figure
00131       // out the right MPI_Datatype for IntType.  Usually it
00132       // is int or long, so these are good enough for now.
00133       template<class IntType> TRILINOS_UNUSED_FUNCTION MPI_Datatype
00134       getMpiDatatype () {
00135         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00136           "Not implemented for IntType != int, long, or long long");
00137       }
00138 
00139       template<> TRILINOS_UNUSED_FUNCTION MPI_Datatype
00140       getMpiDatatype<int> () { return MPI_INT; }
00141 
00142       template<> TRILINOS_UNUSED_FUNCTION MPI_Datatype
00143       getMpiDatatype<long> () { return MPI_LONG; }
00144 
00145       template<> TRILINOS_UNUSED_FUNCTION MPI_Datatype
00146       getMpiDatatype<long long> () { return MPI_LONG_LONG; }
00147 #endif // HAVE_MPI
00148 
00149       template<class IntType>
00150       void
00151       reduceSum (const IntType sendbuf[],
00152                  IntType recvbuf[],
00153                  const int count,
00154                  const int root,
00155                  const Teuchos::RCP<const Teuchos::Comm<int> >& comm)
00156       {
00157 #ifdef HAVE_MPI
00158         using Teuchos::RCP;
00159         using Teuchos::rcp_dynamic_cast;
00160         using Teuchos::MpiComm;
00161 
00162         // Get the raw MPI communicator, so we don't have to wrap MPI_Reduce in Teuchos.
00163         RCP<const MpiComm<int> > mpiComm = rcp_dynamic_cast<const MpiComm<int> > (comm);
00164         if (! mpiComm.is_null ()) {
00165           MPI_Datatype rawMpiType = getMpiDatatype<IntType> ();
00166           MPI_Comm rawMpiComm = * (mpiComm->getRawMpiComm ());
00167           const int err =
00168             MPI_Reduce (reinterpret_cast<void*> (const_cast<IntType*> (sendbuf)),
00169                         reinterpret_cast<void*> (const_cast<IntType*> (recvbuf)),
00170                         count, rawMpiType, MPI_SUM, root, rawMpiComm);
00171           TEUCHOS_TEST_FOR_EXCEPTION(err != MPI_SUCCESS, std::runtime_error, "MPI_Reduce failed");
00172           return;
00173         }
00174 #endif // HAVE_MPI
00175         // Serial communicator case: just copy.  count is the
00176         // amount to receive, so it's the amount to copy.
00177         std::copy (sendbuf, sendbuf + count, recvbuf);
00178       }
00179 
00180     } // namespace (anonymous)
00181 
00245     template<class SparseMatrixType>
00246     class Reader {
00247     public:
00249       typedef SparseMatrixType sparse_matrix_type;
00250       typedef RCP<sparse_matrix_type> sparse_matrix_ptr;
00251 
00254       typedef typename SparseMatrixType::scalar_type scalar_type;
00257       typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
00265       typedef typename SparseMatrixType::global_ordinal_type
00266         global_ordinal_type;
00269       typedef typename SparseMatrixType::node_type node_type;
00270 
00272       typedef MultiVector<scalar_type,
00273                           local_ordinal_type,
00274                           global_ordinal_type,
00275                           node_type> multivector_type;
00276 
00278       typedef Vector<scalar_type,
00279                           local_ordinal_type,
00280                           global_ordinal_type,
00281                           node_type> vector_type;
00282 
00283       typedef RCP<node_type> node_ptr;
00284       typedef Comm<int> comm_type;
00285       typedef RCP<const comm_type> comm_ptr;
00286       typedef Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00287       typedef RCP<const map_type> map_ptr;
00288 
00289     private:
00292       typedef typename ArrayRCP<global_ordinal_type>::size_type size_type;
00293 
00305       static RCP<const map_type>
00306       makeRangeMap (const RCP<const comm_type>& pComm,
00307                     const RCP<node_type>& pNode,
00308                     const global_ordinal_type numRows)
00309       {
00310         // A conventional, uniformly partitioned, contiguous map.
00311         return rcp (new map_type (static_cast<global_size_t> (numRows),
00312                                   static_cast<global_ordinal_type> (0),
00313                                   pComm, GloballyDistributed, pNode));
00314       }
00315 
00344       static RCP<const map_type>
00345       makeRowMap (const RCP<const map_type>& pRowMap,
00346                   const RCP<const comm_type>& pComm,
00347                   const RCP<node_type>& pNode,
00348                   const global_ordinal_type numRows)
00349       {
00350         // If the caller didn't provide a map, return a conventional,
00351         // uniformly partitioned, contiguous map.
00352         if (pRowMap.is_null()) {
00353           return rcp (new map_type (static_cast<global_size_t> (numRows),
00354                                     static_cast<global_ordinal_type> (0),
00355                                     pComm, GloballyDistributed, pNode));
00356         } else {
00357           TEUCHOS_TEST_FOR_EXCEPTION(! pRowMap->isDistributed() && pComm->getSize() > 1,
00358                              std::invalid_argument,
00359                              "The specified row map is not distributed, but "
00360                              "the given communicator includes more than one "
00361                              "rank (in fact, there are " << pComm->getSize()
00362                              << " ranks).");
00363           TEUCHOS_TEST_FOR_EXCEPTION(pRowMap->getComm() != pComm,
00364                              std::invalid_argument,
00365                              "The specified row map's communicator (pRowMap->"
00366                              "getComm()) is different than the given separately "
00367                              "supplied communicator pComm.");
00368           TEUCHOS_TEST_FOR_EXCEPTION(pRowMap->getNode() != pNode,
00369                              std::invalid_argument,
00370                              "The specified row map's node (pRowMap->getNode()) "
00371                              "is different than the given separately supplied "
00372                              "node pNode.");
00373           return pRowMap;
00374         }
00375       }
00376 
00391       static Teuchos::RCP<const map_type>
00392       makeDomainMap (const Teuchos::RCP<const map_type>& pRangeMap,
00393                      const global_ordinal_type numRows,
00394                      const global_ordinal_type numCols)
00395       {
00396         // Abbreviations so that the map creation call isn't too long.
00397         typedef local_ordinal_type LO;
00398         typedef global_ordinal_type GO;
00399         typedef node_type Node;
00400 
00401         if (numRows == numCols) {
00402           return pRangeMap;
00403         } else {
00404           return createUniformContigMapWithNode<LO,GO,Node> (numCols,
00405                                                              pRangeMap->getComm (),
00406                                                              pRangeMap->getNode ());
00407         }
00408       }
00409 
00482       static void
00483       distribute (ArrayRCP<size_t>& myNumEntriesPerRow,
00484                   ArrayRCP<size_t>& myRowPtr,
00485                   ArrayRCP<global_ordinal_type>& myColInd,
00486                   ArrayRCP<scalar_type>& myValues,
00487                   const RCP<const map_type>& pRowMap,
00488                   ArrayRCP<size_t>& numEntriesPerRow,
00489                   ArrayRCP<size_t>& rowPtr,
00490                   ArrayRCP<global_ordinal_type>& colInd,
00491                   ArrayRCP<scalar_type>& values,
00492                   const bool debug=false)
00493       {
00494          using Teuchos::as;
00495          using Teuchos::CommRequest;
00496          using Teuchos::receive;
00497          using Teuchos::send;
00498          using std::cerr;
00499          using std::endl;
00500 
00501          const bool extraDebug = false;
00502          comm_ptr pComm = pRowMap->getComm();
00503          const int numProcs = pComm->getSize ();
00504          const int myRank = pComm->getRank ();
00505          const int rootRank = 0;
00506 
00507          // Type abbreviations to make the code more concise.
00508          typedef global_ordinal_type GO;
00509 
00510          // List of the global indices of my rows.  They may or may
00511          // not be contiguous, and the row map need not be one-to-one.
00512          ArrayView<const GO> myRows = pRowMap->getNodeElementList();
00513          const size_type myNumRows = myRows.size();
00514          TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
00515                             pRowMap->getNodeNumElements(),
00516                             std::logic_error,
00517                             "pRowMap->getNodeElementList().size() = "
00518                             << myNumRows
00519                             << " != pRowMap->getNodeNumElements() = "
00520                             << pRowMap->getNodeNumElements() << ".  "
00521                             "Please report this bug to the Tpetra developers.");
00522          TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
00523                             std::logic_error,
00524                             "On Proc 0: numEntriesPerRow.size() = "
00525                             << numEntriesPerRow.size()
00526                             << " != pRowMap->getNodeElementList().size() = "
00527                             << myNumRows << ".  Please report this bug to the "
00528                             "Tpetra developers.");
00529 
00530          // Space for my proc's number of entries per row.  Will be
00531          // filled in below.
00532          myNumEntriesPerRow = arcp<size_t> (myNumRows);
00533 
00534          if (myRank != rootRank) {
00535            // Tell the root how many rows we have.  If we're sending
00536            // none, then we don't have anything else to send, nor does
00537            // the root have to receive anything else.
00538            send (*pComm, myNumRows, rootRank);
00539            if (myNumRows != 0) {
00540              // Now send my rows' global indices.  Hopefully the cast
00541              // to int doesn't overflow.  This is unlikely, since it
00542              // should fit in a LO, even though it is a GO.
00543              send (*pComm, static_cast<int> (myNumRows),
00544                    myRows.getRawPtr(), rootRank);
00545 
00546              // I (this proc) don't care if my global row indices are
00547              // contiguous, though the root proc does (since otherwise
00548              // it needs to pack noncontiguous data into contiguous
00549              // storage before sending).  That's why we don't check
00550              // for contiguousness here.
00551 
00552              // Ask the root process for my part of the array of the
00553              // number of entries per row.
00554              receive (*pComm, rootRank,
00555                       static_cast<int> (myNumRows),
00556                       myNumEntriesPerRow.getRawPtr());
00557 
00558              // Use the resulting array to figure out how many column
00559              // indices and values I should ask from the root process.
00560              const local_ordinal_type myNumEntries =
00561                std::accumulate (myNumEntriesPerRow.begin(),
00562                                 myNumEntriesPerRow.end(), 0);
00563 
00564              // Make space for my entries of the sparse matrix.  Note
00565              // that they don't have to be sorted by row index.
00566              // Iterating through all my rows requires computing a
00567              // running sum over myNumEntriesPerRow.
00568              myColInd = arcp<GO> (myNumEntries);
00569              myValues = arcp<scalar_type> (myNumEntries);
00570              if (myNumEntries > 0) {
00571                // Ask for that many column indices and values, if
00572                // there are any.
00573                receive (*pComm, rootRank,
00574                         static_cast<int> (myNumEntries),
00575                         myColInd.getRawPtr());
00576                receive (*pComm, rootRank,
00577                         static_cast<int> (myNumEntries),
00578                         myValues.getRawPtr());
00579              }
00580            } // If I own at least one row
00581          } // If I am not the root processor
00582          else { // I _am_ the root processor
00583            if (debug) {
00584              cerr << "-- Proc 0: Copying my data from global arrays" << endl;
00585            }
00586            // Proc 0 still needs to (allocate, if not done already)
00587            // and fill its part of the matrix (my*).
00588            for (size_type k = 0; k < myNumRows; ++k) {
00589              const GO myCurRow = myRows[k];
00590              const local_ordinal_type numEntriesInThisRow = numEntriesPerRow[myCurRow];
00591              myNumEntriesPerRow[k] = numEntriesInThisRow;
00592            }
00593            if (extraDebug && debug) {
00594              cerr << "Proc " << pRowMap->getComm ()->getRank ()
00595                   << ": myNumEntriesPerRow[0.." << (myNumRows-1) << "] = [";
00596              for (size_type k = 0; k < myNumRows; ++k) {
00597                cerr << myNumEntriesPerRow[k];
00598                if (k < myNumRows-1) {
00599                  cerr << " ";
00600                }
00601              }
00602              cerr << "]" << endl;
00603            }
00604            // The total number of matrix entries that my proc owns.
00605            const local_ordinal_type myNumEntries =
00606              std::accumulate (myNumEntriesPerRow.begin(),
00607                               myNumEntriesPerRow.end(), 0);
00608            if (debug) {
00609              cerr << "-- Proc 0: I own " << myNumRows << " rows and "
00610                   << myNumEntries << " entries" << endl;
00611            }
00612            myColInd = arcp<GO> (myNumEntries);
00613            myValues = arcp<scalar_type> (myNumEntries);
00614 
00615            // Copy Proc 0's part of the matrix into the my* arrays.
00616            // It's important that myCurPos be updated _before_ k,
00617            // otherwise myCurPos will get the wrong number of entries
00618            // per row (it should be for the row in the just-completed
00619            // iteration, not for the next iteration's row).
00620            local_ordinal_type myCurPos = 0;
00621            for (size_type k = 0; k < myNumRows;
00622                 myCurPos += myNumEntriesPerRow[k], ++k) {
00623              const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
00624              const GO myRow = myRows[k];
00625              const size_t curPos = rowPtr[myRow];
00626              // Only copy if there are entries to copy, in order not
00627              // to construct empty ranges for the ArrayRCP views.
00628              if (curNumEntries > 0) {
00629                ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
00630                ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
00631                std::copy (colIndView.begin(), colIndView.end(),
00632                           myColIndView.begin());
00633 
00634                ArrayView<scalar_type> valuesView =
00635                  values (curPos, curNumEntries);
00636                ArrayView<scalar_type> myValuesView =
00637                  myValues (myCurPos, curNumEntries);
00638                std::copy (valuesView.begin(), valuesView.end(),
00639                           myValuesView.begin());
00640              }
00641            }
00642 
00643            // Proc 0 processes each other proc p in turn.
00644            for (int p = 1; p < numProcs; ++p) {
00645              if (debug) {
00646                cerr << "-- Proc 0: Processing proc " << p << endl;
00647              }
00648 
00649              size_type theirNumRows = 0;
00650              // Ask Proc p how many rows it has.  If it doesn't
00651              // have any, we can move on to the next proc.  This
00652              // has to be a standard receive so that we can avoid
00653              // the degenerate case of sending zero data.
00654              receive (*pComm, p, &theirNumRows);
00655              if (debug) {
00656                cerr << "-- Proc 0: Proc " << p << " owns "
00657                     << theirNumRows << " rows" << endl;
00658              }
00659              if (theirNumRows != 0) {
00660                // Ask Proc p which rows it owns.  The resulting global
00661                // row indices are not guaranteed to be contiguous or
00662                // sorted.  Global row indices are themselves indices
00663                // into the numEntriesPerRow array.
00664                ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
00665                receive (*pComm, p, as<int> (theirNumRows),
00666                         theirRows.getRawPtr ());
00667                // Extra test to make sure that the rows we received
00668                // are all sensible.  This is a good idea since we are
00669                // going to use the global row indices we've received
00670                // to index into the numEntriesPerRow array.  Better to
00671                // catch any bugs here and print a sensible error
00672                // message, rather than segfault and print a cryptic
00673                // error message.
00674                {
00675                  const global_size_t numRows = pRowMap->getGlobalNumElements ();
00676                  const GO indexBase = pRowMap->getIndexBase ();
00677                  bool theirRowsValid = true;
00678                  for (size_type k = 0; k < theirNumRows; ++k) {
00679                    if (theirRows[k] < indexBase ||
00680                        as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
00681                      theirRowsValid = false;
00682                    }
00683                  }
00684                  if (! theirRowsValid) {
00685                    TEUCHOS_TEST_FOR_EXCEPTION(
00686                      ! theirRowsValid, std::logic_error,
00687                      "Proc " << p << " has at least one invalid row index.  "
00688                      "Here are all of them: " <<
00689                      Teuchos::toString (theirRows ()) << ".  Valid row index "
00690                      "range (zero-based): [0, " << (numRows - 1) << "].");
00691                  }
00692                }
00693 
00694                // Perhaps we could save a little work if we check
00695                // whether Proc p's row indices are contiguous.  That
00696                // would make lookups in the global data arrays
00697                // faster.  For now, we just implement the general
00698                // case and don't prematurely optimize.  (Remember
00699                // that you're making Proc 0 read the whole file, so
00700                // you've already lost scalability.)
00701 
00702                // Compute the number of entries in each of Proc p's
00703                // rows.  (Proc p will compute its row pointer array
00704                // on its own, after it gets the data from Proc 0.)
00705                ArrayRCP<size_t> theirNumEntriesPerRow;
00706                theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
00707                for (size_type k = 0; k < theirNumRows; ++k) {
00708                  theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
00709                }
00710 
00711                // Tell Proc p the number of entries in each of its
00712                // rows.  Hopefully the cast to int doesn't overflow.
00713                // This is unlikely, since it should fit in a LO,
00714                // even though it is a GO.
00715                send (*pComm, static_cast<int> (theirNumRows),
00716                      theirNumEntriesPerRow.getRawPtr(), p);
00717 
00718                // Figure out how many entries Proc p owns.
00719                const local_ordinal_type theirNumEntries =
00720                  std::accumulate (theirNumEntriesPerRow.begin(),
00721                                   theirNumEntriesPerRow.end(), 0);
00722 
00723                if (debug) {
00724                  cerr << "-- Proc 0: Proc " << p << " owns "
00725                       << theirNumEntries << " entries" << endl;
00726                }
00727 
00728                // If there are no entries to send, then we're done
00729                // with Proc p.
00730                if (theirNumEntries == 0) {
00731                  continue;
00732                }
00733 
00734                // Construct (views of) proc p's column indices and
00735                // values.  Later, we might like to optimize for the
00736                // (common) contiguous case, for which we don't need to
00737                // copy data into separate "their*" arrays (we can just
00738                // use contiguous views of the global arrays).
00739                ArrayRCP<GO> theirColInd (theirNumEntries);
00740                ArrayRCP<scalar_type> theirValues (theirNumEntries);
00741                // Copy Proc p's part of the matrix into the their*
00742                // arrays.  It's important that theirCurPos be updated
00743                // _before_ k, otherwise theirCurPos will get the wrong
00744                // number of entries per row (it should be for the row
00745                // in the just-completed iteration, not for the next
00746                // iteration's row).
00747                local_ordinal_type theirCurPos = 0;
00748                for (size_type k = 0; k < theirNumRows;
00749                     theirCurPos += theirNumEntriesPerRow[k], k++) {
00750                  const local_ordinal_type curNumEntries = theirNumEntriesPerRow[k];
00751                  const GO theirRow = theirRows[k];
00752                  const local_ordinal_type curPos = rowPtr[theirRow];
00753 
00754                  // Only copy if there are entries to copy, in order
00755                  // not to construct empty ranges for the ArrayRCP
00756                  // views.
00757                  if (curNumEntries > 0) {
00758                    ArrayView<GO> colIndView =
00759                      colInd (curPos, curNumEntries);
00760                    ArrayView<GO> theirColIndView =
00761                      theirColInd (theirCurPos, curNumEntries);
00762                    std::copy (colIndView.begin(), colIndView.end(),
00763                               theirColIndView.begin());
00764 
00765                    ArrayView<scalar_type> valuesView =
00766                      values (curPos, curNumEntries);
00767                    ArrayView<scalar_type> theirValuesView =
00768                      theirValues (theirCurPos, curNumEntries);
00769                    std::copy (valuesView.begin(), valuesView.end(),
00770                               theirValuesView.begin());
00771                  }
00772                }
00773                // Send Proc p its column indices and values.
00774                // Hopefully the cast to int doesn't overflow.  This
00775                // is unlikely, since it should fit in a LO, even
00776                // though it is a GO.
00777                send (*pComm, static_cast<int> (theirNumEntries),
00778                      theirColInd.getRawPtr(), p);
00779                send (*pComm, static_cast<int> (theirNumEntries),
00780                      theirValues.getRawPtr(), p);
00781 
00782                if (debug) {
00783                  cerr << "-- Proc 0: Finished with proc " << p << endl;
00784                }
00785              } // If proc p owns at least one row
00786            } // For each proc p not the root proc 0
00787          } // If I'm (not) the root proc 0
00788 
00789          // Invalidate the input data to save space, since we don't
00790          // need it anymore.
00791          numEntriesPerRow = null;
00792          rowPtr = null;
00793          colInd = null;
00794          values = null;
00795 
00796          if (debug && myRank == 0) {
00797            cerr << "-- Proc 0: About to fill in myRowPtr" << endl;
00798          }
00799 
00800          // Allocate and fill in myRowPtr (the row pointer array for
00801          // my rank's rows).  We delay this until the end because we
00802          // don't need it to compute anything else in distribute().
00803          // Each proc can do this work for itself, since it only needs
00804          // myNumEntriesPerRow to do so.
00805          myRowPtr = arcp<size_t> (myNumRows+1);
00806          myRowPtr[0] = 0;
00807          for (size_type k = 1; k < myNumRows+1; ++k) {
00808            myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
00809          }
00810          if (extraDebug && debug) {
00811            cerr << "Proc " << Teuchos::rank (*(pRowMap->getComm()))
00812                 << ": myRowPtr[0.." << myNumRows << "] = [";
00813            for (size_type k = 0; k < myNumRows+1; ++k) {
00814              cerr << myRowPtr[k];
00815              if (k < myNumRows) {
00816                cerr << " ";
00817              }
00818            }
00819            cerr << "]" << endl << endl;
00820          }
00821 
00822          if (debug && myRank == 0) {
00823            cerr << "-- Proc 0: Done with distribute" << endl;
00824          }
00825       }
00826 
00840       static sparse_matrix_ptr
00841       makeMatrix (ArrayRCP<size_t>& myNumEntriesPerRow,
00842                   ArrayRCP<size_t>& myRowPtr,
00843                   ArrayRCP<global_ordinal_type>& myColInd,
00844                   ArrayRCP<scalar_type>& myValues,
00845                   const map_ptr& pRowMap,
00846                   const map_ptr& pRangeMap,
00847                   const map_ptr& pDomainMap,
00848                   const bool callFillComplete = true)
00849       {
00850         using std::cerr;
00851         using std::endl;
00852         // Typedef to make certain type declarations shorter.
00853         typedef global_ordinal_type GO;
00854 
00855         // The row pointer array always has at least one entry, even
00856         // if the matrix has zero rows.  myNumEntriesPerRow, myColInd,
00857         // and myValues would all be empty arrays in that degenerate
00858         // case, but the row and domain maps would still be nonnull
00859         // (though they would be trivial maps).
00860         TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
00861           "makeMatrix: myRowPtr array is null.  "
00862           "Please report this bug to the Tpetra developers.");
00863         TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
00864           "makeMatrix: domain map is null.  "
00865           "Please report this bug to the Tpetra developers.");
00866         TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
00867           "makeMatrix: range map is null.  "
00868           "Please report this bug to the Tpetra developers.");
00869         TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
00870           "makeMatrix: row map is null.  "
00871           "Please report this bug to the Tpetra developers.");
00872 
00873         // Construct the CrsMatrix, using the row map, with the
00874         // constructor specifying the number of nonzeros for each row.
00875         // Create with DynamicProfile, so that the fillComplete() can
00876         // do first-touch reallocation (a NUMA (Non-Uniform Memory
00877         // Access) optimization on multicore CPUs).
00878         sparse_matrix_ptr A =
00879           rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
00880                                        DynamicProfile));
00881 
00882         // List of the global indices of my rows.
00883         // They may or may not be contiguous.
00884         ArrayView<const GO> myRows = pRowMap->getNodeElementList ();
00885         const size_type myNumRows = myRows.size ();
00886 
00887         // Add this processor's matrix entries to the CrsMatrix.
00888         const GO indexBase = pRowMap->getIndexBase ();
00889         for (size_type i = 0; i < myNumRows; ++i) {
00890           const size_type myCurPos = myRowPtr[i];
00891           const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
00892           ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
00893           ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
00894 
00895           // Modify the column indices in place to have the right index base.
00896           for (size_type k = 0; k < curNumEntries; ++k) {
00897             curColInd[k] += indexBase;
00898           }
00899           // Avoid constructing empty views of ArrayRCP objects.
00900           if (curNumEntries > 0) {
00901             A->insertGlobalValues (myRows[i], curColInd, curValues);
00902           }
00903         }
00904         // We've entered in all our matrix entries, so we can delete
00905         // the original data.  This will save memory when we call
00906         // fillComplete(), so that we never keep more than two copies
00907         // of the matrix's data in memory at once.
00908         myNumEntriesPerRow = null;
00909         myRowPtr = null;
00910         myColInd = null;
00911         myValues = null;
00912 
00913         if (callFillComplete) {
00914           A->fillComplete (pDomainMap, pRangeMap);
00915         }
00916         return A;
00917       }
00918 
00924       static sparse_matrix_ptr
00925       makeMatrix (ArrayRCP<size_t>& myNumEntriesPerRow,
00926                   ArrayRCP<size_t>& myRowPtr,
00927                   ArrayRCP<global_ordinal_type>& myColInd,
00928                   ArrayRCP<scalar_type>& myValues,
00929                   const map_ptr& pRowMap,
00930                   const map_ptr& pRangeMap,
00931                   const map_ptr& pDomainMap,
00932                   const RCP<Teuchos::ParameterList>& constructorParams,
00933                   const RCP<Teuchos::ParameterList>& fillCompleteParams)
00934       {
00935         using std::cerr;
00936         using std::endl;
00937         // Typedef to make certain type declarations shorter.
00938         typedef global_ordinal_type GO;
00939 
00940         // The row pointer array always has at least one entry, even
00941         // if the matrix has zero rows.  myNumEntriesPerRow, myColInd,
00942         // and myValues would all be empty arrays in that degenerate
00943         // case, but the row and domain maps would still be nonnull
00944         // (though they would be trivial maps).
00945         TEUCHOS_TEST_FOR_EXCEPTION(
00946           myRowPtr.is_null(), std::logic_error,
00947           "makeMatrix: myRowPtr array is null.  "
00948           "Please report this bug to the Tpetra developers.");
00949         TEUCHOS_TEST_FOR_EXCEPTION(
00950           pDomainMap.is_null(), std::logic_error,
00951           "makeMatrix: domain map is null.  "
00952           "Please report this bug to the Tpetra developers.");
00953         TEUCHOS_TEST_FOR_EXCEPTION(
00954           pRangeMap.is_null(), std::logic_error,
00955           "makeMatrix: range map is null.  "
00956           "Please report this bug to the Tpetra developers.");
00957         TEUCHOS_TEST_FOR_EXCEPTION(
00958           pRowMap.is_null(), std::logic_error,
00959           "makeMatrix: row map is null.  "
00960           "Please report this bug to the Tpetra developers.");
00961 
00962         // Construct the CrsMatrix, using the row map, with the
00963         // constructor specifying the number of nonzeros for each row.
00964         // Create with DynamicProfile, so that the fillComplete() can
00965         // do first-touch reallocation (a NUMA (Non-Uniform Memory
00966         // Access) optimization on multicore CPUs).
00967         sparse_matrix_ptr A =
00968           rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
00969                                        DynamicProfile, constructorParams));
00970 
00971         // List of the global indices of my rows.
00972         // They may or may not be contiguous.
00973         ArrayView<const GO> myRows = pRowMap->getNodeElementList();
00974         const size_type myNumRows = myRows.size();
00975 
00976         // Add this processor's matrix entries to the CrsMatrix.
00977         const GO indexBase = pRowMap->getIndexBase ();
00978         for (size_type i = 0; i < myNumRows; ++i) {
00979           const size_type myCurPos = myRowPtr[i];
00980           const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
00981           ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
00982           ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
00983 
00984           // Modify the column indices in place to have the right index base.
00985           for (size_type k = 0; k < curNumEntries; ++k) {
00986             curColInd[k] += indexBase;
00987           }
00988           if (curNumEntries > 0) {
00989             A->insertGlobalValues (myRows[i], curColInd, curValues);
00990           }
00991         }
00992         // We've entered in all our matrix entries, so we can delete
00993         // the original data.  This will save memory when we call
00994         // fillComplete(), so that we never keep more than two copies
00995         // of the matrix's data in memory at once.
00996         myNumEntriesPerRow = null;
00997         myRowPtr = null;
00998         myColInd = null;
00999         myValues = null;
01000 
01001         A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
01002         return A;
01003       }
01004 
01009       static sparse_matrix_ptr
01010       makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
01011                   Teuchos::ArrayRCP<size_t>& myRowPtr,
01012                   Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
01013                   Teuchos::ArrayRCP<scalar_type>& myValues,
01014                   const Teuchos::RCP<const map_type>& rowMap,
01015                   Teuchos::RCP<const map_type>& colMap,
01016                   const Teuchos::RCP<const map_type>& domainMap,
01017                   const Teuchos::RCP<const map_type>& rangeMap,
01018                   const bool callFillComplete = true)
01019       {
01020         using Teuchos::ArrayView;
01021         using Teuchos::as;
01022         using Teuchos::null;
01023         using Teuchos::RCP;
01024         using Teuchos::rcp;
01025         typedef global_ordinal_type GO;
01026         typedef typename ArrayView<const GO>::size_type size_type;
01027 
01028         // Construct the CrsMatrix.
01029         //
01030         // Create with DynamicProfile, so that the fillComplete() can
01031         // do first-touch reallocation.
01032         RCP<sparse_matrix_type> A; // the matrix to return.
01033         if (colMap.is_null ()) { // the user didn't provide a column Map
01034           A = rcp (new sparse_matrix_type (rowMap, myNumEntriesPerRow, DynamicProfile));
01035         } else { // the user provided a column Map
01036           A = rcp (new sparse_matrix_type (rowMap, colMap, myNumEntriesPerRow, DynamicProfile));
01037         }
01038 
01039         // List of the global indices of my rows.
01040         // They may or may not be contiguous.
01041         ArrayView<const GO> myRows = rowMap->getNodeElementList ();
01042         const size_type myNumRows = myRows.size ();
01043 
01044         // Add this process' matrix entries to the CrsMatrix.
01045         const GO indexBase = rowMap->getIndexBase ();
01046         for (size_type i = 0; i < myNumRows; ++i) {
01047           const size_type myCurPos = myRowPtr[i];
01048           const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
01049           ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
01050           ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
01051 
01052           // Modify the column indices in place to have the right index base.
01053           for (size_type k = 0; k < curNumEntries; ++k) {
01054             curColInd[k] += indexBase;
01055           }
01056           if (curNumEntries > 0) {
01057             A->insertGlobalValues (myRows[i], curColInd, curValues);
01058           }
01059         }
01060         // We've entered in all our matrix entries, so we can delete
01061         // the original data.  This will save memory when we call
01062         // fillComplete(), so that we never keep more than two copies
01063         // of the matrix's data in memory at once.
01064         myNumEntriesPerRow = null;
01065         myRowPtr = null;
01066         myColInd = null;
01067         myValues = null;
01068 
01069         if (callFillComplete) {
01070           A->fillComplete (domainMap, rangeMap);
01071           if (colMap.is_null ()) {
01072             colMap = A->getColMap ();
01073           }
01074         }
01075         return A;
01076       }
01077 
01078     private:
01079 
01096       static RCP<const Teuchos::MatrixMarket::Banner>
01097       readBanner (std::istream& in,
01098                   size_t& lineNumber,
01099                   const bool tolerant=false,
01100                   const bool debug=false)
01101       {
01102         using Teuchos::MatrixMarket::Banner;
01103         using std::cerr;
01104         using std::endl;
01105         typedef Teuchos::ScalarTraits<scalar_type> STS;
01106 
01107         RCP<Banner> pBanner; // On output, if successful: the read-in Banner.
01108         std::string line; // If read from stream successful: the Banner line
01109 
01110         // Try to read a line from the input stream.
01111         const bool readFailed = ! getline(in, line);
01112         TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
01113           "Failed to get Matrix Market banner line from input.");
01114 
01115         // We read a line from the input stream.
01116         lineNumber++;
01117 
01118         // Assume that the line we found is the Banner line.
01119         try {
01120           pBanner = rcp (new Banner (line, tolerant));
01121         } catch (std::exception& e) {
01122           TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
01123             "Matrix Market banner line contains syntax error(s): "
01124             << e.what());
01125         }
01126         TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() != "matrix",
01127           std::invalid_argument, "The Matrix Market file does not contain "
01128           "matrix data.  Its Banner (first) line says that its object type is \""
01129           << pBanner->matrixType() << "\", rather than the required \"matrix\".");
01130 
01131         // Validate the data type of the matrix, with respect to the
01132         // Scalar type of the CrsMatrix entries.
01133         TEUCHOS_TEST_FOR_EXCEPTION(
01134           ! STS::isComplex && pBanner->dataType() == "complex",
01135           std::invalid_argument,
01136           "The Matrix Market file contains complex-valued data, but you are "
01137           "trying to read it into a matrix containing entries of the real-"
01138           "valued Scalar type \""
01139           << Teuchos::TypeNameTraits<scalar_type>::name() << "\".");
01140         TEUCHOS_TEST_FOR_EXCEPTION(
01141           pBanner->dataType() != "real" &&
01142           pBanner->dataType() != "complex" &&
01143           pBanner->dataType() != "integer",
01144           std::invalid_argument,
01145           "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
01146           "Matrix Market file may not contain a \"pattern\" matrix.  A "
01147           "pattern matrix is really just a graph with no weights.  It "
01148           "should be stored in a CrsGraph, not a CrsMatrix.");
01149 
01150         return pBanner;
01151       }
01152 
01175       static Tuple<global_ordinal_type, 3>
01176       readCoordDims (std::istream& in,
01177                      size_t& lineNumber,
01178                      const RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
01179                      const comm_ptr& pComm,
01180                      const bool tolerant = false,
01181                      const bool debug = false)
01182       {
01183         using Teuchos::MatrixMarket::readCoordinateDimensions;
01184 
01185         // Packed coordinate matrix dimensions (numRows, numCols,
01186         // numNonzeros); computed on Rank 0 and broadcasted to all MPI
01187         // ranks.
01188         Tuple<global_ordinal_type, 3> dims;
01189 
01190         // Read in the coordinate matrix dimensions from the input
01191         // stream.  "success" tells us whether reading in the
01192         // coordinate matrix dimensions succeeded ("Guilty unless
01193         // proven innocent").
01194         bool success = false;
01195         if (pComm->getRank() == 0) {
01196           TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "coordinate",
01197             std::invalid_argument, "The Tpetra::CrsMatrix Matrix Market reader "
01198             "only accepts \"coordinate\" (sparse) matrix data.");
01199           // Unpacked coordinate matrix dimensions
01200           global_ordinal_type numRows, numCols, numNonzeros;
01201           // Only MPI Rank 0 reads from the input stream
01202           success = readCoordinateDimensions (in, numRows, numCols,
01203                                               numNonzeros, lineNumber,
01204                                               tolerant);
01205           // Pack up the data into a Tuple so we can send them with
01206           // one broadcast instead of three.
01207           dims[0] = numRows;
01208           dims[1] = numCols;
01209           dims[2] = numNonzeros;
01210         }
01211         // Only Rank 0 did the reading, so it decides success.
01212         //
01213         // FIXME (mfh 02 Feb 2011) Teuchos::broadcast doesn't know how
01214         // to send bools.  For now, we convert to/from int instead,
01215         // using the usual "true is 1, false is 0" encoding.
01216         {
01217           int the_success = success ? 1 : 0; // only matters on MPI Rank 0
01218           Teuchos::broadcast (*pComm, 0, &the_success);
01219           success = (the_success == 1);
01220         }
01221         if (success) {
01222           // Broadcast (numRows, numCols, numNonzeros) from Rank 0
01223           // to all the other MPI ranks.
01224           Teuchos::broadcast (*pComm, 0, dims);
01225         }
01226         else {
01227           // Perhaps in tolerant mode, we could set all the
01228           // dimensions to zero for now, and deduce correct
01229           // dimensions by reading all of the file's entries and
01230           // computing the max(row index) and max(column index).
01231           // However, for now we just error out in that case.
01232           TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
01233             "Error reading Matrix Market sparse matrix: failed to read "
01234             "coordinate matrix dimensions.");
01235         }
01236         return dims;
01237       }
01238 
01249       typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
01250 
01276       static RCP<adder_type>
01277       makeAdder (const RCP<const Comm<int> >& pComm,
01278                  RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
01279                  const Tuple<global_ordinal_type, 3>& dims,
01280                  const bool tolerant=false,
01281                  const bool debug=false)
01282       {
01283         if (pComm->getRank() == 0) {
01284           typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> raw_adder_type;
01285           RCP<raw_adder_type> pRaw =
01286             rcp (new raw_adder_type (dims[0], dims[1], dims[2], tolerant, debug));
01287           return rcp (new adder_type (pRaw, pBanner->symmType()));
01288         }
01289         else {
01290           return null;
01291         }
01292       }
01293 
01294     public:
01295 
01320       static sparse_matrix_ptr
01321       readSparseFile (const std::string& filename,
01322                       const RCP<const Comm<int> >& pComm,
01323                       const RCP<node_type>& pNode,
01324                       const bool callFillComplete=true,
01325                       const bool tolerant=false,
01326                       const bool debug=false)
01327       {
01328         const int myRank = Teuchos::rank (*pComm);
01329         std::ifstream in;
01330 
01331         // Only open the file on Rank 0.
01332         if (myRank == 0)
01333           in.open (filename.c_str());
01334         return readSparse (in, pComm, pNode, callFillComplete, tolerant, debug);
01335         // We can rely on the destructor of the input stream to close
01336         // the file on scope exit, even if readSparse() throws an
01337         // exception.
01338       }
01339 
01369       static sparse_matrix_ptr
01370       readSparseFile (const std::string& filename,
01371                       const RCP<const Comm<int> >& pComm,
01372                       const RCP<node_type>& pNode,
01373                       const RCP<Teuchos::ParameterList>& constructorParams,
01374                       const RCP<Teuchos::ParameterList>& fillCompleteParams,
01375                       const bool tolerant=false,
01376                       const bool debug=false)
01377       {
01378         std::ifstream in;
01379         if (pComm->getRank () == 0) { // only open on Process 0
01380           in.open (filename.c_str());
01381         }
01382         return readSparse (in, pComm, pNode, constructorParams,
01383                            fillCompleteParams, tolerant, debug);
01384       }
01385 
01423       static sparse_matrix_ptr
01424       readSparseFile (const std::string& filename,
01425                       const RCP<const map_type>& rowMap,
01426                       RCP<const map_type>& colMap,
01427                       const RCP<const map_type>& domainMap,
01428                       const RCP<const map_type>& rangeMap,
01429                       const bool callFillComplete=true,
01430                       const bool tolerant=false,
01431                       const bool debug=false)
01432       {
01433         using Teuchos::broadcast;
01434         using Teuchos::outArg;
01435         TEUCHOS_TEST_FOR_EXCEPTION(
01436           rowMap.is_null (), std::invalid_argument,
01437           "Row Map must be nonnull.");
01438 
01439         RCP<const Comm<int> > comm = rowMap->getComm ();
01440         RCP<node_type> node = rowMap->getNode ();
01441         const int myRank = comm->getRank ();
01442 
01443         // Only open the file on Process 0.  Test carefully to make
01444         // sure that the file opened successfully (and broadcast that
01445         // result to all processes to prevent a hang on exception
01446         // throw), since it's a common mistake to misspell a filename.
01447         std::ifstream in;
01448         int opened = 0;
01449         if (myRank == 0) {
01450           try {
01451             in.open (filename.c_str ());
01452             opened = 1;
01453           }
01454           catch (...) {
01455             opened = 0;
01456           }
01457         }
01458         broadcast<int, int> (*comm, 0, outArg (opened));
01459         TEUCHOS_TEST_FOR_EXCEPTION(
01460           opened == 0, std::runtime_error,
01461           "readSparseFile: Failed to open file \"" << filename << "\" on "
01462           "Process 0.");
01463         return readSparse (in, rowMap, colMap, domainMap, rangeMap,
01464                            callFillComplete, tolerant, debug);
01465       }
01466 
01493       static sparse_matrix_ptr
01494       readSparse (std::istream& in,
01495                   const RCP<const Comm<int> >& pComm,
01496                   const RCP<node_type>& pNode,
01497                   const bool callFillComplete=true,
01498                   const bool tolerant=false,
01499                   const bool debug=false)
01500       {
01501         using Teuchos::MatrixMarket::Banner;
01502         using Teuchos::broadcast;
01503         using Teuchos::ptr;
01504         using Teuchos::reduceAll;
01505         using std::cerr;
01506         using std::endl;
01507         typedef Teuchos::ScalarTraits<scalar_type> STS;
01508 
01509         const bool extraDebug = false;
01510         const int myRank = pComm->getRank ();
01511         const int rootRank = 0;
01512 
01513         // Current line number in the input stream.  Various calls
01514         // will modify this depending on the number of lines that are
01515         // read from the input stream.  Only Rank 0 modifies this.
01516         size_t lineNumber = 1;
01517 
01518         if (debug && myRank == rootRank) {
01519           cerr << "Matrix Market reader: readSparse:" << endl
01520                << "-- Reading banner line" << endl;
01521         }
01522 
01523         // The "Banner" tells you whether the input stream represents
01524         // a sparse matrix, the symmetry type of the matrix, and the
01525         // type of the data it contains.
01526         //
01527         // pBanner will only be nonnull on MPI Rank 0.  It will be
01528         // null on all other MPI processes.
01529         RCP<const Banner> pBanner;
01530         {
01531           // We read and validate the Banner on Proc 0, but broadcast
01532           // the validation result to all processes.
01533           // Teuchos::broadcast doesn't currently work with bool, so
01534           // we use int (true -> 1, false -> 0).
01535           int bannerIsCorrect = 1;
01536           std::ostringstream errMsg;
01537 
01538           if (myRank == rootRank) {
01539             // Read the Banner line from the input stream.
01540             try {
01541               pBanner = readBanner (in, lineNumber, tolerant, debug);
01542             }
01543             catch (std::exception& e) {
01544               errMsg << "Attempt to read the Matrix Market file's Banner line "
01545                 "threw an exception: " << e.what();
01546               bannerIsCorrect = 0;
01547             }
01548 
01549             if (bannerIsCorrect) {
01550               // Validate the Banner for the case of a sparse matrix.
01551               // We validate on Proc 0, since it reads the Banner.
01552 
01553               // In intolerant mode, the matrix type must be "coordinate".
01554               if (! tolerant && pBanner->matrixType() != "coordinate") {
01555                 bannerIsCorrect = 0;
01556                 errMsg << "The Matrix Market input file must contain a "
01557                   "\"coordinate\"-format sparse matrix in order to create a "
01558                   "Tpetra::CrsMatrix object from it, but the file's matrix "
01559                   "type is \"" << pBanner->matrixType() << "\" instead.";
01560               }
01561               // In tolerant mode, we allow the matrix type to be
01562               // anything other than "array" (which would mean that
01563               // the file contains a dense matrix).
01564               if (tolerant && pBanner->matrixType() == "array") {
01565                 bannerIsCorrect = 0;
01566                 errMsg << "Matrix Market file must contain a \"coordinate\"-"
01567                   "format sparse matrix in order to create a Tpetra::CrsMatrix "
01568                   "object from it, but the file's matrix type is \"array\" "
01569                   "instead.  That probably means the file contains dense matrix "
01570                   "data.";
01571               }
01572             }
01573           } // Proc 0: Done reading the Banner, hopefully successfully.
01574 
01575           // Broadcast from Proc 0 whether the Banner was read correctly.
01576           broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
01577 
01578           // If the Banner is invalid, all processes throw an
01579           // exception.  Only Proc 0 gets the exception message, but
01580           // that's OK, since the main point is to "stop the world"
01581           // (rather than throw an exception on one process and leave
01582           // the others hanging).
01583           TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
01584             std::invalid_argument, errMsg.str ());
01585         } // Done reading the Banner line and broadcasting success.
01586         if (debug && myRank == rootRank) {
01587           cerr << "-- Reading dimensions line" << endl;
01588         }
01589 
01590         // Read the matrix dimensions from the Matrix Market metadata.
01591         // dims = (numRows, numCols, numEntries).  Proc 0 does the
01592         // reading, but it broadcasts the results to all MPI
01593         // processes.  Thus, readCoordDims() is a collective
01594         // operation.  It does a collective check for correctness too.
01595         Tuple<global_ordinal_type, 3> dims =
01596           readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
01597 
01598         if (debug && myRank == rootRank) {
01599           cerr << "-- Making Adder for collecting matrix data" << endl;
01600         }
01601 
01602         // "Adder" object for collecting all the sparse matrix entries
01603         // from the input stream.  This is only nonnull on Proc 0.
01604         RCP<adder_type> pAdder =
01605           makeAdder (pComm, pBanner, dims, tolerant, debug);
01606 
01607         if (debug && myRank == rootRank) {
01608           cerr << "-- Reading matrix data" << endl;
01609         }
01610         //
01611         // Read the matrix entries from the input stream on Proc 0.
01612         //
01613         {
01614           // We use readSuccess to broadcast the results of the read
01615           // (succeeded or not) to all MPI processes.  Since
01616           // Teuchos::broadcast doesn't currently know how to send
01617           // bools, we convert to int (true -> 1, false -> 0).
01618           int readSuccess = 1;
01619           std::ostringstream errMsg; // Exception message (only valid on Proc 0)
01620           if (myRank == rootRank) {
01621             try {
01622               // Reader for "coordinate" format sparse matrix data.
01623               typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
01624                 global_ordinal_type, scalar_type, STS::isComplex> reader_type;
01625               reader_type reader (pAdder);
01626 
01627               // Read the sparse matrix entries.
01628               std::pair<bool, std::vector<size_t> > results =
01629                 reader.read (in, lineNumber, tolerant, debug);
01630               readSuccess = results.first ? 1 : 0;
01631             }
01632             catch (std::exception& e) {
01633               readSuccess = 0;
01634               errMsg << e.what();
01635             }
01636           }
01637           broadcast (*pComm, rootRank, ptr (&readSuccess));
01638 
01639           // It would be nice to add a "verbose" flag, so that in
01640           // tolerant mode, we could log any bad line number(s) on
01641           // Proc 0.  For now, we just throw if the read fails to
01642           // succeed.
01643           //
01644           // Question: If we're in tolerant mode, and if the read did
01645           // not succeed, should we attempt to call fillComplete()?
01646           TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
01647             "Failed to read the Matrix Market sparse matrix file: "
01648             << errMsg.str());
01649         } // Done reading the matrix entries (stored on Proc 0 for now)
01650 
01651         if (debug && myRank == rootRank) {
01652           cerr << "-- Successfully read the Matrix Market data" << endl;
01653         }
01654 
01655         // In tolerant mode, we need to rebroadcast the matrix
01656         // dimensions, since they may be different after reading the
01657         // actual matrix data.  We only need to broadcast the number
01658         // of rows and columns.  Only Rank 0 needs to know the actual
01659         // global number of entries, since (a) we need to merge
01660         // duplicates on Rank 0 first anyway, and (b) when we
01661         // distribute the entries, each rank other than Rank 0 will
01662         // only need to know how many entries it owns, not the total
01663         // number of entries.
01664         if (tolerant) {
01665           if (debug && myRank == rootRank) {
01666             cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
01667                  << endl
01668                  << "----- Dimensions before: "
01669                  << dims[0] << " x " << dims[1]
01670                  << endl;
01671           }
01672           // Packed coordinate matrix dimensions (numRows, numCols).
01673           Tuple<global_ordinal_type, 2> updatedDims;
01674           if (myRank == rootRank) {
01675             // If one or more bottom rows of the matrix contain no
01676             // entries, then the Adder will report that the number
01677             // of rows is less than that specified in the
01678             // metadata.  We allow this case, and favor the
01679             // metadata so that the zero row(s) will be included.
01680             updatedDims[0] =
01681               std::max (dims[0], pAdder->getAdder()->numRows());
01682             updatedDims[1] = pAdder->getAdder()->numCols();
01683           }
01684           broadcast (*pComm, rootRank, updatedDims);
01685           dims[0] = updatedDims[0];
01686           dims[1] = updatedDims[1];
01687           if (debug && myRank == rootRank) {
01688             cerr << "----- Dimensions after: " << dims[0] << " x "
01689                  << dims[1] << endl;
01690           }
01691         }
01692         else {
01693           // In strict mode, we require that the matrix's metadata and
01694           // its actual data agree, at least somewhat.  In particular,
01695           // the number of rows must agree, since otherwise we cannot
01696           // distribute the matrix correctly.
01697 
01698           // Teuchos::broadcast() doesn't know how to broadcast bools,
01699           // so we use an int with the standard 1 == true, 0 == false
01700           // encoding.
01701           int dimsMatch = 1;
01702           if (myRank == rootRank) {
01703             // If one or more bottom rows of the matrix contain no
01704             // entries, then the Adder will report that the number of
01705             // rows is less than that specified in the metadata.  We
01706             // allow this case, and favor the metadata, but do not
01707             // allow the Adder to think there are more rows in the
01708             // matrix than the metadata says.
01709             if (dims[0] < pAdder->getAdder ()->numRows ()) {
01710               dimsMatch = 0;
01711             }
01712           }
01713           broadcast (*pComm, 0, ptr (&dimsMatch));
01714           if (dimsMatch == 0) {
01715             // We're in an error state anyway, so we might as well
01716             // work a little harder to print an informative error
01717             // message.
01718             //
01719             // Broadcast the Adder's idea of the matrix dimensions
01720             // from Proc 0 to all processes.
01721             Tuple<global_ordinal_type, 2> addersDims;
01722             if (myRank == rootRank) {
01723               addersDims[0] = pAdder->getAdder()->numRows();
01724               addersDims[1] = pAdder->getAdder()->numCols();
01725             }
01726             broadcast (*pComm, 0, addersDims);
01727             TEUCHOS_TEST_FOR_EXCEPTION(
01728               dimsMatch == 0, std::runtime_error,
01729               "The matrix metadata says that the matrix is " << dims[0] << " x "
01730               << dims[1] << ", but the actual data says that the matrix is "
01731               << addersDims[0] << " x " << addersDims[1] << ".  That means the "
01732               "data includes more rows than reported in the metadata.  This "
01733               "is not allowed when parsing in strict mode.  Parse the matrix in "
01734               "tolerant mode to ignore the metadata when it disagrees with the "
01735               "data.");
01736           }
01737         } // Matrix dimensions (# rows, # cols, # entries) agree.
01738 
01739         if (debug && myRank == rootRank) {
01740           cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
01741         }
01742 
01743         // Now that we've read in all the matrix entries from the
01744         // input stream into the adder on Proc 0, post-process them
01745         // into CSR format (still on Proc 0).  This will facilitate
01746         // distributing them to all the processors.
01747         //
01748         // These arrays represent the global matrix data as a CSR
01749         // matrix (with numEntriesPerRow as redundant but convenient
01750         // metadata, since it's computable from rowPtr and vice
01751         // versa).  They are valid only on Proc 0.
01752         ArrayRCP<size_t> numEntriesPerRow;
01753         ArrayRCP<size_t> rowPtr;
01754         ArrayRCP<global_ordinal_type> colInd;
01755         ArrayRCP<scalar_type> values;
01756 
01757         // Proc 0 first merges duplicate entries, and then converts
01758         // the coordinate-format matrix data to CSR.
01759         {
01760           int mergeAndConvertSucceeded = 1;
01761           std::ostringstream errMsg;
01762 
01763           if (myRank == rootRank) {
01764             try {
01765               typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
01766                 global_ordinal_type> element_type;
01767 
01768               // Number of rows in the matrix.  If we are in tolerant
01769               // mode, we've already synchronized dims with the actual
01770               // matrix data.  If in strict mode, we should use dims
01771               // (as read from the file's metadata) rather than the
01772               // matrix data to determine the dimensions.  (The matrix
01773               // data will claim fewer rows than the metadata, if one
01774               // or more rows have no entries stored in the file.)
01775               const size_type numRows = dims[0];
01776 
01777               // Additively merge duplicate matrix entries.
01778               pAdder->getAdder()->merge ();
01779 
01780               // Get a temporary const view of the merged matrix entries.
01781               const std::vector<element_type>& entries =
01782                 pAdder->getAdder()->getEntries();
01783 
01784               // Number of matrix entries (after merging).
01785               const size_t numEntries = (size_t)entries.size();
01786 
01787               if (debug) {
01788                 cerr << "----- Proc 0: Matrix has numRows=" << numRows
01789                      << " rows and numEntries=" << numEntries
01790                      << " entries." << endl;
01791               }
01792 
01793               // Make space for the CSR matrix data.  Converting to
01794               // CSR is easier if we fill numEntriesPerRow with zeros
01795               // at first.
01796               numEntriesPerRow = arcp<size_t> (numRows);
01797               std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
01798               rowPtr = arcp<size_t> (numRows+1);
01799               std::fill (rowPtr.begin(), rowPtr.end(), 0);
01800               colInd = arcp<global_ordinal_type> (numEntries);
01801               values = arcp<scalar_type> (numEntries);
01802 
01803               // Convert from array-of-structs coordinate format to CSR
01804               // (compressed sparse row) format.
01805               global_ordinal_type prvRow = 0;
01806               size_t curPos = 0;
01807               rowPtr[0] = 0;
01808               for (curPos = 0; curPos < numEntries; ++curPos) {
01809                 const element_type& curEntry = entries[curPos];
01810                 const global_ordinal_type curRow = curEntry.rowIndex();
01811                 TEUCHOS_TEST_FOR_EXCEPTION(
01812                   curRow < prvRow, std::logic_error,
01813                   "Row indices are out of order, even though they are supposed "
01814                   "to be sorted.  curRow = " << curRow << ", prvRow = "
01815                   << prvRow << ", at curPos = " << curPos << ".  Please report "
01816                   "this bug to the Tpetra developers.");
01817                 if (curRow > prvRow) {
01818                   for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
01819                     rowPtr[r] = curPos;
01820                   }
01821                   prvRow = curRow;
01822                 }
01823                 numEntriesPerRow[curRow]++;
01824                 colInd[curPos] = curEntry.colIndex();
01825                 values[curPos] = curEntry.value();
01826               }
01827               // rowPtr has one more entry than numEntriesPerRow.  The
01828               // last entry of rowPtr is the number of entries in
01829               // colInd and values.
01830               rowPtr[numRows] = numEntries;
01831             } // Finished conversion to CSR format
01832             catch (std::exception& e) {
01833               mergeAndConvertSucceeded = 0;
01834               errMsg << "Failed to merge sparse matrix entries and convert to "
01835                 "CSR format: " << e.what();
01836             }
01837 
01838             if (debug && mergeAndConvertSucceeded) {
01839               // Number of rows in the matrix.
01840               const size_type numRows = dims[0];
01841               const size_type maxToDisplay = 100;
01842 
01843               cerr << "----- Proc 0: numEntriesPerRow[0.."
01844                    << (numEntriesPerRow.size()-1) << "] ";
01845               if (numRows > maxToDisplay) {
01846                 cerr << "(only showing first and last few entries) ";
01847               }
01848               cerr << "= [";
01849               if (numRows > 0) {
01850                 if (numRows > maxToDisplay) {
01851                   for (size_type k = 0; k < 2; ++k) {
01852                     cerr << numEntriesPerRow[k] << " ";
01853                   }
01854                   cerr << "... ";
01855                   for (size_type k = numRows-2; k < numRows-1; ++k) {
01856                     cerr << numEntriesPerRow[k] << " ";
01857                   }
01858                 }
01859                 else {
01860                   for (size_type k = 0; k < numRows-1; ++k) {
01861                     cerr << numEntriesPerRow[k] << " ";
01862                   }
01863                 }
01864                 cerr << numEntriesPerRow[numRows-1];
01865               } // numRows > 0
01866               cerr << "]" << endl;
01867 
01868               cerr << "----- Proc 0: rowPtr ";
01869               if (numRows > maxToDisplay) {
01870                 cerr << "(only showing first and last few entries) ";
01871               }
01872               cerr << "= [";
01873               if (numRows > maxToDisplay) {
01874                 for (size_type k = 0; k < 2; ++k) {
01875                   cerr << rowPtr[k] << " ";
01876                 }
01877                 cerr << "... ";
01878                 for (size_type k = numRows-2; k < numRows; ++k) {
01879                   cerr << rowPtr[k] << " ";
01880                 }
01881               }
01882               else {
01883                 for (size_type k = 0; k < numRows; ++k) {
01884                   cerr << rowPtr[k] << " ";
01885                 }
01886               }
01887               cerr << rowPtr[numRows] << "]" << endl;
01888             }
01889           } // if myRank == rootRank
01890         } // Done converting sparse matrix data to CSR format
01891 
01892         // Now we're done with the Adder, so we can release the
01893         // reference ("free" it) to save space.  This only actually
01894         // does anything on Rank 0, since pAdder is null on all the
01895         // other MPI processes.
01896         pAdder = null;
01897 
01898         if (debug && myRank == rootRank) {
01899           cerr << "-- Making range, domain, and row maps" << endl;
01900         }
01901 
01902         // Make the maps that describe the matrix's range and domain,
01903         // and the distribution of its rows.  Creating a Map is a
01904         // collective operation, so we don't have to do a broadcast of
01905         // a success Boolean.
01906         map_ptr pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
01907         map_ptr pDomainMap = makeDomainMap (pRangeMap, dims[0], dims[1]);
01908         map_ptr pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
01909 
01910         if (debug && myRank == rootRank) {
01911           cerr << "-- Distributing the matrix data" << endl;
01912         }
01913 
01914         // Distribute the matrix data.  Each processor has to add the
01915         // rows that it owns.  If you try to make Proc 0 call
01916         // insertGlobalValues() for _all_ the rows, not just those it
01917         // owns, then fillComplete() will compute the number of
01918         // columns incorrectly.  That's why Proc 0 has to distribute
01919         // the matrix data and why we make all the processors (not
01920         // just Proc 0) call insertGlobalValues() on their own data.
01921         //
01922         // These arrays represent each processor's part of the matrix
01923         // data, in "CSR" format (sort of, since the row indices might
01924         // not be contiguous).
01925         ArrayRCP<size_t> myNumEntriesPerRow;
01926         ArrayRCP<size_t> myRowPtr;
01927         ArrayRCP<global_ordinal_type> myColInd;
01928         ArrayRCP<scalar_type> myValues;
01929         // Distribute the matrix data.  This is a collective operation.
01930         distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
01931                     numEntriesPerRow, rowPtr, colInd, values, debug);
01932 
01933         if (debug && myRank == rootRank) {
01934           cerr << "-- Inserting matrix entries on each processor";
01935           if (callFillComplete) {
01936             cerr << " and calling fillComplete()";
01937           }
01938           cerr << endl;
01939         }
01940         // Each processor inserts its part of the matrix data, and
01941         // then they all call fillComplete().  This method invalidates
01942         // the my* distributed matrix data before calling
01943         // fillComplete(), in order to save space.  In general, we
01944         // never store more than two copies of the matrix's entries in
01945         // memory at once, which is no worse than what Tpetra
01946         // promises.
01947         sparse_matrix_ptr pMatrix =
01948           makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
01949                       pRowMap, pRangeMap, pDomainMap, callFillComplete);
01950         // Only use a reduce-all in debug mode to check if pMatrix is
01951         // null.  Otherwise, just throw an exception.  We never expect
01952         // a null pointer here, so we can save a communication.
01953         if (debug) {
01954           int localIsNull = pMatrix.is_null () ? 1 : 0;
01955           int globalIsNull = 0;
01956           reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
01957           TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
01958             "Reader::makeMatrix() returned a null pointer on at least one "
01959             "process.  Please report this bug to the Tpetra developers.");
01960         }
01961         else {
01962           TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
01963             "Reader::makeMatrix() returned a null pointer.  "
01964             "Please report this bug to the Tpetra developers.");
01965         }
01966 
01967         // We can't get the dimensions of the matrix until after
01968         // fillComplete() is called.  Thus, we can't do the sanity
01969         // check (dimensions read from the Matrix Market data,
01970         // vs. dimensions reported by the CrsMatrix) unless the user
01971         // asked makeMatrix() to call fillComplete().
01972         //
01973         // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
01974         // what one might think it does, so you have to ask the range
01975         // resp. domain map for the number of rows resp. columns.
01976         if (callFillComplete) {
01977           const int numProcs = pComm->getSize ();
01978 
01979           if (extraDebug && debug) {
01980             const global_size_t globalNumRows =
01981               pRangeMap->getGlobalNumElements();
01982             const global_size_t globalNumCols =
01983               pDomainMap->getGlobalNumElements();
01984             if (myRank == rootRank) {
01985               cerr << "-- Matrix is "
01986                    << globalNumRows << " x " << globalNumCols
01987                    << " with " << pMatrix->getGlobalNumEntries()
01988                    << " entries, and index base "
01989                    << pMatrix->getIndexBase() << "." << endl;
01990             }
01991             pComm->barrier ();
01992             for (int p = 0; p < numProcs; ++p) {
01993               if (myRank == p) {
01994                 cerr << "-- Proc " << p << " owns "
01995                      << pMatrix->getNodeNumCols() << " columns, and "
01996                      << pMatrix->getNodeNumEntries() << " entries." << endl;
01997               }
01998               pComm->barrier ();
01999             }
02000           } // if (extraDebug && debug)
02001         } // if (callFillComplete)
02002 
02003         if (debug && myRank == rootRank) {
02004           cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
02005                << endl;
02006         }
02007         return pMatrix;
02008       }
02009 
02010 
02040       static sparse_matrix_ptr
02041       readSparse (std::istream& in,
02042                   const RCP<const Comm<int> >& pComm,
02043                   const RCP<node_type>& pNode,
02044                   const RCP<Teuchos::ParameterList>& constructorParams,
02045                   const RCP<Teuchos::ParameterList>& fillCompleteParams,
02046                   const bool tolerant=false,
02047                   const bool debug=false)
02048       {
02049         using Teuchos::MatrixMarket::Banner;
02050         using Teuchos::broadcast;
02051         using Teuchos::ptr;
02052         using Teuchos::reduceAll;
02053         using std::cerr;
02054         using std::endl;
02055         typedef Teuchos::ScalarTraits<scalar_type> STS;
02056 
02057         const bool extraDebug = false;
02058         const int myRank = pComm->getRank ();
02059         const int rootRank = 0;
02060 
02061         // Current line number in the input stream.  Various calls
02062         // will modify this depending on the number of lines that are
02063         // read from the input stream.  Only Rank 0 modifies this.
02064         size_t lineNumber = 1;
02065 
02066         if (debug && myRank == rootRank) {
02067           cerr << "Matrix Market reader: readSparse:" << endl
02068                << "-- Reading banner line" << endl;
02069         }
02070 
02071         // The "Banner" tells you whether the input stream represents
02072         // a sparse matrix, the symmetry type of the matrix, and the
02073         // type of the data it contains.
02074         //
02075         // pBanner will only be nonnull on MPI Rank 0.  It will be
02076         // null on all other MPI processes.
02077         RCP<const Banner> pBanner;
02078         {
02079           // We read and validate the Banner on Proc 0, but broadcast
02080           // the validation result to all processes.
02081           // Teuchos::broadcast doesn't currently work with bool, so
02082           // we use int (true -> 1, false -> 0).
02083           int bannerIsCorrect = 1;
02084           std::ostringstream errMsg;
02085 
02086           if (myRank == rootRank) {
02087             // Read the Banner line from the input stream.
02088             try {
02089               pBanner = readBanner (in, lineNumber, tolerant, debug);
02090             }
02091             catch (std::exception& e) {
02092               errMsg << "Attempt to read the Matrix Market file's Banner line "
02093                 "threw an exception: " << e.what();
02094               bannerIsCorrect = 0;
02095             }
02096 
02097             if (bannerIsCorrect) {
02098               // Validate the Banner for the case of a sparse matrix.
02099               // We validate on Proc 0, since it reads the Banner.
02100 
02101               // In intolerant mode, the matrix type must be "coordinate".
02102               if (! tolerant && pBanner->matrixType() != "coordinate") {
02103                 bannerIsCorrect = 0;
02104                 errMsg << "The Matrix Market input file must contain a "
02105                   "\"coordinate\"-format sparse matrix in order to create a "
02106                   "Tpetra::CrsMatrix object from it, but the file's matrix "
02107                   "type is \"" << pBanner->matrixType() << "\" instead.";
02108               }
02109               // In tolerant mode, we allow the matrix type to be
02110               // anything other than "array" (which would mean that
02111               // the file contains a dense matrix).
02112               if (tolerant && pBanner->matrixType() == "array") {
02113                 bannerIsCorrect = 0;
02114                 errMsg << "Matrix Market file must contain a \"coordinate\"-"
02115                   "format sparse matrix in order to create a Tpetra::CrsMatrix "
02116                   "object from it, but the file's matrix type is \"array\" "
02117                   "instead.  That probably means the file contains dense matrix "
02118                   "data.";
02119               }
02120             }
02121           } // Proc 0: Done reading the Banner, hopefully successfully.
02122 
02123           // Broadcast from Proc 0 whether the Banner was read correctly.
02124           broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
02125 
02126           // If the Banner is invalid, all processes throw an
02127           // exception.  Only Proc 0 gets the exception message, but
02128           // that's OK, since the main point is to "stop the world"
02129           // (rather than throw an exception on one process and leave
02130           // the others hanging).
02131           TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
02132             std::invalid_argument, errMsg.str ());
02133         } // Done reading the Banner line and broadcasting success.
02134         if (debug && myRank == rootRank) {
02135           cerr << "-- Reading dimensions line" << endl;
02136         }
02137 
02138         // Read the matrix dimensions from the Matrix Market metadata.
02139         // dims = (numRows, numCols, numEntries).  Proc 0 does the
02140         // reading, but it broadcasts the results to all MPI
02141         // processes.  Thus, readCoordDims() is a collective
02142         // operation.  It does a collective check for correctness too.
02143         Tuple<global_ordinal_type, 3> dims =
02144           readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
02145 
02146         if (debug && myRank == rootRank) {
02147           cerr << "-- Making Adder for collecting matrix data" << endl;
02148         }
02149 
02150         // "Adder" object for collecting all the sparse matrix entries
02151         // from the input stream.  This is only nonnull on Proc 0.
02152         RCP<adder_type> pAdder =
02153           makeAdder (pComm, pBanner, dims, tolerant, debug);
02154 
02155         if (debug && myRank == rootRank) {
02156           cerr << "-- Reading matrix data" << endl;
02157         }
02158         //
02159         // Read the matrix entries from the input stream on Proc 0.
02160         //
02161         {
02162           // We use readSuccess to broadcast the results of the read
02163           // (succeeded or not) to all MPI processes.  Since
02164           // Teuchos::broadcast doesn't currently know how to send
02165           // bools, we convert to int (true -> 1, false -> 0).
02166           int readSuccess = 1;
02167           std::ostringstream errMsg; // Exception message (only valid on Proc 0)
02168           if (myRank == rootRank) {
02169             try {
02170               // Reader for "coordinate" format sparse matrix data.
02171               typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
02172                 global_ordinal_type, scalar_type, STS::isComplex> reader_type;
02173               reader_type reader (pAdder);
02174 
02175               // Read the sparse matrix entries.
02176               std::pair<bool, std::vector<size_t> > results =
02177                 reader.read (in, lineNumber, tolerant, debug);
02178               readSuccess = results.first ? 1 : 0;
02179             }
02180             catch (std::exception& e) {
02181               readSuccess = 0;
02182               errMsg << e.what();
02183             }
02184           }
02185           broadcast (*pComm, rootRank, ptr (&readSuccess));
02186 
02187           // It would be nice to add a "verbose" flag, so that in
02188           // tolerant mode, we could log any bad line number(s) on
02189           // Proc 0.  For now, we just throw if the read fails to
02190           // succeed.
02191           //
02192           // Question: If we're in tolerant mode, and if the read did
02193           // not succeed, should we attempt to call fillComplete()?
02194           TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
02195             "Failed to read the Matrix Market sparse matrix file: "
02196             << errMsg.str());
02197         } // Done reading the matrix entries (stored on Proc 0 for now)
02198 
02199         if (debug && myRank == rootRank) {
02200           cerr << "-- Successfully read the Matrix Market data" << endl;
02201         }
02202 
02203         // In tolerant mode, we need to rebroadcast the matrix
02204         // dimensions, since they may be different after reading the
02205         // actual matrix data.  We only need to broadcast the number
02206         // of rows and columns.  Only Rank 0 needs to know the actual
02207         // global number of entries, since (a) we need to merge
02208         // duplicates on Rank 0 first anyway, and (b) when we
02209         // distribute the entries, each rank other than Rank 0 will
02210         // only need to know how many entries it owns, not the total
02211         // number of entries.
02212         if (tolerant) {
02213           if (debug && myRank == rootRank) {
02214             cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
02215                  << endl
02216                  << "----- Dimensions before: "
02217                  << dims[0] << " x " << dims[1]
02218                  << endl;
02219           }
02220           // Packed coordinate matrix dimensions (numRows, numCols).
02221           Tuple<global_ordinal_type, 2> updatedDims;
02222           if (myRank == rootRank) {
02223             // If one or more bottom rows of the matrix contain no
02224             // entries, then the Adder will report that the number
02225             // of rows is less than that specified in the
02226             // metadata.  We allow this case, and favor the
02227             // metadata so that the zero row(s) will be included.
02228             updatedDims[0] =
02229               std::max (dims[0], pAdder->getAdder()->numRows());
02230             updatedDims[1] = pAdder->getAdder()->numCols();
02231           }
02232           broadcast (*pComm, rootRank, updatedDims);
02233           dims[0] = updatedDims[0];
02234           dims[1] = updatedDims[1];
02235           if (debug && myRank == rootRank) {
02236             cerr << "----- Dimensions after: " << dims[0] << " x "
02237                  << dims[1] << endl;
02238           }
02239         }
02240         else {
02241           // In strict mode, we require that the matrix's metadata and
02242           // its actual data agree, at least somewhat.  In particular,
02243           // the number of rows must agree, since otherwise we cannot
02244           // distribute the matrix correctly.
02245 
02246           // Teuchos::broadcast() doesn't know how to broadcast bools,
02247           // so we use an int with the standard 1 == true, 0 == false
02248           // encoding.
02249           int dimsMatch = 1;
02250           if (myRank == rootRank) {
02251             // If one or more bottom rows of the matrix contain no
02252             // entries, then the Adder will report that the number of
02253             // rows is less than that specified in the metadata.  We
02254             // allow this case, and favor the metadata, but do not
02255             // allow the Adder to think there are more rows in the
02256             // matrix than the metadata says.
02257             if (dims[0] < pAdder->getAdder ()->numRows ()) {
02258               dimsMatch = 0;
02259             }
02260           }
02261           broadcast (*pComm, 0, ptr (&dimsMatch));
02262           if (dimsMatch == 0) {
02263             // We're in an error state anyway, so we might as well
02264             // work a little harder to print an informative error
02265             // message.
02266             //
02267             // Broadcast the Adder's idea of the matrix dimensions
02268             // from Proc 0 to all processes.
02269             Tuple<global_ordinal_type, 2> addersDims;
02270             if (myRank == rootRank) {
02271               addersDims[0] = pAdder->getAdder()->numRows();
02272               addersDims[1] = pAdder->getAdder()->numCols();
02273             }
02274             broadcast (*pComm, 0, addersDims);
02275             TEUCHOS_TEST_FOR_EXCEPTION(
02276               dimsMatch == 0, std::runtime_error,
02277               "The matrix metadata says that the matrix is " << dims[0] << " x "
02278               << dims[1] << ", but the actual data says that the matrix is "
02279               << addersDims[0] << " x " << addersDims[1] << ".  That means the "
02280               "data includes more rows than reported in the metadata.  This "
02281               "is not allowed when parsing in strict mode.  Parse the matrix in "
02282               "tolerant mode to ignore the metadata when it disagrees with the "
02283               "data.");
02284           }
02285         } // Matrix dimensions (# rows, # cols, # entries) agree.
02286 
02287         if (debug && myRank == rootRank) {
02288           cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
02289         }
02290 
02291         // Now that we've read in all the matrix entries from the
02292         // input stream into the adder on Proc 0, post-process them
02293         // into CSR format (still on Proc 0).  This will facilitate
02294         // distributing them to all the processors.
02295         //
02296         // These arrays represent the global matrix data as a CSR
02297         // matrix (with numEntriesPerRow as redundant but convenient
02298         // metadata, since it's computable from rowPtr and vice
02299         // versa).  They are valid only on Proc 0.
02300         ArrayRCP<size_t> numEntriesPerRow;
02301         ArrayRCP<size_t> rowPtr;
02302         ArrayRCP<global_ordinal_type> colInd;
02303         ArrayRCP<scalar_type> values;
02304 
02305         // Proc 0 first merges duplicate entries, and then converts
02306         // the coordinate-format matrix data to CSR.
02307         {
02308           int mergeAndConvertSucceeded = 1;
02309           std::ostringstream errMsg;
02310 
02311           if (myRank == rootRank) {
02312             try {
02313               typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
02314                 global_ordinal_type> element_type;
02315 
02316               // Number of rows in the matrix.  If we are in tolerant
02317               // mode, we've already synchronized dims with the actual
02318               // matrix data.  If in strict mode, we should use dims
02319               // (as read from the file's metadata) rather than the
02320               // matrix data to determine the dimensions.  (The matrix
02321               // data will claim fewer rows than the metadata, if one
02322               // or more rows have no entries stored in the file.)
02323               const size_type numRows = dims[0];
02324 
02325               // Additively merge duplicate matrix entries.
02326               pAdder->getAdder()->merge ();
02327 
02328               // Get a temporary const view of the merged matrix entries.
02329               const std::vector<element_type>& entries =
02330                 pAdder->getAdder()->getEntries();
02331 
02332               // Number of matrix entries (after merging).
02333               const size_t numEntries = (size_t)entries.size();
02334 
02335               if (debug) {
02336                 cerr << "----- Proc 0: Matrix has numRows=" << numRows
02337                      << " rows and numEntries=" << numEntries
02338                      << " entries." << endl;
02339               }
02340 
02341               // Make space for the CSR matrix data.  Converting to
02342               // CSR is easier if we fill numEntriesPerRow with zeros
02343               // at first.
02344               numEntriesPerRow = arcp<size_t> (numRows);
02345               std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
02346               rowPtr = arcp<size_t> (numRows+1);
02347               std::fill (rowPtr.begin(), rowPtr.end(), 0);
02348               colInd = arcp<global_ordinal_type> (numEntries);
02349               values = arcp<scalar_type> (numEntries);
02350 
02351               // Convert from array-of-structs coordinate format to CSR
02352               // (compressed sparse row) format.
02353               global_ordinal_type prvRow = 0;
02354               size_t curPos = 0;
02355               rowPtr[0] = 0;
02356               for (curPos = 0; curPos < numEntries; ++curPos) {
02357                 const element_type& curEntry = entries[curPos];
02358                 const global_ordinal_type curRow = curEntry.rowIndex();
02359                 TEUCHOS_TEST_FOR_EXCEPTION(
02360                   curRow < prvRow, std::logic_error,
02361                   "Row indices are out of order, even though they are supposed "
02362                   "to be sorted.  curRow = " << curRow << ", prvRow = "
02363                   << prvRow << ", at curPos = " << curPos << ".  Please report "
02364                   "this bug to the Tpetra developers.");
02365                 if (curRow > prvRow) {
02366                   for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
02367                     rowPtr[r] = curPos;
02368                   }
02369                   prvRow = curRow;
02370                 }
02371                 numEntriesPerRow[curRow]++;
02372                 colInd[curPos] = curEntry.colIndex();
02373                 values[curPos] = curEntry.value();
02374               }
02375               // rowPtr has one more entry than numEntriesPerRow.  The
02376               // last entry of rowPtr is the number of entries in
02377               // colInd and values.
02378               rowPtr[numRows] = numEntries;
02379             } // Finished conversion to CSR format
02380             catch (std::exception& e) {
02381               mergeAndConvertSucceeded = 0;
02382               errMsg << "Failed to merge sparse matrix entries and convert to "
02383                 "CSR format: " << e.what();
02384             }
02385 
02386             if (debug && mergeAndConvertSucceeded) {
02387               // Number of rows in the matrix.
02388               const size_type numRows = dims[0];
02389               const size_type maxToDisplay = 100;
02390 
02391               cerr << "----- Proc 0: numEntriesPerRow[0.."
02392                    << (numEntriesPerRow.size()-1) << "] ";
02393               if (numRows > maxToDisplay) {
02394                 cerr << "(only showing first and last few entries) ";
02395               }
02396               cerr << "= [";
02397               if (numRows > 0) {
02398                 if (numRows > maxToDisplay) {
02399                   for (size_type k = 0; k < 2; ++k) {
02400                     cerr << numEntriesPerRow[k] << " ";
02401                   }
02402                   cerr << "... ";
02403                   for (size_type k = numRows-2; k < numRows-1; ++k) {
02404                     cerr << numEntriesPerRow[k] << " ";
02405                   }
02406                 }
02407                 else {
02408                   for (size_type k = 0; k < numRows-1; ++k) {
02409                     cerr << numEntriesPerRow[k] << " ";
02410                   }
02411                 }
02412                 cerr << numEntriesPerRow[numRows-1];
02413               } // numRows > 0
02414               cerr << "]" << endl;
02415 
02416               cerr << "----- Proc 0: rowPtr ";
02417               if (numRows > maxToDisplay) {
02418                 cerr << "(only showing first and last few entries) ";
02419               }
02420               cerr << "= [";
02421               if (numRows > maxToDisplay) {
02422                 for (size_type k = 0; k < 2; ++k) {
02423                   cerr << rowPtr[k] << " ";
02424                 }
02425                 cerr << "... ";
02426                 for (size_type k = numRows-2; k < numRows; ++k) {
02427                   cerr << rowPtr[k] << " ";
02428                 }
02429               }
02430               else {
02431                 for (size_type k = 0; k < numRows; ++k) {
02432                   cerr << rowPtr[k] << " ";
02433                 }
02434               }
02435               cerr << rowPtr[numRows] << "]" << endl;
02436             }
02437           } // if myRank == rootRank
02438         } // Done converting sparse matrix data to CSR format
02439 
02440         // Now we're done with the Adder, so we can release the
02441         // reference ("free" it) to save space.  This only actually
02442         // does anything on Rank 0, since pAdder is null on all the
02443         // other MPI processes.
02444         pAdder = null;
02445 
02446         if (debug && myRank == rootRank) {
02447           cerr << "-- Making range, domain, and row maps" << endl;
02448         }
02449 
02450         // Make the maps that describe the matrix's range and domain,
02451         // and the distribution of its rows.  Creating a Map is a
02452         // collective operation, so we don't have to do a broadcast of
02453         // a success Boolean.
02454         map_ptr pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
02455         map_ptr pDomainMap = makeDomainMap (pRangeMap, dims[0], dims[1]);
02456         map_ptr pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
02457 
02458         if (debug && myRank == rootRank) {
02459           cerr << "-- Distributing the matrix data" << endl;
02460         }
02461 
02462         // Distribute the matrix data.  Each processor has to add the
02463         // rows that it owns.  If you try to make Proc 0 call
02464         // insertGlobalValues() for _all_ the rows, not just those it
02465         // owns, then fillComplete() will compute the number of
02466         // columns incorrectly.  That's why Proc 0 has to distribute
02467         // the matrix data and why we make all the processors (not
02468         // just Proc 0) call insertGlobalValues() on their own data.
02469         //
02470         // These arrays represent each processor's part of the matrix
02471         // data, in "CSR" format (sort of, since the row indices might
02472         // not be contiguous).
02473         ArrayRCP<size_t> myNumEntriesPerRow;
02474         ArrayRCP<size_t> myRowPtr;
02475         ArrayRCP<global_ordinal_type> myColInd;
02476         ArrayRCP<scalar_type> myValues;
02477         // Distribute the matrix data.  This is a collective operation.
02478         distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
02479                     numEntriesPerRow, rowPtr, colInd, values, debug);
02480 
02481         if (debug && myRank == rootRank) {
02482           cerr << "-- Inserting matrix entries on each process "
02483             "and calling fillComplete()" << endl;
02484         }
02485         // Each processor inserts its part of the matrix data, and
02486         // then they all call fillComplete().  This method invalidates
02487         // the my* distributed matrix data before calling
02488         // fillComplete(), in order to save space.  In general, we
02489         // never store more than two copies of the matrix's entries in
02490         // memory at once, which is no worse than what Tpetra
02491         // promises.
02492         sparse_matrix_ptr pMatrix =
02493           makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
02494                       pRowMap, pRangeMap, pDomainMap, constructorParams,
02495                       fillCompleteParams);
02496         // Only use a reduce-all in debug mode to check if pMatrix is
02497         // null.  Otherwise, just throw an exception.  We never expect
02498         // a null pointer here, so we can save a communication.
02499         if (debug) {
02500           int localIsNull = pMatrix.is_null () ? 1 : 0;
02501           int globalIsNull = 0;
02502           reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
02503           TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
02504             "Reader::makeMatrix() returned a null pointer on at least one "
02505             "process.  Please report this bug to the Tpetra developers.");
02506         }
02507         else {
02508           TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
02509             "Reader::makeMatrix() returned a null pointer.  "
02510             "Please report this bug to the Tpetra developers.");
02511         }
02512 
02513         // Sanity check for dimensions (read from the Matrix Market
02514         // data, vs. dimensions reported by the CrsMatrix).
02515         //
02516         // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
02517         // what one might think it does, so you have to ask the range
02518         // resp. domain map for the number of rows resp. columns.
02519         if (extraDebug && debug) {
02520           const int numProcs = pComm->getSize ();
02521           const global_size_t globalNumRows =
02522             pRangeMap->getGlobalNumElements();
02523           const global_size_t globalNumCols =
02524             pDomainMap->getGlobalNumElements();
02525           if (myRank == rootRank) {
02526             cerr << "-- Matrix is "
02527                  << globalNumRows << " x " << globalNumCols
02528                  << " with " << pMatrix->getGlobalNumEntries()
02529                  << " entries, and index base "
02530                  << pMatrix->getIndexBase() << "." << endl;
02531           }
02532           pComm->barrier ();
02533           for (int p = 0; p < numProcs; ++p) {
02534             if (myRank == p) {
02535               cerr << "-- Proc " << p << " owns "
02536                    << pMatrix->getNodeNumCols() << " columns, and "
02537                    << pMatrix->getNodeNumEntries() << " entries." << endl;
02538             }
02539             pComm->barrier ();
02540           }
02541         } // if (extraDebug && debug)
02542 
02543         if (debug && myRank == rootRank) {
02544           cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
02545                << endl;
02546         }
02547         return pMatrix;
02548       }
02549 
02590       static sparse_matrix_ptr
02591       readSparse (std::istream& in,
02592                   const RCP<const map_type>& rowMap,
02593                   RCP<const map_type>& colMap,
02594                   const RCP<const map_type>& domainMap,
02595                   const RCP<const map_type>& rangeMap,
02596                   const bool callFillComplete=true,
02597                   const bool tolerant=false,
02598                   const bool debug=false)
02599       {
02600         using Teuchos::MatrixMarket::Banner;
02601         using Teuchos::as;
02602         using Teuchos::broadcast;
02603         using Teuchos::ptr;
02604         using Teuchos::reduceAll;
02605         using std::cerr;
02606         using std::endl;
02607         typedef Teuchos::ScalarTraits<scalar_type> STS;
02608 
02609         RCP<const Comm<int> > pComm = rowMap->getComm ();
02610         const int myRank = pComm->getRank ();
02611         const int rootRank = 0;
02612         const bool extraDebug = false;
02613 
02614         // Fast checks for invalid input.  We can't check other
02615         // attributes of the Maps until we've read in the matrix
02616         // dimensions.
02617         TEUCHOS_TEST_FOR_EXCEPTION(
02618           rowMap.is_null (), std::invalid_argument,
02619           "Row Map must be nonnull.");
02620         TEUCHOS_TEST_FOR_EXCEPTION(
02621           rangeMap.is_null (), std::invalid_argument,
02622           "Range Map must be nonnull.");
02623         TEUCHOS_TEST_FOR_EXCEPTION(
02624           domainMap.is_null (), std::invalid_argument,
02625           "Domain Map must be nonnull.");
02626         TEUCHOS_TEST_FOR_EXCEPTION(
02627           rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
02628           std::invalid_argument,
02629           "The specified row Map's communicator (rowMap->getComm())"
02630           "differs from the given separately supplied communicator pComm.");
02631         TEUCHOS_TEST_FOR_EXCEPTION(
02632           domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
02633           std::invalid_argument,
02634           "The specified domain Map's communicator (domainMap->getComm())"
02635           "differs from the given separately supplied communicator pComm.");
02636         TEUCHOS_TEST_FOR_EXCEPTION(
02637           rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
02638           std::invalid_argument,
02639           "The specified range Map's communicator (rangeMap->getComm())"
02640           "differs from the given separately supplied communicator pComm.");
02641 
02642         // Current line number in the input stream.  Various calls
02643         // will modify this depending on the number of lines that are
02644         // read from the input stream.  Only Rank 0 modifies this.
02645         size_t lineNumber = 1;
02646 
02647         if (debug && myRank == rootRank) {
02648           cerr << "Matrix Market reader: readSparse:" << endl
02649                << "-- Reading banner line" << endl;
02650         }
02651 
02652         // The "Banner" tells you whether the input stream represents
02653         // a sparse matrix, the symmetry type of the matrix, and the
02654         // type of the data it contains.
02655         //
02656         // pBanner will only be nonnull on MPI Rank 0.  It will be
02657         // null on all other MPI processes.
02658         RCP<const Banner> pBanner;
02659         {
02660           // We read and validate the Banner on Proc 0, but broadcast
02661           // the validation result to all processes.
02662           // Teuchos::broadcast doesn't currently work with bool, so
02663           // we use int (true -> 1, false -> 0).
02664           int bannerIsCorrect = 1;
02665           std::ostringstream errMsg;
02666 
02667           if (myRank == rootRank) {
02668             // Read the Banner line from the input stream.
02669             try {
02670               pBanner = readBanner (in, lineNumber, tolerant, debug);
02671             }
02672             catch (std::exception& e) {
02673               errMsg << "Attempt to read the Matrix Market file's Banner line "
02674                 "threw an exception: " << e.what();
02675               bannerIsCorrect = 0;
02676             }
02677 
02678             if (bannerIsCorrect) {
02679               // Validate the Banner for the case of a sparse matrix.
02680               // We validate on Proc 0, since it reads the Banner.
02681 
02682               // In intolerant mode, the matrix type must be "coordinate".
02683               if (! tolerant && pBanner->matrixType() != "coordinate") {
02684                 bannerIsCorrect = 0;
02685                 errMsg << "The Matrix Market input file must contain a "
02686                   "\"coordinate\"-format sparse matrix in order to create a "
02687                   "Tpetra::CrsMatrix object from it, but the file's matrix "
02688                   "type is \"" << pBanner->matrixType() << "\" instead.";
02689               }
02690               // In tolerant mode, we allow the matrix type to be
02691               // anything other than "array" (which would mean that
02692               // the file contains a dense matrix).
02693               if (tolerant && pBanner->matrixType() == "array") {
02694                 bannerIsCorrect = 0;
02695                 errMsg << "Matrix Market file must contain a \"coordinate\"-"
02696                   "format sparse matrix in order to create a Tpetra::CrsMatrix "
02697                   "object from it, but the file's matrix type is \"array\" "
02698                   "instead.  That probably means the file contains dense matrix "
02699                   "data.";
02700               }
02701             }
02702           } // Proc 0: Done reading the Banner, hopefully successfully.
02703 
02704           // Broadcast from Proc 0 whether the Banner was read correctly.
02705           broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
02706 
02707           // If the Banner is invalid, all processes throw an
02708           // exception.  Only Proc 0 gets the exception message, but
02709           // that's OK, since the main point is to "stop the world"
02710           // (rather than throw an exception on one process and leave
02711           // the others hanging).
02712           TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
02713             std::invalid_argument, errMsg.str ());
02714         } // Done reading the Banner line and broadcasting success.
02715         if (debug && myRank == rootRank) {
02716           cerr << "-- Reading dimensions line" << endl;
02717         }
02718 
02719         // Read the matrix dimensions from the Matrix Market metadata.
02720         // dims = (numRows, numCols, numEntries).  Proc 0 does the
02721         // reading, but it broadcasts the results to all MPI
02722         // processes.  Thus, readCoordDims() is a collective
02723         // operation.  It does a collective check for correctness too.
02724         Tuple<global_ordinal_type, 3> dims =
02725           readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
02726 
02727         if (debug && myRank == rootRank) {
02728           cerr << "-- Making Adder for collecting matrix data" << endl;
02729         }
02730 
02731         // "Adder" object for collecting all the sparse matrix entries
02732         // from the input stream.  This is only nonnull on Proc 0.
02733         // The Adder internally converts the one-based indices (native
02734         // Matrix Market format) into zero-based indices.
02735         RCP<adder_type> pAdder =
02736           makeAdder (pComm, pBanner, dims, tolerant, debug);
02737 
02738         if (debug && myRank == rootRank) {
02739           cerr << "-- Reading matrix data" << endl;
02740         }
02741         //
02742         // Read the matrix entries from the input stream on Proc 0.
02743         //
02744         {
02745           // We use readSuccess to broadcast the results of the read
02746           // (succeeded or not) to all MPI processes.  Since
02747           // Teuchos::broadcast doesn't currently know how to send
02748           // bools, we convert to int (true -> 1, false -> 0).
02749           int readSuccess = 1;
02750           std::ostringstream errMsg; // Exception message (only valid on Proc 0)
02751           if (myRank == rootRank) {
02752             try {
02753               // Reader for "coordinate" format sparse matrix data.
02754               typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
02755                 global_ordinal_type, scalar_type, STS::isComplex> reader_type;
02756               reader_type reader (pAdder);
02757 
02758               // Read the sparse matrix entries.
02759               std::pair<bool, std::vector<size_t> > results =
02760                 reader.read (in, lineNumber, tolerant, debug);
02761               readSuccess = results.first ? 1 : 0;
02762             }
02763             catch (std::exception& e) {
02764               readSuccess = 0;
02765               errMsg << e.what();
02766             }
02767           }
02768           broadcast (*pComm, rootRank, ptr (&readSuccess));
02769 
02770           // It would be nice to add a "verbose" flag, so that in
02771           // tolerant mode, we could log any bad line number(s) on
02772           // Proc 0.  For now, we just throw if the read fails to
02773           // succeed.
02774           //
02775           // Question: If we're in tolerant mode, and if the read did
02776           // not succeed, should we attempt to call fillComplete()?
02777           TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
02778             "Failed to read the Matrix Market sparse matrix file: "
02779             << errMsg.str());
02780         } // Done reading the matrix entries (stored on Proc 0 for now)
02781 
02782         if (debug && myRank == rootRank) {
02783           cerr << "-- Successfully read the Matrix Market data" << endl;
02784         }
02785 
02786         // In tolerant mode, we need to rebroadcast the matrix
02787         // dimensions, since they may be different after reading the
02788         // actual matrix data.  We only need to broadcast the number
02789         // of rows and columns.  Only Rank 0 needs to know the actual
02790         // global number of entries, since (a) we need to merge
02791         // duplicates on Rank 0 first anyway, and (b) when we
02792         // distribute the entries, each rank other than Rank 0 will
02793         // only need to know how many entries it owns, not the total
02794         // number of entries.
02795         if (tolerant) {
02796           if (debug && myRank == rootRank) {
02797             cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
02798                  << endl
02799                  << "----- Dimensions before: "
02800                  << dims[0] << " x " << dims[1]
02801                  << endl;
02802           }
02803           // Packed coordinate matrix dimensions (numRows, numCols).
02804           Tuple<global_ordinal_type, 2> updatedDims;
02805           if (myRank == rootRank) {
02806             // If one or more bottom rows of the matrix contain no
02807             // entries, then the Adder will report that the number
02808             // of rows is less than that specified in the
02809             // metadata.  We allow this case, and favor the
02810             // metadata so that the zero row(s) will be included.
02811             updatedDims[0] =
02812               std::max (dims[0], pAdder->getAdder()->numRows());
02813             updatedDims[1] = pAdder->getAdder()->numCols();
02814           }
02815           broadcast (*pComm, rootRank, updatedDims);
02816           dims[0] = updatedDims[0];
02817           dims[1] = updatedDims[1];
02818           if (debug && myRank == rootRank) {
02819             cerr << "----- Dimensions after: " << dims[0] << " x "
02820                  << dims[1] << endl;
02821           }
02822         }
02823         else {
02824           // In strict mode, we require that the matrix's metadata and
02825           // its actual data agree, at least somewhat.  In particular,
02826           // the number of rows must agree, since otherwise we cannot
02827           // distribute the matrix correctly.
02828 
02829           // Teuchos::broadcast() doesn't know how to broadcast bools,
02830           // so we use an int with the standard 1 == true, 0 == false
02831           // encoding.
02832           int dimsMatch = 1;
02833           if (myRank == rootRank) {
02834             // If one or more bottom rows of the matrix contain no
02835             // entries, then the Adder will report that the number of
02836             // rows is less than that specified in the metadata.  We
02837             // allow this case, and favor the metadata, but do not
02838             // allow the Adder to think there are more rows in the
02839             // matrix than the metadata says.
02840             if (dims[0] < pAdder->getAdder ()->numRows ()) {
02841               dimsMatch = 0;
02842             }
02843           }
02844           broadcast (*pComm, 0, ptr (&dimsMatch));
02845           if (dimsMatch == 0) {
02846             // We're in an error state anyway, so we might as well
02847             // work a little harder to print an informative error
02848             // message.
02849             //
02850             // Broadcast the Adder's idea of the matrix dimensions
02851             // from Proc 0 to all processes.
02852             Tuple<global_ordinal_type, 2> addersDims;
02853             if (myRank == rootRank) {
02854               addersDims[0] = pAdder->getAdder()->numRows();
02855               addersDims[1] = pAdder->getAdder()->numCols();
02856             }
02857             broadcast (*pComm, 0, addersDims);
02858             TEUCHOS_TEST_FOR_EXCEPTION(
02859               dimsMatch == 0, std::runtime_error,
02860               "The matrix metadata says that the matrix is " << dims[0] << " x "
02861               << dims[1] << ", but the actual data says that the matrix is "
02862               << addersDims[0] << " x " << addersDims[1] << ".  That means the "
02863               "data includes more rows than reported in the metadata.  This "
02864               "is not allowed when parsing in strict mode.  Parse the matrix in "
02865               "tolerant mode to ignore the metadata when it disagrees with the "
02866               "data.");
02867           }
02868         } // Matrix dimensions (# rows, # cols, # entries) agree.
02869 
02870         if (debug && myRank == rootRank) {
02871           cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
02872         }
02873 
02874         // Now that we've read in all the matrix entries from the
02875         // input stream into the adder on Proc 0, post-process them
02876         // into CSR format (still on Proc 0).  This will facilitate
02877         // distributing them to all the processors.
02878         //
02879         // These arrays represent the global matrix data as a CSR
02880         // matrix (with numEntriesPerRow as redundant but convenient
02881         // metadata, since it's computable from rowPtr and vice
02882         // versa).  They are valid only on Proc 0.
02883         ArrayRCP<size_t> numEntriesPerRow;
02884         ArrayRCP<size_t> rowPtr;
02885         ArrayRCP<global_ordinal_type> colInd;
02886         ArrayRCP<scalar_type> values;
02887 
02888         // Proc 0 first merges duplicate entries, and then converts
02889         // the coordinate-format matrix data to CSR.
02890         {
02891           int mergeAndConvertSucceeded = 1;
02892           std::ostringstream errMsg;
02893 
02894           if (myRank == rootRank) {
02895             try {
02896               typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
02897                 global_ordinal_type> element_type;
02898 
02899               // Number of rows in the matrix.  If we are in tolerant
02900               // mode, we've already synchronized dims with the actual
02901               // matrix data.  If in strict mode, we should use dims
02902               // (as read from the file's metadata) rather than the
02903               // matrix data to determine the dimensions.  (The matrix
02904               // data will claim fewer rows than the metadata, if one
02905               // or more rows have no entries stored in the file.)
02906               const size_type numRows = dims[0];
02907 
02908               // Additively merge duplicate matrix entries.
02909               pAdder->getAdder()->merge ();
02910 
02911               // Get a temporary const view of the merged matrix entries.
02912               const std::vector<element_type>& entries =
02913                 pAdder->getAdder()->getEntries();
02914 
02915               // Number of matrix entries (after merging).
02916               const size_t numEntries = (size_t)entries.size();
02917 
02918               if (debug) {
02919                 cerr << "----- Proc 0: Matrix has numRows=" << numRows
02920                      << " rows and numEntries=" << numEntries
02921                      << " entries." << endl;
02922               }
02923 
02924               // Make space for the CSR matrix data.  Converting to
02925               // CSR is easier if we fill numEntriesPerRow with zeros
02926               // at first.
02927               numEntriesPerRow = arcp<size_t> (numRows);
02928               std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
02929               rowPtr = arcp<size_t> (numRows+1);
02930               std::fill (rowPtr.begin(), rowPtr.end(), 0);
02931               colInd = arcp<global_ordinal_type> (numEntries);
02932               values = arcp<scalar_type> (numEntries);
02933 
02934               // Convert from array-of-structs coordinate format to CSR
02935               // (compressed sparse row) format.
02936               global_ordinal_type prvRow = 0;
02937               size_t curPos = 0;
02938               rowPtr[0] = 0;
02939               for (curPos = 0; curPos < numEntries; ++curPos) {
02940                 const element_type& curEntry = entries[curPos];
02941                 const global_ordinal_type curRow = curEntry.rowIndex();
02942                 TEUCHOS_TEST_FOR_EXCEPTION(
02943                   curRow < prvRow, std::logic_error,
02944                   "Row indices are out of order, even though they are supposed "
02945                   "to be sorted.  curRow = " << curRow << ", prvRow = "
02946                   << prvRow << ", at curPos = " << curPos << ".  Please report "
02947                   "this bug to the Tpetra developers.");
02948                 if (curRow > prvRow) {
02949                   for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
02950                     rowPtr[r] = curPos;
02951                   }
02952                   prvRow = curRow;
02953                 }
02954                 numEntriesPerRow[curRow]++;
02955                 colInd[curPos] = curEntry.colIndex();
02956                 values[curPos] = curEntry.value();
02957               }
02958               // rowPtr has one more entry than numEntriesPerRow.  The
02959               // last entry of rowPtr is the number of entries in
02960               // colInd and values.
02961               rowPtr[numRows] = numEntries;
02962             } // Finished conversion to CSR format
02963             catch (std::exception& e) {
02964               mergeAndConvertSucceeded = 0;
02965               errMsg << "Failed to merge sparse matrix entries and convert to "
02966                 "CSR format: " << e.what();
02967             }
02968 
02969             if (debug && mergeAndConvertSucceeded) {
02970               // Number of rows in the matrix.
02971               const size_type numRows = dims[0];
02972               const size_type maxToDisplay = 100;
02973 
02974               cerr << "----- Proc 0: numEntriesPerRow[0.."
02975                    << (numEntriesPerRow.size()-1) << "] ";
02976               if (numRows > maxToDisplay) {
02977                 cerr << "(only showing first and last few entries) ";
02978               }
02979               cerr << "= [";
02980               if (numRows > 0) {
02981                 if (numRows > maxToDisplay) {
02982                   for (size_type k = 0; k < 2; ++k) {
02983                     cerr << numEntriesPerRow[k] << " ";
02984                   }
02985                   cerr << "... ";
02986                   for (size_type k = numRows-2; k < numRows-1; ++k) {
02987                     cerr << numEntriesPerRow[k] << " ";
02988                   }
02989                 }
02990                 else {
02991                   for (size_type k = 0; k < numRows-1; ++k) {
02992                     cerr << numEntriesPerRow[k] << " ";
02993                   }
02994                 }
02995                 cerr << numEntriesPerRow[numRows-1];
02996               } // numRows > 0
02997               cerr << "]" << endl;
02998 
02999               cerr << "----- Proc 0: rowPtr ";
03000               if (numRows > maxToDisplay) {
03001                 cerr << "(only showing first and last few entries) ";
03002               }
03003               cerr << "= [";
03004               if (numRows > maxToDisplay) {
03005                 for (size_type k = 0; k < 2; ++k) {
03006                   cerr << rowPtr[k] << " ";
03007                 }
03008                 cerr << "... ";
03009                 for (size_type k = numRows-2; k < numRows; ++k) {
03010                   cerr << rowPtr[k] << " ";
03011                 }
03012               }
03013               else {
03014                 for (size_type k = 0; k < numRows; ++k) {
03015                   cerr << rowPtr[k] << " ";
03016                 }
03017               }
03018               cerr << rowPtr[numRows] << "]" << endl;
03019 
03020               cerr << "----- Proc 0: colInd = [";
03021               for (size_t k = 0; k < rowPtr[numRows]; ++k) {
03022                 cerr << colInd[k] << " ";
03023               }
03024               cerr << "]" << endl;
03025             }
03026           } // if myRank == rootRank
03027         } // Done converting sparse matrix data to CSR format
03028 
03029         // Now we're done with the Adder, so we can release the
03030         // reference ("free" it) to save space.  This only actually
03031         // does anything on Rank 0, since pAdder is null on all the
03032         // other MPI processes.
03033         pAdder = null;
03034 
03035         // Verify details of the Maps.  Don't count the global number
03036         // of entries in the row Map, since that number doesn't
03037         // correctly count overlap.
03038         if (debug && myRank == rootRank) {
03039           cerr << "-- Verifying Maps" << endl;
03040         }
03041         TEUCHOS_TEST_FOR_EXCEPTION(
03042           as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
03043           std::invalid_argument,
03044           "The range Map has " << rangeMap->getGlobalNumElements ()
03045           << " entries, but the matrix has a global number of rows " << dims[0]
03046           << ".");
03047         TEUCHOS_TEST_FOR_EXCEPTION(
03048           as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
03049           std::invalid_argument,
03050           "The domain Map has " << domainMap->getGlobalNumElements ()
03051           << " entries, but the matrix has a global number of columns "
03052           << dims[1] << ".");
03053 
03054         // Create a row Map which is entirely owned on Proc 0.
03055         RCP<Teuchos::FancyOStream> err = debug ?
03056           Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
03057         RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
03058 
03059         // Create a matrix using this Map, and fill in on Proc 0.  We
03060         // know how many entries there are in each row, so we can use
03061         // static profile.
03062         RCP<sparse_matrix_type> A_proc0 =
03063           rcp (new sparse_matrix_type (gatherRowMap, numEntriesPerRow,
03064                                        Tpetra::StaticProfile));
03065         if (myRank == rootRank) {
03066           if (debug) {
03067             cerr << "-- Proc 0: Filling gather matrix" << endl;
03068           }
03069           ArrayView<const global_ordinal_type> myRows =
03070             gatherRowMap->getNodeElementList ();
03071           if (debug) {
03072             cerr << "---- Rows: " << Teuchos::toString (myRows) << endl;
03073           }
03074           const size_type myNumRows = myRows.size ();
03075 
03076           // Add Proc 0's matrix entries to the CrsMatrix.
03077           const global_ordinal_type indexBase = gatherRowMap->getIndexBase ();
03078           for (size_type i = 0; i < myNumRows; ++i) {
03079             const size_type curPos = as<size_type> (rowPtr[i]);
03080             const local_ordinal_type curNumEntries = numEntriesPerRow[i];
03081             ArrayView<global_ordinal_type> curColInd =
03082               colInd.view (curPos, curNumEntries);
03083             ArrayView<scalar_type> curValues =
03084               values.view (curPos, curNumEntries);
03085 
03086             // Modify the column indices in place to have the right index base.
03087             for (size_type k = 0; k < curNumEntries; ++k) {
03088               curColInd[k] += indexBase;
03089             }
03090             if (debug) {
03091               cerr << "------ Columns: " << Teuchos::toString (curColInd) << endl;
03092               cerr << "------ Values: " << Teuchos::toString (curValues) << endl;
03093             }
03094             // Avoid constructing empty views of ArrayRCP objects.
03095             if (curNumEntries > 0) {
03096               A_proc0->insertGlobalValues (myRows[i], curColInd, curValues);
03097             }
03098           }
03099           // Now we can save space by deallocating numEntriesPerRow,
03100           // rowPtr, colInd, and values, since we've already put those
03101           // data in the matrix.
03102           numEntriesPerRow = null;
03103           rowPtr = null;
03104           colInd = null;
03105           values = null;
03106         } // if myRank == rootRank
03107 
03108         RCP<sparse_matrix_type> A;
03109         if (colMap.is_null ()) {
03110           A = rcp (new sparse_matrix_type (rowMap, 0));
03111         } else {
03112           A = rcp (new sparse_matrix_type (rowMap, colMap, 0));
03113         }
03114         typedef Export<local_ordinal_type, global_ordinal_type, node_type> export_type;
03115         export_type exp (gatherRowMap, rowMap);
03116         A->doExport (*A_proc0, exp, INSERT);
03117 
03118         if (callFillComplete) {
03119           A->fillComplete (domainMap, rangeMap);
03120         }
03121 
03122         // We can't get the dimensions of the matrix until after
03123         // fillComplete() is called.  Thus, we can't do the sanity
03124         // check (dimensions read from the Matrix Market data,
03125         // vs. dimensions reported by the CrsMatrix) unless the user
03126         // asked us to call fillComplete().
03127         //
03128         // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
03129         // what one might think it does, so you have to ask the range
03130         // resp. domain map for the number of rows resp. columns.
03131         if (callFillComplete) {
03132           const int numProcs = pComm->getSize ();
03133 
03134           if (extraDebug && debug) {
03135             const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
03136             const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
03137             if (myRank == rootRank) {
03138               cerr << "-- Matrix is "
03139                    << globalNumRows << " x " << globalNumCols
03140                    << " with " << A->getGlobalNumEntries()
03141                    << " entries, and index base "
03142                    << A->getIndexBase() << "." << endl;
03143             }
03144             pComm->barrier ();
03145             for (int p = 0; p < numProcs; ++p) {
03146               if (myRank == p) {
03147                 cerr << "-- Proc " << p << " owns "
03148                      << A->getNodeNumCols() << " columns, and "
03149                      << A->getNodeNumEntries() << " entries." << endl;
03150               }
03151               pComm->barrier ();
03152             }
03153           } // if (extraDebug && debug)
03154         } // if (callFillComplete)
03155 
03156         if (debug && myRank == rootRank) {
03157           cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
03158                << endl;
03159         }
03160         return A;
03161       }
03162 
03163 
03193       static RCP<multivector_type>
03194       readDenseFile (const std::string& filename,
03195                      const RCP<const comm_type>& comm,
03196                      const RCP<node_type>& node,
03197                      RCP<const map_type>& map,
03198                      const bool tolerant=false,
03199                      const bool debug=false)
03200       {
03201         std::ifstream in;
03202         if (comm->getRank () == 0) { // Only open the file on Proc 0.
03203           in.open (filename.c_str ()); // Destructor closes safely
03204         }
03205         return readDense (in, comm, node, map, tolerant, debug);
03206       }
03207 
03238       static RCP<vector_type>
03239       readVectorFile (const std::string& filename,
03240                      const RCP<const comm_type>& comm,
03241                      const RCP<node_type>& node,
03242                      RCP<const map_type>& map,
03243                      const bool tolerant=false,
03244                      const bool debug=false)
03245       {
03246         std::ifstream in;
03247         if (comm->getRank () == 0) { // Only open the file on Proc 0.
03248           in.open (filename.c_str ()); // Destructor closes safely
03249         }
03250         return readVector (in, comm, node, map, tolerant, debug);
03251       }
03252 
03320       static RCP<multivector_type>
03321       readDense (std::istream& in,
03322                  const RCP<const comm_type>& comm,
03323                  const RCP<node_type>& node,
03324                  RCP<const map_type>& map,
03325                  const bool tolerant=false,
03326                  const bool debug=false)
03327       {
03328         Teuchos::RCP<Teuchos::FancyOStream> err =
03329           Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
03330         return readDenseImpl<scalar_type> (in, comm, node, map,
03331                                            err, tolerant, debug);
03332       }
03333 
03335       static RCP<vector_type>
03336       readVector (std::istream& in,
03337                  const RCP<const comm_type>& comm,
03338                  const RCP<node_type>& node,
03339                  RCP<const map_type>& map,
03340                  const bool tolerant=false,
03341                  const bool debug=false)
03342       {
03343         Teuchos::RCP<Teuchos::FancyOStream> err =
03344           Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
03345         return readVectorImpl<scalar_type> (in, comm, node, map,
03346                                            err, tolerant, debug);
03347       }
03348 
03370       static RCP<const map_type>
03371       readMapFile (const std::string& filename,
03372                    const RCP<const comm_type>& comm,
03373                    const RCP<node_type>& node,
03374                    const bool tolerant=false,
03375                    const bool debug=false)
03376       {
03377         using Teuchos::inOutArg;
03378         using Teuchos::broadcast;
03379         std::ifstream in;
03380 
03381         int success = 1;
03382         if (comm->getRank () == 0) { // Only open the file on Proc 0.
03383           in.open (filename.c_str ()); // Destructor closes safely
03384           if (! in) {
03385             success = 0;
03386           }
03387         }
03388         broadcast<int, int> (*comm, 0, inOutArg (success));
03389         TEUCHOS_TEST_FOR_EXCEPTION(
03390           success == 0, std::runtime_error,
03391           "Tpetra::MatrixMarket::Reader::readMapFile: "
03392           "Failed to read file \"" << filename << "\" on Process 0.");
03393         return readMap (in, comm, node, tolerant, debug);
03394       }
03395 
03396     private:
03397       template<class MultiVectorScalarType>
03398       static RCP<Tpetra::MultiVector<MultiVectorScalarType,
03399                                      local_ordinal_type,
03400                                      global_ordinal_type,
03401                                      node_type> >
03402       readDenseImpl (std::istream& in,
03403                      const RCP<const comm_type>& comm,
03404                      const RCP<node_type>& node,
03405                      RCP<const map_type>& map,
03406                      const Teuchos::RCP<Teuchos::FancyOStream>& err,
03407                      const bool tolerant=false,
03408                      const bool debug=false)
03409       {
03410         using Teuchos::MatrixMarket::Banner;
03411         using Teuchos::MatrixMarket::checkCommentLine;
03412         using Teuchos::as;
03413         using Teuchos::broadcast;
03414         using Teuchos::outArg;
03415         using std::endl;
03416         typedef MultiVectorScalarType ST;
03417         typedef local_ordinal_type LO;
03418         typedef global_ordinal_type GO;
03419         typedef node_type NT;
03420         typedef Teuchos::ScalarTraits<ST> STS;
03421         typedef typename STS::magnitudeType MT;
03422         typedef Teuchos::ScalarTraits<MT> STM;
03423         typedef Tpetra::MultiVector<ST, LO, GO, NT> MV;
03424 
03425         // Rank 0 is the only (MPI) process allowed to read from the
03426         // input stream.
03427         const int myRank = comm->getRank ();
03428 
03429         if (! err.is_null ()) {
03430           err->pushTab ();
03431         }
03432         if (debug) {
03433           *err << myRank << ": readDenseImpl" << endl;
03434         }
03435         if (! err.is_null ()) {
03436           err->pushTab ();
03437         }
03438 
03439         // mfh 17 Feb 2013: It's not strictly necessary that the Comm
03440         // instances be identical and that the Node instances be
03441         // identical.  The essential condition is more complicated to
03442         // test and isn't the same for all Node types.  Thus, we just
03443         // leave it up to the user.
03444 
03445         // // If map is nonnull, check the precondition that its
03446         // // communicator resp. node equal comm resp. node.  Checking
03447         // // now avoids doing a lot of file reading before we detect the
03448         // // violated precondition.
03449         // TEUCHOS_TEST_FOR_EXCEPTION(
03450         //   ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
03451         //   std::invalid_argument, "If you supply a nonnull Map, the Map's "
03452         //   "communicator and node must equal the supplied communicator resp. "
03453         //   "node.");
03454 
03455         // Process 0 will read in the matrix dimensions from the file,
03456         // and broadcast them to all ranks in the given communicator.
03457         // There are only 2 dimensions in the matrix, but we use the
03458         // third element of the Tuple to encode the banner's reported
03459         // data type: "real" == 0, "complex" == 1, and "integer" == 0
03460         // (same as "real").  We don't allow pattern matrices (i.e.,
03461         // graphs) since they only make sense for sparse data.
03462         Tuple<GO, 3> dims;
03463         dims[0] = 0;
03464         dims[1] = 0;
03465 
03466         // Current line number in the input stream.  Only valid on
03467         // Proc 0.  Various calls will modify this depending on the
03468         // number of lines that are read from the input stream.
03469         size_t lineNumber = 1;
03470 
03471         // Capture errors and their messages on Proc 0.
03472         std::ostringstream exMsg;
03473         int localBannerReadSuccess = 1;
03474         int localDimsReadSuccess = 1;
03475 
03476         // Only Proc 0 gets to read matrix data from the input stream.
03477         if (myRank == 0) {
03478           if (debug) {
03479             *err << myRank << ": readDenseImpl: Reading banner line (dense)" << endl;
03480           }
03481 
03482           // The "Banner" tells you whether the input stream
03483           // represents a dense matrix, the symmetry type of the
03484           // matrix, and the type of the data it contains.
03485           RCP<const Banner> pBanner;
03486           try {
03487             pBanner = readBanner (in, lineNumber, tolerant, debug);
03488           } catch (std::exception& e) {
03489             exMsg << e.what ();
03490             localBannerReadSuccess = 0;
03491           }
03492           // Make sure the input stream is the right kind of data.
03493           if (localBannerReadSuccess) {
03494             if (pBanner->matrixType () != "array") {
03495               exMsg << "The Matrix Market file does not contain dense matrix "
03496                 "data.  Its banner (first) line says that its matrix type is \""
03497                 << pBanner->matrixType () << "\", rather that the required "
03498                 "\"array\".";
03499               localBannerReadSuccess = 0;
03500             } else if (pBanner->dataType() == "pattern") {
03501               exMsg << "The Matrix Market file's banner (first) "
03502                 "line claims that the matrix's data type is \"pattern\".  This does "
03503                 "not make sense for a dense matrix, yet the file reports the matrix "
03504                 "as dense.  The only valid data types for a dense matrix are "
03505                 "\"real\", \"complex\", and \"integer\".";
03506               localBannerReadSuccess = 0;
03507             } else {
03508               // Encode the data type reported by the Banner as the
03509               // third element of the dimensions Tuple.
03510               dims[2] = encodeDataType (pBanner->dataType ());
03511             }
03512           } // if we successfully read the banner line
03513 
03514           // At this point, we've successfully read the banner line.
03515           // Now read the dimensions line.
03516           if (localBannerReadSuccess) {
03517             if (debug) {
03518               *err << myRank << ": readDenseImpl: Reading dimensions line (dense)" << endl;
03519             }
03520             // Keep reading lines from the input stream until we find
03521             // a non-comment line, or until we run out of lines.  The
03522             // latter is an error, since every "array" format Matrix
03523             // Market file must have a dimensions line after the
03524             // banner (even if the matrix has zero rows or columns, or
03525             // zero entries).
03526             std::string line;
03527             bool commentLine = true;
03528 
03529             while (commentLine) {
03530               // Test whether it is even valid to read from the input
03531               // stream wrapping the line.
03532               if (in.eof () || in.fail ()) {
03533                 exMsg << "Unable to get array dimensions line (at all) from line "
03534                       << lineNumber << " of input stream.  The input stream "
03535                       << "claims that it is "
03536                       << (in.eof() ? "at end-of-file." : "in a failed state.");
03537                 localDimsReadSuccess = 0;
03538               } else {
03539                 // Try to get the next line from the input stream.
03540                 if (getline (in, line)) {
03541                   ++lineNumber; // We did actually read a line.
03542                 }
03543                 // Is the current line a comment line?  Ignore start
03544                 // and size; they are only useful for reading the
03545                 // actual matrix entries.  (We could use them here as
03546                 // an optimization, but we've chosen not to.)
03547                 size_t start = 0, size = 0;
03548                 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
03549               } // whether we failed to read the line at all
03550             } // while the line we just read is a comment line
03551 
03552             //
03553             // Get <numRows> <numCols> from the line we just read.
03554             //
03555             std::istringstream istr (line);
03556 
03557             // Test whether it is even valid to read from the input
03558             // stream wrapping the line.
03559             if (istr.eof () || istr.fail ()) {
03560               exMsg << "Unable to read any data from line " << lineNumber
03561                     << " of input; the line should contain the matrix dimensions "
03562                     << "\"<numRows> <numCols>\".";
03563               localDimsReadSuccess = 0;
03564             } else { // It's valid to read from the line.
03565               GO theNumRows = 0;
03566               istr >> theNumRows; // Read in the number of rows.
03567               if (istr.fail ()) {
03568                 exMsg << "Failed to get number of rows from line "
03569                       << lineNumber << " of input; the line should contains the "
03570                       << "matrix dimensions \"<numRows> <numCols>\".";
03571                 localDimsReadSuccess = 0;
03572               } else { // We successfully read the number of rows
03573                 dims[0] = theNumRows; // Save the number of rows
03574                 if (istr.eof ()) { // Do we still have data to read?
03575                   exMsg << "No more data after number of rows on line "
03576                         << lineNumber << " of input; the line should contain the "
03577                         << "matrix dimensions \"<numRows> <numCols>\".";
03578                   localDimsReadSuccess = 0;
03579                 } else { // Still data left to read; read in number of columns.
03580                   GO theNumCols = 0;
03581                   istr >> theNumCols; // Read in the number of columns
03582                   if (istr.fail ()) {
03583                     exMsg << "Failed to get number of columns from line "
03584                           << lineNumber << " of input; the line should contain "
03585                           << "the matrix dimensions \"<numRows> <numCols>\".";
03586                     localDimsReadSuccess = 0;
03587                   } else { // We successfully read the number of columns
03588                     dims[1] = theNumCols; // Save the number of columns
03589                   } // if istr.fail ()
03590                 } // if istr.eof ()
03591               } // if we read the number of rows
03592             } // if the input stream wrapping the dims line was (in)valid
03593           } // if we successfully read the banner line
03594         } // if (myRank == 0)
03595 
03596         // Broadcast the matrix dimensions, the encoded data type, and
03597         // whether or not Proc 0 succeeded in reading the banner and
03598         // dimensions.
03599         Tuple<GO, 5> bannerDimsReadResult;
03600         if (myRank == 0) {
03601           bannerDimsReadResult[0] = dims[0]; // numRows
03602           bannerDimsReadResult[1] = dims[1]; // numCols
03603           bannerDimsReadResult[2] = dims[2]; // encoded data type
03604           bannerDimsReadResult[3] = localBannerReadSuccess;
03605           bannerDimsReadResult[4] = localDimsReadSuccess;
03606         }
03607         // Broadcast matrix dimensions and the encoded data type from
03608         // Proc 0 to all the MPI processes.
03609         broadcast (*comm, 0, bannerDimsReadResult);
03610 
03611         TEUCHOS_TEST_FOR_EXCEPTION(
03612           bannerDimsReadResult[3] == 0, std::runtime_error,
03613           "Failed to read banner line: " << exMsg.str ());
03614         TEUCHOS_TEST_FOR_EXCEPTION(
03615           bannerDimsReadResult[4] == 0, std::runtime_error,
03616           "Failed to read matrix dimensions line: " << exMsg.str ());
03617         if (myRank != 0) {
03618           dims[0] = bannerDimsReadResult[0];
03619           dims[1] = bannerDimsReadResult[1];
03620           dims[2] = bannerDimsReadResult[2];
03621         }
03622 
03623         // Tpetra objects want the matrix dimensions in these types.
03624         const global_size_t numRows = static_cast<global_size_t> (dims[0]);
03625         const size_t numCols = static_cast<size_t> (dims[1]);
03626 
03627         // Make a "Proc 0 owns everything" Map that we will use to
03628         // read in the multivector entries in the correct order on
03629         // Proc 0.  This must be a collective
03630         RCP<const map_type> proc0Map; // "Proc 0 owns everything" Map
03631         if (map.is_null ()) {
03632           // The user didn't supply a Map.  Make a contiguous
03633           // distributed Map for them, using the read-in multivector
03634           // dimensions.
03635           map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
03636           const size_t localNumRows = (myRank == 0) ? numRows : 0;
03637           proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows, comm, node);
03638         }
03639         else { // The user supplied a Map.
03640           proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
03641         }
03642 
03643         // Make a multivector X owned entirely by Proc 0.
03644         RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
03645 
03646         //
03647         // On Proc 0, read the Matrix Market data from the input
03648         // stream into the multivector X.
03649         //
03650         int localReadDataSuccess = 1;
03651         if (myRank == 0) {
03652           try {
03653             if (debug) {
03654               *err << myRank << ": readDenseImpl: Reading matrix data (dense)"
03655                    << endl;
03656             }
03657 
03658             // Make sure that we can get a 1-D view of X.
03659             TEUCHOS_TEST_FOR_EXCEPTION(
03660               ! X->isConstantStride (), std::logic_error,
03661               "Can't get a 1-D view of the entries of the MultiVector X on "
03662               "Process 0, because the stride between the columns of X is not "
03663               "constant.  This shouldn't happen because we just created X and "
03664               "haven't filled it in yet.  Please report this bug to the Tpetra "
03665               "developers.");
03666 
03667             // Get a writeable 1-D view of the entries of X.  Rank 0
03668             // owns all of them.  The view will expire at the end of
03669             // scope, so (if necessary) it will be written back to X
03670             // at this time.
03671             ArrayRCP<ST> X_view = X->get1dViewNonConst ();
03672             TEUCHOS_TEST_FOR_EXCEPTION(
03673               as<global_size_t> (X_view.size ()) < numRows * numCols,
03674               std::logic_error,
03675               "The view of X has size " << X_view << " which is not enough to "
03676               "accommodate the expected number of entries numRows*numCols = "
03677               << numRows << "*" << numCols << " = " << numRows*numCols << ".  "
03678               "Please report this bug to the Tpetra developers.");
03679             const size_t stride = X->getStride ();
03680 
03681             // The third element of the dimensions Tuple encodes the data
03682             // type reported by the Banner: "real" == 0, "complex" == 1,
03683             // "integer" == 0 (same as "real"), "pattern" == 2.  We do not
03684             // allow dense matrices to be pattern matrices, so dims[2] ==
03685             // 0 or 1.  We've already checked for this above.
03686             const bool isComplex = (dims[2] == 1);
03687             size_type count = 0, curRow = 0, curCol = 0;
03688 
03689             std::string line;
03690             while (getline (in, line)) {
03691               ++lineNumber;
03692               // Is the current line a comment line?  If it's not,
03693               // line.substr(start,size) contains the data.
03694               size_t start = 0, size = 0;
03695               const bool commentLine =
03696                 checkCommentLine (line, start, size, lineNumber, tolerant);
03697               if (! commentLine) {
03698                 // Make sure we have room in which to put the new matrix
03699                 // entry.  We check this only after checking for a
03700                 // comment line, because there may be one or more
03701                 // comment lines at the end of the file.  In tolerant
03702                 // mode, we simply ignore any extra data.
03703                 if (count >= X_view.size()) {
03704                   if (tolerant) {
03705                     break;
03706                   }
03707                   else {
03708                     TEUCHOS_TEST_FOR_EXCEPTION(
03709                        count >= X_view.size(),
03710                        std::runtime_error,
03711                        "The Matrix Market input stream has more data in it than "
03712                        "its metadata reported.  Current line number is "
03713                        << lineNumber << ".");
03714                   }
03715                 }
03716 
03717                 // mfh 19 Dec 2012: Ignore everything up to the initial
03718                 // colon.  writeDense() has the option to print out the
03719                 // global row index in front of each entry, followed by
03720                 // a colon and space.
03721                 {
03722                   const size_t pos = line.substr (start, size).find (':');
03723                   if (pos != std::string::npos) {
03724                     start = pos+1;
03725                   }
03726                 }
03727                 std::istringstream istr (line.substr (start, size));
03728                 // Does the line contain anything at all?  Can we
03729                 // safely read from the input stream wrapping the
03730                 // line?
03731                 if (istr.eof() || istr.fail()) {
03732                   // In tolerant mode, simply ignore the line.
03733                   if (tolerant) {
03734                     break;
03735                   }
03736                   // We repeat the full test here so the exception
03737                   // message is more informative.
03738                   TEUCHOS_TEST_FOR_EXCEPTION(
03739                     ! tolerant && (istr.eof() || istr.fail()),
03740                     std::runtime_error,
03741                     "Line " << lineNumber << " of the Matrix Market file is "
03742                     "empty, or we cannot read from it for some other reason.");
03743                 }
03744                 // Current matrix entry to read in.
03745                 ST val = STS::zero();
03746                 // Real and imaginary parts of the current matrix entry.
03747                 // The imaginary part is zero if the matrix is real-valued.
03748                 MT real = STM::zero(), imag = STM::zero();
03749 
03750                 // isComplex refers to the input stream's data, not to
03751                 // the scalar type S.  It's OK to read real-valued
03752                 // data into a matrix storing complex-valued data; in
03753                 // that case, all entries' imaginary parts are zero.
03754                 if (isComplex) {
03755                   // STS::real() and STS::imag() return a copy of
03756                   // their respective components, not a writeable
03757                   // reference.  Otherwise we could just assign to
03758                   // them using the istream extraction operator (>>).
03759                   // That's why we have separate magnitude type "real"
03760                   // and "imag" variables.
03761 
03762                   // Attempt to read the real part of the current entry.
03763                   istr >> real;
03764                   if (istr.fail()) {
03765                     TEUCHOS_TEST_FOR_EXCEPTION(
03766                       ! tolerant && istr.eof(), std::runtime_error,
03767                       "Failed to get the real part of a complex-valued matrix "
03768                       "entry from line " << lineNumber << " of the Matrix Market "
03769                       "file.");
03770                     // In tolerant mode, just skip bad lines.
03771                     if (tolerant) {
03772                       break;
03773                     }
03774                   } else if (istr.eof()) {
03775                     TEUCHOS_TEST_FOR_EXCEPTION(
03776                       ! tolerant && istr.eof(), std::runtime_error,
03777                       "Missing imaginary part of a complex-valued matrix entry "
03778                       "on line " << lineNumber << " of the Matrix Market file.");
03779                     // In tolerant mode, let any missing imaginary part be 0.
03780                   } else {
03781                     // Attempt to read the imaginary part of the current
03782                     // matrix entry.
03783                     istr >> imag;
03784                     TEUCHOS_TEST_FOR_EXCEPTION(
03785                       ! tolerant && istr.fail(), std::runtime_error,
03786                       "Failed to get the imaginary part of a complex-valued "
03787                       "matrix entry from line " << lineNumber << " of the "
03788                       "Matrix Market file.");
03789                     // In tolerant mode, let any missing or corrupted
03790                     // imaginary part be 0.
03791                   }
03792                 } else { // Matrix Market file contains real-valued data.
03793                   // Attempt to read the current matrix entry.
03794                   istr >> real;
03795                   TEUCHOS_TEST_FOR_EXCEPTION(
03796                     ! tolerant && istr.fail(), std::runtime_error,
03797                     "Failed to get a real-valued matrix entry from line "
03798                     << lineNumber << " of the Matrix Market file.");
03799                   // In tolerant mode, simply ignore the line if
03800                   // we failed to read a matrix entry.
03801                   if (istr.fail() && tolerant) {
03802                     break;
03803                   }
03804                 }
03805                 // In tolerant mode, we simply let pass through whatever
03806                 // data we got.
03807                 TEUCHOS_TEST_FOR_EXCEPTION(
03808                   ! tolerant && istr.fail(), std::runtime_error,
03809                   "Failed to read matrix data from line " << lineNumber
03810                   << " of the Matrix Market file.");
03811 
03812                 // Assign val = ST(real, imag).
03813                 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
03814 
03815                 curRow = count % numRows;
03816                 curCol = count / numRows;
03817                 X_view[curRow + curCol*stride] = val;
03818                 ++count;
03819               } // if not a comment line
03820             } // while there are still lines in the file, get the next one
03821 
03822             TEUCHOS_TEST_FOR_EXCEPTION(
03823               ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
03824               std::runtime_error,
03825               "The Matrix Market metadata reports that the dense matrix is "
03826               << numRows <<  " x " << numCols << ", and thus has "
03827               << numRows*numCols << " total entries, but we only found " << count
03828               << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
03829           } catch (std::exception& e) {
03830             exMsg << e.what ();
03831             localReadDataSuccess = 0;
03832           }
03833         } // if (myRank == 0)
03834 
03835         if (debug) {
03836           *err << myRank << ": readDenseImpl: done reading data" << endl;
03837         }
03838 
03839         // Synchronize on whether Proc 0 successfully read the data.
03840         int globalReadDataSuccess = localReadDataSuccess;
03841         broadcast (*comm, 0, outArg (globalReadDataSuccess));
03842         TEUCHOS_TEST_FOR_EXCEPTION(
03843           globalReadDataSuccess == 0, std::runtime_error,
03844           "Failed to read the multivector's data: " << exMsg.str ());
03845 
03846         // If there's only one MPI process and the user didn't supply
03847         // a Map (i.e., pMap is null), we're done.  Set pMap to the
03848         // Map used to distribute X, and return X.
03849         if (comm->getSize () == 1 && map.is_null ()) {
03850           map = proc0Map;
03851           if (! err.is_null ()) {
03852             err->popTab ();
03853           }
03854           if (debug) {
03855             *err << myRank << ": readDenseImpl: done" << endl;
03856           }
03857           if (! err.is_null ()) {
03858             err->popTab ();
03859           }
03860           return X;
03861         }
03862 
03863         if (debug) {
03864           *err << myRank << ": readDenseImpl: Creating target MV" << endl;
03865         }
03866 
03867         // Make a multivector Y with the distributed map pMap.
03868         RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
03869 
03870         if (debug) {
03871           *err << myRank << ": readDenseImpl: Creating Export" << endl;
03872         }
03873 
03874         // Make an Export object that will export X to Y.  First
03875         // argument is the source map, second argument is the target
03876         // map.
03877         Export<LO, GO, NT> exporter (proc0Map, map, err);
03878 
03879         if (debug) {
03880           *err << myRank << ": readDenseImpl: Exporting" << endl;
03881         }
03882         // Export X into Y.
03883         Y->doExport (*X, exporter, INSERT);
03884 
03885         if (! err.is_null ()) {
03886           err->popTab ();
03887         }
03888         if (debug) {
03889           *err << myRank << ": readDenseImpl: done" << endl;
03890         }
03891         if (! err.is_null ()) {
03892           err->popTab ();
03893         }
03894 
03895         // Y is distributed over all process(es) in the communicator.
03896         return Y;
03897       }
03898 
03899 
03900       template<class VectorScalarType>
03901       static RCP<Tpetra::Vector<VectorScalarType,
03902                                      local_ordinal_type,
03903                                      global_ordinal_type,
03904                                      node_type> >
03905       readVectorImpl (std::istream& in,
03906                      const RCP<const comm_type>& comm,
03907                      const RCP<node_type>& node,
03908                      RCP<const map_type>& map,
03909                      const Teuchos::RCP<Teuchos::FancyOStream>& err,
03910                      const bool tolerant=false,
03911                      const bool debug=false)
03912       {
03913         using Teuchos::MatrixMarket::Banner;
03914         using Teuchos::MatrixMarket::checkCommentLine;
03915         using Teuchos::as;
03916         using Teuchos::broadcast;
03917         using Teuchos::outArg;
03918         using std::endl;
03919         typedef VectorScalarType ST;
03920         typedef local_ordinal_type LO;
03921         typedef global_ordinal_type GO;
03922         typedef node_type NT;
03923         typedef Teuchos::ScalarTraits<ST> STS;
03924         typedef typename STS::magnitudeType MT;
03925         typedef Teuchos::ScalarTraits<MT> STM;
03926         typedef Tpetra::Vector<ST, LO, GO, NT> MV;
03927 
03928         // Rank 0 is the only (MPI) process allowed to read from the
03929         // input stream.
03930         const int myRank = comm->getRank ();
03931 
03932         if (! err.is_null ()) {
03933           err->pushTab ();
03934         }
03935         if (debug) {
03936           *err << myRank << ": readDenseImpl" << endl;
03937         }
03938         if (! err.is_null ()) {
03939           err->pushTab ();
03940         }
03941 
03942         // mfh 17 Feb 2013: It's not strictly necessary that the Comm
03943         // instances be identical and that the Node instances be
03944         // identical.  The essential condition is more complicated to
03945         // test and isn't the same for all Node types.  Thus, we just
03946         // leave it up to the user.
03947 
03948         // // If map is nonnull, check the precondition that its
03949         // // communicator resp. node equal comm resp. node.  Checking
03950         // // now avoids doing a lot of file reading before we detect the
03951         // // violated precondition.
03952         // TEUCHOS_TEST_FOR_EXCEPTION(
03953         //   ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
03954         //   std::invalid_argument, "If you supply a nonnull Map, the Map's "
03955         //   "communicator and node must equal the supplied communicator resp. "
03956         //   "node.");
03957 
03958         // Process 0 will read in the matrix dimensions from the file,
03959         // and broadcast them to all ranks in the given communicator.
03960         // There are only 2 dimensions in the matrix, but we use the
03961         // third element of the Tuple to encode the banner's reported
03962         // data type: "real" == 0, "complex" == 1, and "integer" == 0
03963         // (same as "real").  We don't allow pattern matrices (i.e.,
03964         // graphs) since they only make sense for sparse data.
03965         Tuple<GO, 3> dims;
03966         dims[0] = 0;
03967         dims[1] = 0;
03968 
03969         // Current line number in the input stream.  Only valid on
03970         // Proc 0.  Various calls will modify this depending on the
03971         // number of lines that are read from the input stream.
03972         size_t lineNumber = 1;
03973 
03974         // Capture errors and their messages on Proc 0.
03975         std::ostringstream exMsg;
03976         int localBannerReadSuccess = 1;
03977         int localDimsReadSuccess = 1;
03978 
03979         // Only Proc 0 gets to read matrix data from the input stream.
03980         if (myRank == 0) {
03981           if (debug) {
03982             *err << myRank << ": readDenseImpl: Reading banner line (dense)" << endl;
03983           }
03984 
03985           // The "Banner" tells you whether the input stream
03986           // represents a dense matrix, the symmetry type of the
03987           // matrix, and the type of the data it contains.
03988           RCP<const Banner> pBanner;
03989           try {
03990             pBanner = readBanner (in, lineNumber, tolerant, debug);
03991           } catch (std::exception& e) {
03992             exMsg << e.what ();
03993             localBannerReadSuccess = 0;
03994           }
03995           // Make sure the input stream is the right kind of data.
03996           if (localBannerReadSuccess) {
03997             if (pBanner->matrixType () != "array") {
03998               exMsg << "The Matrix Market file does not contain dense matrix "
03999                 "data.  Its banner (first) line says that its matrix type is \""
04000                 << pBanner->matrixType () << "\", rather that the required "
04001                 "\"array\".";
04002               localBannerReadSuccess = 0;
04003             } else if (pBanner->dataType() == "pattern") {
04004               exMsg << "The Matrix Market file's banner (first) "
04005                 "line claims that the matrix's data type is \"pattern\".  This does "
04006                 "not make sense for a dense matrix, yet the file reports the matrix "
04007                 "as dense.  The only valid data types for a dense matrix are "
04008                 "\"real\", \"complex\", and \"integer\".";
04009               localBannerReadSuccess = 0;
04010             } else {
04011               // Encode the data type reported by the Banner as the
04012               // third element of the dimensions Tuple.
04013               dims[2] = encodeDataType (pBanner->dataType ());
04014             }
04015           } // if we successfully read the banner line
04016 
04017           // At this point, we've successfully read the banner line.
04018           // Now read the dimensions line.
04019           if (localBannerReadSuccess) {
04020             if (debug) {
04021               *err << myRank << ": readDenseImpl: Reading dimensions line (dense)" << endl;
04022             }
04023             // Keep reading lines from the input stream until we find
04024             // a non-comment line, or until we run out of lines.  The
04025             // latter is an error, since every "array" format Matrix
04026             // Market file must have a dimensions line after the
04027             // banner (even if the matrix has zero rows or columns, or
04028             // zero entries).
04029             std::string line;
04030             bool commentLine = true;
04031 
04032             while (commentLine) {
04033               // Test whether it is even valid to read from the input
04034               // stream wrapping the line.
04035               if (in.eof () || in.fail ()) {
04036                 exMsg << "Unable to get array dimensions line (at all) from line "
04037                       << lineNumber << " of input stream.  The input stream "
04038                       << "claims that it is "
04039                       << (in.eof() ? "at end-of-file." : "in a failed state.");
04040                 localDimsReadSuccess = 0;
04041               } else {
04042                 // Try to get the next line from the input stream.
04043                 if (getline (in, line)) {
04044                   ++lineNumber; // We did actually read a line.
04045                 }
04046                 // Is the current line a comment line?  Ignore start
04047                 // and size; they are only useful for reading the
04048                 // actual matrix entries.  (We could use them here as
04049                 // an optimization, but we've chosen not to.)
04050                 size_t start = 0, size = 0;
04051                 commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
04052               } // whether we failed to read the line at all
04053             } // while the line we just read is a comment line
04054 
04055             //
04056             // Get <numRows> <numCols> from the line we just read.
04057             //
04058             std::istringstream istr (line);
04059 
04060             // Test whether it is even valid to read from the input
04061             // stream wrapping the line.
04062             if (istr.eof () || istr.fail ()) {
04063               exMsg << "Unable to read any data from line " << lineNumber
04064                     << " of input; the line should contain the matrix dimensions "
04065                     << "\"<numRows> <numCols>\".";
04066               localDimsReadSuccess = 0;
04067             } else { // It's valid to read from the line.
04068               GO theNumRows = 0;
04069               istr >> theNumRows; // Read in the number of rows.
04070               if (istr.fail ()) {
04071                 exMsg << "Failed to get number of rows from line "
04072                       << lineNumber << " of input; the line should contains the "
04073                       << "matrix dimensions \"<numRows> <numCols>\".";
04074                 localDimsReadSuccess = 0;
04075               } else { // We successfully read the number of rows
04076                 dims[0] = theNumRows; // Save the number of rows
04077                 if (istr.eof ()) { // Do we still have data to read?
04078                   exMsg << "No more data after number of rows on line "
04079                         << lineNumber << " of input; the line should contain the "
04080                         << "matrix dimensions \"<numRows> <numCols>\".";
04081                   localDimsReadSuccess = 0;
04082                 } else { // Still data left to read; read in number of columns.
04083                   GO theNumCols = 0;
04084                   istr >> theNumCols; // Read in the number of columns
04085                   if (istr.fail ()) {
04086                     exMsg << "Failed to get number of columns from line "
04087                           << lineNumber << " of input; the line should contain "
04088                           << "the matrix dimensions \"<numRows> <numCols>\".";
04089                     localDimsReadSuccess = 0;
04090                   } else { // We successfully read the number of columns
04091                     dims[1] = theNumCols; // Save the number of columns
04092                   } // if istr.fail ()
04093                 } // if istr.eof ()
04094               } // if we read the number of rows
04095             } // if the input stream wrapping the dims line was (in)valid
04096           } // if we successfully read the banner line
04097         } // if (myRank == 0)
04098 
04099         // Check if file has a Vector
04100         if (dims[1]!=1) {
04101           exMsg << "File does not contain a 1D Vector";
04102           localDimsReadSuccess = 0;
04103         }
04104 
04105         // Broadcast the matrix dimensions, the encoded data type, and
04106         // whether or not Proc 0 succeeded in reading the banner and
04107         // dimensions.
04108         Tuple<GO, 5> bannerDimsReadResult;
04109         if (myRank == 0) {
04110           bannerDimsReadResult[0] = dims[0]; // numRows
04111           bannerDimsReadResult[1] = dims[1]; // numCols
04112           bannerDimsReadResult[2] = dims[2]; // encoded data type
04113           bannerDimsReadResult[3] = localBannerReadSuccess;
04114           bannerDimsReadResult[4] = localDimsReadSuccess;
04115         }
04116 
04117         // Broadcast matrix dimensions and the encoded data type from
04118         // Proc 0 to all the MPI processes.
04119         broadcast (*comm, 0, bannerDimsReadResult);
04120 
04121         TEUCHOS_TEST_FOR_EXCEPTION(
04122           bannerDimsReadResult[3] == 0, std::runtime_error,
04123           "Failed to read banner line: " << exMsg.str ());
04124         TEUCHOS_TEST_FOR_EXCEPTION(
04125           bannerDimsReadResult[4] == 0, std::runtime_error,
04126           "Failed to read matrix dimensions line: " << exMsg.str ());
04127         if (myRank != 0) {
04128           dims[0] = bannerDimsReadResult[0];
04129           dims[1] = bannerDimsReadResult[1];
04130           dims[2] = bannerDimsReadResult[2];
04131         }
04132 
04133         // Tpetra objects want the matrix dimensions in these types.
04134         const global_size_t numRows = static_cast<global_size_t> (dims[0]);
04135         const size_t numCols = static_cast<size_t> (dims[1]);
04136 
04137         // Make a "Proc 0 owns everything" Map that we will use to
04138         // read in the multivector entries in the correct order on
04139         // Proc 0.  This must be a collective
04140         RCP<const map_type> proc0Map; // "Proc 0 owns everything" Map
04141         if (map.is_null ()) {
04142           // The user didn't supply a Map.  Make a contiguous
04143           // distributed Map for them, using the read-in multivector
04144           // dimensions.
04145           map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
04146           const size_t localNumRows = (myRank == 0) ? numRows : 0;
04147           proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows, comm, node);
04148         }
04149         else { // The user supplied a Map.
04150           proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
04151         }
04152 
04153         // Make a multivector X owned entirely by Proc 0.
04154         RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
04155 
04156         //
04157         // On Proc 0, read the Matrix Market data from the input
04158         // stream into the multivector X.
04159         //
04160         int localReadDataSuccess = 1;
04161         if (myRank == 0) {
04162           try {
04163             if (debug) {
04164               *err << myRank << ": readDenseImpl: Reading matrix data (dense)"
04165                    << endl;
04166             }
04167 
04168             // Make sure that we can get a 1-D view of X.
04169             TEUCHOS_TEST_FOR_EXCEPTION(
04170               ! X->isConstantStride (), std::logic_error,
04171               "Can't get a 1-D view of the entries of the MultiVector X on "
04172               "Process 0, because the stride between the columns of X is not "
04173               "constant.  This shouldn't happen because we just created X and "
04174               "haven't filled it in yet.  Please report this bug to the Tpetra "
04175               "developers.");
04176 
04177             // Get a writeable 1-D view of the entries of X.  Rank 0
04178             // owns all of them.  The view will expire at the end of
04179             // scope, so (if necessary) it will be written back to X
04180             // at this time.
04181             ArrayRCP<ST> X_view = X->get1dViewNonConst ();
04182             TEUCHOS_TEST_FOR_EXCEPTION(
04183               as<global_size_t> (X_view.size ()) < numRows * numCols,
04184               std::logic_error,
04185               "The view of X has size " << X_view << " which is not enough to "
04186               "accommodate the expected number of entries numRows*numCols = "
04187               << numRows << "*" << numCols << " = " << numRows*numCols << ".  "
04188               "Please report this bug to the Tpetra developers.");
04189             const size_t stride = X->getStride ();
04190 
04191             // The third element of the dimensions Tuple encodes the data
04192             // type reported by the Banner: "real" == 0, "complex" == 1,
04193             // "integer" == 0 (same as "real"), "pattern" == 2.  We do not
04194             // allow dense matrices to be pattern matrices, so dims[2] ==
04195             // 0 or 1.  We've already checked for this above.
04196             const bool isComplex = (dims[2] == 1);
04197             size_type count = 0, curRow = 0, curCol = 0;
04198 
04199             std::string line;
04200             while (getline (in, line)) {
04201               ++lineNumber;
04202               // Is the current line a comment line?  If it's not,
04203               // line.substr(start,size) contains the data.
04204               size_t start = 0, size = 0;
04205               const bool commentLine =
04206                 checkCommentLine (line, start, size, lineNumber, tolerant);
04207               if (! commentLine) {
04208                 // Make sure we have room in which to put the new matrix
04209                 // entry.  We check this only after checking for a
04210                 // comment line, because there may be one or more
04211                 // comment lines at the end of the file.  In tolerant
04212                 // mode, we simply ignore any extra data.
04213                 if (count >= X_view.size()) {
04214                   if (tolerant) {
04215                     break;
04216                   }
04217                   else {
04218                     TEUCHOS_TEST_FOR_EXCEPTION(
04219                        count >= X_view.size(),
04220                        std::runtime_error,
04221                        "The Matrix Market input stream has more data in it than "
04222                        "its metadata reported.  Current line number is "
04223                        << lineNumber << ".");
04224                   }
04225                 }
04226 
04227                 // mfh 19 Dec 2012: Ignore everything up to the initial
04228                 // colon.  writeDense() has the option to print out the
04229                 // global row index in front of each entry, followed by
04230                 // a colon and space.
04231                 {
04232                   const size_t pos = line.substr (start, size).find (':');
04233                   if (pos != std::string::npos) {
04234                     start = pos+1;
04235                   }
04236                 }
04237                 std::istringstream istr (line.substr (start, size));
04238                 // Does the line contain anything at all?  Can we
04239                 // safely read from the input stream wrapping the
04240                 // line?
04241                 if (istr.eof() || istr.fail()) {
04242                   // In tolerant mode, simply ignore the line.
04243                   if (tolerant) {
04244                     break;
04245                   }
04246                   // We repeat the full test here so the exception
04247                   // message is more informative.
04248                   TEUCHOS_TEST_FOR_EXCEPTION(
04249                     ! tolerant && (istr.eof() || istr.fail()),
04250                     std::runtime_error,
04251                     "Line " << lineNumber << " of the Matrix Market file is "
04252                     "empty, or we cannot read from it for some other reason.");
04253                 }
04254                 // Current matrix entry to read in.
04255                 ST val = STS::zero();
04256                 // Real and imaginary parts of the current matrix entry.
04257                 // The imaginary part is zero if the matrix is real-valued.
04258                 MT real = STM::zero(), imag = STM::zero();
04259 
04260                 // isComplex refers to the input stream's data, not to
04261                 // the scalar type S.  It's OK to read real-valued
04262                 // data into a matrix storing complex-valued data; in
04263                 // that case, all entries' imaginary parts are zero.
04264                 if (isComplex) {
04265                   // STS::real() and STS::imag() return a copy of
04266                   // their respective components, not a writeable
04267                   // reference.  Otherwise we could just assign to
04268                   // them using the istream extraction operator (>>).
04269                   // That's why we have separate magnitude type "real"
04270                   // and "imag" variables.
04271 
04272                   // Attempt to read the real part of the current entry.
04273                   istr >> real;
04274                   if (istr.fail()) {
04275                     TEUCHOS_TEST_FOR_EXCEPTION(
04276                       ! tolerant && istr.eof(), std::runtime_error,
04277                       "Failed to get the real part of a complex-valued matrix "
04278                       "entry from line " << lineNumber << " of the Matrix Market "
04279                       "file.");
04280                     // In tolerant mode, just skip bad lines.
04281                     if (tolerant) {
04282                       break;
04283                     }
04284                   } else if (istr.eof()) {
04285                     TEUCHOS_TEST_FOR_EXCEPTION(
04286                       ! tolerant && istr.eof(), std::runtime_error,
04287                       "Missing imaginary part of a complex-valued matrix entry "
04288                       "on line " << lineNumber << " of the Matrix Market file.");
04289                     // In tolerant mode, let any missing imaginary part be 0.
04290                   } else {
04291                     // Attempt to read the imaginary part of the current
04292                     // matrix entry.
04293                     istr >> imag;
04294                     TEUCHOS_TEST_FOR_EXCEPTION(
04295                       ! tolerant && istr.fail(), std::runtime_error,
04296                       "Failed to get the imaginary part of a complex-valued "
04297                       "matrix entry from line " << lineNumber << " of the "
04298                       "Matrix Market file.");
04299                     // In tolerant mode, let any missing or corrupted
04300                     // imaginary part be 0.
04301                   }
04302                 } else { // Matrix Market file contains real-valued data.
04303                   // Attempt to read the current matrix entry.
04304                   istr >> real;
04305                   TEUCHOS_TEST_FOR_EXCEPTION(
04306                     ! tolerant && istr.fail(), std::runtime_error,
04307                     "Failed to get a real-valued matrix entry from line "
04308                     << lineNumber << " of the Matrix Market file.");
04309                   // In tolerant mode, simply ignore the line if
04310                   // we failed to read a matrix entry.
04311                   if (istr.fail() && tolerant) {
04312                     break;
04313                   }
04314                 }
04315                 // In tolerant mode, we simply let pass through whatever
04316                 // data we got.
04317                 TEUCHOS_TEST_FOR_EXCEPTION(
04318                   ! tolerant && istr.fail(), std::runtime_error,
04319                   "Failed to read matrix data from line " << lineNumber
04320                   << " of the Matrix Market file.");
04321 
04322                 // Assign val = ST(real, imag).
04323                 Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
04324 
04325                 curRow = count % numRows;
04326                 curCol = count / numRows;
04327                 X_view[curRow + curCol*stride] = val;
04328                 ++count;
04329               } // if not a comment line
04330             } // while there are still lines in the file, get the next one
04331 
04332             TEUCHOS_TEST_FOR_EXCEPTION(
04333               ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
04334               std::runtime_error,
04335               "The Matrix Market metadata reports that the dense matrix is "
04336               << numRows <<  " x " << numCols << ", and thus has "
04337               << numRows*numCols << " total entries, but we only found " << count
04338               << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
04339           } catch (std::exception& e) {
04340             exMsg << e.what ();
04341             localReadDataSuccess = 0;
04342           }
04343         } // if (myRank == 0)
04344 
04345         if (debug) {
04346           *err << myRank << ": readDenseImpl: done reading data" << endl;
04347         }
04348 
04349         // Synchronize on whether Proc 0 successfully read the data.
04350         int globalReadDataSuccess = localReadDataSuccess;
04351         broadcast (*comm, 0, outArg (globalReadDataSuccess));
04352         TEUCHOS_TEST_FOR_EXCEPTION(
04353           globalReadDataSuccess == 0, std::runtime_error,
04354           "Failed to read the multivector's data: " << exMsg.str ());
04355 
04356         // If there's only one MPI process and the user didn't supply
04357         // a Map (i.e., pMap is null), we're done.  Set pMap to the
04358         // Map used to distribute X, and return X.
04359         if (comm->getSize () == 1 && map.is_null ()) {
04360           map = proc0Map;
04361           if (! err.is_null ()) {
04362             err->popTab ();
04363           }
04364           if (debug) {
04365             *err << myRank << ": readDenseImpl: done" << endl;
04366           }
04367           if (! err.is_null ()) {
04368             err->popTab ();
04369           }
04370           return X;
04371         }
04372 
04373         if (debug) {
04374           *err << myRank << ": readDenseImpl: Creating target MV" << endl;
04375         }
04376 
04377         // Make a multivector Y with the distributed map pMap.
04378         RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
04379 
04380         if (debug) {
04381           *err << myRank << ": readDenseImpl: Creating Export" << endl;
04382         }
04383 
04384         // Make an Export object that will export X to Y.  First
04385         // argument is the source map, second argument is the target
04386         // map.
04387         Export<LO, GO, NT> exporter (proc0Map, map, err);
04388 
04389         if (debug) {
04390           *err << myRank << ": readDenseImpl: Exporting" << endl;
04391         }
04392         // Export X into Y.
04393         Y->doExport (*X, exporter, INSERT);
04394 
04395         if (! err.is_null ()) {
04396           err->popTab ();
04397         }
04398         if (debug) {
04399           *err << myRank << ": readDenseImpl: done" << endl;
04400         }
04401         if (! err.is_null ()) {
04402           err->popTab ();
04403         }
04404 
04405         // Y is distributed over all process(es) in the communicator.
04406         return Y;
04407       }
04408 
04409     public:
04430       static RCP<const map_type>
04431       readMap (std::istream& in,
04432                const RCP<const comm_type>& comm,
04433                const RCP<node_type>& node,
04434                const bool tolerant=false,
04435                const bool debug=false)
04436       {
04437         Teuchos::RCP<Teuchos::FancyOStream> err =
04438           Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
04439         return readMap (in, comm, node, err, tolerant, debug);
04440       }
04441 
04468       static RCP<const map_type>
04469       readMap (std::istream& in,
04470                const RCP<const comm_type>& comm,
04471                const RCP<node_type>& node,
04472                const Teuchos::RCP<Teuchos::FancyOStream>& err,
04473                const bool tolerant=false,
04474                const bool debug=false)
04475       {
04476         using Teuchos::arcp;
04477         using Teuchos::Array;
04478         using Teuchos::ArrayRCP;
04479         using Teuchos::as;
04480         using Teuchos::broadcast;
04481         using Teuchos::Comm;
04482         using Teuchos::CommRequest;
04483         using Teuchos::inOutArg;
04484         using Teuchos::ireceive;
04485         using Teuchos::outArg;
04486         using Teuchos::receive;
04487         using Teuchos::reduceAll;
04488         using Teuchos::REDUCE_MIN;
04489         using Teuchos::isend;
04490         using Teuchos::SerialComm;
04491         using Teuchos::toString;
04492         using Teuchos::wait;
04493         using std::endl;
04494         typedef Tpetra::global_size_t GST;
04495         typedef ptrdiff_t int_type; // Can hold int and GO
04496         typedef local_ordinal_type LO;
04497         typedef global_ordinal_type GO;
04498         typedef node_type NT;
04499         typedef Tpetra::MultiVector<GO, LO, GO, NT> MV;
04500 
04501         const int numProcs = comm->getSize ();
04502         const int myRank = comm->getRank ();
04503 
04504         if (err.is_null ()) {
04505           err->pushTab ();
04506         }
04507         if (debug) {
04508           std::ostringstream os;
04509           os << myRank << ": readMap: " << endl;
04510           *err << os.str ();
04511         }
04512         if (err.is_null ()) {
04513           err->pushTab ();
04514         }
04515 
04516         // Tag for receive-size / send-size messages.  writeMap used
04517         // tags 1337 and 1338; we count up from there.
04518         const int sizeTag = 1339;
04519         // Tag for receive-data / send-data messages.
04520         const int dataTag = 1340;
04521 
04522         // These are for sends on Process 0, and for receives on all
04523         // other processes.  sizeReq is for the {receive,send}-size
04524         // message, and dataReq is for the message containing the
04525         // actual GIDs to belong to the receiving process.
04526         RCP<CommRequest<int> > sizeReq;
04527         RCP<CommRequest<int> > dataReq;
04528 
04529         // Each process will have to receive the number of GIDs to
04530         // expect.  Thus, we can post the receives now, and cancel
04531         // them if something should go wrong in the meantime.
04532         ArrayRCP<int_type> numGidsToRecv (1);
04533         numGidsToRecv[0] = 0;
04534         if (myRank != 0) {
04535           sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
04536         }
04537 
04538         int readSuccess = 1;
04539         std::ostringstream exMsg;
04540         RCP<MV> data; // Will only be valid on Proc 0
04541         if (myRank == 0) {
04542           // If we want to reuse readDenseImpl, we have to make a
04543           // communicator that only contains Proc 0.  Otherwise,
04544           // readDenseImpl will redistribute the data to all
04545           // processes.  While we eventually want that, neither we nor
04546           // readDenseImpl know the correct Map to use at the moment.
04547           // That depends on the second column of the multivector.
04548           RCP<const Comm<int> > proc0Comm (new SerialComm<int> ());
04549           try {
04550             RCP<const map_type> dataMap;
04551             // This is currently the only place where we use the
04552             // 'tolerant' argument.  Later, if we want to be clever,
04553             // we could have tolerant mode allow PIDs out of order.
04554             data = readDenseImpl<GO> (in, proc0Comm, node, dataMap,
04555                                       err, tolerant, debug);
04556             (void) dataMap; // Silence "unused" warnings
04557             if (data.is_null ()) {
04558               readSuccess = 0;
04559               exMsg << "readDenseImpl() returned null." << endl;
04560             }
04561           } catch (std::exception& e) {
04562             readSuccess = 0;
04563             exMsg << e.what () << endl;
04564           }
04565         }
04566 
04567         // Map from PID to all the GIDs for that PID.
04568         // Only populated on Process 0.
04569         std::map<int, Array<GO> > pid2gids;
04570 
04571         // The index base must be the global minimum GID.
04572         // We will compute this on Process 0 and broadcast,
04573         // so that all processes can set up the Map.
04574         int_type globalNumGIDs = 0;
04575 
04576         // The index base must be the global minimum GID.
04577         // We will compute this on Process 0 and broadcast,
04578         // so that all processes can set up the Map.
04579         GO indexBase = 0;
04580 
04581         // Process 0: If the above read of the MultiVector succeeded,
04582         // extract the GIDs and PIDs into pid2gids, and find the
04583         // global min GID.
04584         if (myRank == 0 && readSuccess == 1) {
04585           if (data->getNumVectors () == 2) { // Map format 1.0
04586             ArrayRCP<const GO> GIDs = data->getData (0);
04587             ArrayRCP<const GO> PIDs = data->getData (1); // convert to int
04588             globalNumGIDs = GIDs.size ();
04589 
04590             // Start computing the global min GID, while collecting
04591             // the GIDs for each PID.
04592             if (globalNumGIDs > 0) {
04593               const int pid = static_cast<int> (PIDs[0]);
04594 
04595               if (pid < 0 || pid >= numProcs) {
04596                 readSuccess = 0;
04597                 exMsg << "Tpetra::MatrixMarket::readMap: "
04598                       << "Encountered invalid PID " << pid << "." << endl;
04599               }
04600               else {
04601                 const GO gid = GIDs[0];
04602                 pid2gids[pid].push_back (gid);
04603                 indexBase = gid; // the current min GID
04604               }
04605             }
04606             if (readSuccess == 1) {
04607               // Collect the rest of the GIDs for each PID, and compute
04608               // the global min GID.
04609               for (size_type k = 1; k < globalNumGIDs; ++k) {
04610                 const int pid = static_cast<int> (PIDs[k]);
04611                 if (pid < 0 || pid >= numProcs) {
04612                   readSuccess = 0;
04613                   exMsg << "Tpetra::MatrixMarket::readMap: "
04614                         << "Encountered invalid PID " << pid << "." << endl;
04615                 }
04616                 else {
04617                   const int_type gid = GIDs[k];
04618                   pid2gids[pid].push_back (gid);
04619                   if (gid < indexBase) {
04620                     indexBase = gid; // the current min GID
04621                   }
04622                 }
04623               }
04624             }
04625           }
04626           else if (data->getNumVectors () == 1) { // Map format 2.0
04627             if (data->getGlobalLength () % 2 != static_cast<GST> (0)) {
04628               readSuccess = 0;
04629               exMsg << "Tpetra::MatrixMarket::readMap: Input data has the "
04630                 "wrong format (for Map format 2.0).  The global number of rows "
04631                 "in the MultiVector must be even (divisible by 2)." << endl;
04632             }
04633             else {
04634               ArrayRCP<const GO> theData = data->getData (0);
04635               globalNumGIDs = static_cast<GO> (data->getGlobalLength ()) /
04636                 static_cast<GO> (2);
04637 
04638               // Start computing the global min GID, while
04639               // collecting the GIDs for each PID.
04640               if (globalNumGIDs > 0) {
04641                 const int pid = static_cast<int> (theData[1]);
04642                 if (pid < 0 || pid >= numProcs) {
04643                   readSuccess = 0;
04644                   exMsg << "Tpetra::MatrixMarket::readMap: "
04645                         << "Encountered invalid PID " << pid << "." << endl;
04646                 }
04647                 else {
04648                   const GO gid = theData[0];
04649                   pid2gids[pid].push_back (gid);
04650                   indexBase = gid; // the current min GID
04651                 }
04652               }
04653               // Collect the rest of the GIDs for each PID, and
04654               // compute the global min GID.
04655               for (int_type k = 1; k < globalNumGIDs; ++k) {
04656                 const int pid = static_cast<int> (theData[2*k + 1]);
04657                 if (pid < 0 || pid >= numProcs) {
04658                   readSuccess = 0;
04659                   exMsg << "Tpetra::MatrixMarket::readMap: "
04660                         << "Encountered invalid PID " << pid << "." << endl;
04661                 }
04662                 else {
04663                   const GO gid = theData[2*k];
04664                   pid2gids[pid].push_back (gid);
04665                   if (gid < indexBase) {
04666                     indexBase = gid; // the current min GID
04667                   }
04668                 }
04669               } // for each GID
04670             } // if the amount of data is correct
04671           }
04672           else {
04673             readSuccess = 0;
04674             exMsg << "Tpetra::MatrixMarket::readMap: Input data must have "
04675               "either 1 column (for the new Map format 2.0) or 2 columns (for "
04676               "the old Map format 1.0).";
04677           }
04678         } // myRank is zero
04679 
04680         // Broadcast the indexBase, the global number of GIDs, and the
04681         // current success status.  Use int_type for all of these.
04682         {
04683           int_type readResults[3];
04684           readResults[0] = static_cast<int_type> (indexBase);
04685           readResults[1] = static_cast<int_type> (globalNumGIDs);
04686           readResults[2] = static_cast<int_type> (readSuccess);
04687           broadcast<int, int_type> (*comm, 0, 3, readResults);
04688 
04689           indexBase = static_cast<GO> (readResults[0]);
04690           globalNumGIDs = static_cast<int_type> (readResults[1]);
04691           readSuccess = static_cast<int> (readResults[2]);
04692         }
04693 
04694         // Unwinding the stack will invoke sizeReq's destructor, which
04695         // will cancel the receive-size request on all processes that
04696         // posted it.
04697         TEUCHOS_TEST_FOR_EXCEPTION(
04698           readSuccess != 1, std::runtime_error,
04699           "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
04700           "following exception message: " << exMsg.str ());
04701 
04702         if (myRank == 0) {
04703           // Proc 0: Send each process' number of GIDs to that process.
04704           for (int p = 1; p < numProcs; ++p) {
04705             ArrayRCP<int_type> numGidsToSend (1);
04706 
04707             typename std::map<int, Array<GO> >::const_iterator it = pid2gids.find (p);
04708             if (it == pid2gids.end ()) {
04709               numGidsToSend[0] = 0;
04710             } else {
04711               numGidsToSend[0] = it->second.size ();
04712             }
04713             sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
04714             wait<int> (*comm, outArg (sizeReq));
04715           }
04716         }
04717         else {
04718           // Wait on the receive-size message to finish.
04719           wait<int> (*comm, outArg (sizeReq));
04720         }
04721 
04722         // Allocate / get the array for my GIDs.
04723         // Only Process 0 will have its actual GIDs at this point.
04724         ArrayRCP<GO> myGids;
04725         int_type myNumGids = 0;
04726         if (myRank == 0) {
04727           GO* myGidsRaw = NULL;
04728 
04729           typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
04730           if (it != pid2gids.end ()) {
04731             myGidsRaw = it->second.getRawPtr ();
04732             myNumGids = it->second.size ();
04733             // Nonowning ArrayRCP just views the Array.
04734             myGids = arcp<GO> (myGidsRaw, 0, myNumGids, false);
04735           }
04736         }
04737         else { // myRank != 0
04738           myNumGids = numGidsToRecv[0];
04739           myGids = arcp<GO> (myNumGids);
04740         }
04741 
04742         if (myRank != 0) {
04743           // Post receive for data, now that we know how much data we
04744           // will receive.  Only post receive if my process actually
04745           // has nonzero GIDs.
04746           if (myNumGids > 0) {
04747             dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
04748           }
04749         }
04750 
04751         for (int p = 1; p < numProcs; ++p) {
04752           if (myRank == 0) {
04753             ArrayRCP<GO> sendGids; // to send to Process p
04754             GO* sendGidsRaw = NULL;
04755             int_type numSendGids = 0;
04756 
04757             typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
04758             if (it != pid2gids.end ()) {
04759               numSendGids = it->second.size ();
04760               sendGidsRaw = it->second.getRawPtr ();
04761               sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids, false);
04762             }
04763             // Only send if that process actually has nonzero GIDs.
04764             if (numSendGids > 0) {
04765               dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
04766             }
04767             wait<int> (*comm, outArg (dataReq));
04768           }
04769           else if (myRank == p) {
04770             // Wait on my receive of GIDs to finish.
04771             wait<int> (*comm, outArg (dataReq));
04772           }
04773         } // for each process rank p in 1, 2, ..., numProcs-1
04774 
04775         if (debug) {
04776           std::ostringstream os;
04777           os << myRank << ": readMap: creating Map" << endl;
04778           *err << os.str ();
04779         }
04780         const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
04781         RCP<const map_type> newMap =
04782           rcp (new map_type (INVALID, myGids (), indexBase, comm, node));
04783 
04784         if (err.is_null ()) {
04785           err->popTab ();
04786         }
04787         if (debug) {
04788           std::ostringstream os;
04789           os << myRank << ": readMap: done" << endl;
04790           *err << os.str ();
04791         }
04792         if (err.is_null ()) {
04793           err->popTab ();
04794         }
04795         return newMap;
04796       }
04797 
04798     private:
04799 
04810       static int
04811       encodeDataType (const std::string& dataType)
04812       {
04813         if (dataType == "real" || dataType == "integer") {
04814           return 0;
04815         } else if (dataType == "complex") {
04816           return 1;
04817         } else if (dataType == "pattern") {
04818           return 2;
04819         } else {
04820           // We should never get here, since Banner validates the
04821           // reported data type and ensures it is one of the accepted
04822           // values.
04823           TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
04824             "Unrecognized Matrix Market data type \"" << dataType
04825             << "\".  We should never get here.  "
04826             "Please report this bug to the Tpetra developers.");
04827         }
04828       }
04829     };
04830 
04870     template<class SparseMatrixType>
04871     class Writer {
04872     public:
04873       typedef SparseMatrixType sparse_matrix_type;
04874       typedef RCP<sparse_matrix_type> sparse_matrix_ptr;
04875 
04878       typedef typename SparseMatrixType::scalar_type scalar_type;
04881       typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
04888       typedef typename SparseMatrixType::global_ordinal_type global_ordinal_type;
04891       typedef typename SparseMatrixType::node_type node_type;
04892 
04895       typedef MultiVector<scalar_type,
04896                           local_ordinal_type,
04897                           global_ordinal_type,
04898                           node_type> multivector_type;
04901       typedef Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
04902 
04934       static void
04935       writeSparseFile (const std::string& filename,
04936                        const RCP<const sparse_matrix_type>& pMatrix,
04937                        const std::string& matrixName,
04938                        const std::string& matrixDescription,
04939                        const bool debug=false)
04940       {
04941         const int myRank = Teuchos::rank (*(pMatrix->getComm()));
04942         std::ofstream out;
04943 
04944         // Only open the file on Rank 0.
04945         if (myRank == 0) out.open (filename.c_str());
04946         writeSparse (out, pMatrix, matrixName, matrixDescription, debug);
04947         // We can rely on the destructor of the output stream to close
04948         // the file on scope exit, even if writeSparse() throws an
04949         // exception.
04950       }
04951 
04972       static void
04973       writeSparseFile (const std::string& filename,
04974                        const RCP<const sparse_matrix_type>& pMatrix,
04975                        const bool debug=false)
04976       {
04977         writeSparseFile (filename, pMatrix, "", "", debug);
04978       }
04979 
04980     public:
04981 
05012       static void
05013       writeSparse (std::ostream& out,
05014                    const RCP<const sparse_matrix_type>& pMatrix,
05015                    const std::string& matrixName,
05016                    const std::string& matrixDescription,
05017                    const bool debug=false)
05018       {
05019         using Teuchos::Comm;
05020         using Teuchos::FancyOStream;
05021         using Teuchos::getFancyOStream;
05022         using Teuchos::null;
05023         using Teuchos::RCP;
05024         using Teuchos::rcpFromRef;
05025         using std::cerr;
05026         using std::endl;
05027         typedef scalar_type ST;
05028         typedef local_ordinal_type LO;
05029         typedef global_ordinal_type GO;
05030         typedef typename Teuchos::ScalarTraits<ST> STS;
05031         typedef typename ArrayView<const LO>::const_iterator lo_iter;
05032         typedef typename ArrayView<const GO>::const_iterator go_iter;
05033         typedef typename ArrayView<const ST>::const_iterator st_iter;
05034 
05035         // Make the output stream write floating-point numbers in
05036         // scientific notation.  It will politely put the output
05037         // stream back to its state on input, when this scope
05038         // terminates.
05039         Teuchos::MatrixMarket::details::SetScientific<ST> sci (out);
05040 
05041         RCP<const Comm<int> > comm = pMatrix->getComm();
05042         const int myRank = comm->getRank ();
05043 
05044         RCP<FancyOStream> err =
05045           debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
05046         if (debug) {
05047           std::ostringstream os;
05048           os << myRank << ": writeSparse" << endl;
05049           *err << os.str ();
05050           comm->barrier ();
05051           os << "-- " << myRank << ": past barrier" << endl;
05052           *err << os.str ();
05053         }
05054 
05055         // Whether to print debugging output to stderr.
05056         const bool debugPrint = debug && myRank == 0;
05057 
05058         RCP<const map_type> rowMap = pMatrix->getRowMap ();
05059         RCP<const map_type> colMap = pMatrix->getColMap ();
05060         RCP<const map_type> domainMap = pMatrix->getDomainMap ();
05061         RCP<const map_type> rangeMap = pMatrix->getRangeMap ();
05062 
05063         const global_size_t numRows = rangeMap->getGlobalNumElements ();
05064         const global_size_t numCols = domainMap->getGlobalNumElements ();
05065 
05066         if (debug && myRank == 0) {
05067           std::ostringstream os;
05068           os << "-- Input sparse matrix is:"
05069              << "---- " << numRows << " x " << numCols << " with "
05070              << pMatrix->getGlobalNumEntries() << " entries;" << endl
05071              << "---- "
05072              << (pMatrix->isGloballyIndexed() ? "Globally" : "Locally")
05073              << " indexed." << endl
05074              << "---- Its row map has " << rowMap->getGlobalNumElements ()
05075              << " elements." << endl
05076              << "---- Its col map has " << colMap->getGlobalNumElements ()
05077              << " elements." << endl;
05078           *err << os.str ();
05079         }
05080         // Make the "gather" row map, where Proc 0 owns all rows and
05081         // the other procs own no rows.
05082         const size_t localNumRows = (myRank == 0) ? numRows : 0;
05083         RCP<node_type> node = rowMap->getNode();
05084         if (debug) {
05085           std::ostringstream os;
05086           os << "-- " << myRank << ": making gatherRowMap" << endl;
05087           *err << os.str ();
05088         }
05089         RCP<const map_type> gatherRowMap =
05090           Details::computeGatherMap (rowMap, err, debug);
05091 
05092         // Since the matrix may in general be non-square, we need to
05093         // make a column map as well.  In this case, the column map
05094         // contains all the columns of the original matrix, because we
05095         // are gathering the whole matrix onto Proc 0.  We call
05096         // computeGatherMap to preserve the original order of column
05097         // indices over all the processes.
05098         const size_t localNumCols = (myRank == 0) ? numCols : 0;
05099         RCP<const map_type> gatherColMap =
05100           Details::computeGatherMap (colMap, err, debug);
05101 
05102         // Current map is the source map, gather map is the target map.
05103         typedef Import<LO, GO, node_type> import_type;
05104         import_type importer (rowMap, gatherRowMap);
05105 
05106         // Create a new CrsMatrix to hold the result of the import.
05107         // The constructor needs a column map as well as a row map,
05108         // for the case that the matrix is not square.
05109         RCP<sparse_matrix_type> newMatrix =
05110           rcp (new sparse_matrix_type (gatherRowMap, gatherColMap,
05111                                        static_cast<size_t> (0)));
05112         // Import the sparse matrix onto Proc 0.
05113         newMatrix->doImport (*pMatrix, importer, INSERT);
05114 
05115         // fillComplete() needs the domain and range maps for the case
05116         // that the matrix is not square.
05117         {
05118           RCP<const map_type> gatherDomainMap =
05119             rcp (new map_type (numCols, localNumCols,
05120                                domainMap->getIndexBase (),
05121                                comm, node));
05122           RCP<const map_type> gatherRangeMap =
05123             rcp (new map_type (numRows, localNumRows,
05124                                rangeMap->getIndexBase (),
05125                                comm, node));
05126           newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
05127         }
05128 
05129         if (debugPrint) {
05130           cerr << "-- Output sparse matrix is:"
05131                << "---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
05132                << " x "
05133                << newMatrix->getDomainMap ()->getGlobalNumElements ()
05134                << " with "
05135                << newMatrix->getGlobalNumEntries () << " entries;" << endl
05136                << "---- "
05137                << (newMatrix->isGloballyIndexed () ? "Globally" : "Locally")
05138                << " indexed." << endl
05139                << "---- Its row map has "
05140                << newMatrix->getRowMap ()->getGlobalNumElements ()
05141                << " elements, with index base "
05142                << newMatrix->getRowMap ()->getIndexBase () << "." << endl
05143                << "---- Its col map has "
05144                << newMatrix->getColMap ()->getGlobalNumElements ()
05145                << " elements, with index base "
05146                << newMatrix->getColMap ()->getIndexBase () << "." << endl
05147                << "---- Element count of output matrix's column Map may differ "
05148                << "from that of the input matrix's column Map, if some columns "
05149                << "of the matrix contain no entries." << endl;
05150         }
05151 
05152         //
05153         // Print the metadata and the matrix entries on Rank 0.
05154         //
05155         if (myRank == 0) {
05156           // Print the Matrix Market banner line.  CrsMatrix stores
05157           // data nonsymmetrically ("general").  This implies that
05158           // readSparse() on a symmetrically stored input file,
05159           // followed by writeSparse() on the resulting sparse matrix,
05160           // will result in an output file with a different banner
05161           // line than the original input file.
05162           out << "%%MatrixMarket matrix coordinate "
05163               << (STS::isComplex ? "complex" : "real")
05164               << " general" << endl;
05165 
05166           // Print comments (the matrix name and / or description).
05167           if (matrixName != "") {
05168             printAsComment (out, matrixName);
05169           }
05170           if (matrixDescription != "") {
05171             printAsComment (out, matrixDescription);
05172           }
05173 
05174           // Print the Matrix Market header (# rows, # columns, #
05175           // nonzeros).  Use the range resp. domain map for the number
05176           // of rows resp. columns, since Tpetra::CrsMatrix uses the
05177           // column map for the number of columns.  That only
05178           // corresponds to the "linear-algebraic" number of columns
05179           // when the column map is uniquely owned (a.k.a. one-to-one),
05180           // which only happens if the matrix is (block) diagonal.
05181           out << newMatrix->getRangeMap ()->getGlobalNumElements () << " "
05182               << newMatrix->getDomainMap ()->getGlobalNumElements () << " "
05183               << newMatrix->getGlobalNumEntries () << endl;
05184 
05185           // The Matrix Market format expects one-based row and column
05186           // indices.  We'll convert the indices on output from
05187           // whatever index base they use to one-based indices.
05188           const GO rowIndexBase = gatherRowMap->getIndexBase ();
05189           const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
05190           //
05191           // Print the entries of the matrix.
05192           //
05193           // newMatrix can never be globally indexed, since we called
05194           // fillComplete() on it.  We include code for both cases
05195           // (globally or locally indexed) just in case that ever
05196           // changes.
05197           if (newMatrix->isGloballyIndexed()) {
05198             // We know that the "gather" row Map is contiguous, so we
05199             // don't need to get the list of GIDs.
05200             const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
05201             const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
05202             for (GO globalRowIndex = minAllGlobalIndex;
05203                  globalRowIndex <= maxAllGlobalIndex; // inclusive range
05204                  ++globalRowIndex) {
05205               ArrayView<const GO> ind;
05206               ArrayView<const ST> val;
05207               newMatrix->getGlobalRowView (globalRowIndex, ind, val);
05208               go_iter indIter = ind.begin ();
05209               st_iter valIter = val.begin ();
05210               for (; indIter != ind.end() && valIter != val.end();
05211                    ++indIter, ++valIter) {
05212                 const GO globalColIndex = *indIter;
05213                 // Convert row and column indices to 1-based.
05214                 // This works because the global index type is signed.
05215                 out << (globalRowIndex + 1 - rowIndexBase) << " "
05216                     << (globalColIndex + 1 - colIndexBase) << " ";
05217                 if (STS::isComplex) {
05218                   out << STS::real (*valIter) << " " << STS::imag (*valIter);
05219                 } else {
05220                   out << *valIter;
05221                 }
05222                 out << endl;
05223               } // For each entry in the current row
05224             } // For each row of the "gather" matrix
05225           } else { // newMatrix is locally indexed
05226             typedef OrdinalTraits<GO> OTG;
05227             for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
05228                  localRowIndex <= gatherRowMap->getMaxLocalIndex();
05229                  ++localRowIndex) {
05230               // Convert from local to global row index.
05231               const GO globalRowIndex =
05232                 gatherRowMap->getGlobalElement (localRowIndex);
05233               TEUCHOS_TEST_FOR_EXCEPTION(
05234                 globalRowIndex == OTG::invalid(), std::logic_error,
05235                 "Failed to convert the supposed local row index "
05236                 << localRowIndex << " into a global row index.  "
05237                 "Please report this bug to the Tpetra developers.");
05238               ArrayView<const LO> ind;
05239               ArrayView<const ST> val;
05240               newMatrix->getLocalRowView (localRowIndex, ind, val);
05241               lo_iter indIter = ind.begin ();
05242               st_iter valIter = val.begin ();
05243               for (; indIter != ind.end() && valIter != val.end();
05244                    ++indIter, ++valIter) {
05245                 // Convert the column index from local to global.
05246                 const GO globalColIndex =
05247                   newMatrix->getColMap()->getGlobalElement (*indIter);
05248                 TEUCHOS_TEST_FOR_EXCEPTION(
05249                   globalColIndex == OTG::invalid(), std::logic_error,
05250                   "On local row " << localRowIndex << " of the sparse matrix: "
05251                   "Failed to convert the supposed local column index "
05252                   << *indIter << " into a global column index.  Please report "
05253                   "this bug to the Tpetra developers.");
05254                 // Convert row and column indices to 1-based.
05255                 // This works because the global index type is signed.
05256                 out << (globalRowIndex + 1 - rowIndexBase) << " "
05257                     << (globalColIndex + 1 - colIndexBase) << " ";
05258                 if (STS::isComplex) {
05259                   out << STS::real (*valIter) << " " << STS::imag (*valIter);
05260                 } else {
05261                   out << *valIter;
05262                 }
05263                 out << endl;
05264               } // For each entry in the current row
05265             } // For each row of the "gather" matrix
05266           } // Whether the "gather" matrix is locally or globally indexed
05267         } // If my process' rank is 0
05268       }
05269 
05292       static void
05293       writeSparse (std::ostream& out,
05294                    const RCP<const sparse_matrix_type>& pMatrix,
05295                    const bool debug=false)
05296       {
05297         writeSparse (out, pMatrix, "", "", debug);
05298       }
05299 
05326       static void
05327       writeDenseFile (const std::string& filename,
05328                       const RCP<const multivector_type>& X,
05329                       const std::string& matrixName,
05330                       const std::string& matrixDescription)
05331       {
05332         const int myRank = Teuchos::rank (*(X->getMap()->getComm()));
05333         std::ofstream out;
05334 
05335         if (myRank == 0) // Only open the file on Rank 0.
05336           out.open (filename.c_str());
05337 
05338         writeDense (out, X, matrixName, matrixDescription);
05339         // We can rely on the destructor of the output stream to close
05340         // the file on scope exit, even if writeDense() throws an
05341         // exception.
05342       }
05343 
05362       static void
05363       writeDenseFile (const std::string& filename,
05364                       const RCP<const multivector_type>& X)
05365       {
05366         writeDenseFile (filename, X, "", "");
05367       }
05368 
05395       static void
05396       writeDense (std::ostream& out,
05397                   const RCP<const multivector_type>& X,
05398                   const std::string& matrixName,
05399                   const std::string& matrixDescription)
05400       {
05401         writeDenseImpl<scalar_type> (out, *X, matrixName, matrixDescription);
05402       }
05403 
05422       static void
05423       writeDense (std::ostream& out,
05424                   const RCP<const multivector_type>& X)
05425       {
05426         writeDense (out, X, "", "");
05427       }
05428 
05448       static void
05449       writeMap (std::ostream& out, const map_type& map, const bool debug=false)
05450       {
05451         Teuchos::RCP<Teuchos::FancyOStream> err =
05452           Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
05453         writeMap (out, map, err, debug);
05454       }
05455 
05464       static void
05465       writeMap (std::ostream& out,
05466                 const map_type& map,
05467                 const Teuchos::RCP<Teuchos::FancyOStream>& err,
05468                 const bool debug=false)
05469       {
05470         using Teuchos::ArrayRCP;
05471         using Teuchos::ArrayView;
05472         using Teuchos::CommRequest;
05473         using Teuchos::ireceive;
05474         using Teuchos::isend;
05475         using Teuchos::RCP;
05476         using Teuchos::TypeNameTraits;
05477         using Teuchos::wait;
05478         using std::endl;
05479         typedef global_ordinal_type GO;
05480         typedef int pid_type;
05481 
05482         // Treat the Map as a 1-column "multivector."  This differs
05483         // from the previous two-column format, in which column 0 held
05484         // the GIDs, and column 1 held the corresponding PIDs.  It
05485         // differs because printing that format requires knowing the
05486         // entire first column -- that is, all the GIDs -- in advance.
05487         // Sending messages from each process one at a time saves
05488         // memory, but it means that Process 0 doesn't ever have all
05489         // the GIDs at once.
05490         //
05491         // We pack the entries as ptrdiff_t, since this should be the
05492         // biggest signed built-in integer type that can hold any GO
05493         // or pid_type (= int) quantity without overflow.  Test this
05494         // assumption at run time.
05495         typedef ptrdiff_t int_type;
05496         TEUCHOS_TEST_FOR_EXCEPTION(
05497           sizeof (GO) > sizeof (int_type), std::logic_error,
05498           "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
05499           << " is too big for ptrdiff_t.  sizeof(GO) = " << sizeof (GO)
05500           << " > sizeof(ptrdiff_t) = " << sizeof (ptrdiff_t) << ".");
05501         TEUCHOS_TEST_FOR_EXCEPTION(
05502           sizeof (pid_type) > sizeof (int_type), std::logic_error,
05503           "The (MPI) process rank type pid_type=" <<
05504           TypeNameTraits<pid_type>::name () << " is too big for ptrdiff_t.  "
05505           "sizeof(pid_type) = " << sizeof (pid_type) << " > sizeof(ptrdiff_t)"
05506           " = " << sizeof (ptrdiff_t) << ".");
05507 
05508         const Comm<int>& comm = * (map.getComm ());
05509         const int myRank = comm.getRank ();
05510         const int numProcs = comm.getSize ();
05511 
05512         if (! err.is_null ()) {
05513           err->pushTab ();
05514         }
05515         if (debug) {
05516           std::ostringstream os;
05517           os << myRank << ": writeMap" << endl;
05518           *err << os.str ();
05519         }
05520         if (! err.is_null ()) {
05521           err->pushTab ();
05522         }
05523 
05524         const size_t myNumRows = map.getNodeNumElements ();
05525         // Use a different tag for the "size" messages than for the
05526         // "data" messages, in order to help us debug any mix-ups.
05527         const int sizeTag = 1337;
05528         const int dataTag = 1338;
05529 
05530         // Process 0 pipelines nonblocking receives with file output.
05531         //
05532         // Constraints:
05533         //   - Process 0 can't post a receive for another process'
05534         //     actual data, until it posts and waits on the receive
05535         //     from that process with the amount of data to receive.
05536         //     (We could just post receives with a max data size, but
05537         //     I feel uncomfortable about that.)
05538         //   - The C++ standard library doesn't allow nonblocking
05539         //     output to an std::ostream.
05540         //
05541         // Process 0: Post receive-size receives from Processes 1 and 2.
05542         // Process 1: Post send-size send to Process 0.
05543         // Process 2: Post send-size send to Process 0.
05544         //
05545         // All processes: Pack my GIDs and PIDs.
05546         //
05547         // Process 1:
05548         //   - Post send-data send to Process 0.
05549         //   - Wait on my send-size send to Process 0.
05550         //
05551         // Process 0:
05552         //   - Print MatrixMarket header.
05553         //   - Print my GIDs and PIDs.
05554         //   - Wait on receive-size receive from Process 1.
05555         //   - Post receive-data receive from Process 1.
05556         //
05557         // For each process p = 1, 2, ... numProcs-1:
05558         //   If I am Process 0:
05559         //     - Post receive-size receive from Process p + 2
05560         //     - Wait on receive-size receive from Process p + 1
05561         //     - Post receive-data receive from Process p + 1
05562         //     - Wait on receive-data receive from Process p
05563         //     - Write data from Process p.
05564         //   Else if I am Process p:
05565         //     - Wait on my send-data send.
05566         //   Else if I am Process p+1:
05567         //     - Post send-data send to Process 0.
05568         //     - Wait on my send-size send.
05569         //   Else if I am Process p+2:
05570         //     - Post send-size send to Process 0.
05571         //
05572         // Pipelining has three goals here:
05573         //   1. Overlap communication (the receives) with file I/O
05574         //   2. Give Process 0 a chance to prepost some receives,
05575         //      before sends show up, by packing local data before
05576         //      posting sends
05577         //   3. Don't post _all_ receives or _all_ sends, because that
05578         //      wouldn't be memory scalable.  (Just because we can't
05579         //      see how much memory MPI consumes, doesn't mean that it
05580         //      doesn't consume any!)
05581 
05582         // These are used on every process.  sendReqSize[0] holds the
05583         // number of rows on this process, and sendReqBuf holds this
05584         // process' data.  Process 0 packs into sendReqBuf, but
05585         // doesn't send; it only uses that for printing.  All other
05586         // processes send both of these to Process 0.
05587         RCP<CommRequest<int> > sendReqSize, sendReqData;
05588 
05589         // These are used only on Process 0, for received data.  Keep
05590         // 3 of each, and treat the arrays as circular buffers.  When
05591         // receiving from Process p, the corresponding array index
05592         // here is p % 3.
05593         Array<ArrayRCP<int_type> > recvSizeBufs (3);
05594         Array<ArrayRCP<int_type> > recvDataBufs (3);
05595         Array<RCP<CommRequest<int> > > recvSizeReqs (3);
05596         Array<RCP<CommRequest<int> > > recvDataReqs (3);
05597 
05598         // Buffer for nonblocking send of the "send size."
05599         ArrayRCP<int_type> sendDataSize (1);
05600         sendDataSize[0] = myNumRows;
05601 
05602         if (myRank == 0) {
05603           if (debug) {
05604             std::ostringstream os;
05605             os << myRank << ": Post receive-size receives from "
05606               "Procs 1 and 2: tag = " << sizeTag << endl;
05607             *err << os.str ();
05608           }
05609           // Process 0: Post receive-size receives from Processes 1 and 2.
05610           recvSizeBufs[0].resize (1);
05611           (recvSizeBufs[0])[0] = -1; // error flag
05612           recvSizeBufs[1].resize (1);
05613           (recvSizeBufs[1])[0] = -1; // error flag
05614           recvSizeBufs[2].resize (1);
05615           (recvSizeBufs[2])[0] = -1; // error flag
05616           if (numProcs > 1) {
05617             recvSizeReqs[1] =
05618               ireceive<int, int_type> (recvSizeBufs[1], 1, sizeTag, comm);
05619           }
05620           if (numProcs > 2) {
05621             recvSizeReqs[2] =
05622               ireceive<int, int_type> (recvSizeBufs[2], 2, sizeTag, comm);
05623           }
05624         }
05625         else if (myRank == 1 || myRank == 2) {
05626           if (debug) {
05627             std::ostringstream os;
05628             os << myRank << ": Post send-size send: size = "
05629                << sendDataSize[0] << ", tag = " << sizeTag << endl;
05630             *err << os.str ();
05631           }
05632           // Prime the pipeline by having Processes 1 and 2 start
05633           // their send-size sends.  We don't want _all_ the processes
05634           // to start their send-size sends, because that wouldn't be
05635           // memory scalable.
05636           sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
05637         }
05638         else {
05639           if (debug) {
05640             std::ostringstream os;
05641             os << myRank << ": Not posting my send-size send yet" << endl;
05642             *err << os.str ();
05643           }
05644         }
05645 
05646         //
05647         // Pack my GIDs and PIDs.  Each (GID,PID) pair gets packed
05648         // consecutively, for better locality.
05649         //
05650 
05651         if (debug) {
05652           std::ostringstream os;
05653           os << myRank << ": Pack my GIDs and PIDs" << endl;
05654           *err << os.str ();
05655         }
05656 
05657         ArrayRCP<int_type> sendDataBuf (myNumRows * 2);
05658 
05659         if (map.isContiguous ()) {
05660           const int_type myMinGblIdx =
05661             static_cast<int_type> (map.getMinGlobalIndex ());
05662           for (size_t k = 0; k < myNumRows; ++k) {
05663             const int_type gid = myMinGblIdx + static_cast<int_type> (k);
05664             const int_type pid = static_cast<int_type> (myRank);
05665             sendDataBuf[2*k] = gid;
05666             sendDataBuf[2*k+1] = pid;
05667           }
05668         }
05669         else {
05670           ArrayView<const GO> myGblInds = map.getNodeElementList ();
05671           for (size_t k = 0; k < myNumRows; ++k) {
05672             const int_type gid = static_cast<int_type> (myGblInds[k]);
05673             const int_type pid = static_cast<int_type> (myRank);
05674             sendDataBuf[2*k] = gid;
05675             sendDataBuf[2*k+1] = pid;
05676           }
05677         }
05678 
05679         if (debug) {
05680           std::ostringstream os;
05681           os << myRank << ": Done packing my GIDs and PIDs" << endl;
05682           *err << os.str ();
05683         }
05684 
05685         if (myRank == 1) {
05686           // Process 1: post send-data send to Process 0.
05687           if (debug) {
05688             *err << myRank << ": Post send-data send: tag = " << dataTag
05689                  << endl;
05690           }
05691           sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
05692         }
05693 
05694         if (myRank == 0) {
05695           if (debug) {
05696             *err << myRank << ": Write MatrixMarket header" << endl;
05697           }
05698 
05699           // Process 0: Write the MatrixMarket header.
05700           // Description section explains each column.
05701           std::ostringstream hdr;
05702 
05703           // Print the Matrix Market header.  MultiVector stores data
05704           // nonsymmetrically, hence "general" in the banner line.
05705           hdr << "%%MatrixMarket matrix array integer general" << endl
05706               << "% Format: Version 2.0" << endl
05707               << "%" << endl
05708               << "% This file encodes a Tpetra::Map." << endl
05709               << "% It is stored as a dense vector, with twice as many " << endl
05710               << "% entries as the global number of GIDs (global indices)." << endl
05711               << "% (GID, PID) pairs are stored contiguously, where the PID " << endl
05712               << "% is the rank of the process owning that GID." << endl
05713               << (2 * map.getGlobalNumElements ()) << " " << 1 << endl;
05714           out << hdr.str ();
05715 
05716           if (debug) {
05717             std::ostringstream os;
05718             os << myRank << ": Write my GIDs and PIDs" << endl;
05719             *err << os.str ();
05720           }
05721 
05722           // Write Process 0's data to the output stream.
05723           // Matrix Market prints dense matrices in column-major order.
05724           const int_type printNumRows = myNumRows;
05725           ArrayView<const int_type> printData = sendDataBuf ();
05726           for (int_type k = 0; k < printNumRows; ++k) {
05727             const int_type gid = printData[2*k];
05728             const int_type pid = printData[2*k+1];
05729             out << gid << endl << pid << endl;
05730           }
05731         }
05732 
05733         if (myRank == 0) {
05734           // Wait on receive-size receive from Process 1.
05735           const int recvRank = 1;
05736           const int circBufInd = recvRank % 3;
05737           if (debug) {
05738             std::ostringstream os;
05739             os << myRank << ": Wait on receive-size receive from Process "
05740                << recvRank << endl;
05741             *err << os.str ();
05742           }
05743           if (numProcs > 1) {
05744             wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
05745 
05746             // We received the number of rows of data.  (The data
05747             // come in two columns.)
05748             const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
05749             if (debug && recvNumRows == -1) {
05750               std::ostringstream os;
05751               os << myRank << ": Result of receive-size receive from Process "
05752                  << recvRank << " is -1.  This should never happen, and "
05753                 "suggests that the receive never got posted.  Please report "
05754                 "this bug to the Tpetra developers." << endl;
05755               *err << os.str ();
05756             }
05757 
05758             // Post receive-data receive from Process 1.
05759             recvDataBufs[circBufInd].resize (recvNumRows * 2);
05760             if (debug) {
05761               std::ostringstream os;
05762               os << myRank << ": Post receive-data receive from Process "
05763                  << recvRank << ": tag = " << dataTag << ", buffer size = "
05764                  << recvDataBufs[circBufInd].size () << endl;
05765               *err << os.str ();
05766             }
05767             if (! recvSizeReqs[circBufInd].is_null ()) {
05768               std::ostringstream os;
05769               os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
05770                 "null, before posting the receive-data receive from Process "
05771                  << recvRank << ".  This should never happen.  Please report "
05772                 "this bug to the Tpetra developers." << endl;
05773               *err << os.str ();
05774             }
05775             recvDataReqs[circBufInd] =
05776               ireceive<int, int_type> (recvDataBufs[circBufInd],
05777                                        recvRank, dataTag, comm);
05778           } // numProcs > 1
05779         }
05780         else if (myRank == 1) {
05781           // Wait on my send-size send.
05782           if (debug) {
05783             std::ostringstream os;
05784             os << myRank << ": Wait on my send-size send" << endl;
05785             *err << os.str ();
05786           }
05787           wait<int> (comm, outArg (sendReqSize));
05788         }
05789 
05790         //
05791         // Pipeline loop
05792         //
05793         for (int p = 1; p < numProcs; ++p) {
05794           if (myRank == 0) {
05795             if (p + 2 < numProcs) {
05796               // Post receive-size receive from Process p + 2.
05797               const int recvRank = p + 2;
05798               const int circBufInd = recvRank % 3;
05799               if (debug) {
05800                 std::ostringstream os;
05801                 os << myRank << ": Post receive-size receive from Process "
05802                    << recvRank << ": tag = " << sizeTag << endl;
05803                 *err << os.str ();
05804               }
05805               if (! recvSizeReqs[circBufInd].is_null ()) {
05806                 std::ostringstream os;
05807                 os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
05808                    << "null, for the receive-size receive from Process "
05809                    << recvRank << "!  This may mean that this process never "
05810                    << "finished waiting for the receive from Process "
05811                    << (recvRank - 3) << "." << endl;
05812                 *err << os.str ();
05813               }
05814               recvSizeReqs[circBufInd] =
05815                 ireceive<int, int_type> (recvSizeBufs[circBufInd],
05816                                          recvRank, sizeTag, comm);
05817             }
05818 
05819             if (p + 1 < numProcs) {
05820               const int recvRank = p + 1;
05821               const int circBufInd = recvRank % 3;
05822 
05823               // Wait on receive-size receive from Process p + 1.
05824               if (debug) {
05825                 std::ostringstream os;
05826                 os << myRank << ": Wait on receive-size receive from Process "
05827                    << recvRank << endl;
05828                 *err << os.str ();
05829               }
05830               wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
05831 
05832               // We received the number of rows of data.  (The data
05833               // come in two columns.)
05834               const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
05835               if (debug && recvNumRows == -1) {
05836                 std::ostringstream os;
05837                 os << myRank << ": Result of receive-size receive from Process "
05838                    << recvRank << " is -1.  This should never happen, and "
05839                   "suggests that the receive never got posted.  Please report "
05840                   "this bug to the Tpetra developers." << endl;
05841                 *err << os.str ();
05842               }
05843 
05844               // Post receive-data receive from Process p + 1.
05845               recvDataBufs[circBufInd].resize (recvNumRows * 2);
05846               if (debug) {
05847                 std::ostringstream os;
05848                 os << myRank << ": Post receive-data receive from Process "
05849                    << recvRank << ": tag = " << dataTag << ", buffer size = "
05850                    << recvDataBufs[circBufInd].size () << endl;
05851                 *err << os.str ();
05852               }
05853               if (! recvDataReqs[circBufInd].is_null ()) {
05854                 std::ostringstream os;
05855                 os << myRank << ": recvDataReqs[" << circBufInd << "] is not "
05856                    << "null, for the receive-data receive from Process "
05857                    << recvRank << "!  This may mean that this process never "
05858                    << "finished waiting for the receive from Process "
05859                    << (recvRank - 3) << "." << endl;
05860                 *err << os.str ();
05861               }
05862               recvDataReqs[circBufInd] =
05863                 ireceive<int, int_type> (recvDataBufs[circBufInd],
05864                                          recvRank, dataTag, comm);
05865             }
05866 
05867             // Wait on receive-data receive from Process p.
05868             const int recvRank = p;
05869             const int circBufInd = recvRank % 3;
05870             if (debug) {
05871               std::ostringstream os;
05872               os << myRank << ": Wait on receive-data receive from Process "
05873                  << recvRank << endl;
05874               *err << os.str ();
05875             }
05876             wait<int> (comm, outArg (recvDataReqs[circBufInd]));
05877 
05878             // Write Process p's data.  Number of rows lives in
05879             // recvSizeBufs[circBufInd], and the actual data live in
05880             // recvDataBufs[circBufInd].  Do this after posting receives,
05881             // in order to expose overlap of comm. with file I/O.
05882             if (debug) {
05883               std::ostringstream os;
05884               os << myRank << ": Write GIDs and PIDs from Process "
05885                  << recvRank << endl;
05886               *err << os.str () << endl;
05887             }
05888             const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
05889             if (debug && printNumRows == -1) {
05890               std::ostringstream os;
05891               os << myRank << ": Result of receive-size receive from Process "
05892                  << recvRank << " was -1.  This should never happen, and "
05893                 "suggests that its receive-size receive was never posted.  "
05894                 "Please report this bug to the Tpetra developers." << endl;
05895               *err << os.str ();
05896             }
05897             if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
05898               std::ostringstream os;
05899               os << myRank << ": Result of receive-size receive from Proc "
05900                  << recvRank << " was " << printNumRows << " > 0, but "
05901                 "recvDataBufs[" << circBufInd << "] is null.  This should "
05902                 "never happen.  Please report this bug to the Tpetra "
05903                 "developers." << endl;
05904               *err << os.str ();
05905             }
05906             ArrayView<const int_type> printData = (recvDataBufs[circBufInd]) ();
05907             for (int_type k = 0; k < printNumRows; ++k) {
05908               const int_type gid = printData[2*k];
05909               const int_type pid = printData[2*k+1];
05910               out << gid << endl << pid << endl;
05911             }
05912           }
05913           else if (myRank == p) { // Process p
05914             // Wait on my send-data send.
05915             if (debug) {
05916               std::ostringstream os;
05917               os << myRank << ": Wait on my send-data send" << endl;
05918               *err << os.str ();
05919             }
05920             wait<int> (comm, outArg (sendReqData));
05921           }
05922           else if (myRank == p + 1) { // Process p + 1
05923             // Post send-data send to Process 0.
05924             if (debug) {
05925               std::ostringstream os;
05926               os << myRank << ": Post send-data send: tag = " << dataTag
05927                  << endl;
05928               *err << os.str ();
05929             }
05930             sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
05931             // Wait on my send-size send.
05932             if (debug) {
05933               std::ostringstream os;
05934               os << myRank << ": Wait on my send-size send" << endl;
05935               *err << os.str ();
05936             }
05937             wait<int> (comm, outArg (sendReqSize));
05938           }
05939           else if (myRank == p + 2) { // Process p + 2
05940             // Post send-size send to Process 0.
05941             if (debug) {
05942               std::ostringstream os;
05943               os << myRank << ": Post send-size send: size = "
05944                  << sendDataSize[0] << ", tag = " << sizeTag << endl;
05945               *err << os.str ();
05946             }
05947             sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
05948           }
05949         }
05950 
05951         if (! err.is_null ()) {
05952           err->popTab ();
05953         }
05954         if (debug) {
05955           *err << myRank << ": writeMap: Done" << endl;
05956         }
05957         if (! err.is_null ()) {
05958           err->popTab ();
05959         }
05960       }
05961 
05963       static void
05964       writeMapFile (const std::string& filename,
05965                     const map_type& map)
05966       {
05967         const int myRank = map.getComm ()->getRank ();
05968         std::ofstream out;
05969         if (myRank == 0) { // Only open the file on Proc 0.
05970           out.open (filename.c_str());
05971         }
05972         writeMap (out, map);
05973         // We can rely on the destructor of the output stream to close
05974         // the file on scope exit, even if writeDense() throws an
05975         // exception.
05976       }
05977 
05978     private:
06010       template<class MultiVectorScalarType>
06011       static void
06012       writeDenseImpl (std::ostream& out,
06013                       const Tpetra::MultiVector<MultiVectorScalarType, local_ordinal_type, global_ordinal_type, node_type>& X,
06014                       const std::string& matrixName,
06015                       const std::string& matrixDescription,
06016                       const Teuchos::RCP<Teuchos::FancyOStream>& err=Teuchos::null,
06017                       const bool debug=false)
06018       {
06019         using Tpetra::createOneToOne;
06020         using Tpetra::createNonContigMapWithNode;
06021         using Teuchos::as;
06022         using Teuchos::rcp;
06023         using Teuchos::rcp_const_cast;
06024         using Teuchos::rcp_dynamic_cast;
06025         using std::endl;
06026 
06027         typedef MultiVectorScalarType ST;
06028         typedef local_ordinal_type LO;
06029         typedef global_ordinal_type GO;
06030         typedef node_type NT;
06031 
06032         typedef typename Teuchos::ScalarTraits<ST> STS;
06033         typedef Tpetra::MultiVector<ST, LO, GO, NT> MV;
06034 
06035         RCP<const map_type> map = X.getMap ();
06036         RCP<const Comm<int> > comm = map->getComm ();
06037         const int myRank = comm->getRank ();
06038 
06039         if (! err.is_null ()) {
06040           err->pushTab ();
06041         }
06042         if (debug) {
06043           std::ostringstream os;
06044           os << myRank << ": writeDenseImpl" << endl;
06045           *err << os.str ();
06046         }
06047         if (! err.is_null ()) {
06048           err->pushTab ();
06049         }
06050 
06051         // Make the output stream write floating-point numbers in
06052         // scientific notation.  It will politely put the output
06053         // stream back to its state on input, when this scope
06054         // terminates.
06055         Teuchos::MatrixMarket::details::SetScientific<ST> sci (out);
06056 
06057         const global_size_t numRows = map->getGlobalNumElements ();
06058         // Promote to global_size_t.
06059         const global_size_t numCols = as<global_size_t> (X.getNumVectors ());
06060 
06061         if (debug) {
06062           std::ostringstream os;
06063           os << myRank << ": writeDenseImpl: Computing gather Map" << endl;
06064           *err << os.str ();
06065         }
06066 
06067         // Make the "gather" map, where Proc 0 owns all rows of X, and
06068         // the other procs own no rows.
06069         RCP<NT> node = map->getNode();
06070         RCP<const map_type> gatherMap =
06071           Details::computeGatherMap<map_type> (map, err, debug);
06072 
06073         if (debug) {
06074           std::ostringstream os;
06075           os << myRank << ": writeDenseImpl: Computing Import" << endl;
06076           *err << os.str ();
06077         }
06078 
06079         // Create an Import object to import X's data into a
06080         // multivector Y owned entirely by Proc 0.  In the Import
06081         // constructor, X's map is the source map, and the "gather
06082         // map" is the target map.
06083         Import<LO, GO, NT> importer (map, gatherMap, err);
06084 
06085         if (debug) {
06086           std::ostringstream os;
06087           os << myRank << ": writeDenseImpl: Creating target MultiVector" << endl;
06088           *err << os.str ();
06089         }
06090 
06091         // Create a new multivector Y to hold the result of the import.
06092         RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (gatherMap, numCols);
06093 
06094         if (debug) {
06095           std::ostringstream os;
06096           os << myRank << ": writeDenseImpl: Doing Import" << endl;
06097           *err << os.str ();
06098         }
06099 
06100         // Import the multivector onto Proc 0.
06101         Y->doImport (X, importer, INSERT);
06102 
06103         if (debug) {
06104           std::ostringstream os;
06105           os << myRank << ": writeDenseImpl: Writing data to out" << endl;
06106           *err << os.str ();
06107         }
06108 
06109         //
06110         // Print the matrix in Matrix Market format on Proc 0.
06111         //
06112         std::ostringstream exMsg; // for exception messages on Proc 0
06113         int localWriteSuccess = 1; // whether the write succeeded on Proc 0
06114         if (myRank == 0) {
06115           try {
06116             std::string dataType;
06117             if (STS::isComplex) {
06118               dataType = "complex";
06119             } else if (STS::isOrdinal) {
06120               dataType = "integer";
06121             } else {
06122               dataType = "real";
06123             }
06124 
06125             // Print the Matrix Market banner line.  MultiVector
06126             // stores data nonsymmetrically, hence "general".
06127             out << "%%MatrixMarket matrix array "
06128                 << dataType << " general" << endl;
06129 
06130             // Print comments (the matrix name and / or description).
06131             if (matrixName != "") {
06132               printAsComment (out, matrixName);
06133             }
06134             if (matrixDescription != "") {
06135               printAsComment (out, matrixDescription);
06136             }
06137             // Print the Matrix Market dimensions header for dense matrices.
06138             out << numRows << " " << numCols << endl;
06139 
06140             // Get a read-only view of the entries of Y.
06141             // Rank 0 owns all of the entries of Y.
06142             ArrayRCP<ArrayRCP<const ST> > Y_view = Y->get2dView ();
06143             TEUCHOS_TEST_FOR_EXCEPTION(
06144               as<global_size_t> (Y_view.size ()) < numCols, std::logic_error,
06145               "Y_view has size " << Y_view.size () << " < numCols = "
06146               << numCols << ".");
06147 
06148             // Print the entries of the matrix, in column-major order.
06149             for (global_size_t j = 0; j < numCols; ++j) {
06150               ArrayRCP<const ST> Y_j = Y_view[j];
06151               TEUCHOS_TEST_FOR_EXCEPTION(
06152                 as<global_size_t> (Y_j.size ()) < numRows, std::logic_error,
06153                 "The current column's data has length " << Y_j.size ()
06154                 << " < numRows=" << numRows << ".");
06155               for (global_size_t i = 0; i < numRows; ++i) {
06156                 const ST Y_ij = Y_j[i];
06157                 if (STS::isComplex) {
06158                   out << STS::real (Y_ij) << " " << STS::imag (Y_ij) << endl;
06159                 } else {
06160                   out << Y_ij << endl;
06161                 }
06162               }
06163             }
06164           } catch (std::exception& e) {
06165             exMsg << e.what ();
06166             localWriteSuccess = 0;
06167           }
06168         } // if (myRank == 0)
06169 
06170         // Synchronize on error.  This is not necessary here, but it
06171         // will prevent later code from hanging in Proc 0 threw an
06172         // exception above.
06173         int globalWriteSuccess = localWriteSuccess;
06174         Teuchos::broadcast (*comm, 0, outArg (globalWriteSuccess));
06175         TEUCHOS_TEST_FOR_EXCEPTION(
06176           globalWriteSuccess == 0, std::runtime_error,
06177           "Failed to write data: " << exMsg.str ());
06178 
06179         if (! err.is_null ()) {
06180           err->popTab ();
06181         }
06182         if (debug) {
06183           *err << myRank << ": writeDenseImpl: done" << endl;
06184         }
06185         if (! err.is_null ()) {
06186           err->popTab ();
06187         }
06188       }
06189 
06213       static void
06214       printAsComment (std::ostream& out, const std::string& str)
06215       {
06216         using std::endl;
06217         std::istringstream inpstream (str);
06218         std::string line;
06219 
06220         while (getline (inpstream, line)) {
06221           if (! line.empty()) {
06222             // Note that getline() doesn't store '\n', so we have to
06223             // append the endline ourselves.
06224             if (line[0] == '%') { // Line starts with a comment character.
06225               out << line << endl;
06226             }
06227             else { // Line doesn't start with a comment character.
06228               out << "%% " << line << endl;
06229             }
06230           }
06231         }
06232       }
06233     }; // class Writer
06234 
06235   } // namespace MatrixMarket
06236 } // namespace Tpetra
06237 
06238 #endif // __MatrixMarket_Tpetra_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines