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