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_CrsMatrix.hpp"
00057 #include "Teuchos_MatrixMarket_Raw_Adder.hpp"
00058 #include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp"
00059 #include "Teuchos_MatrixMarket_assignScalar.hpp"
00060 #include "Teuchos_MatrixMarket_Banner.hpp"
00061 #include "Teuchos_MatrixMarket_CoordDataReader.hpp"
00062 #include "Teuchos_MatrixMarket_SetScientific.hpp"
00063 
00064 #include <algorithm>
00065 #include <fstream>
00066 #include <iostream>
00067 #include <iterator>
00068 #include <vector>
00069 #include <stdexcept>
00070 
00071 namespace Tpetra {
00099   namespace MatrixMarket {
00100 
00164     template<class SparseMatrixType>
00165     class Reader {
00166     public:
00169       typedef SparseMatrixType sparse_matrix_type;
00170       typedef RCP<sparse_matrix_type> sparse_matrix_ptr;
00171 
00174       typedef typename SparseMatrixType::scalar_type scalar_type;
00177       typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
00184       typedef typename SparseMatrixType::global_ordinal_type
00185         global_ordinal_type;
00186 
00189       typedef typename SparseMatrixType::node_type node_type;
00190 
00193       typedef MultiVector<scalar_type,
00194                           local_ordinal_type,
00195                           global_ordinal_type,
00196                           node_type> multivector_type;
00197 
00198       typedef RCP<node_type> node_ptr;
00199       typedef Comm<int> comm_type;
00200       typedef RCP<const comm_type> comm_ptr;
00201       typedef Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
00202       typedef RCP<const map_type> map_ptr;
00203 
00204     private:
00205 
00208       typedef typename ArrayRCP<global_ordinal_type>::size_type size_type;
00209 
00222       static RCP<const map_type>
00223       makeRangeMap (const RCP<const comm_type>& pComm,
00224                     const RCP<node_type>& pNode,
00225                     const global_ordinal_type numRows)
00226       {
00227         // A conventional, uniformly partitioned, contiguous map.
00228         return rcp (new map_type (static_cast<global_size_t> (numRows),
00229                                   static_cast<global_ordinal_type> (0),
00230                                   pComm, GloballyDistributed, pNode));
00231       }
00232 
00265       static RCP<const map_type>
00266       makeRowMap (const RCP<const map_type>& pRowMap,
00267                   const RCP<const comm_type>& pComm,
00268                   const RCP<node_type>& pNode,
00269                   const global_ordinal_type numRows)
00270       {
00271         // If the caller didn't provide a map, return a conventional,
00272         // uniformly partitioned, contiguous map.
00273         if (pRowMap.is_null()) {
00274           return rcp (new map_type (static_cast<global_size_t> (numRows),
00275                                     static_cast<global_ordinal_type> (0),
00276                                     pComm, GloballyDistributed, pNode));
00277         } else {
00278           TEUCHOS_TEST_FOR_EXCEPTION(! pRowMap->isDistributed() && pComm->getSize() > 1,
00279                              std::invalid_argument,
00280                              "The specified row map is not distributed, but "
00281                              "the given communicator includes more than one "
00282                              "rank (in fact, there are " << pComm->getSize()
00283                              << " ranks).");
00284           TEUCHOS_TEST_FOR_EXCEPTION(pRowMap->getComm() != pComm,
00285                              std::invalid_argument,
00286                              "The specified row map's communicator (pRowMap->"
00287                              "getComm()) is different than the given separately "
00288                              "supplied communicator pComm.");
00289           TEUCHOS_TEST_FOR_EXCEPTION(pRowMap->getNode() != pNode,
00290                              std::invalid_argument,
00291                              "The specified row map's node (pRowMap->getNode()) "
00292                              "is different than the given separately supplied "
00293                              "node pNode.");
00294           return pRowMap;
00295         }
00296       }
00297 
00312       static map_ptr
00313       makeDomainMap (const map_ptr& pRangeMap,
00314                      const global_ordinal_type numRows,
00315                      const global_ordinal_type numCols)
00316       {
00317         // Abbreviations so that the map creation call isn't too long.
00318         typedef local_ordinal_type LO;
00319         typedef global_ordinal_type GO;
00320         typedef node_type Node;
00321 
00322         if (numRows == numCols) {
00323           return pRangeMap;
00324         } else {
00325           comm_ptr pComm = pRangeMap->getComm();
00326           node_ptr pNode = pRangeMap->getNode();
00327           return createUniformContigMapWithNode<LO,GO,Node> (numCols,
00328                                                              pComm,
00329                                                              pNode);
00330         }
00331       }
00332 
00405       static void
00406       distribute (ArrayRCP<size_t>& myNumEntriesPerRow,
00407                   ArrayRCP<size_t>& myRowPtr,
00408                   ArrayRCP<global_ordinal_type>& myColInd,
00409                   ArrayRCP<scalar_type>& myValues,
00410                   const RCP<const map_type>& pRowMap,
00411                   ArrayRCP<size_t>& numEntriesPerRow,
00412                   ArrayRCP<size_t>& rowPtr,
00413                   ArrayRCP<global_ordinal_type>& colInd,
00414                   ArrayRCP<scalar_type>& values,
00415                   const bool debug=false)
00416       {
00417          using Teuchos::CommRequest;
00418          using Teuchos::receive;
00419          using Teuchos::send;
00420          using std::cerr;
00421          using std::endl;
00422 
00423          const bool extraDebug = false;
00424          comm_ptr pComm = pRowMap->getComm();
00425          const int numProcs = Teuchos::size (*pComm);
00426          const int myRank = Teuchos::rank (*pComm);
00427          const int rootRank = 0;
00428 
00429          // Type abbreviations to make the code more concise.
00430          typedef global_ordinal_type GO;
00431          typedef local_ordinal_type LO;
00432 
00433          // List of the global indices of my rows.  They may or may
00434          // not be contiguous, and the row map need not be one-to-one.
00435          ArrayView<const GO> myRows = pRowMap->getNodeElementList();
00436          const size_type myNumRows = myRows.size();
00437          TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
00438                             pRowMap->getNodeNumElements(),
00439                             std::logic_error,
00440                             "pRowMap->getNodeElementList().size() = "
00441                             << myNumRows
00442                             << " != pRowMap->getNodeNumElements() = "
00443                             << pRowMap->getNodeNumElements() << ".  "
00444                             "Please report this bug to the Tpetra developers.");
00445          TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
00446                             std::logic_error,
00447                             "On Proc 0: numEntriesPerRow.size() = "
00448                             << numEntriesPerRow.size()
00449                             << " != pRowMap->getNodeElementList().size() = "
00450                             << myNumRows << ".  Please report this bug to the "
00451                             "Tpetra developers.");
00452 
00453          // Space for my proc's number of entries per row.
00454          // Will be filled in below.
00455          myNumEntriesPerRow = arcp<size_t> (myNumRows);
00456 
00457          // Teuchos::receive() returns an int; here is space for it.
00458          int recvResult = 0;
00459 
00460          if (myRank != rootRank)
00461            {
00462              // Tell the root how many rows we have.  If we're sending
00463              // none, then we don't have anything else to send, nor
00464              // does the root have to receive anything else.
00465              send (*pComm, myNumRows, rootRank);
00466              if (myNumRows != 0)
00467                {
00468                  // Now send my rows' global indices.  Hopefully the
00469                  // cast to int doesn't overflow.  This is unlikely,
00470                  // since it should fit in a LO, even
00471                  // though it is a GO.
00472                  send (*pComm, static_cast<int> (myNumRows),
00473                        myRows.getRawPtr(), rootRank);
00474 
00475                  // I (this proc) don't care if my global row indices
00476                  // are contiguous, though the root proc does (since
00477                  // otherwise it needs to pack noncontiguous data into
00478                  // contiguous storage before sending).  That's why we
00479                  // don't check for contiguousness here.
00480 
00481                  // Ask the root processor for my part of the array of the
00482                  // number of entries per row.
00483                  recvResult = receive (*pComm, rootRank,
00484                                        static_cast<int> (myNumRows),
00485                                        myNumEntriesPerRow.getRawPtr());
00486 
00487                  // Use the resulting array to figure out how many column
00488                  // indices and values for which I should ask from the root
00489                  // processor.
00490                  const local_ordinal_type myNumEntries =
00491                    std::accumulate (myNumEntriesPerRow.begin(),
00492                                     myNumEntriesPerRow.end(),
00493                                     0);
00494 
00495                  // Make space for my entries of the sparse matrix.
00496                  // Note that they don't have to be sorted by row
00497                  // index.  Iterating through all my rows requires
00498                  // computing a running sum over myNumEntriesPerRow.
00499                  myColInd = arcp<GO> (myNumEntries);
00500                  myValues = arcp<scalar_type> (myNumEntries);
00501                  if (myNumEntries > 0)
00502                    { // Ask for that many column indices and values,
00503                      // if there are any.
00504                      recvResult = receive (*pComm, rootRank,
00505                                            static_cast<int> (myNumEntries),
00506                                            myColInd.getRawPtr());
00507                      recvResult = receive (*pComm, rootRank,
00508                                            static_cast<int> (myNumEntries),
00509                                            myValues.getRawPtr());
00510                    }
00511                } // If I own at least one row
00512            } // If I am not the root processor
00513          else // I _am_ the root processor
00514            {
00515              if (debug)
00516                cerr << "-- Proc 0: Copying my data from global arrays" << endl;
00517 
00518              // Proc 0 (the root processor) still needs to (allocate,
00519              // if not done already) and fill its part of the matrix
00520              // (my*).
00521              for (size_type k = 0; k < myNumRows; ++k)
00522                {
00523                  const GO myCurRow = myRows[k];
00524                  const local_ordinal_type numEntriesInThisRow = numEntriesPerRow[myCurRow];
00525                  //myNumEntriesPerRow[k] = numEntriesPerRow[myCurRow];
00526                  myNumEntriesPerRow[k] = numEntriesInThisRow;
00527                }
00528              if (extraDebug && debug)
00529                {
00530                  cerr << "Proc " << Teuchos::rank (*(pRowMap->getComm()))
00531                       << ": myNumEntriesPerRow[0.." << (myNumRows-1) << "] = [";
00532                  for (size_type k = 0; k < myNumRows; ++k)
00533                    {
00534                      cerr << myNumEntriesPerRow[k];
00535                      if (k < myNumRows-1)
00536                        cerr << " ";
00537                    }
00538                  cerr << "]" << endl;
00539                }
00540              // The total number of matrix entries that my proc owns.
00541              const local_ordinal_type myNumEntries =
00542                std::accumulate (myNumEntriesPerRow.begin(),
00543                                 myNumEntriesPerRow.end(),
00544                                 0);
00545              if (debug)
00546                cerr << "-- Proc 0: I own " << myNumRows << " rows and "
00547                     << myNumEntries << " entries" << endl;
00548 
00549              myColInd = arcp<GO> (myNumEntries);
00550              myValues = arcp<scalar_type> (myNumEntries);
00551 
00552              // Copy Proc 0's part of the matrix into the my* arrays.
00553              // It's important that myCurPos be updated _before_ k,
00554              // otherwise myCurPos will get the wrong number of entries
00555              // per row (it should be for the row in the just-completed
00556              // iteration, not for the next iteration's row).
00557              local_ordinal_type myCurPos = 0;
00558              for (size_type k = 0; k < myNumRows;
00559                   myCurPos += myNumEntriesPerRow[k], ++k)
00560                {
00561                  const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
00562                  const GO myRow = myRows[k];
00563                  const size_t curPos = rowPtr[myRow];
00564                  if (extraDebug && debug)
00565                    {
00566                      cerr << "k = " << k << ", myRow = " << myRow << ": colInd("
00567                           << curPos << "," << curNumEntries << "), myColInd("
00568                           << myCurPos << "," << curNumEntries << ")" << endl;
00569                    }
00570                  // Only copy if there are entries to copy, in order
00571                  // not to construct empty ranges for the ArrayRCP
00572                  // views.
00573                  if (curNumEntries > 0)
00574                    {
00575                      ArrayView<GO> colIndView = colInd(curPos, curNumEntries);
00576                      ArrayView<GO> myColIndView =
00577                        myColInd(myCurPos, curNumEntries);
00578                      std::copy (colIndView.begin(), colIndView.end(),
00579                                 myColIndView.begin());
00580 
00581                      ArrayView<scalar_type> valuesView =
00582                        values(curPos, curNumEntries);
00583                      ArrayView<scalar_type> myValuesView =
00584                        myValues(myCurPos, curNumEntries);
00585                      std::copy (valuesView.begin(), valuesView.end(),
00586                                 myValuesView.begin());
00587                    }
00588                }
00589 
00590              // Proc 0 processes each other proc p in turn.
00591              for (int p = 1; p < numProcs; ++p)
00592                {
00593                  if (debug)
00594                    cerr << "-- Proc 0: Processing proc " << p << endl;
00595 
00596                  size_type theirNumRows = 0;
00597                  // Ask Proc p how many rows it has.  If it doesn't
00598                  // have any, we can move on to the next proc.  This
00599                  // has to be a standard receive so that we can avoid
00600                  // the degenerate case of sending zero data.
00601                  recvResult = receive (*pComm, p, &theirNumRows);
00602                  if (debug)
00603                    cerr << "-- Proc 0: Proc " << p << " owns "
00604                         << theirNumRows << " rows" << endl;
00605                  if (theirNumRows != 0)
00606                    { // Ask Proc p which rows it owns.  The resulting
00607                      // global row indices are not guaranteed to be
00608                      // contiguous or sorted.  Global row indices are
00609                      // themselves indices into the numEntriesPerRow
00610                      // array.
00611                      ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
00612                      recvResult = receive (*pComm, p,
00613                                            static_cast<int> (theirNumRows),
00614                                            theirRows.getRawPtr());
00615                      // Extra test to make sure that the rows we
00616                      // received are all sensible.  This is a good
00617                      // idea since we are going to use the global row
00618                      // indices we've received to index into the
00619                      // numEntriesPerRow array.  Better to catch any
00620                      // bugs here and print a sensible error message,
00621                      // rather than segfault and print a cryptic error
00622                      // message.
00623                      {
00624                        const global_size_t numRows =
00625                          pRowMap->getGlobalNumElements();
00626                        bool theirRowsValid = true;
00627                        for (size_type k = 0; k < theirNumRows; ++k)
00628                          {
00629                            // global_ordinal_t is generally a signed type.
00630                            if (theirRows[k] < 0)
00631                              theirRowsValid = false;
00632                            // Signed to unsigned cast should not
00633                            // overflow.  We cast in order to avoid
00634                            // compiler warnings about comparing signed
00635                            // and unsigned types, in case
00636                            // global_size_t is unsigned.
00637                            else if (static_cast<global_size_t> (theirRows[k]) >= numRows)
00638                              theirRowsValid = false;
00639                          }
00640                        if (! theirRowsValid)
00641                          {
00642                            std::ostringstream os;
00643                            std::copy (theirRows.begin(), theirRows.end(),
00644                                       std::ostream_iterator<GO>(os, " "));
00645                            TEUCHOS_TEST_FOR_EXCEPTION(! theirRowsValid,
00646                                               std::logic_error,
00647                                               "Proc " << p << " has at least "
00648                                               "one invalid row index.  Here are"
00649                                               " all of them: [ " << os.str()
00650                                               << "]");
00651                          }
00652                      }
00653 
00654                      // Perhaps we could save a little work if we
00655                      // check whether Proc p's row indices are
00656                      // contiguous.  That would make lookups in the
00657                      // global data arrays faster.  For now, we just
00658                      // implement the general case and don't
00659                      // prematurely optimize.  (Remember that you're
00660                      // making Proc 0 read the whole file, so you've
00661                      // already lost scalability.)
00662 
00663                      // Compute the number of entries in each of Proc
00664                      // p's rows.  (Proc p will compute its row
00665                      // pointer array on its own, after it gets the
00666                      // data from Proc 0.)
00667                      ArrayRCP<size_t> theirNumEntriesPerRow;
00668                      theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
00669                      for (size_type k = 0; k < theirNumRows; ++k)
00670                        theirNumEntriesPerRow[k] =
00671                          numEntriesPerRow[theirRows[k]];
00672 
00673                      // Tell Proc p the number of entries in each of
00674                      // its rows.  Hopefully the cast to int doesn't
00675                      // overflow.  This is unlikely, since it should
00676                      // fit in a LO, even though it is
00677                      // a GO.
00678                      send (*pComm, static_cast<int> (theirNumRows),
00679                            theirNumEntriesPerRow.getRawPtr(), p);
00680 
00681                      // Figure out how many entries Proc p owns.
00682                      const local_ordinal_type theirNumEntries =
00683                        std::accumulate (theirNumEntriesPerRow.begin(),
00684                                         theirNumEntriesPerRow.end(),
00685                                         0);
00686 
00687                      if (debug)
00688                        cerr << "-- Proc 0: Proc " << p << " owns "
00689                             << theirNumEntries << " entries" << endl;
00690 
00691                      // If there are no entries to send, then we're
00692                      // done with Proc p.
00693                      if (theirNumEntries == 0)
00694                        continue;
00695 
00696                      // Construct (views of) proc p's column indices
00697                      // and values.  Later, we might like to optimize
00698                      // for the (common) contiguous case, for which we
00699                      // don't need to copy data into separate "their*"
00700                      // arrays (we can just use contiguous views of
00701                      // the global arrays).
00702                      ArrayRCP<GO> theirColInd = arcp<GO> (theirNumEntries);
00703                      ArrayRCP<scalar_type> theirValues =
00704                        arcp<scalar_type> (theirNumEntries);
00705                      // Copy Proc p's part of the matrix into the
00706                      // their* arrays.  It's important that
00707                      // theirCurPos be updated _before_ k, otherwise
00708                      // theirCurPos will get the wrong number of
00709                      // entries per row (it should be for the row in
00710                      // the just-completed iteration, not for the next
00711                      // iteration's row).
00712                      local_ordinal_type theirCurPos = 0;
00713                      for (size_type k = 0; k < theirNumRows;
00714                           theirCurPos += theirNumEntriesPerRow[k], k++)
00715                        {
00716                          const local_ordinal_type curNumEntries = theirNumEntriesPerRow[k];
00717                          const GO theirRow = theirRows[k];
00718                          const local_ordinal_type curPos = rowPtr[theirRow];
00719 
00720                          // Only copy if there are entries to copy, in
00721                          // order not to construct empty ranges for
00722                          // the ArrayRCP views.
00723                          if (curNumEntries > 0)
00724                            {
00725                              ArrayView<GO> colIndView =
00726                                colInd(curPos, curNumEntries);
00727                              ArrayView<GO> theirColIndView =
00728                                theirColInd(theirCurPos, curNumEntries);
00729                              std::copy (colIndView.begin(), colIndView.end(),
00730                                         theirColIndView.begin());
00731 
00732                              ArrayView<scalar_type> valuesView =
00733                                values(curPos, curNumEntries);
00734                              ArrayView<scalar_type> theirValuesView =
00735                                theirValues(theirCurPos, curNumEntries);
00736                              std::copy (valuesView.begin(), valuesView.end(),
00737                                         theirValuesView.begin());
00738                            }
00739                        }
00740                      // Send Proc p its column indices and values.
00741                      // Hopefully the cast to int doesn't overflow.
00742                      // This is unlikely, since it should fit in a LO,
00743                      // even though it is a GO.
00744                      send (*pComm, static_cast<int> (theirNumEntries),
00745                            theirColInd.getRawPtr(), p);
00746                      send (*pComm, static_cast<int> (theirNumEntries),
00747                            theirValues.getRawPtr(), p);
00748 
00749                      if (debug)
00750                        cerr << "-- Proc 0: Finished with proc " << p << endl;
00751                    } // If proc p owns at least one row
00752                } // For each proc p not the root proc 0
00753            } // If I'm (not) the root proc 0
00754 
00755          // Invalidate the input data to save space, since we don't
00756          // need it anymore.
00757          numEntriesPerRow = null;
00758          rowPtr = null;
00759          colInd = null;
00760          values = null;
00761 
00762          if (debug && myRank == 0)
00763            cerr << "-- Proc 0: About to fill in myRowPtr" << endl;
00764 
00765          // Allocate and fill in myRowPtr (the row pointer array for
00766          // my rank's rows).  We delay this until the end because we
00767          // don't need it to compute anything else in distribute().
00768          // Each proc can do this work for itself, since it only needs
00769          // myNumEntriesPerRow to do so.
00770          myRowPtr = arcp<size_t> (myNumRows+1);
00771          myRowPtr[0] = 0;
00772          for (size_type k = 1; k < myNumRows+1; ++k) {
00773            myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
00774          }
00775          if (extraDebug && debug)
00776            {
00777              cerr << "Proc " << Teuchos::rank (*(pRowMap->getComm()))
00778                   << ": myRowPtr[0.." << myNumRows << "] = [";
00779              for (size_type k = 0; k < myNumRows+1; ++k)
00780                {
00781                  cerr << myRowPtr[k];
00782                  if (k < myNumRows)
00783                    cerr << " ";
00784                }
00785              cerr << "]" << endl << endl;
00786            }
00787 
00788          if (debug && myRank == 0)
00789            cerr << "-- Proc 0: Done with distribute" << endl;
00790       }
00791 
00802       static sparse_matrix_ptr
00803       makeMatrix (ArrayRCP<size_t>& myNumEntriesPerRow,
00804                   ArrayRCP<size_t>& myRowPtr,
00805                   ArrayRCP<global_ordinal_type>& myColInd,
00806                   ArrayRCP<scalar_type>& myValues,
00807                   const map_ptr& pRowMap,
00808                   const map_ptr& pRangeMap,
00809                   const map_ptr& pDomainMap,
00810                   const bool callFillComplete = true)
00811       {
00812         using std::cerr;
00813         using std::endl;
00814         // Typedef to make certain type declarations shorter.
00815         typedef global_ordinal_type GO;
00816 
00817         const bool extraDebug = false;
00818         const bool debug = false;
00819 
00820         // The row pointer array always has at least one entry, even
00821         // if the matrix has zero rows.  myNumEntriesPerRow, myColInd,
00822         // and myValues would all be empty arrays in that degenerate
00823         // case, but the row and domain maps would still be nonnull
00824         // (though they would be trivial maps).
00825         TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
00826                            "makeMatrix: myRowPtr array is null.  "
00827                            "Please report this bug to the Tpetra developers.");
00828         TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
00829                            "makeMatrix: domain map is null.  "
00830                            "Please report this bug to the Tpetra developers.");
00831         TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
00832                            "makeMatrix: range map is null.  "
00833                            "Please report this bug to the Tpetra developers.");
00834         TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
00835                            "makeMatrix: row map is null.  "
00836                            "Please report this bug to the Tpetra developers.");
00837 
00838         // Handy for debugging output; not needed otherwise.
00839         const int myRank = Teuchos::rank (*(pRangeMap->getComm()));
00840 
00841         if (extraDebug && debug)
00842           {
00843             cerr << "Proc " << myRank << ":" << endl
00844                  << "-- myRowPtr = [ ";
00845             std::copy (myRowPtr.begin(), myRowPtr.end(),
00846                        std::ostream_iterator<size_type>(cerr, " "));
00847             cerr << "]" << endl << "-- myColInd = [ ";
00848             std::copy (myColInd.begin(), myColInd.end(),
00849                        std::ostream_iterator<size_type>(cerr, " "));
00850             cerr << "]" << endl << endl;
00851           }
00852 
00853         // Go through all of my columns, and see if any are not in the
00854         // domain map.  This is possible if numProcs > 1, otherwise
00855         // not.
00856         if (extraDebug && debug)
00857           {
00858             size_type numRemote = 0;
00859             std::vector<GO> remoteGIDs;
00860 
00861             typedef typename ArrayRCP<GO>::const_iterator iter_type;
00862             for (iter_type it = myColInd.begin(); it != myColInd.end(); ++it)
00863               {
00864                 if (! pDomainMap->isNodeGlobalElement (*it))
00865                   {
00866                     numRemote++;
00867                     remoteGIDs.push_back (*it);
00868                   }
00869               }
00870             if (numRemote > 0)
00871               {
00872                 cerr << "Proc " << myRank << ": " << numRemote
00873                      << " remote GIDs = [ " << endl;
00874                 std::copy (remoteGIDs.begin(), remoteGIDs.end(),
00875                            std::ostream_iterator<GO>(cerr, " "));
00876                 cerr << "]" << endl;
00877               }
00878           }
00879 
00880         // Construct the CrsMatrix, using the row map, with the
00881         // constructor specifying the number of nonzeros for each row.
00882         // Create with DynamicProfile, so that the fillComplete() can
00883         // do first-touch reallocation (a NUMA (Non-Uniform Memory
00884         // Access) optimization on multicore CPUs).
00885         sparse_matrix_ptr A =
00886           rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
00887                                        DynamicProfile));
00888         TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), std::logic_error,
00889                            "makeMatrix: Initial allocation of CrsMatrix failed"
00890                            ".  Please report this bug to the Tpetra developers"
00891                            ".");
00892         // List of the global indices of my rows.
00893         // They may or may not be contiguous.
00894         ArrayView<const GO> myRows = pRowMap->getNodeElementList();
00895         const size_type myNumRows = myRows.size();
00896 
00897         // Add this processor's matrix entries to the CrsMatrix.
00898         for (size_type k = 0; k < myNumRows; ++k)
00899           {
00900             const size_type myCurPos = myRowPtr[k];
00901             const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
00902 
00903             if (extraDebug && debug)
00904               {
00905                 cerr << "Proc " << myRank << ": k = " << k
00906                      << ", myCurPos = " << myCurPos
00907                      << ", curNumEntries = " << curNumEntries
00908                      << endl;
00909               }
00910             // Avoid constructing empty views of ArrayRCP objects.
00911             if (curNumEntries > 0)
00912               A->insertGlobalValues (myRows[k],
00913                                      myColInd(myCurPos, curNumEntries),
00914                                      myValues(myCurPos, curNumEntries));
00915           }
00916         // We've entered in all our matrix entries, so we can delete
00917         // the original data.  This will save memory when we call
00918         // fillComplete(), so that we never keep more than two copies
00919         // of the matrix's data in memory at once.
00920         myNumEntriesPerRow = null;
00921         myRowPtr = null;
00922         myColInd = null;
00923         myValues = null;
00924 
00925         if (callFillComplete)
00926           A->fillComplete (pDomainMap, pRangeMap);
00927         return A;
00928       }
00929 
00930     private:
00931 
00948       static RCP<const Teuchos::MatrixMarket::Banner>
00949       readBanner (std::istream& in,
00950                   size_t& lineNumber,
00951                   const bool tolerant=false,
00952                   const bool debug=false)
00953       {
00954         using Teuchos::MatrixMarket::Banner;
00955         using std::cerr;
00956         using std::endl;
00957         typedef Teuchos::ScalarTraits<scalar_type> STS;
00958 
00959         RCP<Banner> pBanner; // On output, if successful: the read-in Banner.
00960         std::string line; // If read from stream successful: the Banner line
00961 
00962         // Try to read a line from the input stream.
00963         const bool readFailed = ! getline(in, line);
00964         TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
00965           "Failed to get Matrix Market banner line from input.");
00966 
00967         // We read a line from the input stream.
00968         lineNumber++;
00969 
00970         // Assume that the line we found is the Banner line.
00971         try {
00972           pBanner = rcp (new Banner (line, tolerant));
00973         } catch (std::exception& e) {
00974           TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
00975             "Matrix Market banner line contains syntax error(s): "
00976             << e.what());
00977         }
00978         TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() != "matrix",
00979           std::invalid_argument, "The Matrix Market file does not contain "
00980           "matrix data.  Its Banner (first) line says that its object type is \""
00981           << pBanner->matrixType() << "\", rather than the required \"matrix\".");
00982 
00983         // Validate the data type of the matrix, with respect to the
00984         // Scalar type of the CrsMatrix entries.
00985         TEUCHOS_TEST_FOR_EXCEPTION(
00986           ! STS::isComplex && pBanner->dataType() == "complex",
00987           std::invalid_argument,
00988           "The Matrix Market file contains complex-valued data, but you are "
00989           "trying to read it into a matrix containing entries of the real-"
00990           "valued Scalar type \""
00991           << Teuchos::TypeNameTraits<scalar_type>::name() << "\".");
00992         TEUCHOS_TEST_FOR_EXCEPTION(
00993           pBanner->dataType() != "real" &&
00994           pBanner->dataType() != "complex" &&
00995           pBanner->dataType() != "integer",
00996           std::invalid_argument,
00997           "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
00998           "Matrix Market file may not contain a \"pattern\" matrix.  A "
00999           "pattern matrix is really just a graph with no weights.  It "
01000           "should be stored in a CrsGraph, not a CrsMatrix.");
01001 
01002         return pBanner;
01003       }
01004 
01027       static Tuple<global_ordinal_type, 3>
01028       readCoordDims (std::istream& in,
01029                      size_t& lineNumber,
01030                      const RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
01031                      const comm_ptr& pComm,
01032                      const bool tolerant = false,
01033                      const bool debug = false)
01034       {
01035         using Teuchos::MatrixMarket::readCoordinateDimensions;
01036 
01037         // Packed coordinate matrix dimensions (numRows, numCols,
01038         // numNonzeros); computed on Rank 0 and broadcasted to all MPI
01039         // ranks.
01040         Tuple<global_ordinal_type, 3> dims;
01041 
01042         // Read in the coordinate matrix dimensions from the input
01043         // stream.  "success" tells us whether reading in the
01044         // coordinate matrix dimensions succeeded ("Guilty unless
01045         // proven innocent").
01046         bool success = false;
01047         if (pComm->getRank() == 0) {
01048           TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "coordinate",
01049             std::invalid_argument, "The Tpetra::CrsMatrix Matrix Market reader "
01050             "only accepts \"coordinate\" (sparse) matrix data.");
01051           // Unpacked coordinate matrix dimensions
01052           global_ordinal_type numRows, numCols, numNonzeros;
01053           // Only MPI Rank 0 reads from the input stream
01054           success = readCoordinateDimensions (in, numRows, numCols,
01055                                               numNonzeros, lineNumber,
01056                                               tolerant);
01057           // Pack up the data into a Tuple so we can send them with
01058           // one broadcast instead of three.
01059           dims[0] = numRows;
01060           dims[1] = numCols;
01061           dims[2] = numNonzeros;
01062         }
01063         // Only Rank 0 did the reading, so it decides success.
01064         //
01065         // FIXME (mfh 02 Feb 2011) Teuchos::broadcast doesn't know how
01066         // to send bools.  For now, we convert to/from int instead,
01067         // using the usual "true is 1, false is 0" encoding.
01068         {
01069           int the_success = success ? 1 : 0; // only matters on MPI Rank 0
01070           Teuchos::broadcast (*pComm, 0, &the_success);
01071           success = (the_success == 1);
01072         }
01073         if (success) {
01074           // Broadcast (numRows, numCols, numNonzeros) from Rank 0
01075           // to all the other MPI ranks.
01076           Teuchos::broadcast (*pComm, 0, dims);
01077         }
01078         else {
01079           // Perhaps in tolerant mode, we could set all the
01080           // dimensions to zero for now, and deduce correct
01081           // dimensions by reading all of the file's entries and
01082           // computing the max(row index) and max(column index).
01083           // However, for now we just error out in that case.
01084           TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
01085             "Error reading Matrix Market sparse matrix: failed to read "
01086             "coordinate matrix dimensions.");
01087         }
01088         return dims;
01089       }
01090 
01101       typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
01102 
01128       static RCP<adder_type>
01129       makeAdder (const RCP<const Comm<int> >& pComm,
01130                  RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
01131                  const Tuple<global_ordinal_type, 3>& dims,
01132                  const bool tolerant=false,
01133                  const bool debug=false)
01134       {
01135         if (pComm->getRank() == 0) {
01136           typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> raw_adder_type;
01137           RCP<raw_adder_type> pRaw =
01138             rcp (new raw_adder_type (dims[0], dims[1], dims[2], tolerant, debug));
01139           return rcp (new adder_type (pRaw, pBanner->symmType()));
01140         }
01141         else {
01142           return null;
01143         }
01144       }
01145 
01146     public:
01147 
01172       static sparse_matrix_ptr
01173       readSparseFile (const std::string& filename,
01174                       const RCP<const Comm<int> >& pComm,
01175                       const RCP<node_type>& pNode,
01176                       const bool callFillComplete=true,
01177                       const bool tolerant=false,
01178                       const bool debug=false)
01179       {
01180         const int myRank = Teuchos::rank (*pComm);
01181         std::ifstream in;
01182 
01183         // Only open the file on Rank 0.
01184         if (myRank == 0)
01185           in.open (filename.c_str());
01186         return readSparse (in, pComm, pNode, callFillComplete, tolerant, debug);
01187         // We can rely on the destructor of the input stream to close
01188         // the file on scope exit, even if readSparse() throws an
01189         // exception.
01190       }
01191 
01218       static sparse_matrix_ptr
01219       readSparse (std::istream& in,
01220                   const RCP<const Comm<int> >& pComm,
01221                   const RCP<node_type>& pNode,
01222                   const bool callFillComplete=true,
01223                   const bool tolerant=false,
01224                   const bool debug=false)
01225       {
01226         using Teuchos::MatrixMarket::Banner;
01227         using Teuchos::broadcast;
01228         using Teuchos::ptr;
01229         using Teuchos::reduceAll;
01230         using std::cerr;
01231         using std::endl;
01232         typedef Teuchos::ScalarTraits<scalar_type> STS;
01233 
01234         const bool extraDebug = false;
01235         const int myRank = pComm->getRank ();
01236         const int rootRank = 0;
01237 
01238         // Current line number in the input stream.  Various calls
01239         // will modify this depending on the number of lines that are
01240         // read from the input stream.  Only Rank 0 modifies this.
01241         size_t lineNumber = 1;
01242 
01243         if (debug && myRank == rootRank) {
01244           cerr << "Matrix Market reader: readSparse:" << endl
01245                << "-- Reading banner line" << endl;
01246         }
01247 
01248         // The "Banner" tells you whether the input stream represents
01249         // a sparse matrix, the symmetry type of the matrix, and the
01250         // type of the data it contains.
01251         //
01252         // pBanner will only be nonnull on MPI Rank 0.  It will be
01253         // null on all other MPI processes.
01254         RCP<const Banner> pBanner;
01255         {
01256           // We read and validate the Banner on Proc 0, but broadcast
01257           // the validation result to all processes.
01258           // Teuchos::broadcast doesn't currently work with bool, so
01259           // we use int (true -> 1, false -> 0).
01260           int bannerIsCorrect = 1;
01261           std::ostringstream errMsg;
01262 
01263           if (myRank == rootRank) {
01264             // Read the Banner line from the input stream.
01265             try {
01266               pBanner = readBanner (in, lineNumber, tolerant, debug);
01267             }
01268             catch (std::exception& e) {
01269               errMsg << "Attempt to read the Matrix Market file's Banner line "
01270                 "threw an exception: " << e.what();
01271               bannerIsCorrect = 0;
01272             }
01273 
01274             if (bannerIsCorrect) {
01275               // Validate the Banner for the case of a sparse matrix.
01276               // We validate on Proc 0, since it reads the Banner.
01277 
01278               // In intolerant mode, the matrix type must be "coordinate".
01279               if (! tolerant && pBanner->matrixType() != "coordinate") {
01280                 bannerIsCorrect = 0;
01281                 errMsg << "The Matrix Market input file must contain a "
01282                   "\"coordinate\"-format sparse matrix in order to create a "
01283                   "Tpetra::CrsMatrix object from it, but the file's matrix "
01284                   "type is \"" << pBanner->matrixType() << "\" instead.";
01285               }
01286               // In tolerant mode, we allow the matrix type to be
01287               // anything other than "array" (which would mean that
01288               // the file contains a dense matrix).
01289               if (tolerant && pBanner->matrixType() == "array") {
01290                 bannerIsCorrect = 0;
01291                 errMsg << "Matrix Market file must contain a \"coordinate\"-"
01292                   "format sparse matrix in order to create a Tpetra::CrsMatrix "
01293                   "object from it, but the file's matrix type is \"array\" "
01294                   "instead.  That probably means the file contains dense matrix "
01295                   "data.";
01296               }
01297             }
01298           } // Proc 0: Done reading the Banner, hopefully successfully.
01299 
01300           // Broadcast from Proc 0 whether the Banner was read correctly.
01301           broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
01302 
01303           // If the Banner is invalid, all processes throw an
01304           // exception.  Only Proc 0 gets the exception message, but
01305           // that's OK, since the main point is to "stop the world"
01306           // (rather than throw an exception on one process and leave
01307           // the others hanging).
01308           TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
01309             std::invalid_argument, errMsg.str ());
01310         } // Done reading the Banner line and broadcasting success.
01311         if (debug && myRank == rootRank) {
01312           cerr << "-- Reading dimensions line" << endl;
01313         }
01314 
01315         // Read the matrix dimensions from the Matrix Market metadata.
01316         // dims = (numRows, numCols, numEntries).  Proc 0 does the
01317         // reading, but it broadcasts the results to all MPI
01318         // processes.  Thus, readCoordDims() is a collective
01319         // operation.  It does a collective check for correctness too.
01320         Tuple<global_ordinal_type, 3> dims =
01321           readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
01322 
01323         if (debug && myRank == rootRank) {
01324           cerr << "-- Making Adder for collecting matrix data" << endl;
01325         }
01326 
01327         // "Adder" object for collecting all the sparse matrix entries
01328         // from the input stream.  This is only nonnull on Proc 0.
01329         RCP<adder_type> pAdder =
01330           makeAdder (pComm, pBanner, dims, tolerant, debug);
01331 
01332         if (debug && myRank == rootRank) {
01333           cerr << "-- Reading matrix data" << endl;
01334         }
01335         //
01336         // Read the matrix entries from the input stream on Proc 0.
01337         //
01338         {
01339           // We use readSuccess to broadcast the results of the read
01340           // (succeeded or not) to all MPI processes.  Since
01341           // Teuchos::broadcast doesn't currently know how to send
01342           // bools, we convert to int (true -> 1, false -> 0).
01343           int readSuccess = 1;
01344           std::ostringstream errMsg; // Exception message (only valid on Proc 0)
01345           if (myRank == rootRank) {
01346             try {
01347               // Reader for "coordinate" format sparse matrix data.
01348               typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
01349                 global_ordinal_type, scalar_type, STS::isComplex> reader_type;
01350               reader_type reader (pAdder);
01351 
01352               // Read the sparse matrix entries.
01353               std::pair<bool, std::vector<size_t> > results =
01354                 reader.read (in, lineNumber, tolerant, debug);
01355               readSuccess = results.first ? 1 : 0;
01356             }
01357             catch (std::exception& e) {
01358               readSuccess = 0;
01359               errMsg << e.what();
01360             }
01361           }
01362           broadcast (*pComm, rootRank, ptr (&readSuccess));
01363 
01364           // It would be nice to add a "verbose" flag, so that in
01365           // tolerant mode, we could log any bad line number(s) on
01366           // Proc 0.  For now, we just throw if the read fails to
01367           // succeed.
01368           //
01369           // Question: If we're in tolerant mode, and if the read did
01370           // not succeed, should we attempt to call fillComplete()?
01371           TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
01372             "Failed to read the Matrix Market sparse matrix file: "
01373             << errMsg.str());
01374         } // Done reading the matrix entries (stored on Proc 0 for now)
01375 
01376         if (debug && myRank == rootRank) {
01377           cerr << "-- Successfully read the Matrix Market data" << endl;
01378         }
01379 
01380         // In tolerant mode, we need to rebroadcast the matrix
01381         // dimensions, since they may be different after reading the
01382         // actual matrix data.  We only need to broadcast the number
01383         // of rows and columns.  Only Rank 0 needs to know the actual
01384         // global number of entries, since (a) we need to merge
01385         // duplicates on Rank 0 first anyway, and (b) when we
01386         // distribute the entries, each rank other than Rank 0 will
01387         // only need to know how many entries it owns, not the total
01388         // number of entries.
01389         if (tolerant) {
01390           if (debug && myRank == rootRank) {
01391             cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
01392                  << endl
01393                  << "----- Dimensions before: "
01394                  << dims[0] << " x " << dims[1]
01395                  << endl;
01396           }
01397           // Packed coordinate matrix dimensions (numRows, numCols).
01398           Tuple<global_ordinal_type, 2> updatedDims;
01399           if (myRank == rootRank) {
01400             // If one or more bottom rows of the matrix contain no
01401             // entries, then the Adder will report that the number
01402             // of rows is less than that specified in the
01403             // metadata.  We allow this case, and favor the
01404             // metadata so that the zero row(s) will be included.
01405             updatedDims[0] =
01406               std::max (dims[0], pAdder->getAdder()->numRows());
01407             updatedDims[1] = pAdder->getAdder()->numCols();
01408           }
01409           broadcast (*pComm, rootRank, updatedDims);
01410           dims[0] = updatedDims[0];
01411           dims[1] = updatedDims[1];
01412           if (debug && myRank == rootRank) {
01413             cerr << "----- Dimensions after: " << dims[0] << " x "
01414                  << dims[1] << endl;
01415           }
01416         }
01417         else {
01418           // In strict mode, we require that the matrix's metadata and
01419           // its actual data agree, at least somewhat.  In particular,
01420           // the number of rows must agree, since otherwise we cannot
01421           // distribute the matrix correctly.
01422 
01423           // Teuchos::broadcast() doesn't know how to broadcast bools,
01424           // so we use an int with the standard 1 == true, 0 == false
01425           // encoding.
01426           int dimsMatch = 1;
01427           if (myRank == rootRank) {
01428             // If one or more bottom rows of the matrix contain no
01429             // entries, then the Adder will report that the number of
01430             // rows is less than that specified in the metadata.  We
01431             // allow this case, and favor the metadata, but do not
01432             // allow the Adder to think there are more rows in the
01433             // matrix than the metadata says.
01434             if (dims[0] < pAdder->getAdder ()->numRows ()) {
01435               dimsMatch = 0;
01436             }
01437           }
01438           broadcast (*pComm, 0, ptr (&dimsMatch));
01439           if (dimsMatch == 0) {
01440             // We're in an error state anyway, so we might as well
01441             // work a little harder to print an informative error
01442             // message.
01443             //
01444             // Broadcast the Adder's idea of the matrix dimensions
01445             // from Proc 0 to all processes.
01446             Tuple<global_ordinal_type, 2> addersDims;
01447             if (myRank == rootRank) {
01448               addersDims[0] = pAdder->getAdder()->numRows();
01449               addersDims[1] = pAdder->getAdder()->numCols();
01450             }
01451             broadcast (*pComm, 0, addersDims);
01452             TEUCHOS_TEST_FOR_EXCEPTION(
01453               dimsMatch == 0, std::runtime_error,
01454               "The matrix metadata says that the matrix is " << dims[0] << " x "
01455               << dims[1] << ", but the actual data says that the matrix is "
01456               << addersDims[0] << " x " << addersDims[1] << ".  That means the "
01457               "data includes more rows than reported in the metadata.  This "
01458               "is not allowed when parsing in strict mode.  Parse the matrix in "
01459               "tolerant mode to ignore the metadata when it disagrees with the "
01460               "data.");
01461           }
01462         } // Matrix dimensions (# rows, # cols, # entries) agree.
01463 
01464         if (debug && myRank == rootRank) {
01465           cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
01466         }
01467 
01468         // Now that we've read in all the matrix entries from the
01469         // input stream into the adder on Proc 0, post-process them
01470         // into CSR format (still on Proc 0).  This will facilitate
01471         // distributing them to all the processors.
01472         //
01473         // These arrays represent the global matrix data as a CSR
01474         // matrix (with numEntriesPerRow as redundant but convenient
01475         // metadata, since it's computable from rowPtr and vice
01476         // versa).  They are valid only on Proc 0.
01477         ArrayRCP<size_t> numEntriesPerRow;
01478         ArrayRCP<size_t> rowPtr;
01479         ArrayRCP<global_ordinal_type> colInd;
01480         ArrayRCP<scalar_type> values;
01481 
01482         // Proc 0 first merges duplicate entries, and then converts
01483         // the coordinate-format matrix data to CSR.
01484         {
01485           int mergeAndConvertSucceeded = 1;
01486           std::ostringstream errMsg;
01487 
01488           if (myRank == rootRank) {
01489             try {
01490               typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
01491                 global_ordinal_type> element_type;
01492               typedef typename std::vector<element_type>::const_iterator iter_type;
01493 
01494               // Number of rows in the matrix.  If we are in tolerant
01495               // mode, we've already synchronized dims with the actual
01496               // matrix data.  If in strict mode, we should use dims
01497               // (as read from the file's metadata) rather than the
01498               // matrix data to determine the dimensions.  (The matrix
01499               // data will claim fewer rows than the metadata, if one
01500               // or more rows have no entries stored in the file.)
01501               const size_type numRows = dims[0];
01502 
01503               // Additively merge duplicate matrix entries.
01504               pAdder->getAdder()->merge ();
01505 
01506               // Get a temporary const view of the merged matrix entries.
01507               const std::vector<element_type>& entries =
01508                 pAdder->getAdder()->getEntries();
01509 
01510               // Number of matrix entries (after merging).
01511               const size_t numEntries = (size_t)entries.size();
01512 
01513               if (debug) {
01514                 cerr << "----- Proc 0: Matrix has numRows=" << numRows
01515                      << " rows and numEntries=" << numEntries
01516                      << " entries." << endl;
01517               }
01518 
01519               // Make space for the CSR matrix data.  Converting to
01520               // CSR is easier if we fill numEntriesPerRow with zeros
01521               // at first.
01522               numEntriesPerRow = arcp<size_t> (numRows);
01523               std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
01524               rowPtr = arcp<size_t> (numRows+1);
01525               std::fill (rowPtr.begin(), rowPtr.end(), 0);
01526               colInd = arcp<global_ordinal_type> (numEntries);
01527               values = arcp<scalar_type> (numEntries);
01528 
01529               // Convert from array-of-structs coordinate format to CSR
01530               // (compressed sparse row) format.
01531               global_ordinal_type prvRow = 0;
01532               size_t curPos = 0;
01533               rowPtr[0] = 0;
01534               for (curPos = 0; curPos < numEntries; ++curPos) {
01535                 const element_type& curEntry = entries[curPos];
01536                 const global_ordinal_type curRow = curEntry.rowIndex();
01537                 TEUCHOS_TEST_FOR_EXCEPTION(
01538                   curRow < prvRow, std::logic_error,
01539                   "Row indices are out of order, even though they are supposed "
01540                   "to be sorted.  curRow = " << curRow << ", prvRow = "
01541                   << prvRow << ", at curPos = " << curPos << ".  Please report "
01542                   "this bug to the Tpetra developers.");
01543                 if (curRow > prvRow) {
01544                   for (size_type r = prvRow+1; r <= curRow; ++r) {
01545                     rowPtr[r] = curPos;
01546                   }
01547                   prvRow = curRow;
01548                 }
01549                 numEntriesPerRow[curRow]++;
01550                 colInd[curPos] = curEntry.colIndex();
01551                 values[curPos] = curEntry.value();
01552               }
01553               // rowPtr has one more entry than numEntriesPerRow.  The
01554               // last entry of rowPtr is the number of entries in
01555               // colInd and values.
01556               rowPtr[numRows] = numEntries;
01557             } // Finished conversion to CSR format
01558             catch (std::exception& e) {
01559               mergeAndConvertSucceeded = 0;
01560               errMsg << "Failed to merge sparse matrix entries and convert to "
01561                 "CSR format: " << e.what();
01562             }
01563 
01564             if (debug && mergeAndConvertSucceeded) {
01565               // Number of rows in the matrix.
01566               const size_type numRows = dims[0];
01567               const size_type maxToDisplay = 100;
01568 
01569               cerr << "----- Proc 0: numEntriesPerRow[0.."
01570                    << (numEntriesPerRow.size()-1) << "] ";
01571               if (numRows > maxToDisplay) {
01572                 cerr << "(only showing first and last few entries) ";
01573               }
01574               cerr << "= [";
01575               if (numRows > 0) {
01576                 if (numRows > maxToDisplay) {
01577                   for (size_type k = 0; k < 2; ++k) {
01578                     cerr << numEntriesPerRow[k] << " ";
01579                   }
01580                   cerr << "... ";
01581                   for (size_type k = numRows-2; k < numRows-1; ++k) {
01582                     cerr << numEntriesPerRow[k] << " ";
01583                   }
01584                 }
01585                 else {
01586                   for (size_type k = 0; k < numRows-1; ++k) {
01587                     cerr << numEntriesPerRow[k] << " ";
01588                   }
01589                 }
01590                 cerr << numEntriesPerRow[numRows-1];
01591               } // numRows > 0
01592               cerr << "]" << endl;
01593 
01594               cerr << "----- Proc 0: rowPtr ";
01595               if (numRows > maxToDisplay) {
01596                 cerr << "(only showing first and last few entries) ";
01597               }
01598               cerr << "= [";
01599               if (numRows > maxToDisplay) {
01600                 for (size_type k = 0; k < 2; ++k) {
01601                   cerr << rowPtr[k] << " ";
01602                 }
01603                 cerr << "... ";
01604                 for (size_type k = numRows-2; k < numRows; ++k) {
01605                   cerr << rowPtr[k] << " ";
01606                 }
01607               }
01608               else {
01609                 for (size_type k = 0; k < numRows; ++k) {
01610                   cerr << rowPtr[k] << " ";
01611                 }
01612               }
01613               cerr << rowPtr[numRows] << "]" << endl;
01614             }
01615           } // if myRank == rootRank
01616         } // Done converting sparse matrix data to CSR format
01617 
01618         // Now we're done with the Adder, so we can release the
01619         // reference ("free" it) to save space.  This only actually
01620         // does anything on Rank 0, since pAdder is null on all the
01621         // other MPI processes.
01622         pAdder = null;
01623 
01624         if (debug && myRank == rootRank) {
01625           cerr << "-- Making range, domain, and row maps" << endl;
01626         }
01627 
01628         // Make the maps that describe the matrix's range and domain,
01629         // and the distribution of its rows.  Creating a Map is a
01630         // collective operation, so we don't have to do a broadcast of
01631         // a success Boolean.
01632         map_ptr pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
01633         map_ptr pDomainMap = makeDomainMap (pRangeMap, dims[0], dims[1]);
01634         map_ptr pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
01635 
01636         if (debug && myRank == rootRank) {
01637           cerr << "-- Distributing the matrix data" << endl;
01638         }
01639 
01640         // Distribute the matrix data.  Each processor has to add the
01641         // rows that it owns.  If you try to make Proc 0 call
01642         // insertGlobalValues() for _all_ the rows, not just those it
01643         // owns, then fillComplete() will compute the number of
01644         // columns incorrectly.  That's why Proc 0 has to distribute
01645         // the matrix data and why we make all the processors (not
01646         // just Proc 0) call insertGlobalValues() on their own data.
01647         //
01648         // These arrays represent each processor's part of the matrix
01649         // data, in "CSR" format (sort of, since the row indices might
01650         // not be contiguous).
01651         ArrayRCP<size_t> myNumEntriesPerRow;
01652         ArrayRCP<size_t> myRowPtr;
01653         ArrayRCP<global_ordinal_type> myColInd;
01654         ArrayRCP<scalar_type> myValues;
01655         // Distribute the matrix data.  This is a collective operation.
01656         distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
01657                     numEntriesPerRow, rowPtr, colInd, values, debug);
01658 
01659         if (debug && myRank == rootRank) {
01660           cerr << "-- Inserting matrix entries on each processor";
01661           if (callFillComplete) {
01662             cerr << " and calling fillComplete()";
01663           }
01664           cerr << endl;
01665         }
01666         // Each processor inserts its part of the matrix data, and
01667         // then they all call fillComplete().  This method invalidates
01668         // the my* distributed matrix data before calling
01669         // fillComplete(), in order to save space.  In general, we
01670         // never store more than two copies of the matrix's entries in
01671         // memory at once, which is no worse than what Tpetra
01672         // promises.
01673         sparse_matrix_ptr pMatrix =
01674           makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
01675                       pRowMap, pRangeMap, pDomainMap, callFillComplete);
01676         // Only use a reduce-all in debug mode to check if pMatrix is
01677         // null.  Otherwise, just throw an exception.  We never expect
01678         // a null pointer here, so we can save a communication.
01679         if (debug) {
01680           int localIsNull = pMatrix.is_null () ? 1 : 0;
01681           int globalIsNull = 0;
01682           reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
01683           TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
01684             "Reader::makeMatrix() returned a null pointer on at least one "
01685             "process.  Please report this bug to the Tpetra developers.");
01686         }
01687         else {
01688           TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
01689             "Reader::makeMatrix() returned a null pointer.  "
01690             "Please report this bug to the Tpetra developers.");
01691         }
01692 
01693         // We can't get the dimensions of the matrix until after
01694         // fillComplete() is called.  Thus, we can't do the sanity
01695         // check (dimensions read from the Matrix Market data,
01696         // vs. dimensions reported by the CrsMatrix) unless the user
01697         // asked makeMatrix() to call fillComplete().
01698         //
01699         // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
01700         // what one might think it does, so you have to ask the range
01701         // resp. domain map for the number of rows resp. columns.
01702         if (callFillComplete) {
01703           const int numProcs = pComm->getSize ();
01704 
01705           if (extraDebug && debug) {
01706             const global_size_t globalNumRows =
01707               pRangeMap->getGlobalNumElements();
01708             const global_size_t globalNumCols =
01709               pDomainMap->getGlobalNumElements();
01710             if (myRank == rootRank) {
01711               cerr << "-- Matrix is "
01712                    << globalNumRows << " x " << globalNumCols
01713                    << " with " << pMatrix->getGlobalNumEntries()
01714                    << " entries, and index base "
01715                    << pMatrix->getIndexBase() << "." << endl;
01716             }
01717             pComm->barrier ();
01718             for (int p = 0; p < numProcs; ++p) {
01719               if (myRank == p) {
01720                 cerr << "-- Proc " << p << " owns "
01721                      << pMatrix->getNodeNumCols() << " columns, and "
01722                      << pMatrix->getNodeNumEntries() << " entries." << endl;
01723               }
01724               pComm->barrier ();
01725             }
01726           } // if (extraDebug && debug)
01727         } // if (callFillComplete)
01728 
01729         if (debug && myRank == rootRank) {
01730           cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
01731                << endl;
01732         }
01733         return pMatrix;
01734       }
01735 
01765       static RCP<multivector_type>
01766       readDenseFile (const std::string& filename,
01767                      const RCP<const comm_type>& pComm,
01768                      const RCP<node_type>& pNode,
01769                      RCP<const map_type>& pMap,
01770                      const bool tolerant=false,
01771                      const bool debug=false)
01772       {
01773         const int myRank = Teuchos::rank (*pComm);
01774         std::ifstream in;
01775 
01776         // Only open the file on Rank 0.
01777         if (myRank == 0)
01778           in.open (filename.c_str());
01779         return readDense (in, pComm, pNode, pMap, tolerant, debug);
01780         // We can rely on the destructor of the input stream to close
01781         // the file on scope exit, even if readSparse() throws an
01782         // exception.
01783       }
01784 
01849       static RCP<multivector_type>
01850       readDense (std::istream& in,
01851                  const RCP<const comm_type>& pComm,
01852                  const RCP<node_type>& pNode,
01853                  RCP<const map_type>& pMap,
01854                  const bool tolerant=false,
01855                  const bool debug=false)
01856       {
01857         using Teuchos::MatrixMarket::Banner;
01858         using Teuchos::MatrixMarket::checkCommentLine;
01859         using std::cerr;
01860         using std::endl;
01861 
01862         // Abbreviations for typedefs, to make the code more concise.
01863         typedef scalar_type S;
01864         typedef local_ordinal_type LO;
01865         typedef global_ordinal_type GO;
01866         typedef node_type Node;
01867         typedef Teuchos::ScalarTraits<S> STS;
01868         typedef typename STS::magnitudeType M;
01869         typedef Teuchos::ScalarTraits<M> STM;
01870 
01871         // Rank 0 is the only (MPI) process allowed to read from the
01872         // input stream.
01873         const int myRank = pComm->getRank ();
01874 
01875         if (debug && myRank == 0) {
01876           cerr << "Matrix Market reader: readDense:" << endl;
01877         }
01878 
01879         // If pMap is nonnull, check the precondition that its
01880         // communicator resp. node equal pComm resp. pNode.  Checking
01881         // now avoids doing a lot of file reading before we detect the
01882         // violated precondition.
01883         TEUCHOS_TEST_FOR_EXCEPTION(! pMap.is_null() &&
01884                                    (pMap->getComm() != pComm ||
01885                                     pMap->getNode() != pNode),
01886           std::invalid_argument, "If you supply a nonnull Map, the Map's "
01887           "communicator and node must equal the supplied communicator resp. "
01888           "node.");
01889 
01890         // Rank 0 will read in the matrix dimensions from the file,
01891         // and broadcast them to all ranks in the given communicator.
01892         // There are only 2 dimensions in the matrix, but we use the
01893         // third element of the Tuple to encode the banner's reported
01894         // data type: "real" == 0, "complex" == 1, and "integer" == 0
01895         // (same as "real").  We don't allow pattern matrices (i.e.,
01896         // graphs) since they only make sense for sparse data.
01897         Tuple<GO, 3> dims;
01898         dims[0] = 0;
01899         dims[1] = 0;
01900 
01901         // Current line number in the input stream.  Only valid on
01902         // Rank 0.  Various calls will modify this depending on the
01903         // number of lines that are read from the input stream.
01904         size_t lineNumber = 1;
01905 
01906         // Only Rank 0 gets to read matrix data from the input stream.
01907         if (myRank == 0) {
01908           if (debug && myRank == 0) {
01909             cerr << "-- Reading banner line (dense)" << endl;
01910           }
01911 
01912           // The "Banner" tells you whether the input stream
01913           // represents a dense matrix, the symmetry type of the
01914           // matrix, and the type of the data it contains.
01915           RCP<const Banner> pBanner =
01916             readBanner (in, lineNumber, tolerant, debug);
01917           TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "array",
01918             std::invalid_argument, "The Matrix Market file does not contain "
01919             "dense matrix data.  Its banner (first) line says that its matrix "
01920             "type is \"" << pBanner->matrixType() << "\", rather than the "
01921             "required \"array\".");
01922           TEUCHOS_TEST_FOR_EXCEPTION(pBanner->dataType() == "pattern",
01923             std::invalid_argument, "The Matrix Market file's banner (first) "
01924             "line claims that the matrix's data type is \"pattern\".  This does "
01925             "not make sense for a dense matrix, yet the file reports the matrix "
01926             "as dense.  The only valid data types for a dense matrix are "
01927             "\"real\", \"complex\", and \"integer\".");
01928           // Encode the data type reported by the Banner as the third
01929           // element of the dimensions Tuple.
01930           dims[2] = encodeDataType (pBanner->dataType());
01931 
01932           if (debug && myRank == 0) {
01933             cerr << "-- Reading dimensions line (dense)" << endl;
01934           }
01935           // Keep reading lines from the input stream until we find a
01936           // non-comment line, or until we run out of lines.  The
01937           // latter is an error, since every "array" format Matrix
01938           // Market file must have a dimensions line after the banner
01939           // (even if the matrix has zero rows or columns, or zero
01940           // entries).
01941           std::string line;
01942           bool commentLine = true;
01943           while (commentLine) {
01944             // Test whether it is even valid to read from the
01945             // input stream wrapping the line.
01946             TEUCHOS_TEST_FOR_EXCEPTION(in.eof() || in.fail(), std::runtime_error,
01947               "Unable to get array dimensions line (at all) from line "
01948               << lineNumber << " of input stream.  The input stream claims that "
01949               "it is " << (in.eof() ? "at end-of-file." : "in a failed state."));
01950             // Try to get the next line from the input stream.
01951             if (getline(in, line)) {
01952               ++lineNumber; // We did actually read a line.
01953             }
01954             // Is the current line a comment line?  Ignore start and
01955             // size; they are only useful for reading the actual
01956             // matrix entries.  (We could use them here as an
01957             // optimization, but we've chosen not to.)
01958             size_t start = 0, size = 0;
01959             commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
01960           }
01961           //
01962           // Read in <numRows> <numCols> from input line.
01963           //
01964           std::istringstream istr (line);
01965 
01966           // Test whether it is even valid to read from the input
01967           // stream wrapping the line.
01968           TEUCHOS_TEST_FOR_EXCEPTION(istr.eof() || istr.fail(), std::runtime_error,
01969                                      "Unable to read any data from line " << lineNumber
01970                                      << " of input; the line should contain the matrix "
01971                                      << "dimensions \"<numRows> <numCols>\".");
01972           // Read in <numRows>.
01973           {
01974             GO theNumRows = 0;
01975             istr >> theNumRows;
01976             TEUCHOS_TEST_FOR_EXCEPTION(istr.fail(), std::runtime_error,
01977                                        "Failed to get number of rows from line "
01978                                        << lineNumber << " of input; the line should "
01979                                        "contain the matrix dimensions \"<numRows> "
01980                                        "<numCols>\".");
01981             // Capture the validly read result before checking for eof.
01982             dims[0] = theNumRows;
01983           }
01984           // There should be one more thing to read.
01985           TEUCHOS_TEST_FOR_EXCEPTION(istr.eof(), std::runtime_error,
01986                                      "No more data after number of rows on line "
01987                                      << lineNumber << " of input; the line should "
01988                                      "contain the matrix dimensions \"<numRows> "
01989                                      "<numCols>\".");
01990           // Read in <numCols>
01991           {
01992             GO theNumCols = 0;
01993             istr >> theNumCols;
01994             TEUCHOS_TEST_FOR_EXCEPTION(istr.fail(), std::runtime_error,
01995                                        "Failed to get number of columns from line "
01996                                        << lineNumber << " of input; the line should "
01997                                        "contain the matrix dimensions \"<numRows> "
01998                                        "<numCols>\".");
01999             // Capture the validly read result.
02000             dims[1] = theNumCols;
02001           }
02002         } // if (myRank == 0)
02003 
02004         // Broadcast matrix dimensions and the encoded data type from
02005         // MPI Rank 0 to all the MPI processes.
02006         Teuchos::broadcast (*pComm, 0, dims);
02007 
02008         // Tpetra objects want the matrix dimensions in these types.
02009         const global_size_t numRows = static_cast<global_size_t> (dims[0]);
02010         const size_t numCols = static_cast<size_t> (dims[1]);
02011 
02012         // Make an "everything on Rank 0" map pRank0Map.  We'll use
02013         // this to construct a "distributed" vector X owned entirely
02014         // by Rank 0.  Rank 0 will then read all the matrix entries
02015         // and put them in X.
02016         RCP<const map_type> pRank0Map =
02017           createContigMapWithNode<LO, GO, Node> (numRows,
02018                                                  (myRank == 0 ? numRows : 0),
02019                                                  pComm, pNode);
02020 
02021         // Make a multivector X owned entirely by Rank 0.
02022         RCP<multivector_type> X =
02023           createMultiVector<S, LO, GO, Node> (pRank0Map, numCols);
02024 
02025         //
02026         // On Rank 0, read the Matrix Market data from the input
02027         // stream into the multivector X.
02028         //
02029         if (myRank == 0) {
02030           if (debug && myRank == 0) {
02031             cerr << "-- Reading matrix data (dense)" << endl;
02032           }
02033 
02034           // Make sure that we can get a 1-D view of X.
02035           TEUCHOS_TEST_FOR_EXCEPTION(! X->isConstantStride (), std::logic_error,
02036                                      "Can't get a 1-D view of the entries of the "
02037                                      "MultiVector X on Rank 0, because the stride "
02038                                      "between the columns of X is not constant.  "
02039                                      "This shouldn't happen because we just created "
02040                                      "X and haven't filled it in yet.  Please report "
02041                                      "this bug to the Tpetra developers.");
02042 
02043           // Get a writeable 1-D view of the entries of X.  Rank 0
02044           // owns all of them.  The view will expire at the end of
02045           // scope, so (if necessary) it will be written back to X
02046           // at this time.
02047           ArrayRCP<S> X_view = X->get1dViewNonConst ();
02048           TEUCHOS_TEST_FOR_EXCEPTION(
02049             static_cast<global_size_t> (X_view.size()) < numRows * numCols,
02050             std::logic_error,
02051             "The view of X has size " << X_view << " which is not enough to "
02052             "accommodate the expected number of entries numRows*numCols = "
02053             << numRows << "*" << numCols << " = " << numRows*numCols << ".  "
02054             "Please report this bug to the Tpetra developers.");
02055           const size_t stride = X->getStride ();
02056 
02057           // The third element of the dimensions Tuple encodes the data
02058           // type reported by the Banner: "real" == 0, "complex" == 1,
02059           // "integer" == 0 (same as "real"), "pattern" == 2.  We do not
02060           // allow dense matrices to be pattern matrices, so dims[2] ==
02061           // 0 or 1.  We've already checked for this above.
02062           const bool isComplex = (dims[2] == 1);
02063           typename ArrayRCP<S>::size_type count = 0, curRow = 0, curCol = 0;
02064 
02065           std::string line;
02066           while (getline (in, line)) {
02067             ++lineNumber;
02068             // Is the current line a comment line?  If it's not,
02069             // line.substr(start,size) contains the data.
02070             size_t start = 0, size = 0;
02071             const bool commentLine =
02072               checkCommentLine (line, start, size, lineNumber, tolerant);
02073             if (! commentLine) {
02074               // Make sure we have room in which to put the new matrix
02075               // entry.  We check this only after checking for a
02076               // comment line, because there may be one or more
02077               // comment lines at the end of the file.  In tolerant
02078               // mode, we simply ignore any extra data.
02079               if (count >= X_view.size()) {
02080                 if (tolerant) {
02081                   break;
02082                 }
02083                 else {
02084                   TEUCHOS_TEST_FOR_EXCEPTION(
02085                      count >= X_view.size(),
02086                      std::runtime_error,
02087                      "The Matrix Market input stream has more data in it than "
02088                      "its metadata reported.  Current line number is "
02089                      << lineNumber << ".");
02090                 }
02091               }
02092               std::istringstream istr (line.substr (start, size));
02093               // Does the line contain anything at all?  Can we
02094               // safely read from the input stream wrapping the
02095               // line?
02096               if (istr.eof() || istr.fail()) {
02097                 // In tolerant mode, simply ignore the line.
02098                 if (tolerant) {
02099                   break;
02100                 }
02101                 // We repeat the full test here so the exception
02102                 // message is more informative.
02103                 TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && (istr.eof() || istr.fail()),
02104                   std::runtime_error, "Line " << lineNumber << " of the Matrix "
02105                   "Market file is empty, or we cannot read from it for some "
02106                   "other reason.");
02107               }
02108               // Current matrix entry to read in.
02109               S val = STS::zero();
02110               // Real and imaginary parts of the current matrix entry.
02111               // The imaginary part is zero if the matrix is
02112               // real-valued.
02113               M real = STM::zero(), imag = STM::zero();
02114 
02115               // isComplex refers to the input stream's data, not to
02116               // the scalar type S.  It's OK to read real-valued data
02117               // into a matrix storing complex-valued data; in that
02118               // case, all entries' imaginary parts are zero.
02119               if (isComplex) {
02120                 // STS::real() and STS::imag() return a copy of their
02121                 // respective components, not a writeable reference.
02122                 // Otherwise we could just assign to them using the
02123                 // istream extraction operator (>>).  That's why we
02124                 // have separate magnitude type "real" and "imag"
02125                 // variables.
02126 
02127                 // Attempt to read the real part of the current matrix
02128                 // entry.
02129                 istr >> real;
02130                 if (istr.fail()) {
02131                   TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
02132                                              std::runtime_error,
02133                                              "Failed to get real part of a "
02134                                              "complex-valued matrix entry "
02135                                              "from line " << lineNumber
02136                                              << " of Matrix Market file.");
02137                   // In tolerant mode, just skip bad lines.
02138                   if (tolerant) {
02139                     break;
02140                   }
02141                 }
02142                 else if (istr.eof()) {
02143                   TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.eof(),
02144                                              std::runtime_error,
02145                                              "Missing imaginary part of a "
02146                                              "complex-valued matrix entry "
02147                                              "on line " << lineNumber
02148                                              << " of Matrix Market file.");
02149                   // In tolerant mode, let any missing imaginary part
02150                   // be 0.
02151                 }
02152                 else {
02153                   // Attempt to read the imaginary part of the current
02154                   // matrix entry.
02155                   istr >> imag;
02156                   TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
02157                                              std::runtime_error,
02158                                              "Failed to get imaginary part "
02159                                              "of a complex-valued matrix "
02160                                              "entry from line " << lineNumber
02161                                              << " of Matrix Market file.");
02162                   // In tolerant mode, let any missing
02163                   // or corrupted imaginary part be 0.
02164                 }
02165               }
02166               else { // Matrix Market file contains real-valued data.
02167                 // Attempt to read the current matrix entry.
02168                 istr >> real;
02169                 TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
02170                                            std::runtime_error,
02171                                            "Failed to get a real-valued matrix "
02172                                            "entry from line " << lineNumber
02173                                            << " of Matrix Market file.");
02174                 // In tolerant mode, simply ignore the line if
02175                 // we failed to read a matrix entry.
02176                 if (istr.fail() && tolerant) {
02177                   break;
02178                 }
02179               }
02180 
02181               // In tolerant mode, we simply let pass through whatever
02182               // data we got.
02183               TEUCHOS_TEST_FOR_EXCEPTION(! tolerant && istr.fail(),
02184                                          std::runtime_error,
02185                                          "Failed to read matrix data from line "
02186                                          << lineNumber << " of the Matrix Market "
02187                                          "file.");
02188 
02189               // val = S(real, imag), without potential badness
02190               // if S is a real type.
02191               Teuchos::MatrixMarket::details::assignScalar<S> (val, real, imag);
02192 
02193               curRow = count % numRows;
02194               curCol = count / numRows;
02195               X_view[curRow + curCol*stride] = val;
02196               ++count;
02197             } // if not a comment line
02198           } // while there are still lines in the file, get the next one
02199 
02200           TEUCHOS_TEST_FOR_EXCEPTION(! tolerant &&
02201                                      static_cast<global_size_t> (count) < numRows * numCols,
02202                                      std::runtime_error,
02203                                      "The Matrix Market metadata reports that the "
02204                                      "dense matrix is " << numRows <<  " x "
02205                                      << numCols << ", and thus has "
02206                                      << numRows*numCols << " total entries, but we "
02207                                      "only managed to find " << count << " entr"
02208                                      << (count == 1 ? "y" : "ies")
02209                                      << " in the file.");
02210         } // if (myRank == 0)
02211 
02212         // If there's only one MPI process and the user didn't supply
02213         // a Map (i.e., pMap is null), we're done.  Set pMap to the
02214         // Map used to distribute X, and return X.
02215         if (pComm->getSize() == 1 && pMap.is_null()) {
02216           pMap = pRank0Map;
02217           if (debug && myRank == 0) {
02218             cerr << "-- Done reading multivector" << endl;
02219           }
02220           return X;
02221         }
02222 
02223         if (debug && myRank == 0) {
02224           cerr << "-- Creating distributed Map and target MultiVector" << endl;
02225         }
02226 
02227         // If pMap is null, make a distributed map.  We've already
02228         // checked the preconditions above, in the case that pMap was
02229         // not null on input to this method.
02230         if (pMap.is_null()) {
02231           pMap = createUniformContigMapWithNode<LO, GO, Node> (numRows, pComm, pNode);
02232         }
02233         // Make a multivector Y with the distributed map pMap.
02234         RCP<multivector_type> Y = createMultiVector<S, LO, GO, Node> (pMap, numCols);
02235 
02236         // Make an Export object that will export X to Y.  First
02237         // argument is the source map, second argument is the target
02238         // map.
02239         Export<LO, GO, Node> exporter (pRank0Map, pMap);
02240 
02241         if (debug && myRank == 0) {
02242           cerr << "-- Exporting from X (owned by Rank 0) to globally owned Y"
02243                << endl;
02244         }
02245         // Export X into Y.
02246         Y->doExport (*X, exporter, INSERT);
02247 
02248         if (debug && myRank == 0) {
02249           cerr << "-- Done reading multivector" << endl;
02250         }
02251 
02252         // Y is distributed over all process(es) in the communicator.
02253         return Y;
02254       }
02255 
02256     private:
02257 
02268       static int
02269       encodeDataType (const std::string& dataType)
02270       {
02271         if (dataType == "real" || dataType == "integer") {
02272           return 0;
02273         }
02274         else if (dataType == "complex") {
02275           return 1;
02276         }
02277         else if (dataType == "pattern") {
02278           return 2;
02279         }
02280         else {
02281           // We should never get here, since Banner validates the
02282           // reported data type and ensures it is one of the accepted
02283           // values.
02284           TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
02285             "Unrecognized Matrix Market data type \"" << dataType << "\".");
02286         }
02287       }
02288     };
02289 
02329     template<class SparseMatrixType>
02330     class Writer {
02331     public:
02332       typedef SparseMatrixType sparse_matrix_type;
02333       typedef RCP<sparse_matrix_type> sparse_matrix_ptr;
02334 
02337       typedef typename SparseMatrixType::scalar_type scalar_type;
02340       typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
02347       typedef typename SparseMatrixType::global_ordinal_type global_ordinal_type;
02350       typedef typename SparseMatrixType::node_type node_type;
02351 
02354       typedef MultiVector<scalar_type,
02355                           local_ordinal_type,
02356                           global_ordinal_type,
02357                           node_type> multivector_type;
02360       typedef Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
02361 
02393       static void
02394       writeSparseFile (const std::string& filename,
02395                        const RCP<const sparse_matrix_type>& pMatrix,
02396                        const std::string& matrixName,
02397                        const std::string& matrixDescription,
02398                        const bool debug=false)
02399       {
02400         const int myRank = Teuchos::rank (*(pMatrix->getComm()));
02401         std::ofstream out;
02402 
02403         // Only open the file on Rank 0.
02404         if (myRank == 0) out.open (filename.c_str());
02405         writeSparse (out, pMatrix, matrixName, matrixDescription, debug);
02406         // We can rely on the destructor of the output stream to close
02407         // the file on scope exit, even if writeSparse() throws an
02408         // exception.
02409       }
02410 
02431       static void
02432       writeSparseFile (const std::string& filename,
02433                        const RCP<const sparse_matrix_type>& pMatrix,
02434                        const bool debug=false)
02435       {
02436         writeSparseFile (filename, pMatrix, "", "", debug);
02437       }
02438 
02439     public:
02440 
02471       static void
02472       writeSparse (std::ostream& out,
02473                    const RCP<const sparse_matrix_type>& pMatrix,
02474                    const std::string& matrixName,
02475                    const std::string& matrixDescription,
02476                    const bool debug=false)
02477       {
02478         using std::cerr;
02479         using std::endl;
02480         typedef typename Teuchos::ScalarTraits<scalar_type> STS;
02481         typedef typename STS::magnitudeType magnitude_type;
02482         typedef typename Teuchos::ScalarTraits<magnitude_type> STM;
02483         typedef typename ArrayView<scalar_type>::size_type size_type;
02484 
02485         // Convenient abbreviations.
02486         typedef local_ordinal_type LO;
02487         typedef global_ordinal_type GO;
02488 
02489         // Make the output stream write floating-point numbers in
02490         // scientific notation.  It will politely put the output
02491         // stream back to its state on input, when this scope
02492         // terminates.
02493         Teuchos::MatrixMarket::details::SetScientific<scalar_type> sci (out);
02494 
02495         RCP<const Comm<int> > pComm = pMatrix->getComm();
02496         const int myRank = Teuchos::rank (*pComm);
02497         // Whether to print debugging output to stderr.
02498         const bool debugPrint = debug && myRank == 0;
02499 
02500         const global_size_t numRows =
02501           pMatrix->getRangeMap()->getGlobalNumElements();
02502         const global_size_t numCols =
02503           pMatrix->getDomainMap()->getGlobalNumElements();
02504         if (debugPrint) {
02505           cerr << "writeSparse:" << endl
02506                << "-- Input sparse matrix is:"
02507                << "---- " << numRows << " x " << numCols << " with "
02508                << pMatrix->getGlobalNumEntries() << " entries;" << endl
02509                << "---- "
02510                << (pMatrix->isGloballyIndexed() ? "Globally" : "Locally")
02511                << " indexed." << endl;
02512         }
02513         RCP<const map_type> pRowMap = pMatrix->getRowMap();
02514         if (debugPrint) {
02515           RCP<const map_type> pColMap = pMatrix->getColMap();
02516           cerr << "---- Its row map has " << pRowMap->getGlobalNumElements()
02517                << " elements." << endl
02518                << "---- Its col map has " << pColMap->getGlobalNumElements()
02519                << " elements." << endl;
02520         }
02521         // Make the "gather" row map, where Proc 0 owns all rows and
02522         // the other procs own no rows.
02523         const size_t localNumRows = (myRank == 0) ? numRows : 0;
02524         RCP<node_type> pNode = pRowMap->getNode();
02525         RCP<const map_type> pGatherRowMap =
02526           createContigMapWithNode<LO, GO, node_type> (numRows, localNumRows,
02527                                                       pComm, pNode);
02528         // Since the matrix may in general be non-square, we need to
02529         // make a column map as well.  In this case, the column map
02530         // contains all the columns of the original matrix, because we
02531         // are gathering the whole matrix onto Proc 0.
02532         const size_t localNumCols = (myRank == 0) ? numCols : 0;
02533         RCP<const map_type> pGatherColMap =
02534           createContigMapWithNode<LO, GO, node_type> (numCols, localNumCols,
02535                                                       pComm, pNode);
02536 
02537         // Current map is the source map, gather map is the target map.
02538         typedef Import<LO, GO, node_type> import_type;
02539         import_type importer (pRowMap, pGatherRowMap);
02540 
02541         // Create a new CrsMatrix to hold the result of the import.
02542         // The constructor needs a column map as well as a row map,
02543         // for the case that the matrix is not square.
02544         RCP<sparse_matrix_type> newMatrix =
02545           rcp (new sparse_matrix_type (pGatherRowMap,
02546                                        pGatherColMap,
02547                                        size_t(0)));
02548         // Import the sparse matrix onto Proc 0.
02549         newMatrix->doImport (*pMatrix, importer, INSERT);
02550 
02551         // fillComplete() needs the domain and range maps for the case
02552         // that the matrix is not square.
02553         {
02554           RCP<const map_type> pGatherDomainMap =
02555             createContigMapWithNode<LO, GO, node_type> (numCols, localNumCols,
02556                                                         pComm, pNode);
02557           RCP<const map_type> pGatherRangeMap =
02558             createContigMapWithNode<LO, GO, node_type> (numRows, localNumRows,
02559                                                         pComm, pNode);
02560           newMatrix->fillComplete (pGatherDomainMap, pGatherRangeMap);
02561         }
02562 
02563         if (debugPrint) {
02564           cerr << "-- Output sparse matrix is:"
02565                << "---- " << newMatrix->getRangeMap()->getGlobalNumElements()
02566                << " x "
02567                << newMatrix->getDomainMap()->getGlobalNumElements()
02568                << " with "
02569                << newMatrix->getGlobalNumEntries() << " entries;" << endl
02570                << "---- "
02571                << (newMatrix->isGloballyIndexed() ? "Globally" : "Locally")
02572                << " indexed." << endl
02573                << "---- Its row map has "
02574                << newMatrix->getRowMap()->getGlobalNumElements()
02575                << " elements." << endl
02576                << "---- Its col map has "
02577                << newMatrix->getColMap()->getGlobalNumElements()
02578                << " elements (may be different than the input matrix's column"
02579                << " map, if some columns of the matrix contain no entries)."
02580                << endl;
02581         }
02582 
02583         //
02584         // Print the metadata and the matrix entries on Rank 0.
02585         //
02586         if (myRank == 0) {
02587           // Print the Matrix Market banner line.  CrsMatrix stores
02588           // data nonsymmetrically ("general").  This implies that
02589           // readSparse() on a symmetrically stored input file,
02590           // followed by writeSparse() on the resulting sparse matrix,
02591           // will result in an output file with a different banner
02592           // line than the original input file.
02593           out << "%%MatrixMarket matrix coordinate "
02594               << (STS::isComplex ? "complex" : "real")
02595               << " general" << endl;
02596 
02597           // Print comments (the matrix name and / or description).
02598           if (matrixName != "")
02599             printAsComment (out, matrixName);
02600           if (matrixDescription != "")
02601             printAsComment (out, matrixDescription);
02602 
02603           // Print the Matrix Market header (# rows, # columns, #
02604           // nonzeros).  Use the range resp. domain map for the
02605           // number of rows resp. columns, since Tpetra::CrsMatrix
02606           // uses the column map for the number of columns.  That
02607           // only corresponds to the "linear-algebraic" number of
02608           // columns when the column map is 1-1, which only happens
02609           // if the matrix is (block) diagonal.
02610           out << newMatrix->getRangeMap()->getGlobalNumElements()
02611               << " "
02612               << newMatrix->getDomainMap()->getGlobalNumElements()
02613               << " "
02614               << newMatrix->getGlobalNumEntries()
02615               << endl;
02616 
02617           // Index base (0-based or 1-based) for the row map.
02618           const GO rowIndexBase = pGatherRowMap->getIndexBase();
02619           // Index base (0-based or 1-based) for the column map.
02620           const GO colIndexBase = newMatrix->getColMap()->getIndexBase();
02621 
02622           //
02623           // Print the entries of the matrix.
02624           //
02625           // newMatrix can never be globally indexed, since we
02626           // called fillComplete() on it.  We include code for both
02627           // cases (globally or locally indexed) just in case that
02628           // ever changes.
02629           if (newMatrix->isGloballyIndexed()) {
02630             for (GO globalRowIndex = pGatherRowMap->getMinAllGlobalIndex();
02631                  globalRowIndex <= pGatherRowMap->getMaxAllGlobalIndex();
02632                  ++globalRowIndex)
02633               {
02634                 ArrayView<const GO> ind;
02635                 ArrayView<const scalar_type> val;
02636 
02637                 newMatrix->getGlobalRowView (globalRowIndex, ind, val);
02638                 typename ArrayView<const GO>::const_iterator indIter =
02639                   ind.begin();
02640                 typename ArrayView<const scalar_type>::const_iterator valIter =
02641                   val.begin();
02642                 for (; indIter != ind.end() && valIter != val.end();
02643                      ++indIter, ++valIter)
02644                   {
02645                     const GO globalColIndex = *indIter;
02646 
02647                     // Matrix Market files use 1-based indices.
02648                     out << (globalRowIndex + 1 - rowIndexBase) << " "
02649                         << (globalColIndex + 1 - colIndexBase) << " ";
02650                     if (STS::isComplex)
02651                       out << STS::real(*valIter) << " " << STS::imag(*valIter);
02652                     else
02653                       out << *valIter;
02654                     out << endl;
02655                   }
02656               }
02657           }
02658           else // newMatrix is locally indexed
02659             {
02660               typedef OrdinalTraits<GO> OTG;
02661               RCP<const map_type> pColMap = newMatrix->getColMap ();
02662 
02663               for (LO localRowIndex = pGatherRowMap->getMinLocalIndex();
02664                    localRowIndex <= pGatherRowMap->getMaxLocalIndex();
02665                    ++localRowIndex)
02666                 {
02667                   // Convert from local to global row index.
02668                   const GO globalRowIndex =
02669                     pGatherRowMap->getGlobalElement (localRowIndex);
02670                   TEUCHOS_TEST_FOR_EXCEPTION(globalRowIndex == OTG::invalid(),
02671                                      std::logic_error,
02672                                      "Failed to convert \"local\" row index "
02673                                      << localRowIndex << " into a global row "
02674                                      "index.  Please report this bug to the "
02675                                      "Tpetra developers.");
02676                   ArrayView<const LO> ind;
02677                   ArrayView<const scalar_type> val;
02678 
02679                   newMatrix->getLocalRowView (localRowIndex, ind, val);
02680                   typename ArrayView<const LO>::const_iterator indIter =
02681                     ind.begin();
02682                   typename ArrayView<const scalar_type>::const_iterator
02683                     valIter = val.begin();
02684                   for (; indIter != ind.end() && valIter != val.end();
02685                        ++indIter, ++valIter)
02686                     {
02687                       // Convert from local to global index.
02688                       const GO globalColIndex =
02689                         pColMap->getGlobalElement (*indIter);
02690                       TEUCHOS_TEST_FOR_EXCEPTION(globalColIndex == OTG::invalid(),
02691                                          std::logic_error,
02692                                          "Failed to convert \"local\" column "
02693                                          "index " << *indIter << " into a "
02694                                          "global column index.  Please report "
02695                                          "this bug to the Tpetra developers.");
02696                       // Matrix Market files use 1-based indices.
02697                       out << (globalRowIndex + 1 - rowIndexBase) << " "
02698                           << (globalColIndex + 1 - colIndexBase) << " ";
02699                       if (STS::isComplex)
02700                         out << STS::real(*valIter) << " "
02701                             << STS::imag(*valIter);
02702                       else
02703                         out << *valIter;
02704                       out << endl;
02705                     }
02706                 }
02707             }
02708         }
02709       }
02710 
02733       static void
02734       writeSparse (std::ostream& out,
02735                    const RCP<const sparse_matrix_type>& pMatrix,
02736                    const bool debug=false)
02737       {
02738         writeSparse (out, pMatrix, "", "", debug);
02739       }
02740 
02767       static void
02768       writeDenseFile (const std::string& filename,
02769                       const RCP<const multivector_type>& X,
02770                       const std::string& matrixName,
02771                       const std::string& matrixDescription)
02772       {
02773         const int myRank = Teuchos::rank (*(X->getMap()->getComm()));
02774         std::ofstream out;
02775 
02776         if (myRank == 0) // Only open the file on Rank 0.
02777           out.open (filename.c_str());
02778 
02779         writeDense (out, X, matrixName, matrixDescription);
02780         // We can rely on the destructor of the output stream to close
02781         // the file on scope exit, even if writeDense() throws an
02782         // exception.
02783       }
02784 
02803       static void
02804       writeDenseFile (const std::string& filename,
02805                       const RCP<const multivector_type>& X)
02806       {
02807         writeDenseFile (filename, X, "", "");
02808       }
02809 
02836       static void
02837       writeDense (std::ostream& out,
02838                   const RCP<const multivector_type>& X,
02839                   const std::string& matrixName,
02840                   const std::string& matrixDescription)
02841       {
02842         using Teuchos::rcp;
02843         using Teuchos::rcp_const_cast;
02844         using std::endl;
02845 
02846         typedef typename Teuchos::ScalarTraits<scalar_type> STS;
02847         typedef typename STS::magnitudeType magnitude_type;
02848         typedef typename Teuchos::ScalarTraits<magnitude_type> STM;
02849         typedef typename ArrayView<scalar_type>::size_type size_type;
02850 
02851         // Convenient abbreviations.
02852         typedef local_ordinal_type LO;
02853         typedef global_ordinal_type GO;
02854         typedef node_type NT;
02855 
02856         // If true, when making the "gather" Map, use the same index
02857         // base as X's Map.  Otherwise, just use the default index
02858         // base for the gather Map.
02859         const bool keepTheSameIndexBase = false;
02860 
02861         // Make the output stream write floating-point numbers in
02862         // scientific notation.  It will politely put the output
02863         // stream back to its state on input, when this scope
02864         // terminates.
02865         Teuchos::MatrixMarket::details::SetScientific<scalar_type> sci (out);
02866 
02867         RCP<const Comm<int> > comm = X->getMap()->getComm();
02868         const int myRank = comm->getRank ();
02869         RCP<const map_type> map = X->getMap();
02870         const global_size_t numRows = map->getGlobalNumElements();
02871         // Promote to global_size_t.
02872         const global_size_t numCols = X->getNumVectors();
02873 
02874         // Make the "gather" map, where Proc 0 owns all rows of X, and
02875         // the other procs own no rows.
02876         const size_t localNumRows = (myRank == 0) ? numRows : 0;
02877         RCP<NT> node = map->getNode();
02878         RCP<const map_type> gatherMap = keepTheSameIndexBase ?
02879           rcp_const_cast<const map_type> (rcp (new map_type (numRows, localNumRows, map->getIndexBase(), comm, node))) :
02880           createContigMapWithNode<LO, GO, NT> (numRows, localNumRows, comm, node);
02881 
02882         // Create an Import object to import X's data into a
02883         // multivector Y owned entirely by Proc 0.  In the Import
02884         // constructor, X's map is the source map, and the "gather
02885         // map" is the target map.
02886         Import<LO, GO, NT> importer (map, gatherMap);
02887 
02888         // Create a new multivector Y to hold the result of the import.
02889         RCP<multivector_type> Y =
02890           createMultiVector<scalar_type, LO, GO, NT> (gatherMap, numCols);
02891 
02892         // Import the multivector onto Proc 0.
02893         Y->doImport (*X, importer, INSERT);
02894 
02895         //
02896         // Print the matrix in Matrix Market format on Rank 0.
02897         //
02898         if (myRank == 0) {
02899           // Print the Matrix Market banner line.  MultiVector
02900           // stores data nonsymmetrically, hence "general".
02901           out << "%%MatrixMarket matrix array "
02902               << (STS::isComplex ? "complex" : "real")
02903               << " general" << endl;
02904 
02905           // Print comments (the matrix name and / or description).
02906           if (matrixName != "") {
02907             printAsComment (out, matrixName);
02908           }
02909           if (matrixDescription != "") {
02910             printAsComment (out, matrixDescription);
02911           }
02912           // Print the Matrix Market dimensions header for dense matrices.
02913           out << numRows << " " << numCols << endl;
02914           //
02915           // Get a read-only view of the entries of Y.
02916           // Rank 0 owns all of the entries of Y.
02917           //
02918           // Y must have constant stride, since multivectors have
02919           // constant stride if created with a map and number of
02920           // vectors.  This should be the case even if X does not
02921           // have constant stride.  However, we check Y, just to
02922           // make sure.
02923           TEUCHOS_TEST_FOR_EXCEPTION(! Y->isConstantStride (), std::logic_error,
02924             "The multivector Y imported onto Rank 0 does not have constant "
02925             "stride.  Please report this bug to the Tpetra developers.");
02926           ArrayRCP<const scalar_type> Y_view = Y->get1dView();
02927           //
02928           // Print the entries of the matrix, in column-major order.
02929           //
02930           const global_size_t stride = Y->getStride ();
02931           for (global_size_t j = 0; j < numCols; ++j) {
02932             for (global_size_t i = 0; i < numRows; ++i) {
02933               const scalar_type Y_ij = Y_view[i + j*stride];
02934               if (STS::isComplex) {
02935                 out << STS::real(Y_ij) << " " << STS::imag(Y_ij) << endl;
02936               }
02937               else {
02938                 out << Y_ij << endl;
02939               }
02940             }
02941           }
02942         } // if (myRank == 0)
02943       }
02944 
02963       static void
02964       writeDense (std::ostream& out,
02965                   const RCP<const multivector_type>& X)
02966       {
02967         writeDense (out, X, "", "");
02968       }
02969 
02970     private:
02994       static void
02995       printAsComment (std::ostream& out, const std::string& str)
02996       {
02997         using std::endl;
02998         std::istringstream istream (str);
02999         std::string line;
03000 
03001         while (getline (istream, line)) {
03002           if (! line.empty()) {
03003             // Note that getline() doesn't store '\n', so we have to
03004             // append the endline ourselves.
03005             if (line[0] == '%') { // Line starts with a comment character.
03006               out << line << endl;
03007             }
03008             else { // Line doesn't start with a comment character.
03009               out << "%% " << line << endl;
03010             }
03011           }
03012         }
03013       }
03014     }; // class Writer
03015 
03016   } // namespace MatrixMarket
03017 } // namespace Tpetra
03018 
03019 #endif // __MatrixMarket_Tpetra_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines