Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MatrixMarket_Tpetra.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 
42 #ifndef __MatrixMarket_Tpetra_hpp
43 #define __MatrixMarket_Tpetra_hpp
44 
57 #include "Tpetra_CrsMatrix.hpp"
58 #include "Tpetra_Operator.hpp"
59 #include "Tpetra_Vector.hpp"
61 #include "Teuchos_MatrixMarket_Raw_Adder.hpp"
62 #include "Teuchos_MatrixMarket_Raw_Graph_Adder.hpp"
63 #include "Teuchos_MatrixMarket_SymmetrizingAdder.hpp"
64 #include "Teuchos_MatrixMarket_SymmetrizingGraphAdder.hpp"
65 #include "Teuchos_MatrixMarket_assignScalar.hpp"
66 #include "Teuchos_MatrixMarket_Banner.hpp"
67 #include "Teuchos_MatrixMarket_CoordDataReader.hpp"
68 #include "Teuchos_SetScientific.hpp"
69 #include "Teuchos_TimeMonitor.hpp"
70 
71 extern "C" {
72 #include "mmio_Tpetra.h"
73 }
74 #include "Tpetra_Distribution.hpp"
75 
76 
77 #include <algorithm>
78 #include <fstream>
79 #include <iostream>
80 #include <iterator>
81 #include <vector>
82 #include <stdexcept>
83 #include <numeric>
84 
85 namespace Tpetra {
115  namespace MatrixMarket {
171  template<class SparseMatrixType>
172  class Reader {
173  public:
175  typedef SparseMatrixType sparse_matrix_type;
176  typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
177 
180  typedef typename SparseMatrixType::scalar_type scalar_type;
183  typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
191  typedef typename SparseMatrixType::global_ordinal_type
194  typedef typename SparseMatrixType::node_type node_type;
195 
200 
202  typedef MultiVector<scalar_type,
206 
208  typedef Vector<scalar_type,
212 
214 
216  using trcp_tcomm_t = Teuchos::RCP<const Teuchos::Comm<int>>;
217 
218  private:
224  typedef Teuchos::ArrayRCP<int>::size_type size_type;
225 
236  static Teuchos::RCP<const map_type>
237  makeRangeMap (const trcp_tcomm_t& pComm,
238  const global_ordinal_type numRows)
239  {
240  // Return a conventional, uniformly partitioned, contiguous map.
241  return rcp (new map_type (static_cast<global_size_t> (numRows),
242  static_cast<global_ordinal_type> (0),
243  pComm, GloballyDistributed));
244  }
245 
273  static Teuchos::RCP<const map_type>
274  makeRowMap (const Teuchos::RCP<const map_type>& pRowMap,
275  const trcp_tcomm_t& pComm,
276  const global_ordinal_type numRows)
277  {
278  // If the caller didn't provide a map, return a conventional,
279  // uniformly partitioned, contiguous map.
280  if (pRowMap.is_null ()) {
281  return rcp (new map_type (static_cast<global_size_t> (numRows),
282  static_cast<global_ordinal_type> (0),
283  pComm, GloballyDistributed));
284  }
285  else {
286  TEUCHOS_TEST_FOR_EXCEPTION
287  (! pRowMap->isDistributed () && pComm->getSize () > 1,
288  std::invalid_argument, "The specified row map is not distributed, "
289  "but the given communicator includes more than one process (in "
290  "fact, there are " << pComm->getSize () << " processes).");
291  TEUCHOS_TEST_FOR_EXCEPTION
292  (pRowMap->getComm () != pComm, std::invalid_argument,
293  "The specified row Map's communicator (pRowMap->getComm()) "
294  "differs from the given separately supplied communicator pComm.");
295  return pRowMap;
296  }
297  }
298 
313  static Teuchos::RCP<const map_type>
314  makeDomainMap (const Teuchos::RCP<const map_type>& pRangeMap,
315  const global_ordinal_type numRows,
316  const global_ordinal_type numCols)
317  {
318  // Abbreviations so that the map creation call isn't too long.
319  typedef local_ordinal_type LO;
320  typedef global_ordinal_type GO;
321  typedef node_type NT;
322 
323  if (numRows == numCols) {
324  return pRangeMap;
325  } else {
326  return createUniformContigMapWithNode<LO,GO,NT> (numCols,
327  pRangeMap->getComm ());
328  }
329  }
330 
403  static void
404  distribute (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
405  Teuchos::ArrayRCP<size_t>& myRowPtr,
406  Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
407  Teuchos::ArrayRCP<scalar_type>& myValues,
408  const Teuchos::RCP<const map_type>& pRowMap,
409  Teuchos::ArrayRCP<size_t>& numEntriesPerRow,
410  Teuchos::ArrayRCP<size_t>& rowPtr,
411  Teuchos::ArrayRCP<global_ordinal_type>& colInd,
412  Teuchos::ArrayRCP<scalar_type>& values,
413  const bool debug=false)
414  {
415  using Teuchos::arcp;
416  using Teuchos::ArrayRCP;
417  using Teuchos::ArrayView;
418  using Teuchos::as;
419  using Teuchos::Comm;
420  using Teuchos::CommRequest;
421  using Teuchos::null;
422  using Teuchos::RCP;
423  using Teuchos::receive;
424  using Teuchos::send;
425  using std::cerr;
426  using std::endl;
427 
428  const bool extraDebug = false;
429  trcp_tcomm_t pComm = pRowMap->getComm ();
430  const int numProcs = pComm->getSize ();
431  const int myRank = pComm->getRank ();
432  const int rootRank = 0;
433 
434  // Type abbreviations to make the code more concise.
435  typedef global_ordinal_type GO;
436 
437  // List of the global indices of my rows. They may or may
438  // not be contiguous, and the row map need not be one-to-one.
439  ArrayView<const GO> myRows = pRowMap->getLocalElementList();
440  const size_type myNumRows = myRows.size();
441  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
442  pRowMap->getLocalNumElements(),
443  std::logic_error,
444  "pRowMap->getLocalElementList().size() = "
445  << myNumRows
446  << " != pRowMap->getLocalNumElements() = "
447  << pRowMap->getLocalNumElements() << ". "
448  "Please report this bug to the Tpetra developers.");
449  TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
450  std::logic_error,
451  "On Proc 0: numEntriesPerRow.size() = "
452  << numEntriesPerRow.size()
453  << " != pRowMap->getLocalElementList().size() = "
454  << myNumRows << ". Please report this bug to the "
455  "Tpetra developers.");
456 
457  // Space for my proc's number of entries per row. Will be
458  // filled in below.
459  myNumEntriesPerRow = arcp<size_t> (myNumRows);
460 
461  if (myRank != rootRank) {
462  // Tell the root how many rows we have. If we're sending
463  // none, then we don't have anything else to send, nor does
464  // the root have to receive anything else.
465  send (*pComm, myNumRows, rootRank);
466  if (myNumRows != 0) {
467  // Now send my rows' global indices. Hopefully the cast
468  // to int doesn't overflow. This is unlikely, since it
469  // should fit in a LO, even though it is a GO.
470  send (*pComm, static_cast<int> (myNumRows),
471  myRows.getRawPtr(), rootRank);
472 
473  // I (this proc) don't care if my global row indices are
474  // contiguous, though the root proc does (since otherwise
475  // it needs to pack noncontiguous data into contiguous
476  // storage before sending). That's why we don't check
477  // for contiguousness here.
478 
479  // Ask the root process for my part of the array of the
480  // number of entries per row.
481  receive (*pComm, rootRank,
482  static_cast<int> (myNumRows),
483  myNumEntriesPerRow.getRawPtr());
484 
485  // Use the resulting array to figure out how many column
486  // indices and values I should ask from the root process.
487  const local_ordinal_type myNumEntries =
488  std::accumulate (myNumEntriesPerRow.begin(),
489  myNumEntriesPerRow.end(), 0);
490 
491  // Make space for my entries of the sparse matrix. Note
492  // that they don't have to be sorted by row index.
493  // Iterating through all my rows requires computing a
494  // running sum over myNumEntriesPerRow.
495  myColInd = arcp<GO> (myNumEntries);
496  myValues = arcp<scalar_type> (myNumEntries);
497  if (myNumEntries > 0) {
498  // Ask for that many column indices and values, if
499  // there are any.
500  receive (*pComm, rootRank,
501  static_cast<int> (myNumEntries),
502  myColInd.getRawPtr());
503  receive (*pComm, rootRank,
504  static_cast<int> (myNumEntries),
505  myValues.getRawPtr());
506  }
507  } // If I own at least one row
508  } // If I am not the root processor
509  else { // I _am_ the root processor
510  if (debug) {
511  cerr << "-- Proc 0: Copying my data from global arrays" << endl;
512  }
513  // Proc 0 still needs to (allocate, if not done already)
514  // and fill its part of the matrix (my*).
515  for (size_type k = 0; k < myNumRows; ++k) {
516  const GO myCurRow = myRows[k];
517  const local_ordinal_type numEntriesInThisRow = numEntriesPerRow[myCurRow];
518  myNumEntriesPerRow[k] = numEntriesInThisRow;
519  }
520  if (extraDebug && debug) {
521  cerr << "Proc " << pRowMap->getComm ()->getRank ()
522  << ": myNumEntriesPerRow[0.." << (myNumRows-1) << "] = [";
523  for (size_type k = 0; k < myNumRows; ++k) {
524  cerr << myNumEntriesPerRow[k];
525  if (k < myNumRows-1) {
526  cerr << " ";
527  }
528  }
529  cerr << "]" << endl;
530  }
531  // The total number of matrix entries that my proc owns.
532  const local_ordinal_type myNumEntries =
533  std::accumulate (myNumEntriesPerRow.begin(),
534  myNumEntriesPerRow.end(), 0);
535  if (debug) {
536  cerr << "-- Proc 0: I own " << myNumRows << " rows and "
537  << myNumEntries << " entries" << endl;
538  }
539  myColInd = arcp<GO> (myNumEntries);
540  myValues = arcp<scalar_type> (myNumEntries);
541 
542  // Copy Proc 0's part of the matrix into the my* arrays.
543  // It's important that myCurPos be updated _before_ k,
544  // otherwise myCurPos will get the wrong number of entries
545  // per row (it should be for the row in the just-completed
546  // iteration, not for the next iteration's row).
547  local_ordinal_type myCurPos = 0;
548  for (size_type k = 0; k < myNumRows;
549  myCurPos += myNumEntriesPerRow[k], ++k) {
550  const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
551  const GO myRow = myRows[k];
552  const size_t curPos = rowPtr[myRow];
553  // Only copy if there are entries to copy, in order not
554  // to construct empty ranges for the ArrayRCP views.
555  if (curNumEntries > 0) {
556  ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
557  ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
558  std::copy (colIndView.begin(), colIndView.end(),
559  myColIndView.begin());
560 
561  ArrayView<scalar_type> valuesView =
562  values (curPos, curNumEntries);
563  ArrayView<scalar_type> myValuesView =
564  myValues (myCurPos, curNumEntries);
565  std::copy (valuesView.begin(), valuesView.end(),
566  myValuesView.begin());
567  }
568  }
569 
570  // Proc 0 processes each other proc p in turn.
571  for (int p = 1; p < numProcs; ++p) {
572  if (debug) {
573  cerr << "-- Proc 0: Processing proc " << p << endl;
574  }
575 
576  size_type theirNumRows = 0;
577  // Ask Proc p how many rows it has. If it doesn't
578  // have any, we can move on to the next proc. This
579  // has to be a standard receive so that we can avoid
580  // the degenerate case of sending zero data.
581  receive (*pComm, p, &theirNumRows);
582  if (debug) {
583  cerr << "-- Proc 0: Proc " << p << " owns "
584  << theirNumRows << " rows" << endl;
585  }
586  if (theirNumRows != 0) {
587  // Ask Proc p which rows it owns. The resulting global
588  // row indices are not guaranteed to be contiguous or
589  // sorted. Global row indices are themselves indices
590  // into the numEntriesPerRow array.
591  ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
592  receive (*pComm, p, as<int> (theirNumRows),
593  theirRows.getRawPtr ());
594  // Extra test to make sure that the rows we received
595  // are all sensible. This is a good idea since we are
596  // going to use the global row indices we've received
597  // to index into the numEntriesPerRow array. Better to
598  // catch any bugs here and print a sensible error
599  // message, rather than segfault and print a cryptic
600  // error message.
601  {
602  const global_size_t numRows = pRowMap->getGlobalNumElements ();
603  const GO indexBase = pRowMap->getIndexBase ();
604  bool theirRowsValid = true;
605  for (size_type k = 0; k < theirNumRows; ++k) {
606  if (theirRows[k] < indexBase ||
607  as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
608  theirRowsValid = false;
609  }
610  }
611  if (! theirRowsValid) {
612  TEUCHOS_TEST_FOR_EXCEPTION(
613  ! theirRowsValid, std::logic_error,
614  "Proc " << p << " has at least one invalid row index. "
615  "Here are all of them: " <<
616  Teuchos::toString (theirRows ()) << ". Valid row index "
617  "range (zero-based): [0, " << (numRows - 1) << "].");
618  }
619  }
620 
621  // Perhaps we could save a little work if we check
622  // whether Proc p's row indices are contiguous. That
623  // would make lookups in the global data arrays
624  // faster. For now, we just implement the general
625  // case and don't prematurely optimize. (Remember
626  // that you're making Proc 0 read the whole file, so
627  // you've already lost scalability.)
628 
629  // Compute the number of entries in each of Proc p's
630  // rows. (Proc p will compute its row pointer array
631  // on its own, after it gets the data from Proc 0.)
632  ArrayRCP<size_t> theirNumEntriesPerRow;
633  theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
634  for (size_type k = 0; k < theirNumRows; ++k) {
635  theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
636  }
637 
638  // Tell Proc p the number of entries in each of its
639  // rows. Hopefully the cast to int doesn't overflow.
640  // This is unlikely, since it should fit in a LO,
641  // even though it is a GO.
642  send (*pComm, static_cast<int> (theirNumRows),
643  theirNumEntriesPerRow.getRawPtr(), p);
644 
645  // Figure out how many entries Proc p owns.
646  const local_ordinal_type theirNumEntries =
647  std::accumulate (theirNumEntriesPerRow.begin(),
648  theirNumEntriesPerRow.end(), 0);
649 
650  if (debug) {
651  cerr << "-- Proc 0: Proc " << p << " owns "
652  << theirNumEntries << " entries" << endl;
653  }
654 
655  // If there are no entries to send, then we're done
656  // with Proc p.
657  if (theirNumEntries == 0) {
658  continue;
659  }
660 
661  // Construct (views of) proc p's column indices and
662  // values. Later, we might like to optimize for the
663  // (common) contiguous case, for which we don't need to
664  // copy data into separate "their*" arrays (we can just
665  // use contiguous views of the global arrays).
666  ArrayRCP<GO> theirColInd (theirNumEntries);
667  ArrayRCP<scalar_type> theirValues (theirNumEntries);
668  // Copy Proc p's part of the matrix into the their*
669  // arrays. It's important that theirCurPos be updated
670  // _before_ k, otherwise theirCurPos will get the wrong
671  // number of entries per row (it should be for the row
672  // in the just-completed iteration, not for the next
673  // iteration's row).
674  local_ordinal_type theirCurPos = 0;
675  for (size_type k = 0; k < theirNumRows;
676  theirCurPos += theirNumEntriesPerRow[k], k++) {
677  const local_ordinal_type curNumEntries = theirNumEntriesPerRow[k];
678  const GO theirRow = theirRows[k];
679  const local_ordinal_type curPos = rowPtr[theirRow];
680 
681  // Only copy if there are entries to copy, in order
682  // not to construct empty ranges for the ArrayRCP
683  // views.
684  if (curNumEntries > 0) {
685  ArrayView<GO> colIndView =
686  colInd (curPos, curNumEntries);
687  ArrayView<GO> theirColIndView =
688  theirColInd (theirCurPos, curNumEntries);
689  std::copy (colIndView.begin(), colIndView.end(),
690  theirColIndView.begin());
691 
692  ArrayView<scalar_type> valuesView =
693  values (curPos, curNumEntries);
694  ArrayView<scalar_type> theirValuesView =
695  theirValues (theirCurPos, curNumEntries);
696  std::copy (valuesView.begin(), valuesView.end(),
697  theirValuesView.begin());
698  }
699  }
700  // Send Proc p its column indices and values.
701  // Hopefully the cast to int doesn't overflow. This
702  // is unlikely, since it should fit in a LO, even
703  // though it is a GO.
704  send (*pComm, static_cast<int> (theirNumEntries),
705  theirColInd.getRawPtr(), p);
706  send (*pComm, static_cast<int> (theirNumEntries),
707  theirValues.getRawPtr(), p);
708 
709  if (debug) {
710  cerr << "-- Proc 0: Finished with proc " << p << endl;
711  }
712  } // If proc p owns at least one row
713  } // For each proc p not the root proc 0
714  } // If I'm (not) the root proc 0
715 
716  // Invalidate the input data to save space, since we don't
717  // need it anymore.
718  numEntriesPerRow = null;
719  rowPtr = null;
720  colInd = null;
721  values = null;
722 
723  if (debug && myRank == 0) {
724  cerr << "-- Proc 0: About to fill in myRowPtr" << endl;
725  }
726 
727  // Allocate and fill in myRowPtr (the row pointer array for
728  // my rank's rows). We delay this until the end because we
729  // don't need it to compute anything else in distribute().
730  // Each proc can do this work for itself, since it only needs
731  // myNumEntriesPerRow to do so.
732  myRowPtr = arcp<size_t> (myNumRows+1);
733  myRowPtr[0] = 0;
734  for (size_type k = 1; k < myNumRows+1; ++k) {
735  myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
736  }
737  if (extraDebug && debug) {
738  cerr << "Proc " << Teuchos::rank (*(pRowMap->getComm()))
739  << ": myRowPtr[0.." << myNumRows << "] = [";
740  for (size_type k = 0; k < myNumRows+1; ++k) {
741  cerr << myRowPtr[k];
742  if (k < myNumRows) {
743  cerr << " ";
744  }
745  }
746  cerr << "]" << endl << endl;
747  }
748 
749  if (debug && myRank == 0) {
750  cerr << "-- Proc 0: Done with distribute" << endl;
751  }
752  }
753 
767  static Teuchos::RCP<sparse_matrix_type>
768  makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
769  Teuchos::ArrayRCP<size_t>& myRowPtr,
770  Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
771  Teuchos::ArrayRCP<scalar_type>& myValues,
772  const Teuchos::RCP<const map_type>& pRowMap,
773  const Teuchos::RCP<const map_type>& pRangeMap,
774  const Teuchos::RCP<const map_type>& pDomainMap,
775  const bool callFillComplete = true)
776  {
777  using Teuchos::ArrayView;
778  using Teuchos::null;
779  using Teuchos::RCP;
780  using Teuchos::rcp;
781  using std::cerr;
782  using std::endl;
783  // Typedef to make certain type declarations shorter.
784  typedef global_ordinal_type GO;
785 
786  // The row pointer array always has at least one entry, even
787  // if the matrix has zero rows. myNumEntriesPerRow, myColInd,
788  // and myValues would all be empty arrays in that degenerate
789  // case, but the row and domain maps would still be nonnull
790  // (though they would be trivial maps).
791  TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
792  "makeMatrix: myRowPtr array is null. "
793  "Please report this bug to the Tpetra developers.");
794  TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
795  "makeMatrix: domain map is null. "
796  "Please report this bug to the Tpetra developers.");
797  TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
798  "makeMatrix: range map is null. "
799  "Please report this bug to the Tpetra developers.");
800  TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
801  "makeMatrix: row map is null. "
802  "Please report this bug to the Tpetra developers.");
803 
804  // Construct the CrsMatrix, using the row map, with the
805  // constructor specifying the number of nonzeros for each row.
806  RCP<sparse_matrix_type> A =
807  rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow ()));
808 
809  // List of the global indices of my rows.
810  // They may or may not be contiguous.
811  ArrayView<const GO> myRows = pRowMap->getLocalElementList ();
812  const size_type myNumRows = myRows.size ();
813 
814  // Add this processor's matrix entries to the CrsMatrix.
815  const GO indexBase = pRowMap->getIndexBase ();
816  for (size_type i = 0; i < myNumRows; ++i) {
817  const size_type myCurPos = myRowPtr[i];
818  const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
819  ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
820  ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
821 
822  // Modify the column indices in place to have the right index base.
823  for (size_type k = 0; k < curNumEntries; ++k) {
824  curColInd[k] += indexBase;
825  }
826  // Avoid constructing empty views of ArrayRCP objects.
827  if (curNumEntries > 0) {
828  A->insertGlobalValues (myRows[i], curColInd, curValues);
829  }
830  }
831  // We've entered in all our matrix entries, so we can delete
832  // the original data. This will save memory when we call
833  // fillComplete(), so that we never keep more than two copies
834  // of the matrix's data in memory at once.
835  myNumEntriesPerRow = null;
836  myRowPtr = null;
837  myColInd = null;
838  myValues = null;
839 
840  if (callFillComplete) {
841  A->fillComplete (pDomainMap, pRangeMap);
842  }
843  return A;
844  }
845 
851  static Teuchos::RCP<sparse_matrix_type>
852  makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
853  Teuchos::ArrayRCP<size_t>& myRowPtr,
854  Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
855  Teuchos::ArrayRCP<scalar_type>& myValues,
856  const Teuchos::RCP<const map_type>& pRowMap,
857  const Teuchos::RCP<const map_type>& pRangeMap,
858  const Teuchos::RCP<const map_type>& pDomainMap,
859  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
860  const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams)
861  {
862  using Teuchos::ArrayView;
863  using Teuchos::null;
864  using Teuchos::RCP;
865  using Teuchos::rcp;
866  using std::cerr;
867  using std::endl;
868  // Typedef to make certain type declarations shorter.
869  typedef global_ordinal_type GO;
870 
871  // The row pointer array always has at least one entry, even
872  // if the matrix has zero rows. myNumEntriesPerRow, myColInd,
873  // and myValues would all be empty arrays in that degenerate
874  // case, but the row and domain maps would still be nonnull
875  // (though they would be trivial maps).
876  TEUCHOS_TEST_FOR_EXCEPTION(
877  myRowPtr.is_null(), std::logic_error,
878  "makeMatrix: myRowPtr array is null. "
879  "Please report this bug to the Tpetra developers.");
880  TEUCHOS_TEST_FOR_EXCEPTION(
881  pDomainMap.is_null(), std::logic_error,
882  "makeMatrix: domain map is null. "
883  "Please report this bug to the Tpetra developers.");
884  TEUCHOS_TEST_FOR_EXCEPTION(
885  pRangeMap.is_null(), std::logic_error,
886  "makeMatrix: range map is null. "
887  "Please report this bug to the Tpetra developers.");
888  TEUCHOS_TEST_FOR_EXCEPTION(
889  pRowMap.is_null(), std::logic_error,
890  "makeMatrix: row map is null. "
891  "Please report this bug to the Tpetra developers.");
892 
893  // Construct the CrsMatrix, using the row map, with the
894  // constructor specifying the number of nonzeros for each row.
895  RCP<sparse_matrix_type> A =
896  rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow(),
897  constructorParams));
898 
899  // List of the global indices of my rows.
900  // They may or may not be contiguous.
901  ArrayView<const GO> myRows = pRowMap->getLocalElementList();
902  const size_type myNumRows = myRows.size();
903 
904  // Add this processor's matrix entries to the CrsMatrix.
905  const GO indexBase = pRowMap->getIndexBase ();
906  for (size_type i = 0; i < myNumRows; ++i) {
907  const size_type myCurPos = myRowPtr[i];
908  const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
909  ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
910  ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
911 
912  // Modify the column indices in place to have the right index base.
913  for (size_type k = 0; k < curNumEntries; ++k) {
914  curColInd[k] += indexBase;
915  }
916  if (curNumEntries > 0) {
917  A->insertGlobalValues (myRows[i], curColInd, curValues);
918  }
919  }
920  // We've entered in all our matrix entries, so we can delete
921  // the original data. This will save memory when we call
922  // fillComplete(), so that we never keep more than two copies
923  // of the matrix's data in memory at once.
924  myNumEntriesPerRow = null;
925  myRowPtr = null;
926  myColInd = null;
927  myValues = null;
928 
929  A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
930  return A;
931  }
932 
937  static Teuchos::RCP<sparse_matrix_type>
938  makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
939  Teuchos::ArrayRCP<size_t>& myRowPtr,
940  Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
941  Teuchos::ArrayRCP<scalar_type>& myValues,
942  const Teuchos::RCP<const map_type>& rowMap,
943  Teuchos::RCP<const map_type>& colMap,
944  const Teuchos::RCP<const map_type>& domainMap,
945  const Teuchos::RCP<const map_type>& rangeMap,
946  const bool callFillComplete = true)
947  {
948  using Teuchos::ArrayView;
949  using Teuchos::as;
950  using Teuchos::null;
951  using Teuchos::RCP;
952  using Teuchos::rcp;
953  typedef global_ordinal_type GO;
954 
955  // Construct the CrsMatrix.
956 
957  RCP<sparse_matrix_type> A; // the matrix to return.
958  if (colMap.is_null ()) { // the user didn't provide a column Map
959  A = rcp (new sparse_matrix_type (rowMap, myNumEntriesPerRow));
960  } else { // the user provided a column Map
961  A = rcp (new sparse_matrix_type (rowMap, colMap, myNumEntriesPerRow));
962  }
963 
964  // List of the global indices of my rows.
965  // They may or may not be contiguous.
966  ArrayView<const GO> myRows = rowMap->getLocalElementList ();
967  const size_type myNumRows = myRows.size ();
968 
969  // Add this process' matrix entries to the CrsMatrix.
970  const GO indexBase = rowMap->getIndexBase ();
971  for (size_type i = 0; i < myNumRows; ++i) {
972  const size_type myCurPos = myRowPtr[i];
973  const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
974  ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
975  ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
976 
977  // Modify the column indices in place to have the right index base.
978  for (size_type k = 0; k < curNumEntries; ++k) {
979  curColInd[k] += indexBase;
980  }
981  if (curNumEntries > 0) {
982  A->insertGlobalValues (myRows[i], curColInd, curValues);
983  }
984  }
985  // We've entered in all our matrix entries, so we can delete
986  // the original data. This will save memory when we call
987  // fillComplete(), so that we never keep more than two copies
988  // of the matrix's data in memory at once.
989  myNumEntriesPerRow = null;
990  myRowPtr = null;
991  myColInd = null;
992  myValues = null;
993 
994  if (callFillComplete) {
995  A->fillComplete (domainMap, rangeMap);
996  if (colMap.is_null ()) {
997  colMap = A->getColMap ();
998  }
999  }
1000  return A;
1001  }
1002 
1003  private:
1004 
1021  static Teuchos::RCP<const Teuchos::MatrixMarket::Banner>
1022  readBanner (std::istream& in,
1023  size_t& lineNumber,
1024  const bool tolerant=false,
1025  const bool /* debug */=false,
1026  const bool isGraph=false)
1027  {
1028  using Teuchos::MatrixMarket::Banner;
1029  using Teuchos::RCP;
1030  using Teuchos::rcp;
1031  using std::cerr;
1032  using std::endl;
1033  typedef Teuchos::ScalarTraits<scalar_type> STS;
1034 
1035  RCP<Banner> pBanner; // On output, if successful: the read-in Banner.
1036  std::string line; // If read from stream successful: the Banner line
1037 
1038  // Try to read a line from the input stream.
1039  const bool readFailed = ! getline(in, line);
1040  TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1041  "Failed to get Matrix Market banner line from input.");
1042 
1043  // We read a line from the input stream.
1044  lineNumber++;
1045 
1046  // Assume that the line we found is the Banner line.
1047  try {
1048  pBanner = rcp (new Banner (line, tolerant));
1049  } catch (std::exception& e) {
1050 
1051  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
1052  "Matrix Market banner line contains syntax error(s): "
1053  << e.what());
1054  }
1055 
1056  TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() != "matrix",
1057  std::invalid_argument, "The Matrix Market file does not contain "
1058  "matrix data. Its Banner (first) line says that its object type is \""
1059  << pBanner->matrixType() << "\", rather than the required \"matrix\".");
1060 
1061  // Validate the data type of the matrix, with respect to the
1062  // Scalar type of the CrsMatrix entries.
1063  TEUCHOS_TEST_FOR_EXCEPTION(
1064  ! STS::isComplex && pBanner->dataType() == "complex",
1065  std::invalid_argument,
1066  "The Matrix Market file contains complex-valued data, but you are "
1067  "trying to read it into a matrix containing entries of the real-"
1068  "valued Scalar type \""
1069  << Teuchos::TypeNameTraits<scalar_type>::name() << "\".");
1070  TEUCHOS_TEST_FOR_EXCEPTION(
1071  !isGraph &&
1072  pBanner->dataType() != "real" &&
1073  pBanner->dataType() != "complex" &&
1074  pBanner->dataType() != "integer",
1075  std::invalid_argument,
1076  "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
1077  "Matrix Market file may not contain a \"pattern\" matrix. A "
1078  "pattern matrix is really just a graph with no weights. It "
1079  "should be stored in a CrsGraph, not a CrsMatrix.");
1080 
1081  TEUCHOS_TEST_FOR_EXCEPTION(
1082  isGraph &&
1083  pBanner->dataType() != "pattern",
1084  std::invalid_argument,
1085  "When reading Matrix Market data into a Tpetra::CrsGraph, the "
1086  "Matrix Market file must contain a \"pattern\" matrix.");
1087 
1088  return pBanner;
1089  }
1090 
1113  static Teuchos::Tuple<global_ordinal_type, 3>
1114  readCoordDims (std::istream& in,
1115  size_t& lineNumber,
1116  const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1117  const trcp_tcomm_t& pComm,
1118  const bool tolerant = false,
1119  const bool /* debug */ = false)
1120  {
1121  using Teuchos::MatrixMarket::readCoordinateDimensions;
1122  using Teuchos::Tuple;
1123 
1124  // Packed coordinate matrix dimensions (numRows, numCols,
1125  // numNonzeros); computed on Rank 0 and broadcasted to all MPI
1126  // ranks.
1127  Tuple<global_ordinal_type, 3> dims;
1128 
1129  // Read in the coordinate matrix dimensions from the input
1130  // stream. "success" tells us whether reading in the
1131  // coordinate matrix dimensions succeeded ("Guilty unless
1132  // proven innocent").
1133  bool success = false;
1134  if (pComm->getRank() == 0) {
1135  TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "coordinate",
1136  std::invalid_argument, "The Tpetra::CrsMatrix Matrix Market reader "
1137  "only accepts \"coordinate\" (sparse) matrix data.");
1138  // Unpacked coordinate matrix dimensions
1139  global_ordinal_type numRows, numCols, numNonzeros;
1140  // Only MPI Rank 0 reads from the input stream
1141  success = readCoordinateDimensions (in, numRows, numCols,
1142  numNonzeros, lineNumber,
1143  tolerant);
1144  // Pack up the data into a Tuple so we can send them with
1145  // one broadcast instead of three.
1146  dims[0] = numRows;
1147  dims[1] = numCols;
1148  dims[2] = numNonzeros;
1149  }
1150  // Only Rank 0 did the reading, so it decides success.
1151  //
1152  // FIXME (mfh 02 Feb 2011) Teuchos::broadcast doesn't know how
1153  // to send bools. For now, we convert to/from int instead,
1154  // using the usual "true is 1, false is 0" encoding.
1155  {
1156  int the_success = success ? 1 : 0; // only matters on MPI Rank 0
1157  Teuchos::broadcast (*pComm, 0, &the_success);
1158  success = (the_success == 1);
1159  }
1160  if (success) {
1161  // Broadcast (numRows, numCols, numNonzeros) from Rank 0
1162  // to all the other MPI ranks.
1163  Teuchos::broadcast (*pComm, 0, dims);
1164  }
1165  else {
1166  // Perhaps in tolerant mode, we could set all the
1167  // dimensions to zero for now, and deduce correct
1168  // dimensions by reading all of the file's entries and
1169  // computing the max(row index) and max(column index).
1170  // However, for now we just error out in that case.
1171  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
1172  "Error reading Matrix Market sparse matrix: failed to read "
1173  "coordinate matrix dimensions.");
1174  }
1175  return dims;
1176  }
1177 
1188  typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1189 
1190  typedef Teuchos::MatrixMarket::SymmetrizingGraphAdder<Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> > graph_adder_type;
1191 
1217  static Teuchos::RCP<adder_type>
1218  makeAdder (const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1219  Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1220  const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1221  const bool tolerant=false,
1222  const bool debug=false)
1223  {
1224  if (pComm->getRank () == 0) {
1225  typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type,
1227  raw_adder_type;
1228  Teuchos::RCP<raw_adder_type> pRaw =
1229  Teuchos::rcp (new raw_adder_type (dims[0], dims[1], dims[2],
1230  tolerant, debug));
1231  return Teuchos::rcp (new adder_type (pRaw, pBanner->symmType ()));
1232  }
1233  else {
1234  return Teuchos::null;
1235  }
1236  }
1237 
1263  static Teuchos::RCP<graph_adder_type>
1264  makeGraphAdder (const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1265  Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1266  const Teuchos::Tuple<global_ordinal_type, 3>& dims,
1267  const bool tolerant=false,
1268  const bool debug=false)
1269  {
1270  if (pComm->getRank () == 0) {
1271  typedef Teuchos::MatrixMarket::Raw::GraphAdder<global_ordinal_type> raw_adder_type;
1272  Teuchos::RCP<raw_adder_type> pRaw =
1273  Teuchos::rcp (new raw_adder_type (dims[0], dims[1], dims[2],
1274  tolerant, debug));
1275  return Teuchos::rcp (new graph_adder_type (pRaw, pBanner->symmType ()));
1276  }
1277  else {
1278  return Teuchos::null;
1279  }
1280  }
1281 
1283  static Teuchos::RCP<sparse_graph_type>
1284  readSparseGraphHelper (std::istream& in,
1285  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1286  const Teuchos::RCP<const map_type>& rowMap,
1287  Teuchos::RCP<const map_type>& colMap,
1288  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1289  const bool tolerant,
1290  const bool debug)
1291  {
1292  using Teuchos::MatrixMarket::Banner;
1293  using Teuchos::RCP;
1294  using Teuchos::ptr;
1295  using Teuchos::Tuple;
1296  using std::cerr;
1297  using std::endl;
1298 
1299  const int myRank = pComm->getRank ();
1300  const int rootRank = 0;
1301 
1302  // Current line number in the input stream. Various calls
1303  // will modify this depending on the number of lines that are
1304  // read from the input stream. Only Rank 0 modifies this.
1305  size_t lineNumber = 1;
1306 
1307  if (debug && myRank == rootRank) {
1308  cerr << "Matrix Market reader: readGraph:" << endl
1309  << "-- Reading banner line" << endl;
1310  }
1311 
1312  // The "Banner" tells you whether the input stream represents
1313  // a sparse matrix, the symmetry type of the matrix, and the
1314  // type of the data it contains.
1315  //
1316  // pBanner will only be nonnull on MPI Rank 0. It will be
1317  // null on all other MPI processes.
1318  RCP<const Banner> pBanner;
1319  {
1320  // We read and validate the Banner on Proc 0, but broadcast
1321  // the validation result to all processes.
1322  // Teuchos::broadcast doesn't currently work with bool, so
1323  // we use int (true -> 1, false -> 0).
1324  int bannerIsCorrect = 1;
1325  std::ostringstream errMsg;
1326 
1327  if (myRank == rootRank) {
1328  // Read the Banner line from the input stream.
1329  try {
1330  pBanner = readBanner (in, lineNumber, tolerant, debug, true);
1331  }
1332  catch (std::exception& e) {
1333  errMsg << "Attempt to read the Matrix Market file's Banner line "
1334  "threw an exception: " << e.what();
1335  bannerIsCorrect = 0;
1336  }
1337 
1338  if (bannerIsCorrect) {
1339  // Validate the Banner for the case of a sparse graph.
1340  // We validate on Proc 0, since it reads the Banner.
1341 
1342  // In intolerant mode, the matrix type must be "coordinate".
1343  if (! tolerant && pBanner->matrixType() != "coordinate") {
1344  bannerIsCorrect = 0;
1345  errMsg << "The Matrix Market input file must contain a "
1346  "\"coordinate\"-format sparse graph in order to create a "
1347  "Tpetra::CrsGraph object from it, but the file's matrix "
1348  "type is \"" << pBanner->matrixType() << "\" instead.";
1349  }
1350  // In tolerant mode, we allow the matrix type to be
1351  // anything other than "array" (which would mean that
1352  // the file contains a dense matrix).
1353  if (tolerant && pBanner->matrixType() == "array") {
1354  bannerIsCorrect = 0;
1355  errMsg << "Matrix Market file must contain a \"coordinate\"-"
1356  "format sparse graph in order to create a Tpetra::CrsGraph "
1357  "object from it, but the file's matrix type is \"array\" "
1358  "instead. That probably means the file contains dense matrix "
1359  "data.";
1360  }
1361  }
1362  } // Proc 0: Done reading the Banner, hopefully successfully.
1363 
1364  // Broadcast from Proc 0 whether the Banner was read correctly.
1365  broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1366 
1367  // If the Banner is invalid, all processes throw an
1368  // exception. Only Proc 0 gets the exception message, but
1369  // that's OK, since the main point is to "stop the world"
1370  // (rather than throw an exception on one process and leave
1371  // the others hanging).
1372  TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1373  std::invalid_argument, errMsg.str ());
1374  } // Done reading the Banner line and broadcasting success.
1375  if (debug && myRank == rootRank) {
1376  cerr << "-- Reading dimensions line" << endl;
1377  }
1378 
1379  // Read the graph dimensions from the Matrix Market metadata.
1380  // dims = (numRows, numCols, numEntries). Proc 0 does the
1381  // reading, but it broadcasts the results to all MPI
1382  // processes. Thus, readCoordDims() is a collective
1383  // operation. It does a collective check for correctness too.
1384  Tuple<global_ordinal_type, 3> dims =
1385  readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1386 
1387  if (debug && myRank == rootRank) {
1388  cerr << "-- Making Adder for collecting graph data" << endl;
1389  }
1390 
1391  // "Adder" object for collecting all the sparse graph entries
1392  // from the input stream. This is only nonnull on Proc 0.
1393  // The Adder internally converts the one-based indices (native
1394  // Matrix Market format) into zero-based indices.
1395  RCP<graph_adder_type> pAdder =
1396  makeGraphAdder (pComm, pBanner, dims, tolerant, debug);
1397 
1398  if (debug && myRank == rootRank) {
1399  cerr << "-- Reading graph data" << endl;
1400  }
1401  //
1402  // Read the graph entries from the input stream on Proc 0.
1403  //
1404  {
1405  // We use readSuccess to broadcast the results of the read
1406  // (succeeded or not) to all MPI processes. Since
1407  // Teuchos::broadcast doesn't currently know how to send
1408  // bools, we convert to int (true -> 1, false -> 0).
1409  int readSuccess = 1;
1410  std::ostringstream errMsg; // Exception message (only valid on Proc 0)
1411  if (myRank == rootRank) {
1412  try {
1413  // Reader for "coordinate" format sparse graph data.
1414  typedef Teuchos::MatrixMarket::CoordPatternReader<graph_adder_type,
1415  global_ordinal_type> reader_type;
1416  reader_type reader (pAdder);
1417 
1418  // Read the sparse graph entries.
1419  std::pair<bool, std::vector<size_t> > results =
1420  reader.read (in, lineNumber, tolerant, debug);
1421  readSuccess = results.first ? 1 : 0;
1422  }
1423  catch (std::exception& e) {
1424  readSuccess = 0;
1425  errMsg << e.what();
1426  }
1427  }
1428  broadcast (*pComm, rootRank, ptr (&readSuccess));
1429 
1430  // It would be nice to add a "verbose" flag, so that in
1431  // tolerant mode, we could log any bad line number(s) on
1432  // Proc 0. For now, we just throw if the read fails to
1433  // succeed.
1434  //
1435  // Question: If we're in tolerant mode, and if the read did
1436  // not succeed, should we attempt to call fillComplete()?
1437  TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1438  "Failed to read the Matrix Market sparse graph file: "
1439  << errMsg.str());
1440  } // Done reading the graph entries (stored on Proc 0 for now)
1441 
1442  if (debug && myRank == rootRank) {
1443  cerr << "-- Successfully read the Matrix Market data" << endl;
1444  }
1445 
1446  // In tolerant mode, we need to rebroadcast the graph
1447  // dimensions, since they may be different after reading the
1448  // actual graph data. We only need to broadcast the number
1449  // of rows and columns. Only Rank 0 needs to know the actual
1450  // global number of entries, since (a) we need to merge
1451  // duplicates on Rank 0 first anyway, and (b) when we
1452  // distribute the entries, each rank other than Rank 0 will
1453  // only need to know how many entries it owns, not the total
1454  // number of entries.
1455  if (tolerant) {
1456  if (debug && myRank == rootRank) {
1457  cerr << "-- Tolerant mode: rebroadcasting graph dimensions"
1458  << endl
1459  << "----- Dimensions before: "
1460  << dims[0] << " x " << dims[1]
1461  << endl;
1462  }
1463  // Packed coordinate graph dimensions (numRows, numCols).
1464  Tuple<global_ordinal_type, 2> updatedDims;
1465  if (myRank == rootRank) {
1466  // If one or more bottom rows of the graph contain no
1467  // entries, then the Adder will report that the number
1468  // of rows is less than that specified in the
1469  // metadata. We allow this case, and favor the
1470  // metadata so that the zero row(s) will be included.
1471  updatedDims[0] =
1472  std::max (dims[0], pAdder->getAdder()->numRows());
1473  updatedDims[1] = pAdder->getAdder()->numCols();
1474  }
1475  broadcast (*pComm, rootRank, updatedDims);
1476  dims[0] = updatedDims[0];
1477  dims[1] = updatedDims[1];
1478  if (debug && myRank == rootRank) {
1479  cerr << "----- Dimensions after: " << dims[0] << " x "
1480  << dims[1] << endl;
1481  }
1482  }
1483  else {
1484  // In strict mode, we require that the graph's metadata and
1485  // its actual data agree, at least somewhat. In particular,
1486  // the number of rows must agree, since otherwise we cannot
1487  // distribute the graph correctly.
1488 
1489  // Teuchos::broadcast() doesn't know how to broadcast bools,
1490  // so we use an int with the standard 1 == true, 0 == false
1491  // encoding.
1492  int dimsMatch = 1;
1493  if (myRank == rootRank) {
1494  // If one or more bottom rows of the graph contain no
1495  // entries, then the Adder will report that the number of
1496  // rows is less than that specified in the metadata. We
1497  // allow this case, and favor the metadata, but do not
1498  // allow the Adder to think there are more rows in the
1499  // graph than the metadata says.
1500  if (dims[0] < pAdder->getAdder ()->numRows ()) {
1501  dimsMatch = 0;
1502  }
1503  }
1504  broadcast (*pComm, 0, ptr (&dimsMatch));
1505  if (dimsMatch == 0) {
1506  // We're in an error state anyway, so we might as well
1507  // work a little harder to print an informative error
1508  // message.
1509  //
1510  // Broadcast the Adder's idea of the graph dimensions
1511  // from Proc 0 to all processes.
1512  Tuple<global_ordinal_type, 2> addersDims;
1513  if (myRank == rootRank) {
1514  addersDims[0] = pAdder->getAdder()->numRows();
1515  addersDims[1] = pAdder->getAdder()->numCols();
1516  }
1517  broadcast (*pComm, 0, addersDims);
1518  TEUCHOS_TEST_FOR_EXCEPTION(
1519  dimsMatch == 0, std::runtime_error,
1520  "The graph metadata says that the graph is " << dims[0] << " x "
1521  << dims[1] << ", but the actual data says that the graph is "
1522  << addersDims[0] << " x " << addersDims[1] << ". That means the "
1523  "data includes more rows than reported in the metadata. This "
1524  "is not allowed when parsing in strict mode. Parse the graph in "
1525  "tolerant mode to ignore the metadata when it disagrees with the "
1526  "data.");
1527  }
1528  } // Matrix dimensions (# rows, # cols, # entries) agree.
1529 
1530  // Create a map describing a distribution where the root owns EVERYTHING
1531  RCP<map_type> proc0Map;
1532  global_ordinal_type indexBase;
1533  if(Teuchos::is_null(rowMap)) {
1534  indexBase = 0;
1535  }
1536  else {
1537  indexBase = rowMap->getIndexBase();
1538  }
1539  if(myRank == rootRank) {
1540  proc0Map = rcp(new map_type(dims[0],dims[0],indexBase,pComm));
1541  }
1542  else {
1543  proc0Map = rcp(new map_type(dims[0],0,indexBase,pComm));
1544  }
1545 
1546  // Create the graph where the root owns EVERYTHING
1547  std::map<global_ordinal_type, size_t> numEntriesPerRow_map;
1548  if (myRank == rootRank) {
1549  const auto& entries = pAdder()->getAdder()->getEntries();
1550  // This will count duplicates, but it's better than dense.
1551  // An even better approach would use a classic algorithm,
1552  // likely in Saad's old textbook, for converting COO (entries)
1553  // to CSR (the local part of the sparse matrix data structure).
1554  for (const auto& entry : entries) {
1555  const global_ordinal_type gblRow = entry.rowIndex () + indexBase;
1556  ++numEntriesPerRow_map[gblRow];
1557  }
1558  }
1559 
1560  Teuchos::Array<size_t> numEntriesPerRow (proc0Map->getLocalNumElements ());
1561  for (const auto& ent : numEntriesPerRow_map) {
1562  const local_ordinal_type lclRow = proc0Map->getLocalElement (ent.first);
1563  numEntriesPerRow[lclRow] = ent.second;
1564  }
1565  // Free anything we don't need before allocating the graph.
1566  // Swapping with an empty data structure is the standard idiom
1567  // for freeing memory used by Standard Library containers.
1568  // (Just resizing to 0 doesn't promise to free memory.)
1569  {
1570  std::map<global_ordinal_type, size_t> empty_map;
1571  std::swap (numEntriesPerRow_map, empty_map);
1572  }
1573 
1574  RCP<sparse_graph_type> proc0Graph =
1575  rcp(new sparse_graph_type(proc0Map,numEntriesPerRow (),
1576  constructorParams));
1577  if(myRank == rootRank) {
1578  typedef Teuchos::MatrixMarket::Raw::GraphElement<global_ordinal_type> element_type;
1579 
1580  // Get the entries
1581  const std::vector<element_type>& entries =
1582  pAdder->getAdder()->getEntries();
1583 
1584  // Insert them one at a time
1585  for(size_t curPos=0; curPos<entries.size(); curPos++) {
1586  const element_type& curEntry = entries[curPos];
1587  const global_ordinal_type curRow = curEntry.rowIndex()+indexBase;
1588  const global_ordinal_type curCol = curEntry.colIndex()+indexBase;
1589  Teuchos::ArrayView<const global_ordinal_type> colView(&curCol,1);
1590  proc0Graph->insertGlobalIndices(curRow,colView);
1591  }
1592  }
1593  proc0Graph->fillComplete();
1594 
1595  RCP<sparse_graph_type> distGraph;
1596  if(Teuchos::is_null(rowMap))
1597  {
1598  // Create a map describing the distribution we actually want
1599  RCP<map_type> distMap =
1600  rcp(new map_type(dims[0],0,pComm,GloballyDistributed));
1601 
1602  // Create the graph with that distribution too
1603  distGraph = rcp(new sparse_graph_type(distMap,colMap,0,constructorParams));
1604 
1605  // Create an importer/exporter/vandelay to redistribute the graph
1606  typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1607  import_type importer (proc0Map, distMap);
1608 
1609  // Import the data
1610  distGraph->doImport(*proc0Graph,importer,INSERT);
1611  }
1612  else {
1613  distGraph = rcp(new sparse_graph_type(rowMap,colMap,0,constructorParams));
1614 
1615  // Create an importer/exporter/vandelay to redistribute the graph
1616  typedef Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
1617  import_type importer (proc0Map, rowMap);
1618 
1619  // Import the data
1620  distGraph->doImport(*proc0Graph,importer,INSERT);
1621  }
1622 
1623  return distGraph;
1624  }
1625 
1626  public:
1650  static Teuchos::RCP<sparse_graph_type>
1651  readSparseGraphFile (const std::string& filename,
1652  const trcp_tcomm_t& comm,
1653  const bool callFillComplete=true,
1654  const bool tolerant=false,
1655  const bool debug=false)
1656  {
1657  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true);
1658 
1659  return readSparseGraph (in, comm,
1660  callFillComplete,
1661  tolerant, debug);
1662  // We can rely on the destructor of the input stream to close
1663  // the file on scope exit, even if readSparseGraph() throws an
1664  // exception.
1665  }
1666 
1667 
1696  static Teuchos::RCP<sparse_graph_type>
1697  readSparseGraphFile (const std::string& filename,
1698  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1699  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1700  const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1701  const bool tolerant=false,
1702  const bool debug=false)
1703  {
1704  std::ifstream in = Reader::openInFileOnRankZero(pComm, filename, true);
1705 
1706  return readSparseGraph (in, pComm,
1707  constructorParams,
1708  fillCompleteParams, tolerant, debug);
1709  // We can rely on the destructor of the input stream to close
1710  // the file on scope exit, even if readSparseGraph() throws an
1711  // exception.
1712  }
1713 
1714 
1752  static Teuchos::RCP<sparse_graph_type>
1753  readSparseGraphFile (const std::string& filename,
1754  const Teuchos::RCP<const map_type>& rowMap,
1755  Teuchos::RCP<const map_type>& colMap,
1756  const Teuchos::RCP<const map_type>& domainMap,
1757  const Teuchos::RCP<const map_type>& rangeMap,
1758  const bool callFillComplete=true,
1759  const bool tolerant=false,
1760  const bool debug=false)
1761  {
1762  TEUCHOS_TEST_FOR_EXCEPTION
1763  (rowMap.is_null (), std::invalid_argument,
1764  "Input rowMap must be nonnull.");
1765  trcp_tcomm_t comm = rowMap->getComm ();
1766  if (comm.is_null ()) {
1767  // If the input communicator is null on some process, then
1768  // that process does not participate in the collective.
1769  return Teuchos::null;
1770  }
1771 
1772  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true);
1773 
1774  return readSparseGraph (in, rowMap, colMap, domainMap, rangeMap,
1775  callFillComplete, tolerant, debug);
1776  }
1777 
1803  static Teuchos::RCP<sparse_graph_type>
1804  readSparseGraph (std::istream& in,
1805  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1806  const bool callFillComplete=true,
1807  const bool tolerant=false,
1808  const bool debug=false)
1809  {
1810  Teuchos::RCP<const map_type> fakeRowMap;
1811  Teuchos::RCP<const map_type> fakeColMap;
1812  Teuchos::RCP<Teuchos::ParameterList> fakeCtorParams;
1813 
1814  Teuchos::RCP<sparse_graph_type> graph =
1815  readSparseGraphHelper (in, pComm,
1816  fakeRowMap, fakeColMap,
1817  fakeCtorParams, tolerant, debug);
1818  if (callFillComplete) {
1819  graph->fillComplete ();
1820  }
1821  return graph;
1822  }
1823 
1824 
1854  static Teuchos::RCP<sparse_graph_type>
1855  readSparseGraph (std::istream& in,
1856  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1857  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1858  const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1859  const bool tolerant=false,
1860  const bool debug=false)
1861  {
1862  Teuchos::RCP<const map_type> fakeRowMap;
1863  Teuchos::RCP<const map_type> fakeColMap;
1864  Teuchos::RCP<sparse_graph_type> graph =
1865  readSparseGraphHelper (in, pComm,
1866  fakeRowMap, fakeColMap,
1867  constructorParams, tolerant, debug);
1868  graph->fillComplete (fillCompleteParams);
1869  return graph;
1870  }
1871 
1872 
1913  static Teuchos::RCP<sparse_graph_type>
1914  readSparseGraph (std::istream& in,
1915  const Teuchos::RCP<const map_type>& rowMap,
1916  Teuchos::RCP<const map_type>& colMap,
1917  const Teuchos::RCP<const map_type>& domainMap,
1918  const Teuchos::RCP<const map_type>& rangeMap,
1919  const bool callFillComplete=true,
1920  const bool tolerant=false,
1921  const bool debug=false)
1922  {
1923  Teuchos::RCP<sparse_graph_type> graph =
1924  readSparseGraphHelper (in, rowMap->getComm (),
1925  rowMap, colMap, Teuchos::null, tolerant,
1926  debug);
1927  if (callFillComplete) {
1928  graph->fillComplete (domainMap, rangeMap);
1929  }
1930  return graph;
1931  }
1932 
1933 #include "MatrixMarket_TpetraNew.hpp"
1934 
1958  static Teuchos::RCP<sparse_matrix_type>
1959  readSparseFile (const std::string& filename,
1960  const trcp_tcomm_t& comm,
1961  const bool callFillComplete=true,
1962  const bool tolerant=false,
1963  const bool debug=false)
1964  {
1965  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true);
1966 
1967  return readSparse (in, comm, callFillComplete, tolerant, debug);
1968  // We can rely on the destructor of the input stream to close
1969  // the file on scope exit, even if readSparse() throws an
1970  // exception.
1971  }
1972 
1973 
2002  static Teuchos::RCP<sparse_matrix_type>
2003  readSparseFile (const std::string& filename,
2004  const trcp_tcomm_t& comm,
2005  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2006  const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2007  const bool tolerant=false,
2008  const bool debug=false)
2009  {
2010  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true);
2011 
2012  return readSparse (in, comm, constructorParams,
2013  fillCompleteParams, tolerant, debug);
2014  }
2015 
2016 
2054  static Teuchos::RCP<sparse_matrix_type>
2055  readSparseFile (const std::string& filename,
2056  const Teuchos::RCP<const map_type>& rowMap,
2057  Teuchos::RCP<const map_type>& colMap,
2058  const Teuchos::RCP<const map_type>& domainMap,
2059  const Teuchos::RCP<const map_type>& rangeMap,
2060  const bool callFillComplete=true,
2061  const bool tolerant=false,
2062  const bool debug=false)
2063  {
2064  TEUCHOS_TEST_FOR_EXCEPTION(
2065  rowMap.is_null (), std::invalid_argument,
2066  "Row Map must be nonnull.");
2067 
2068  trcp_tcomm_t comm = rowMap->getComm ();
2069 
2070  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true);
2071 
2072  return readSparse (in, rowMap, colMap, domainMap, rangeMap,
2073  callFillComplete, tolerant, debug);
2074  }
2075 
2101  static Teuchos::RCP<sparse_matrix_type>
2102  readSparse (std::istream& in,
2103  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
2104  const bool callFillComplete=true,
2105  const bool tolerant=false,
2106  const bool debug=false)
2107  {
2108  using Teuchos::MatrixMarket::Banner;
2109  using Teuchos::arcp;
2110  using Teuchos::ArrayRCP;
2111  using Teuchos::broadcast;
2112  using Teuchos::null;
2113  using Teuchos::ptr;
2114  using Teuchos::RCP;
2115  using Teuchos::REDUCE_MAX;
2116  using Teuchos::reduceAll;
2117  using Teuchos::Tuple;
2118  using std::cerr;
2119  using std::endl;
2120  typedef Teuchos::ScalarTraits<scalar_type> STS;
2121 
2122  const bool extraDebug = false;
2123  const int myRank = pComm->getRank ();
2124  const int rootRank = 0;
2125 
2126  // Current line number in the input stream. Various calls
2127  // will modify this depending on the number of lines that are
2128  // read from the input stream. Only Rank 0 modifies this.
2129  size_t lineNumber = 1;
2130 
2131  if (debug && myRank == rootRank) {
2132  cerr << "Matrix Market reader: readSparse:" << endl
2133  << "-- Reading banner line" << endl;
2134  }
2135 
2136  // The "Banner" tells you whether the input stream represents
2137  // a sparse matrix, the symmetry type of the matrix, and the
2138  // type of the data it contains.
2139  //
2140  // pBanner will only be nonnull on MPI Rank 0. It will be
2141  // null on all other MPI processes.
2142  RCP<const Banner> pBanner;
2143  {
2144  // We read and validate the Banner on Proc 0, but broadcast
2145  // the validation result to all processes.
2146  // Teuchos::broadcast doesn't currently work with bool, so
2147  // we use int (true -> 1, false -> 0).
2148  int bannerIsCorrect = 1;
2149  std::ostringstream errMsg;
2150 
2151  if (myRank == rootRank) {
2152  // Read the Banner line from the input stream.
2153  try {
2154  pBanner = readBanner (in, lineNumber, tolerant, debug);
2155  }
2156  catch (std::exception& e) {
2157  errMsg << "Attempt to read the Matrix Market file's Banner line "
2158  "threw an exception: " << e.what();
2159  bannerIsCorrect = 0;
2160  }
2161 
2162  if (bannerIsCorrect) {
2163  // Validate the Banner for the case of a sparse matrix.
2164  // We validate on Proc 0, since it reads the Banner.
2165 
2166  // In intolerant mode, the matrix type must be "coordinate".
2167  if (! tolerant && pBanner->matrixType() != "coordinate") {
2168  bannerIsCorrect = 0;
2169  errMsg << "The Matrix Market input file must contain a "
2170  "\"coordinate\"-format sparse matrix in order to create a "
2171  "Tpetra::CrsMatrix object from it, but the file's matrix "
2172  "type is \"" << pBanner->matrixType() << "\" instead.";
2173  }
2174  // In tolerant mode, we allow the matrix type to be
2175  // anything other than "array" (which would mean that
2176  // the file contains a dense matrix).
2177  if (tolerant && pBanner->matrixType() == "array") {
2178  bannerIsCorrect = 0;
2179  errMsg << "Matrix Market file must contain a \"coordinate\"-"
2180  "format sparse matrix in order to create a Tpetra::CrsMatrix "
2181  "object from it, but the file's matrix type is \"array\" "
2182  "instead. That probably means the file contains dense matrix "
2183  "data.";
2184  }
2185  }
2186  } // Proc 0: Done reading the Banner, hopefully successfully.
2187 
2188  // Broadcast from Proc 0 whether the Banner was read correctly.
2189  broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2190 
2191  // If the Banner is invalid, all processes throw an
2192  // exception. Only Proc 0 gets the exception message, but
2193  // that's OK, since the main point is to "stop the world"
2194  // (rather than throw an exception on one process and leave
2195  // the others hanging).
2196  TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2197  std::invalid_argument, errMsg.str ());
2198  } // Done reading the Banner line and broadcasting success.
2199  if (debug && myRank == rootRank) {
2200  cerr << "-- Reading dimensions line" << endl;
2201  }
2202 
2203  // Read the matrix dimensions from the Matrix Market metadata.
2204  // dims = (numRows, numCols, numEntries). Proc 0 does the
2205  // reading, but it broadcasts the results to all MPI
2206  // processes. Thus, readCoordDims() is a collective
2207  // operation. It does a collective check for correctness too.
2208  Tuple<global_ordinal_type, 3> dims =
2209  readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2210 
2211  if (debug && myRank == rootRank) {
2212  cerr << "-- Making Adder for collecting matrix data" << endl;
2213  }
2214 
2215  // "Adder" object for collecting all the sparse matrix entries
2216  // from the input stream. This is only nonnull on Proc 0.
2217  RCP<adder_type> pAdder =
2218  makeAdder (pComm, pBanner, dims, tolerant, debug);
2219 
2220  if (debug && myRank == rootRank) {
2221  cerr << "-- Reading matrix data" << endl;
2222  }
2223  //
2224  // Read the matrix entries from the input stream on Proc 0.
2225  //
2226  {
2227  // We use readSuccess to broadcast the results of the read
2228  // (succeeded or not) to all MPI processes. Since
2229  // Teuchos::broadcast doesn't currently know how to send
2230  // bools, we convert to int (true -> 1, false -> 0).
2231  int readSuccess = 1;
2232  std::ostringstream errMsg; // Exception message (only valid on Proc 0)
2233  if (myRank == rootRank) {
2234  try {
2235  // Reader for "coordinate" format sparse matrix data.
2236  typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2237  global_ordinal_type, scalar_type, STS::isComplex> reader_type;
2238  reader_type reader (pAdder);
2239 
2240  // Read the sparse matrix entries.
2241  std::pair<bool, std::vector<size_t> > results =
2242  reader.read (in, lineNumber, tolerant, debug);
2243  readSuccess = results.first ? 1 : 0;
2244  }
2245  catch (std::exception& e) {
2246  readSuccess = 0;
2247  errMsg << e.what();
2248  }
2249  }
2250  broadcast (*pComm, rootRank, ptr (&readSuccess));
2251 
2252  // It would be nice to add a "verbose" flag, so that in
2253  // tolerant mode, we could log any bad line number(s) on
2254  // Proc 0. For now, we just throw if the read fails to
2255  // succeed.
2256  //
2257  // Question: If we're in tolerant mode, and if the read did
2258  // not succeed, should we attempt to call fillComplete()?
2259  TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2260  "Failed to read the Matrix Market sparse matrix file: "
2261  << errMsg.str());
2262  } // Done reading the matrix entries (stored on Proc 0 for now)
2263 
2264  if (debug && myRank == rootRank) {
2265  cerr << "-- Successfully read the Matrix Market data" << endl;
2266  }
2267 
2268  // In tolerant mode, we need to rebroadcast the matrix
2269  // dimensions, since they may be different after reading the
2270  // actual matrix data. We only need to broadcast the number
2271  // of rows and columns. Only Rank 0 needs to know the actual
2272  // global number of entries, since (a) we need to merge
2273  // duplicates on Rank 0 first anyway, and (b) when we
2274  // distribute the entries, each rank other than Rank 0 will
2275  // only need to know how many entries it owns, not the total
2276  // number of entries.
2277  if (tolerant) {
2278  if (debug && myRank == rootRank) {
2279  cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
2280  << endl
2281  << "----- Dimensions before: "
2282  << dims[0] << " x " << dims[1]
2283  << endl;
2284  }
2285  // Packed coordinate matrix dimensions (numRows, numCols).
2286  Tuple<global_ordinal_type, 2> updatedDims;
2287  if (myRank == rootRank) {
2288  // If one or more bottom rows of the matrix contain no
2289  // entries, then the Adder will report that the number
2290  // of rows is less than that specified in the
2291  // metadata. We allow this case, and favor the
2292  // metadata so that the zero row(s) will be included.
2293  updatedDims[0] =
2294  std::max (dims[0], pAdder->getAdder()->numRows());
2295  updatedDims[1] = pAdder->getAdder()->numCols();
2296  }
2297  broadcast (*pComm, rootRank, updatedDims);
2298  dims[0] = updatedDims[0];
2299  dims[1] = updatedDims[1];
2300  if (debug && myRank == rootRank) {
2301  cerr << "----- Dimensions after: " << dims[0] << " x "
2302  << dims[1] << endl;
2303  }
2304  }
2305  else {
2306  // In strict mode, we require that the matrix's metadata and
2307  // its actual data agree, at least somewhat. In particular,
2308  // the number of rows must agree, since otherwise we cannot
2309  // distribute the matrix correctly.
2310 
2311  // Teuchos::broadcast() doesn't know how to broadcast bools,
2312  // so we use an int with the standard 1 == true, 0 == false
2313  // encoding.
2314  int dimsMatch = 1;
2315  if (myRank == rootRank) {
2316  // If one or more bottom rows of the matrix contain no
2317  // entries, then the Adder will report that the number of
2318  // rows is less than that specified in the metadata. We
2319  // allow this case, and favor the metadata, but do not
2320  // allow the Adder to think there are more rows in the
2321  // matrix than the metadata says.
2322  if (dims[0] < pAdder->getAdder ()->numRows ()) {
2323  dimsMatch = 0;
2324  }
2325  }
2326  broadcast (*pComm, 0, ptr (&dimsMatch));
2327  if (dimsMatch == 0) {
2328  // We're in an error state anyway, so we might as well
2329  // work a little harder to print an informative error
2330  // message.
2331  //
2332  // Broadcast the Adder's idea of the matrix dimensions
2333  // from Proc 0 to all processes.
2334  Tuple<global_ordinal_type, 2> addersDims;
2335  if (myRank == rootRank) {
2336  addersDims[0] = pAdder->getAdder()->numRows();
2337  addersDims[1] = pAdder->getAdder()->numCols();
2338  }
2339  broadcast (*pComm, 0, addersDims);
2340  TEUCHOS_TEST_FOR_EXCEPTION(
2341  dimsMatch == 0, std::runtime_error,
2342  "The matrix metadata says that the matrix is " << dims[0] << " x "
2343  << dims[1] << ", but the actual data says that the matrix is "
2344  << addersDims[0] << " x " << addersDims[1] << ". That means the "
2345  "data includes more rows than reported in the metadata. This "
2346  "is not allowed when parsing in strict mode. Parse the matrix in "
2347  "tolerant mode to ignore the metadata when it disagrees with the "
2348  "data.");
2349  }
2350  } // Matrix dimensions (# rows, # cols, # entries) agree.
2351 
2352  if (debug && myRank == rootRank) {
2353  cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
2354  }
2355 
2356  // Now that we've read in all the matrix entries from the
2357  // input stream into the adder on Proc 0, post-process them
2358  // into CSR format (still on Proc 0). This will facilitate
2359  // distributing them to all the processors.
2360  //
2361  // These arrays represent the global matrix data as a CSR
2362  // matrix (with numEntriesPerRow as redundant but convenient
2363  // metadata, since it's computable from rowPtr and vice
2364  // versa). They are valid only on Proc 0.
2365  ArrayRCP<size_t> numEntriesPerRow;
2366  ArrayRCP<size_t> rowPtr;
2367  ArrayRCP<global_ordinal_type> colInd;
2368  ArrayRCP<scalar_type> values;
2369 
2370  // Proc 0 first merges duplicate entries, and then converts
2371  // the coordinate-format matrix data to CSR.
2372  {
2373  int mergeAndConvertSucceeded = 1;
2374  std::ostringstream errMsg;
2375 
2376  if (myRank == rootRank) {
2377  try {
2378  typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
2379  global_ordinal_type> element_type;
2380 
2381  // Number of rows in the matrix. If we are in tolerant
2382  // mode, we've already synchronized dims with the actual
2383  // matrix data. If in strict mode, we should use dims
2384  // (as read from the file's metadata) rather than the
2385  // matrix data to determine the dimensions. (The matrix
2386  // data will claim fewer rows than the metadata, if one
2387  // or more rows have no entries stored in the file.)
2388  const size_type numRows = dims[0];
2389 
2390  // Additively merge duplicate matrix entries.
2391  pAdder->getAdder()->merge ();
2392 
2393  // Get a temporary const view of the merged matrix entries.
2394  const std::vector<element_type>& entries =
2395  pAdder->getAdder()->getEntries();
2396 
2397  // Number of matrix entries (after merging).
2398  const size_t numEntries = (size_t)entries.size();
2399 
2400  if (debug) {
2401  cerr << "----- Proc 0: Matrix has numRows=" << numRows
2402  << " rows and numEntries=" << numEntries
2403  << " entries." << endl;
2404  }
2405 
2406  // Make space for the CSR matrix data. Converting to
2407  // CSR is easier if we fill numEntriesPerRow with zeros
2408  // at first.
2409  numEntriesPerRow = arcp<size_t> (numRows);
2410  std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2411  rowPtr = arcp<size_t> (numRows+1);
2412  std::fill (rowPtr.begin(), rowPtr.end(), 0);
2413  colInd = arcp<global_ordinal_type> (numEntries);
2414  values = arcp<scalar_type> (numEntries);
2415 
2416  // Convert from array-of-structs coordinate format to CSR
2417  // (compressed sparse row) format.
2418  global_ordinal_type prvRow = 0;
2419  size_t curPos = 0;
2420  rowPtr[0] = 0;
2421  for (curPos = 0; curPos < numEntries; ++curPos) {
2422  const element_type& curEntry = entries[curPos];
2423  const global_ordinal_type curRow = curEntry.rowIndex();
2424  TEUCHOS_TEST_FOR_EXCEPTION(
2425  curRow < prvRow, std::logic_error,
2426  "Row indices are out of order, even though they are supposed "
2427  "to be sorted. curRow = " << curRow << ", prvRow = "
2428  << prvRow << ", at curPos = " << curPos << ". Please report "
2429  "this bug to the Tpetra developers.");
2430  if (curRow > prvRow) {
2431  for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
2432  rowPtr[r] = curPos;
2433  }
2434  prvRow = curRow;
2435  }
2436  numEntriesPerRow[curRow]++;
2437  colInd[curPos] = curEntry.colIndex();
2438  values[curPos] = curEntry.value();
2439  }
2440  // rowPtr has one more entry than numEntriesPerRow. The
2441  // last entry of rowPtr is the number of entries in
2442  // colInd and values.
2443  rowPtr[numRows] = numEntries;
2444  } // Finished conversion to CSR format
2445  catch (std::exception& e) {
2446  mergeAndConvertSucceeded = 0;
2447  errMsg << "Failed to merge sparse matrix entries and convert to "
2448  "CSR format: " << e.what();
2449  }
2450 
2451  if (debug && mergeAndConvertSucceeded) {
2452  // Number of rows in the matrix.
2453  const size_type numRows = dims[0];
2454  const size_type maxToDisplay = 100;
2455 
2456  cerr << "----- Proc 0: numEntriesPerRow[0.."
2457  << (numEntriesPerRow.size()-1) << "] ";
2458  if (numRows > maxToDisplay) {
2459  cerr << "(only showing first and last few entries) ";
2460  }
2461  cerr << "= [";
2462  if (numRows > 0) {
2463  if (numRows > maxToDisplay) {
2464  for (size_type k = 0; k < 2; ++k) {
2465  cerr << numEntriesPerRow[k] << " ";
2466  }
2467  cerr << "... ";
2468  for (size_type k = numRows-2; k < numRows-1; ++k) {
2469  cerr << numEntriesPerRow[k] << " ";
2470  }
2471  }
2472  else {
2473  for (size_type k = 0; k < numRows-1; ++k) {
2474  cerr << numEntriesPerRow[k] << " ";
2475  }
2476  }
2477  cerr << numEntriesPerRow[numRows-1];
2478  } // numRows > 0
2479  cerr << "]" << endl;
2480 
2481  cerr << "----- Proc 0: rowPtr ";
2482  if (numRows > maxToDisplay) {
2483  cerr << "(only showing first and last few entries) ";
2484  }
2485  cerr << "= [";
2486  if (numRows > maxToDisplay) {
2487  for (size_type k = 0; k < 2; ++k) {
2488  cerr << rowPtr[k] << " ";
2489  }
2490  cerr << "... ";
2491  for (size_type k = numRows-2; k < numRows; ++k) {
2492  cerr << rowPtr[k] << " ";
2493  }
2494  }
2495  else {
2496  for (size_type k = 0; k < numRows; ++k) {
2497  cerr << rowPtr[k] << " ";
2498  }
2499  }
2500  cerr << rowPtr[numRows] << "]" << endl;
2501  }
2502  } // if myRank == rootRank
2503  } // Done converting sparse matrix data to CSR format
2504 
2505  // Now we're done with the Adder, so we can release the
2506  // reference ("free" it) to save space. This only actually
2507  // does anything on Rank 0, since pAdder is null on all the
2508  // other MPI processes.
2509  pAdder = null;
2510 
2511  if (debug && myRank == rootRank) {
2512  cerr << "-- Making range, domain, and row maps" << endl;
2513  }
2514 
2515  // Make the maps that describe the matrix's range and domain,
2516  // and the distribution of its rows. Creating a Map is a
2517  // collective operation, so we don't have to do a broadcast of
2518  // a success Boolean.
2519  RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
2520  RCP<const map_type> pDomainMap =
2521  makeDomainMap (pRangeMap, dims[0], dims[1]);
2522  RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
2523 
2524  if (debug && myRank == rootRank) {
2525  cerr << "-- Distributing the matrix data" << endl;
2526  }
2527 
2528  // Distribute the matrix data. Each processor has to add the
2529  // rows that it owns. If you try to make Proc 0 call
2530  // insertGlobalValues() for _all_ the rows, not just those it
2531  // owns, then fillComplete() will compute the number of
2532  // columns incorrectly. That's why Proc 0 has to distribute
2533  // the matrix data and why we make all the processors (not
2534  // just Proc 0) call insertGlobalValues() on their own data.
2535  //
2536  // These arrays represent each processor's part of the matrix
2537  // data, in "CSR" format (sort of, since the row indices might
2538  // not be contiguous).
2539  ArrayRCP<size_t> myNumEntriesPerRow;
2540  ArrayRCP<size_t> myRowPtr;
2541  ArrayRCP<global_ordinal_type> myColInd;
2542  ArrayRCP<scalar_type> myValues;
2543  // Distribute the matrix data. This is a collective operation.
2544  distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
2545  numEntriesPerRow, rowPtr, colInd, values, debug);
2546 
2547  if (debug && myRank == rootRank) {
2548  cerr << "-- Inserting matrix entries on each processor";
2549  if (callFillComplete) {
2550  cerr << " and calling fillComplete()";
2551  }
2552  cerr << endl;
2553  }
2554  // Each processor inserts its part of the matrix data, and
2555  // then they all call fillComplete(). This method invalidates
2556  // the my* distributed matrix data before calling
2557  // fillComplete(), in order to save space. In general, we
2558  // never store more than two copies of the matrix's entries in
2559  // memory at once, which is no worse than what Tpetra
2560  // promises.
2561  RCP<sparse_matrix_type> pMatrix =
2562  makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
2563  pRowMap, pRangeMap, pDomainMap, callFillComplete);
2564  // Only use a reduce-all in debug mode to check if pMatrix is
2565  // null. Otherwise, just throw an exception. We never expect
2566  // a null pointer here, so we can save a communication.
2567  if (debug) {
2568  int localIsNull = pMatrix.is_null () ? 1 : 0;
2569  int globalIsNull = 0;
2570  reduceAll (*pComm, REDUCE_MAX, localIsNull, ptr (&globalIsNull));
2571  TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
2572  "Reader::makeMatrix() returned a null pointer on at least one "
2573  "process. Please report this bug to the Tpetra developers.");
2574  }
2575  else {
2576  TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
2577  "Reader::makeMatrix() returned a null pointer. "
2578  "Please report this bug to the Tpetra developers.");
2579  }
2580 
2581  // We can't get the dimensions of the matrix until after
2582  // fillComplete() is called. Thus, we can't do the sanity
2583  // check (dimensions read from the Matrix Market data,
2584  // vs. dimensions reported by the CrsMatrix) unless the user
2585  // asked makeMatrix() to call fillComplete().
2586  //
2587  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
2588  // what one might think it does, so you have to ask the range
2589  // resp. domain map for the number of rows resp. columns.
2590  if (callFillComplete) {
2591  const int numProcs = pComm->getSize ();
2592 
2593  if (extraDebug && debug) {
2594  const global_size_t globalNumRows =
2595  pRangeMap->getGlobalNumElements ();
2596  const global_size_t globalNumCols =
2597  pDomainMap->getGlobalNumElements ();
2598  if (myRank == rootRank) {
2599  cerr << "-- Matrix is "
2600  << globalNumRows << " x " << globalNumCols
2601  << " with " << pMatrix->getGlobalNumEntries()
2602  << " entries, and index base "
2603  << pMatrix->getIndexBase() << "." << endl;
2604  }
2605  pComm->barrier ();
2606  for (int p = 0; p < numProcs; ++p) {
2607  if (myRank == p) {
2608  cerr << "-- Proc " << p << " owns "
2609  << pMatrix->getLocalNumCols() << " columns, and "
2610  << pMatrix->getLocalNumEntries() << " entries." << endl;
2611  }
2612  pComm->barrier ();
2613  }
2614  } // if (extraDebug && debug)
2615  } // if (callFillComplete)
2616 
2617  if (debug && myRank == rootRank) {
2618  cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
2619  << endl;
2620  }
2621  return pMatrix;
2622  }
2623 
2624 
2653  static Teuchos::RCP<sparse_matrix_type>
2654  readSparse (std::istream& in,
2655  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
2656  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
2657  const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
2658  const bool tolerant=false,
2659  const bool debug=false)
2660  {
2661  using Teuchos::MatrixMarket::Banner;
2662  using Teuchos::arcp;
2663  using Teuchos::ArrayRCP;
2664  using Teuchos::broadcast;
2665  using Teuchos::null;
2666  using Teuchos::ptr;
2667  using Teuchos::RCP;
2668  using Teuchos::reduceAll;
2669  using Teuchos::Tuple;
2670  using std::cerr;
2671  using std::endl;
2672  typedef Teuchos::ScalarTraits<scalar_type> STS;
2673 
2674  const bool extraDebug = false;
2675  const int myRank = pComm->getRank ();
2676  const int rootRank = 0;
2677 
2678  // Current line number in the input stream. Various calls
2679  // will modify this depending on the number of lines that are
2680  // read from the input stream. Only Rank 0 modifies this.
2681  size_t lineNumber = 1;
2682 
2683  if (debug && myRank == rootRank) {
2684  cerr << "Matrix Market reader: readSparse:" << endl
2685  << "-- Reading banner line" << endl;
2686  }
2687 
2688  // The "Banner" tells you whether the input stream represents
2689  // a sparse matrix, the symmetry type of the matrix, and the
2690  // type of the data it contains.
2691  //
2692  // pBanner will only be nonnull on MPI Rank 0. It will be
2693  // null on all other MPI processes.
2694  RCP<const Banner> pBanner;
2695  {
2696  // We read and validate the Banner on Proc 0, but broadcast
2697  // the validation result to all processes.
2698  // Teuchos::broadcast doesn't currently work with bool, so
2699  // we use int (true -> 1, false -> 0).
2700  int bannerIsCorrect = 1;
2701  std::ostringstream errMsg;
2702 
2703  if (myRank == rootRank) {
2704  // Read the Banner line from the input stream.
2705  try {
2706  pBanner = readBanner (in, lineNumber, tolerant, debug);
2707  }
2708  catch (std::exception& e) {
2709  errMsg << "Attempt to read the Matrix Market file's Banner line "
2710  "threw an exception: " << e.what();
2711  bannerIsCorrect = 0;
2712  }
2713 
2714  if (bannerIsCorrect) {
2715  // Validate the Banner for the case of a sparse matrix.
2716  // We validate on Proc 0, since it reads the Banner.
2717 
2718  // In intolerant mode, the matrix type must be "coordinate".
2719  if (! tolerant && pBanner->matrixType() != "coordinate") {
2720  bannerIsCorrect = 0;
2721  errMsg << "The Matrix Market input file must contain a "
2722  "\"coordinate\"-format sparse matrix in order to create a "
2723  "Tpetra::CrsMatrix object from it, but the file's matrix "
2724  "type is \"" << pBanner->matrixType() << "\" instead.";
2725  }
2726  // In tolerant mode, we allow the matrix type to be
2727  // anything other than "array" (which would mean that
2728  // the file contains a dense matrix).
2729  if (tolerant && pBanner->matrixType() == "array") {
2730  bannerIsCorrect = 0;
2731  errMsg << "Matrix Market file must contain a \"coordinate\"-"
2732  "format sparse matrix in order to create a Tpetra::CrsMatrix "
2733  "object from it, but the file's matrix type is \"array\" "
2734  "instead. That probably means the file contains dense matrix "
2735  "data.";
2736  }
2737  }
2738  } // Proc 0: Done reading the Banner, hopefully successfully.
2739 
2740  // Broadcast from Proc 0 whether the Banner was read correctly.
2741  broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2742 
2743  // If the Banner is invalid, all processes throw an
2744  // exception. Only Proc 0 gets the exception message, but
2745  // that's OK, since the main point is to "stop the world"
2746  // (rather than throw an exception on one process and leave
2747  // the others hanging).
2748  TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2749  std::invalid_argument, errMsg.str ());
2750  } // Done reading the Banner line and broadcasting success.
2751  if (debug && myRank == rootRank) {
2752  cerr << "-- Reading dimensions line" << endl;
2753  }
2754 
2755  // Read the matrix dimensions from the Matrix Market metadata.
2756  // dims = (numRows, numCols, numEntries). Proc 0 does the
2757  // reading, but it broadcasts the results to all MPI
2758  // processes. Thus, readCoordDims() is a collective
2759  // operation. It does a collective check for correctness too.
2760  Tuple<global_ordinal_type, 3> dims =
2761  readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2762 
2763  if (debug && myRank == rootRank) {
2764  cerr << "-- Making Adder for collecting matrix data" << endl;
2765  }
2766 
2767  // "Adder" object for collecting all the sparse matrix entries
2768  // from the input stream. This is only nonnull on Proc 0.
2769  RCP<adder_type> pAdder =
2770  makeAdder (pComm, pBanner, dims, tolerant, debug);
2771 
2772  if (debug && myRank == rootRank) {
2773  cerr << "-- Reading matrix data" << endl;
2774  }
2775  //
2776  // Read the matrix entries from the input stream on Proc 0.
2777  //
2778  {
2779  // We use readSuccess to broadcast the results of the read
2780  // (succeeded or not) to all MPI processes. Since
2781  // Teuchos::broadcast doesn't currently know how to send
2782  // bools, we convert to int (true -> 1, false -> 0).
2783  int readSuccess = 1;
2784  std::ostringstream errMsg; // Exception message (only valid on Proc 0)
2785  if (myRank == rootRank) {
2786  try {
2787  // Reader for "coordinate" format sparse matrix data.
2788  typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2789  global_ordinal_type, scalar_type, STS::isComplex> reader_type;
2790  reader_type reader (pAdder);
2791 
2792  // Read the sparse matrix entries.
2793  std::pair<bool, std::vector<size_t> > results =
2794  reader.read (in, lineNumber, tolerant, debug);
2795  readSuccess = results.first ? 1 : 0;
2796  }
2797  catch (std::exception& e) {
2798  readSuccess = 0;
2799  errMsg << e.what();
2800  }
2801  }
2802  broadcast (*pComm, rootRank, ptr (&readSuccess));
2803 
2804  // It would be nice to add a "verbose" flag, so that in
2805  // tolerant mode, we could log any bad line number(s) on
2806  // Proc 0. For now, we just throw if the read fails to
2807  // succeed.
2808  //
2809  // Question: If we're in tolerant mode, and if the read did
2810  // not succeed, should we attempt to call fillComplete()?
2811  TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2812  "Failed to read the Matrix Market sparse matrix file: "
2813  << errMsg.str());
2814  } // Done reading the matrix entries (stored on Proc 0 for now)
2815 
2816  if (debug && myRank == rootRank) {
2817  cerr << "-- Successfully read the Matrix Market data" << endl;
2818  }
2819 
2820  // In tolerant mode, we need to rebroadcast the matrix
2821  // dimensions, since they may be different after reading the
2822  // actual matrix data. We only need to broadcast the number
2823  // of rows and columns. Only Rank 0 needs to know the actual
2824  // global number of entries, since (a) we need to merge
2825  // duplicates on Rank 0 first anyway, and (b) when we
2826  // distribute the entries, each rank other than Rank 0 will
2827  // only need to know how many entries it owns, not the total
2828  // number of entries.
2829  if (tolerant) {
2830  if (debug && myRank == rootRank) {
2831  cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
2832  << endl
2833  << "----- Dimensions before: "
2834  << dims[0] << " x " << dims[1]
2835  << endl;
2836  }
2837  // Packed coordinate matrix dimensions (numRows, numCols).
2838  Tuple<global_ordinal_type, 2> updatedDims;
2839  if (myRank == rootRank) {
2840  // If one or more bottom rows of the matrix contain no
2841  // entries, then the Adder will report that the number
2842  // of rows is less than that specified in the
2843  // metadata. We allow this case, and favor the
2844  // metadata so that the zero row(s) will be included.
2845  updatedDims[0] =
2846  std::max (dims[0], pAdder->getAdder()->numRows());
2847  updatedDims[1] = pAdder->getAdder()->numCols();
2848  }
2849  broadcast (*pComm, rootRank, updatedDims);
2850  dims[0] = updatedDims[0];
2851  dims[1] = updatedDims[1];
2852  if (debug && myRank == rootRank) {
2853  cerr << "----- Dimensions after: " << dims[0] << " x "
2854  << dims[1] << endl;
2855  }
2856  }
2857  else {
2858  // In strict mode, we require that the matrix's metadata and
2859  // its actual data agree, at least somewhat. In particular,
2860  // the number of rows must agree, since otherwise we cannot
2861  // distribute the matrix correctly.
2862 
2863  // Teuchos::broadcast() doesn't know how to broadcast bools,
2864  // so we use an int with the standard 1 == true, 0 == false
2865  // encoding.
2866  int dimsMatch = 1;
2867  if (myRank == rootRank) {
2868  // If one or more bottom rows of the matrix contain no
2869  // entries, then the Adder will report that the number of
2870  // rows is less than that specified in the metadata. We
2871  // allow this case, and favor the metadata, but do not
2872  // allow the Adder to think there are more rows in the
2873  // matrix than the metadata says.
2874  if (dims[0] < pAdder->getAdder ()->numRows ()) {
2875  dimsMatch = 0;
2876  }
2877  }
2878  broadcast (*pComm, 0, ptr (&dimsMatch));
2879  if (dimsMatch == 0) {
2880  // We're in an error state anyway, so we might as well
2881  // work a little harder to print an informative error
2882  // message.
2883  //
2884  // Broadcast the Adder's idea of the matrix dimensions
2885  // from Proc 0 to all processes.
2886  Tuple<global_ordinal_type, 2> addersDims;
2887  if (myRank == rootRank) {
2888  addersDims[0] = pAdder->getAdder()->numRows();
2889  addersDims[1] = pAdder->getAdder()->numCols();
2890  }
2891  broadcast (*pComm, 0, addersDims);
2892  TEUCHOS_TEST_FOR_EXCEPTION(
2893  dimsMatch == 0, std::runtime_error,
2894  "The matrix metadata says that the matrix is " << dims[0] << " x "
2895  << dims[1] << ", but the actual data says that the matrix is "
2896  << addersDims[0] << " x " << addersDims[1] << ". That means the "
2897  "data includes more rows than reported in the metadata. This "
2898  "is not allowed when parsing in strict mode. Parse the matrix in "
2899  "tolerant mode to ignore the metadata when it disagrees with the "
2900  "data.");
2901  }
2902  } // Matrix dimensions (# rows, # cols, # entries) agree.
2903 
2904  if (debug && myRank == rootRank) {
2905  cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
2906  }
2907 
2908  // Now that we've read in all the matrix entries from the
2909  // input stream into the adder on Proc 0, post-process them
2910  // into CSR format (still on Proc 0). This will facilitate
2911  // distributing them to all the processors.
2912  //
2913  // These arrays represent the global matrix data as a CSR
2914  // matrix (with numEntriesPerRow as redundant but convenient
2915  // metadata, since it's computable from rowPtr and vice
2916  // versa). They are valid only on Proc 0.
2917  ArrayRCP<size_t> numEntriesPerRow;
2918  ArrayRCP<size_t> rowPtr;
2919  ArrayRCP<global_ordinal_type> colInd;
2920  ArrayRCP<scalar_type> values;
2921 
2922  // Proc 0 first merges duplicate entries, and then converts
2923  // the coordinate-format matrix data to CSR.
2924  {
2925  int mergeAndConvertSucceeded = 1;
2926  std::ostringstream errMsg;
2927 
2928  if (myRank == rootRank) {
2929  try {
2930  typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
2931  global_ordinal_type> element_type;
2932 
2933  // Number of rows in the matrix. If we are in tolerant
2934  // mode, we've already synchronized dims with the actual
2935  // matrix data. If in strict mode, we should use dims
2936  // (as read from the file's metadata) rather than the
2937  // matrix data to determine the dimensions. (The matrix
2938  // data will claim fewer rows than the metadata, if one
2939  // or more rows have no entries stored in the file.)
2940  const size_type numRows = dims[0];
2941 
2942  // Additively merge duplicate matrix entries.
2943  pAdder->getAdder()->merge ();
2944 
2945  // Get a temporary const view of the merged matrix entries.
2946  const std::vector<element_type>& entries =
2947  pAdder->getAdder()->getEntries();
2948 
2949  // Number of matrix entries (after merging).
2950  const size_t numEntries = (size_t)entries.size();
2951 
2952  if (debug) {
2953  cerr << "----- Proc 0: Matrix has numRows=" << numRows
2954  << " rows and numEntries=" << numEntries
2955  << " entries." << endl;
2956  }
2957 
2958  // Make space for the CSR matrix data. Converting to
2959  // CSR is easier if we fill numEntriesPerRow with zeros
2960  // at first.
2961  numEntriesPerRow = arcp<size_t> (numRows);
2962  std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2963  rowPtr = arcp<size_t> (numRows+1);
2964  std::fill (rowPtr.begin(), rowPtr.end(), 0);
2965  colInd = arcp<global_ordinal_type> (numEntries);
2966  values = arcp<scalar_type> (numEntries);
2967 
2968  // Convert from array-of-structs coordinate format to CSR
2969  // (compressed sparse row) format.
2970  global_ordinal_type prvRow = 0;
2971  size_t curPos = 0;
2972  rowPtr[0] = 0;
2973  for (curPos = 0; curPos < numEntries; ++curPos) {
2974  const element_type& curEntry = entries[curPos];
2975  const global_ordinal_type curRow = curEntry.rowIndex();
2976  TEUCHOS_TEST_FOR_EXCEPTION(
2977  curRow < prvRow, std::logic_error,
2978  "Row indices are out of order, even though they are supposed "
2979  "to be sorted. curRow = " << curRow << ", prvRow = "
2980  << prvRow << ", at curPos = " << curPos << ". Please report "
2981  "this bug to the Tpetra developers.");
2982  if (curRow > prvRow) {
2983  for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
2984  rowPtr[r] = curPos;
2985  }
2986  prvRow = curRow;
2987  }
2988  numEntriesPerRow[curRow]++;
2989  colInd[curPos] = curEntry.colIndex();
2990  values[curPos] = curEntry.value();
2991  }
2992  // rowPtr has one more entry than numEntriesPerRow. The
2993  // last entry of rowPtr is the number of entries in
2994  // colInd and values.
2995  rowPtr[numRows] = numEntries;
2996  } // Finished conversion to CSR format
2997  catch (std::exception& e) {
2998  mergeAndConvertSucceeded = 0;
2999  errMsg << "Failed to merge sparse matrix entries and convert to "
3000  "CSR format: " << e.what();
3001  }
3002 
3003  if (debug && mergeAndConvertSucceeded) {
3004  // Number of rows in the matrix.
3005  const size_type numRows = dims[0];
3006  const size_type maxToDisplay = 100;
3007 
3008  cerr << "----- Proc 0: numEntriesPerRow[0.."
3009  << (numEntriesPerRow.size()-1) << "] ";
3010  if (numRows > maxToDisplay) {
3011  cerr << "(only showing first and last few entries) ";
3012  }
3013  cerr << "= [";
3014  if (numRows > 0) {
3015  if (numRows > maxToDisplay) {
3016  for (size_type k = 0; k < 2; ++k) {
3017  cerr << numEntriesPerRow[k] << " ";
3018  }
3019  cerr << "... ";
3020  for (size_type k = numRows-2; k < numRows-1; ++k) {
3021  cerr << numEntriesPerRow[k] << " ";
3022  }
3023  }
3024  else {
3025  for (size_type k = 0; k < numRows-1; ++k) {
3026  cerr << numEntriesPerRow[k] << " ";
3027  }
3028  }
3029  cerr << numEntriesPerRow[numRows-1];
3030  } // numRows > 0
3031  cerr << "]" << endl;
3032 
3033  cerr << "----- Proc 0: rowPtr ";
3034  if (numRows > maxToDisplay) {
3035  cerr << "(only showing first and last few entries) ";
3036  }
3037  cerr << "= [";
3038  if (numRows > maxToDisplay) {
3039  for (size_type k = 0; k < 2; ++k) {
3040  cerr << rowPtr[k] << " ";
3041  }
3042  cerr << "... ";
3043  for (size_type k = numRows-2; k < numRows; ++k) {
3044  cerr << rowPtr[k] << " ";
3045  }
3046  }
3047  else {
3048  for (size_type k = 0; k < numRows; ++k) {
3049  cerr << rowPtr[k] << " ";
3050  }
3051  }
3052  cerr << rowPtr[numRows] << "]" << endl;
3053  }
3054  } // if myRank == rootRank
3055  } // Done converting sparse matrix data to CSR format
3056 
3057  // Now we're done with the Adder, so we can release the
3058  // reference ("free" it) to save space. This only actually
3059  // does anything on Rank 0, since pAdder is null on all the
3060  // other MPI processes.
3061  pAdder = null;
3062 
3063  if (debug && myRank == rootRank) {
3064  cerr << "-- Making range, domain, and row maps" << endl;
3065  }
3066 
3067  // Make the maps that describe the matrix's range and domain,
3068  // and the distribution of its rows. Creating a Map is a
3069  // collective operation, so we don't have to do a broadcast of
3070  // a success Boolean.
3071  RCP<const map_type> pRangeMap = makeRangeMap (pComm, dims[0]);
3072  RCP<const map_type> pDomainMap =
3073  makeDomainMap (pRangeMap, dims[0], dims[1]);
3074  RCP<const map_type> pRowMap = makeRowMap (null, pComm, dims[0]);
3075 
3076  if (debug && myRank == rootRank) {
3077  cerr << "-- Distributing the matrix data" << endl;
3078  }
3079 
3080  // Distribute the matrix data. Each processor has to add the
3081  // rows that it owns. If you try to make Proc 0 call
3082  // insertGlobalValues() for _all_ the rows, not just those it
3083  // owns, then fillComplete() will compute the number of
3084  // columns incorrectly. That's why Proc 0 has to distribute
3085  // the matrix data and why we make all the processors (not
3086  // just Proc 0) call insertGlobalValues() on their own data.
3087  //
3088  // These arrays represent each processor's part of the matrix
3089  // data, in "CSR" format (sort of, since the row indices might
3090  // not be contiguous).
3091  ArrayRCP<size_t> myNumEntriesPerRow;
3092  ArrayRCP<size_t> myRowPtr;
3093  ArrayRCP<global_ordinal_type> myColInd;
3094  ArrayRCP<scalar_type> myValues;
3095  // Distribute the matrix data. This is a collective operation.
3096  distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
3097  numEntriesPerRow, rowPtr, colInd, values, debug);
3098 
3099  if (debug && myRank == rootRank) {
3100  cerr << "-- Inserting matrix entries on each process "
3101  "and calling fillComplete()" << endl;
3102  }
3103  // Each processor inserts its part of the matrix data, and
3104  // then they all call fillComplete(). This method invalidates
3105  // the my* distributed matrix data before calling
3106  // fillComplete(), in order to save space. In general, we
3107  // never store more than two copies of the matrix's entries in
3108  // memory at once, which is no worse than what Tpetra
3109  // promises.
3110  Teuchos::RCP<sparse_matrix_type> pMatrix =
3111  makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
3112  pRowMap, pRangeMap, pDomainMap, constructorParams,
3113  fillCompleteParams);
3114  // Only use a reduce-all in debug mode to check if pMatrix is
3115  // null. Otherwise, just throw an exception. We never expect
3116  // a null pointer here, so we can save a communication.
3117  if (debug) {
3118  int localIsNull = pMatrix.is_null () ? 1 : 0;
3119  int globalIsNull = 0;
3120  reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
3121  TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
3122  "Reader::makeMatrix() returned a null pointer on at least one "
3123  "process. Please report this bug to the Tpetra developers.");
3124  }
3125  else {
3126  TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
3127  "Reader::makeMatrix() returned a null pointer. "
3128  "Please report this bug to the Tpetra developers.");
3129  }
3130 
3131  // Sanity check for dimensions (read from the Matrix Market
3132  // data, vs. dimensions reported by the CrsMatrix).
3133  //
3134  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
3135  // what one might think it does, so you have to ask the range
3136  // resp. domain map for the number of rows resp. columns.
3137  if (extraDebug && debug) {
3138  const int numProcs = pComm->getSize ();
3139  const global_size_t globalNumRows =
3140  pRangeMap->getGlobalNumElements();
3141  const global_size_t globalNumCols =
3142  pDomainMap->getGlobalNumElements();
3143  if (myRank == rootRank) {
3144  cerr << "-- Matrix is "
3145  << globalNumRows << " x " << globalNumCols
3146  << " with " << pMatrix->getGlobalNumEntries()
3147  << " entries, and index base "
3148  << pMatrix->getIndexBase() << "." << endl;
3149  }
3150  pComm->barrier ();
3151  for (int p = 0; p < numProcs; ++p) {
3152  if (myRank == p) {
3153  cerr << "-- Proc " << p << " owns "
3154  << pMatrix->getLocalNumCols() << " columns, and "
3155  << pMatrix->getLocalNumEntries() << " entries." << endl;
3156  }
3157  pComm->barrier ();
3158  }
3159  } // if (extraDebug && debug)
3160 
3161  if (debug && myRank == rootRank) {
3162  cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
3163  << endl;
3164  }
3165  return pMatrix;
3166  }
3167 
3168 
3209  static Teuchos::RCP<sparse_matrix_type>
3210  readSparse (std::istream& in,
3211  const Teuchos::RCP<const map_type>& rowMap,
3212  Teuchos::RCP<const map_type>& colMap,
3213  const Teuchos::RCP<const map_type>& domainMap,
3214  const Teuchos::RCP<const map_type>& rangeMap,
3215  const bool callFillComplete=true,
3216  const bool tolerant=false,
3217  const bool debug=false)
3218  {
3219  using Teuchos::MatrixMarket::Banner;
3220  using Teuchos::arcp;
3221  using Teuchos::ArrayRCP;
3222  using Teuchos::ArrayView;
3223  using Teuchos::as;
3224  using Teuchos::broadcast;
3225  using Teuchos::Comm;
3226  using Teuchos::null;
3227  using Teuchos::ptr;
3228  using Teuchos::RCP;
3229  using Teuchos::reduceAll;
3230  using Teuchos::Tuple;
3231  using std::cerr;
3232  using std::endl;
3233  typedef Teuchos::ScalarTraits<scalar_type> STS;
3234 
3235  RCP<const Comm<int> > pComm = rowMap->getComm ();
3236  const int myRank = pComm->getRank ();
3237  const int rootRank = 0;
3238  const bool extraDebug = false;
3239 
3240  // Fast checks for invalid input. We can't check other
3241  // attributes of the Maps until we've read in the matrix
3242  // dimensions.
3243  TEUCHOS_TEST_FOR_EXCEPTION(
3244  rowMap.is_null (), std::invalid_argument,
3245  "Row Map must be nonnull.");
3246  TEUCHOS_TEST_FOR_EXCEPTION(
3247  rangeMap.is_null (), std::invalid_argument,
3248  "Range Map must be nonnull.");
3249  TEUCHOS_TEST_FOR_EXCEPTION(
3250  domainMap.is_null (), std::invalid_argument,
3251  "Domain Map must be nonnull.");
3252  TEUCHOS_TEST_FOR_EXCEPTION(
3253  rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
3254  std::invalid_argument,
3255  "The specified row Map's communicator (rowMap->getComm())"
3256  "differs from the given separately supplied communicator pComm.");
3257  TEUCHOS_TEST_FOR_EXCEPTION(
3258  domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
3259  std::invalid_argument,
3260  "The specified domain Map's communicator (domainMap->getComm())"
3261  "differs from the given separately supplied communicator pComm.");
3262  TEUCHOS_TEST_FOR_EXCEPTION(
3263  rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
3264  std::invalid_argument,
3265  "The specified range Map's communicator (rangeMap->getComm())"
3266  "differs from the given separately supplied communicator pComm.");
3267 
3268  // Current line number in the input stream. Various calls
3269  // will modify this depending on the number of lines that are
3270  // read from the input stream. Only Rank 0 modifies this.
3271  size_t lineNumber = 1;
3272 
3273  if (debug && myRank == rootRank) {
3274  cerr << "Matrix Market reader: readSparse:" << endl
3275  << "-- Reading banner line" << endl;
3276  }
3277 
3278  // The "Banner" tells you whether the input stream represents
3279  // a sparse matrix, the symmetry type of the matrix, and the
3280  // type of the data it contains.
3281  //
3282  // pBanner will only be nonnull on MPI Rank 0. It will be
3283  // null on all other MPI processes.
3284  RCP<const Banner> pBanner;
3285  {
3286  // We read and validate the Banner on Proc 0, but broadcast
3287  // the validation result to all processes.
3288  // Teuchos::broadcast doesn't currently work with bool, so
3289  // we use int (true -> 1, false -> 0).
3290  int bannerIsCorrect = 1;
3291  std::ostringstream errMsg;
3292 
3293  if (myRank == rootRank) {
3294  // Read the Banner line from the input stream.
3295  try {
3296  pBanner = readBanner (in, lineNumber, tolerant, debug);
3297  }
3298  catch (std::exception& e) {
3299  errMsg << "Attempt to read the Matrix Market file's Banner line "
3300  "threw an exception: " << e.what();
3301  bannerIsCorrect = 0;
3302  }
3303 
3304  if (bannerIsCorrect) {
3305  // Validate the Banner for the case of a sparse matrix.
3306  // We validate on Proc 0, since it reads the Banner.
3307 
3308  // In intolerant mode, the matrix type must be "coordinate".
3309  if (! tolerant && pBanner->matrixType() != "coordinate") {
3310  bannerIsCorrect = 0;
3311  errMsg << "The Matrix Market input file must contain a "
3312  "\"coordinate\"-format sparse matrix in order to create a "
3313  "Tpetra::CrsMatrix object from it, but the file's matrix "
3314  "type is \"" << pBanner->matrixType() << "\" instead.";
3315  }
3316  // In tolerant mode, we allow the matrix type to be
3317  // anything other than "array" (which would mean that
3318  // the file contains a dense matrix).
3319  if (tolerant && pBanner->matrixType() == "array") {
3320  bannerIsCorrect = 0;
3321  errMsg << "Matrix Market file must contain a \"coordinate\"-"
3322  "format sparse matrix in order to create a Tpetra::CrsMatrix "
3323  "object from it, but the file's matrix type is \"array\" "
3324  "instead. That probably means the file contains dense matrix "
3325  "data.";
3326  }
3327  }
3328  } // Proc 0: Done reading the Banner, hopefully successfully.
3329 
3330  // Broadcast from Proc 0 whether the Banner was read correctly.
3331  broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
3332 
3333  // If the Banner is invalid, all processes throw an
3334  // exception. Only Proc 0 gets the exception message, but
3335  // that's OK, since the main point is to "stop the world"
3336  // (rather than throw an exception on one process and leave
3337  // the others hanging).
3338  TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
3339  std::invalid_argument, errMsg.str ());
3340  } // Done reading the Banner line and broadcasting success.
3341  if (debug && myRank == rootRank) {
3342  cerr << "-- Reading dimensions line" << endl;
3343  }
3344 
3345  // Read the matrix dimensions from the Matrix Market metadata.
3346  // dims = (numRows, numCols, numEntries). Proc 0 does the
3347  // reading, but it broadcasts the results to all MPI
3348  // processes. Thus, readCoordDims() is a collective
3349  // operation. It does a collective check for correctness too.
3350  Tuple<global_ordinal_type, 3> dims =
3351  readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
3352 
3353  if (debug && myRank == rootRank) {
3354  cerr << "-- Making Adder for collecting matrix data" << endl;
3355  }
3356 
3357  // "Adder" object for collecting all the sparse matrix entries
3358  // from the input stream. This is only nonnull on Proc 0.
3359  // The Adder internally converts the one-based indices (native
3360  // Matrix Market format) into zero-based indices.
3361  RCP<adder_type> pAdder =
3362  makeAdder (pComm, pBanner, dims, tolerant, debug);
3363 
3364  if (debug && myRank == rootRank) {
3365  cerr << "-- Reading matrix data" << endl;
3366  }
3367  //
3368  // Read the matrix entries from the input stream on Proc 0.
3369  //
3370  {
3371  // We use readSuccess to broadcast the results of the read
3372  // (succeeded or not) to all MPI processes. Since
3373  // Teuchos::broadcast doesn't currently know how to send
3374  // bools, we convert to int (true -> 1, false -> 0).
3375  int readSuccess = 1;
3376  std::ostringstream errMsg; // Exception message (only valid on Proc 0)
3377  if (myRank == rootRank) {
3378  try {
3379  // Reader for "coordinate" format sparse matrix data.
3380  typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
3381  global_ordinal_type, scalar_type, STS::isComplex> reader_type;
3382  reader_type reader (pAdder);
3383 
3384  // Read the sparse matrix entries.
3385  std::pair<bool, std::vector<size_t> > results =
3386  reader.read (in, lineNumber, tolerant, debug);
3387  readSuccess = results.first ? 1 : 0;
3388  }
3389  catch (std::exception& e) {
3390  readSuccess = 0;
3391  errMsg << e.what();
3392  }
3393  }
3394  broadcast (*pComm, rootRank, ptr (&readSuccess));
3395 
3396  // It would be nice to add a "verbose" flag, so that in
3397  // tolerant mode, we could log any bad line number(s) on
3398  // Proc 0. For now, we just throw if the read fails to
3399  // succeed.
3400  //
3401  // Question: If we're in tolerant mode, and if the read did
3402  // not succeed, should we attempt to call fillComplete()?
3403  TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
3404  "Failed to read the Matrix Market sparse matrix file: "
3405  << errMsg.str());
3406  } // Done reading the matrix entries (stored on Proc 0 for now)
3407 
3408  if (debug && myRank == rootRank) {
3409  cerr << "-- Successfully read the Matrix Market data" << endl;
3410  }
3411 
3412  // In tolerant mode, we need to rebroadcast the matrix
3413  // dimensions, since they may be different after reading the
3414  // actual matrix data. We only need to broadcast the number
3415  // of rows and columns. Only Rank 0 needs to know the actual
3416  // global number of entries, since (a) we need to merge
3417  // duplicates on Rank 0 first anyway, and (b) when we
3418  // distribute the entries, each rank other than Rank 0 will
3419  // only need to know how many entries it owns, not the total
3420  // number of entries.
3421  if (tolerant) {
3422  if (debug && myRank == rootRank) {
3423  cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
3424  << endl
3425  << "----- Dimensions before: "
3426  << dims[0] << " x " << dims[1]
3427  << endl;
3428  }
3429  // Packed coordinate matrix dimensions (numRows, numCols).
3430  Tuple<global_ordinal_type, 2> updatedDims;
3431  if (myRank == rootRank) {
3432  // If one or more bottom rows of the matrix contain no
3433  // entries, then the Adder will report that the number
3434  // of rows is less than that specified in the
3435  // metadata. We allow this case, and favor the
3436  // metadata so that the zero row(s) will be included.
3437  updatedDims[0] =
3438  std::max (dims[0], pAdder->getAdder()->numRows());
3439  updatedDims[1] = pAdder->getAdder()->numCols();
3440  }
3441  broadcast (*pComm, rootRank, updatedDims);
3442  dims[0] = updatedDims[0];
3443  dims[1] = updatedDims[1];
3444  if (debug && myRank == rootRank) {
3445  cerr << "----- Dimensions after: " << dims[0] << " x "
3446  << dims[1] << endl;
3447  }
3448  }
3449  else {
3450  // In strict mode, we require that the matrix's metadata and
3451  // its actual data agree, at least somewhat. In particular,
3452  // the number of rows must agree, since otherwise we cannot
3453  // distribute the matrix correctly.
3454 
3455  // Teuchos::broadcast() doesn't know how to broadcast bools,
3456  // so we use an int with the standard 1 == true, 0 == false
3457  // encoding.
3458  int dimsMatch = 1;
3459  if (myRank == rootRank) {
3460  // If one or more bottom rows of the matrix contain no
3461  // entries, then the Adder will report that the number of
3462  // rows is less than that specified in the metadata. We
3463  // allow this case, and favor the metadata, but do not
3464  // allow the Adder to think there are more rows in the
3465  // matrix than the metadata says.
3466  if (dims[0] < pAdder->getAdder ()->numRows ()) {
3467  dimsMatch = 0;
3468  }
3469  }
3470  broadcast (*pComm, 0, ptr (&dimsMatch));
3471  if (dimsMatch == 0) {
3472  // We're in an error state anyway, so we might as well
3473  // work a little harder to print an informative error
3474  // message.
3475  //
3476  // Broadcast the Adder's idea of the matrix dimensions
3477  // from Proc 0 to all processes.
3478  Tuple<global_ordinal_type, 2> addersDims;
3479  if (myRank == rootRank) {
3480  addersDims[0] = pAdder->getAdder()->numRows();
3481  addersDims[1] = pAdder->getAdder()->numCols();
3482  }
3483  broadcast (*pComm, 0, addersDims);
3484  TEUCHOS_TEST_FOR_EXCEPTION(
3485  dimsMatch == 0, std::runtime_error,
3486  "The matrix metadata says that the matrix is " << dims[0] << " x "
3487  << dims[1] << ", but the actual data says that the matrix is "
3488  << addersDims[0] << " x " << addersDims[1] << ". That means the "
3489  "data includes more rows than reported in the metadata. This "
3490  "is not allowed when parsing in strict mode. Parse the matrix in "
3491  "tolerant mode to ignore the metadata when it disagrees with the "
3492  "data.");
3493  }
3494  } // Matrix dimensions (# rows, # cols, # entries) agree.
3495 
3496  if (debug && myRank == rootRank) {
3497  cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
3498  }
3499 
3500  // Now that we've read in all the matrix entries from the
3501  // input stream into the adder on Proc 0, post-process them
3502  // into CSR format (still on Proc 0). This will facilitate
3503  // distributing them to all the processors.
3504  //
3505  // These arrays represent the global matrix data as a CSR
3506  // matrix (with numEntriesPerRow as redundant but convenient
3507  // metadata, since it's computable from rowPtr and vice
3508  // versa). They are valid only on Proc 0.
3509  ArrayRCP<size_t> numEntriesPerRow;
3510  ArrayRCP<size_t> rowPtr;
3511  ArrayRCP<global_ordinal_type> colInd;
3512  ArrayRCP<scalar_type> values;
3513 
3514  // Proc 0 first merges duplicate entries, and then converts
3515  // the coordinate-format matrix data to CSR.
3516  {
3517  int mergeAndConvertSucceeded = 1;
3518  std::ostringstream errMsg;
3519 
3520  if (myRank == rootRank) {
3521  try {
3522  typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
3523  global_ordinal_type> element_type;
3524 
3525  // Number of rows in the matrix. If we are in tolerant
3526  // mode, we've already synchronized dims with the actual
3527  // matrix data. If in strict mode, we should use dims
3528  // (as read from the file's metadata) rather than the
3529  // matrix data to determine the dimensions. (The matrix
3530  // data will claim fewer rows than the metadata, if one
3531  // or more rows have no entries stored in the file.)
3532  const size_type numRows = dims[0];
3533 
3534  // Additively merge duplicate matrix entries.
3535  pAdder->getAdder()->merge ();
3536 
3537  // Get a temporary const view of the merged matrix entries.
3538  const std::vector<element_type>& entries =
3539  pAdder->getAdder()->getEntries();
3540 
3541  // Number of matrix entries (after merging).
3542  const size_t numEntries = (size_t)entries.size();
3543 
3544  if (debug) {
3545  cerr << "----- Proc 0: Matrix has numRows=" << numRows
3546  << " rows and numEntries=" << numEntries
3547  << " entries." << endl;
3548  }
3549 
3550  // Make space for the CSR matrix data. Converting to
3551  // CSR is easier if we fill numEntriesPerRow with zeros
3552  // at first.
3553  numEntriesPerRow = arcp<size_t> (numRows);
3554  std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
3555  rowPtr = arcp<size_t> (numRows+1);
3556  std::fill (rowPtr.begin(), rowPtr.end(), 0);
3557  colInd = arcp<global_ordinal_type> (numEntries);
3558  values = arcp<scalar_type> (numEntries);
3559 
3560  // Convert from array-of-structs coordinate format to CSR
3561  // (compressed sparse row) format.
3562  global_ordinal_type prvRow = 0;
3563  size_t curPos = 0;
3564  rowPtr[0] = 0;
3565  for (curPos = 0; curPos < numEntries; ++curPos) {
3566  const element_type& curEntry = entries[curPos];
3567  const global_ordinal_type curRow = curEntry.rowIndex();
3568  TEUCHOS_TEST_FOR_EXCEPTION(
3569  curRow < prvRow, std::logic_error,
3570  "Row indices are out of order, even though they are supposed "
3571  "to be sorted. curRow = " << curRow << ", prvRow = "
3572  << prvRow << ", at curPos = " << curPos << ". Please report "
3573  "this bug to the Tpetra developers.");
3574  if (curRow > prvRow) {
3575  prvRow = curRow;
3576  }
3577  numEntriesPerRow[curRow]++;
3578  colInd[curPos] = curEntry.colIndex();
3579  values[curPos] = curEntry.value();
3580  }
3581 
3582  rowPtr[0] = 0;
3583  for (global_ordinal_type row = 1; row <= numRows; ++row) {
3584  rowPtr[row] = numEntriesPerRow[row-1] + rowPtr[row-1];
3585  }
3586  } // Finished conversion to CSR format
3587  catch (std::exception& e) {
3588  mergeAndConvertSucceeded = 0;
3589  errMsg << "Failed to merge sparse matrix entries and convert to "
3590  "CSR format: " << e.what();
3591  }
3592 
3593  if (debug && mergeAndConvertSucceeded) {
3594  // Number of rows in the matrix.
3595  const size_type numRows = dims[0];
3596  const size_type maxToDisplay = 100;
3597 
3598  cerr << "----- Proc 0: numEntriesPerRow[0.."
3599  << (numEntriesPerRow.size()-1) << "] ";
3600  if (numRows > maxToDisplay) {
3601  cerr << "(only showing first and last few entries) ";
3602  }
3603  cerr << "= [";
3604  if (numRows > 0) {
3605  if (numRows > maxToDisplay) {
3606  for (size_type k = 0; k < 2; ++k) {
3607  cerr << numEntriesPerRow[k] << " ";
3608  }
3609  cerr << "... ";
3610  for (size_type k = numRows-2; k < numRows-1; ++k) {
3611  cerr << numEntriesPerRow[k] << " ";
3612  }
3613  }
3614  else {
3615  for (size_type k = 0; k < numRows-1; ++k) {
3616  cerr << numEntriesPerRow[k] << " ";
3617  }
3618  }
3619  cerr << numEntriesPerRow[numRows-1];
3620  } // numRows > 0
3621  cerr << "]" << endl;
3622 
3623  cerr << "----- Proc 0: rowPtr ";
3624  if (numRows > maxToDisplay) {
3625  cerr << "(only showing first and last few entries) ";
3626  }
3627  cerr << "= [";
3628  if (numRows > maxToDisplay) {
3629  for (size_type k = 0; k < 2; ++k) {
3630  cerr << rowPtr[k] << " ";
3631  }
3632  cerr << "... ";
3633  for (size_type k = numRows-2; k < numRows; ++k) {
3634  cerr << rowPtr[k] << " ";
3635  }
3636  }
3637  else {
3638  for (size_type k = 0; k < numRows; ++k) {
3639  cerr << rowPtr[k] << " ";
3640  }
3641  }
3642  cerr << rowPtr[numRows] << "]" << endl;
3643 
3644  cerr << "----- Proc 0: colInd = [";
3645  for (size_t k = 0; k < rowPtr[numRows]; ++k) {
3646  cerr << colInd[k] << " ";
3647  }
3648  cerr << "]" << endl;
3649  }
3650  } // if myRank == rootRank
3651  } // Done converting sparse matrix data to CSR format
3652 
3653  // Now we're done with the Adder, so we can release the
3654  // reference ("free" it) to save space. This only actually
3655  // does anything on Rank 0, since pAdder is null on all the
3656  // other MPI processes.
3657  pAdder = null;
3658 
3659  // Verify details of the Maps. Don't count the global number
3660  // of entries in the row Map, since that number doesn't
3661  // correctly count overlap.
3662  if (debug && myRank == rootRank) {
3663  cerr << "-- Verifying Maps" << endl;
3664  }
3665  TEUCHOS_TEST_FOR_EXCEPTION(
3666  as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
3667  std::invalid_argument,
3668  "The range Map has " << rangeMap->getGlobalNumElements ()
3669  << " entries, but the matrix has a global number of rows " << dims[0]
3670  << ".");
3671  TEUCHOS_TEST_FOR_EXCEPTION(
3672  as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
3673  std::invalid_argument,
3674  "The domain Map has " << domainMap->getGlobalNumElements ()
3675  << " entries, but the matrix has a global number of columns "
3676  << dims[1] << ".");
3677 
3678  // Create a row Map which is entirely owned on Proc 0.
3679  RCP<Teuchos::FancyOStream> err = debug ?
3680  Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
3681 
3682  RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
3683  ArrayView<const global_ordinal_type> myRows =
3684  gatherRowMap->getLocalElementList ();
3685  const size_type myNumRows = myRows.size ();
3686  const global_ordinal_type indexBase = gatherRowMap->getIndexBase ();
3687 
3688  ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
3689  for (size_type i_ = 0; i_ < myNumRows; i_++) {
3690  gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_]-indexBase];
3691  }
3692 
3693  // Create a matrix using this Map, and fill in on Proc 0. We
3694  // know how many entries there are in each row, so we can use
3695  // static profile.
3696  RCP<sparse_matrix_type> A_proc0 =
3697  rcp (new sparse_matrix_type (gatherRowMap, gatherNumEntriesPerRow()));
3698  if (myRank == rootRank) {
3699  if (debug) {
3700  cerr << "-- Proc 0: Filling gather matrix" << endl;
3701  }
3702  if (debug) {
3703  cerr << "---- Rows: " << Teuchos::toString (myRows) << endl;
3704  }
3705 
3706  // Add Proc 0's matrix entries to the CrsMatrix.
3707  for (size_type i_ = 0; i_ < myNumRows; ++i_) {
3708  size_type i = myRows[i_] - indexBase;
3709 
3710  const size_type curPos = as<size_type> (rowPtr[i]);
3711  const local_ordinal_type curNumEntries = numEntriesPerRow[i];
3712  ArrayView<global_ordinal_type> curColInd =
3713  colInd.view (curPos, curNumEntries);
3714  ArrayView<scalar_type> curValues =
3715  values.view (curPos, curNumEntries);
3716 
3717  // Modify the column indices in place to have the right index base.
3718  for (size_type k = 0; k < curNumEntries; ++k) {
3719  curColInd[k] += indexBase;
3720  }
3721  if (debug) {
3722  cerr << "------ Columns: " << Teuchos::toString (curColInd) << endl;
3723  cerr << "------ Values: " << Teuchos::toString (curValues) << endl;
3724  }
3725  // Avoid constructing empty views of ArrayRCP objects.
3726  if (curNumEntries > 0) {
3727  A_proc0->insertGlobalValues (myRows[i_], curColInd, curValues);
3728  }
3729  }
3730  // Now we can save space by deallocating numEntriesPerRow,
3731  // rowPtr, colInd, and values, since we've already put those
3732  // data in the matrix.
3733  numEntriesPerRow = null;
3734  rowPtr = null;
3735  colInd = null;
3736  values = null;
3737  } // if myRank == rootRank
3738 
3739  RCP<sparse_matrix_type> A;
3741  export_type exp (gatherRowMap, rowMap);
3742 
3743  // Communicate the precise number of nonzeros per row, which was already
3744  // calculated above.
3745  typedef local_ordinal_type LO;
3746  typedef global_ordinal_type GO;
3747  typedef Tpetra::MultiVector<GO, LO, GO, node_type> mv_type_go;
3748  mv_type_go target_nnzPerRow(rowMap,1);
3749  mv_type_go source_nnzPerRow(gatherRowMap,1);
3750  Teuchos::ArrayRCP<GO> srcData = source_nnzPerRow.getDataNonConst(0);
3751  for (int i=0; i<myNumRows; i++)
3752  srcData[i] = gatherNumEntriesPerRow[i];
3753  srcData = Teuchos::null;
3754  target_nnzPerRow.doExport(source_nnzPerRow,exp,Tpetra::INSERT);
3755  Teuchos::ArrayRCP<GO> targetData = target_nnzPerRow.getDataNonConst(0);
3756  ArrayRCP<size_t> targetData_size_t = arcp<size_t>(targetData.size());
3757  for (int i=0; i<targetData.size(); i++)
3758  targetData_size_t[i] = targetData[i];
3759 
3760  if (colMap.is_null ()) {
3761  A = rcp (new sparse_matrix_type (rowMap, targetData_size_t()));
3762  } else {
3763  A = rcp (new sparse_matrix_type (rowMap, colMap, targetData_size_t()));
3764  }
3765  A->doExport (*A_proc0, exp, INSERT);
3766  if (callFillComplete) {
3767  A->fillComplete (domainMap, rangeMap);
3768  }
3769 
3770  // We can't get the dimensions of the matrix until after
3771  // fillComplete() is called. Thus, we can't do the sanity
3772  // check (dimensions read from the Matrix Market data,
3773  // vs. dimensions reported by the CrsMatrix) unless the user
3774  // asked us to call fillComplete().
3775  //
3776  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
3777  // what one might think it does, so you have to ask the range
3778  // resp. domain map for the number of rows resp. columns.
3779  if (callFillComplete) {
3780  const int numProcs = pComm->getSize ();
3781 
3782  if (extraDebug && debug) {
3783  const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
3784  const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
3785  if (myRank == rootRank) {
3786  cerr << "-- Matrix is "
3787  << globalNumRows << " x " << globalNumCols
3788  << " with " << A->getGlobalNumEntries()
3789  << " entries, and index base "
3790  << A->getIndexBase() << "." << endl;
3791  }
3792  pComm->barrier ();
3793  for (int p = 0; p < numProcs; ++p) {
3794  if (myRank == p) {
3795  cerr << "-- Proc " << p << " owns "
3796  << A->getLocalNumCols() << " columns, and "
3797  << A->getLocalNumEntries() << " entries." << endl;
3798  }
3799  pComm->barrier ();
3800  }
3801  } // if (extraDebug && debug)
3802  } // if (callFillComplete)
3803 
3804  if (debug && myRank == rootRank) {
3805  cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
3806  << endl;
3807  }
3808  return A;
3809  }
3810 
3840  static Teuchos::RCP<multivector_type>
3841  readDenseFile (const std::string& filename,
3842  const trcp_tcomm_t& comm,
3843  Teuchos::RCP<const map_type>& map,
3844  const bool tolerant=false,
3845  const bool debug=false,
3846  const bool binary=false)
3847  {
3848  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true, binary ? std::ios::binary : std::ios::in);
3849 
3850  return readDense (in, comm, map, tolerant, debug, binary);
3851  }
3852 
3863  static std::ifstream openInFileOnRankZero(
3864  const trcp_tcomm_t& comm,
3865  const std::string& filename, const bool safe = true,
3866  std::ios_base::openmode mode = std::ios_base::in
3867  ){
3868  // Input stream.
3869  std::ifstream in;
3870 
3871  // Placeholder for broadcasting in-stream state.
3872  int all_should_stop = false;
3873 
3874  // Try to open the in-stream on root rank.
3875  if (comm->getRank() == 0) {
3876  in.open(filename, mode);
3877  all_should_stop = !in && safe;
3878  }
3879 
3880  // Broadcast state and possibly throw.
3881  if(comm) Teuchos::broadcast(*comm, 0, &all_should_stop);
3882 
3883  TEUCHOS_TEST_FOR_EXCEPTION(
3884  all_should_stop,
3885  std::runtime_error,
3886  "Could not open input file '" + filename + "' on root rank 0."
3887  );
3888 
3889  return in;
3890  }
3891 
3892 
3922  static Teuchos::RCP<vector_type>
3923  readVectorFile (const std::string& filename,
3924  const trcp_tcomm_t& comm,
3925  Teuchos::RCP<const map_type>& map,
3926  const bool tolerant=false,
3927  const bool debug=false)
3928  {
3929  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true);
3930 
3931  return readVector (in, comm, map, tolerant, debug);
3932  }
3933 
3934 
4002 
4003  static Teuchos::RCP<multivector_type>
4004  readDense (std::istream& in,
4005  const trcp_tcomm_t& comm,
4006  Teuchos::RCP<const map_type>& map,
4007  const bool tolerant=false,
4008  const bool debug=false,
4009  const bool binary=false)
4010  {
4011  Teuchos::RCP<Teuchos::FancyOStream> err =
4012  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4013  return readDenseImpl<scalar_type> (in, comm, map, err, tolerant, debug, binary);
4014  }
4015 
4016 
4018  static Teuchos::RCP<vector_type>
4019  readVector (std::istream& in,
4020  const trcp_tcomm_t& comm,
4021  Teuchos::RCP<const map_type>& map,
4022  const bool tolerant=false,
4023  const bool debug=false)
4024  {
4025  Teuchos::RCP<Teuchos::FancyOStream> err =
4026  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4027  return readVectorImpl<scalar_type> (in, comm, map, err, tolerant, debug);
4028  }
4029 
4030 
4052  static Teuchos::RCP<const map_type>
4053  readMapFile (const std::string& filename,
4054  const trcp_tcomm_t& comm,
4055  const bool tolerant=false,
4056  const bool debug=false,
4057  const bool binary=false)
4058  {
4059  std::ifstream in = Reader::openInFileOnRankZero(comm, filename, true, binary ? std::ios::binary : std::ios::in);
4060 
4061  return readMap (in, comm, tolerant, debug, binary);
4062  }
4063 
4064 
4065  private:
4066  template<class MultiVectorScalarType>
4067  static Teuchos::RCP<Tpetra::MultiVector<MultiVectorScalarType,
4070  node_type> >
4071  readDenseImpl (std::istream& in,
4072  const trcp_tcomm_t& comm,
4073  Teuchos::RCP<const map_type>& map,
4074  const Teuchos::RCP<Teuchos::FancyOStream>& err,
4075  const bool tolerant=false,
4076  const bool debug=false,
4077  const bool binary=false)
4078  {
4079  using Teuchos::MatrixMarket::Banner;
4080  using Teuchos::MatrixMarket::checkCommentLine;
4081  using Teuchos::ArrayRCP;
4082  using Teuchos::as;
4083  using Teuchos::broadcast;
4084  using Teuchos::outArg;
4085  using Teuchos::RCP;
4086  using Teuchos::Tuple;
4087  using std::endl;
4088  typedef MultiVectorScalarType ST;
4089  typedef local_ordinal_type LO;
4090  typedef global_ordinal_type GO;
4091  typedef node_type NT;
4092  typedef Teuchos::ScalarTraits<ST> STS;
4093  typedef typename STS::magnitudeType MT;
4094  typedef Teuchos::ScalarTraits<MT> STM;
4096 
4097  // Rank 0 is the only (MPI) process allowed to read from the
4098  // input stream.
4099  const int myRank = comm->getRank ();
4100 
4101  if (! err.is_null ()) {
4102  err->pushTab ();
4103  }
4104  if (debug) {
4105  *err << myRank << ": readDenseImpl" << endl;
4106  }
4107  if (! err.is_null ()) {
4108  err->pushTab ();
4109  }
4110 
4111  // mfh 17 Feb 2013: It's not strictly necessary that the Comm
4112  // instances be identical and that the Node instances be
4113  // identical. The essential condition is more complicated to
4114  // test and isn't the same for all Node types. Thus, we just
4115  // leave it up to the user.
4116 
4117  // // If map is nonnull, check the precondition that its
4118  // // communicator resp. node equal comm resp. node. Checking
4119  // // now avoids doing a lot of file reading before we detect the
4120  // // violated precondition.
4121  // TEUCHOS_TEST_FOR_EXCEPTION(
4122  // ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
4123  // std::invalid_argument, "If you supply a nonnull Map, the Map's "
4124  // "communicator and node must equal the supplied communicator resp. "
4125  // "node.");
4126 
4127  // Process 0 will read in the matrix dimensions from the file,
4128  // and broadcast them to all ranks in the given communicator.
4129  // There are only 2 dimensions in the matrix, but we use the
4130  // third element of the Tuple to encode the banner's reported
4131  // data type: "real" == 0, "complex" == 1, and "integer" == 0
4132  // (same as "real"). We don't allow pattern matrices (i.e.,
4133  // graphs) since they only make sense for sparse data.
4134  Tuple<GO, 3> dims;
4135  dims[0] = 0;
4136  dims[1] = 0;
4137  dims[2] = 0;
4138 
4139  // Current line number in the input stream. Only valid on
4140  // Proc 0. Various calls will modify this depending on the
4141  // number of lines that are read from the input stream.
4142  size_t lineNumber = 1;
4143 
4144  // Capture errors and their messages on Proc 0.
4145  std::ostringstream exMsg;
4146  int localBannerReadSuccess = 1;
4147  int localDimsReadSuccess = 1;
4148 
4149  // Only Proc 0 gets to read matrix data from the input stream.
4150  if (myRank == 0) {
4151  if (! binary) {
4152  if (debug) {
4153  *err << myRank << ": readDenseImpl: Reading banner line (dense)" << endl;
4154  }
4155 
4156  // The "Banner" tells you whether the input stream
4157  // represents a dense matrix, the symmetry type of the
4158  // matrix, and the type of the data it contains.
4159  RCP<const Banner> pBanner;
4160  try {
4161  pBanner = readBanner (in, lineNumber, tolerant, debug);
4162  } catch (std::exception& e) {
4163  exMsg << e.what ();
4164  localBannerReadSuccess = 0;
4165  }
4166  // Make sure the input stream is the right kind of data.
4167  if (localBannerReadSuccess) {
4168  if (pBanner->matrixType () != "array") {
4169  exMsg << "The Matrix Market file does not contain dense matrix "
4170  "data. Its banner (first) line says that its matrix type is \""
4171  << pBanner->matrixType () << "\", rather that the required "
4172  "\"array\".";
4173  localBannerReadSuccess = 0;
4174  } else if (pBanner->dataType() == "pattern") {
4175  exMsg << "The Matrix Market file's banner (first) "
4176  "line claims that the matrix's data type is \"pattern\". This does "
4177  "not make sense for a dense matrix, yet the file reports the matrix "
4178  "as dense. The only valid data types for a dense matrix are "
4179  "\"real\", \"complex\", and \"integer\".";
4180  localBannerReadSuccess = 0;
4181  } else {
4182  // Encode the data type reported by the Banner as the
4183  // third element of the dimensions Tuple.
4184  dims[2] = encodeDataType (pBanner->dataType ());
4185  }
4186  } // if we successfully read the banner line
4187 
4188  // At this point, we've successfully read the banner line.
4189  // Now read the dimensions line.
4190  if (localBannerReadSuccess) {
4191  if (debug) {
4192  *err << myRank << ": readDenseImpl: Reading dimensions line (dense)" << endl;
4193  }
4194  // Keep reading lines from the input stream until we find
4195  // a non-comment line, or until we run out of lines. The
4196  // latter is an error, since every "array" format Matrix
4197  // Market file must have a dimensions line after the
4198  // banner (even if the matrix has zero rows or columns, or
4199  // zero entries).
4200  std::string line;
4201  bool commentLine = true;
4202 
4203  while (commentLine) {
4204  // Test whether it is even valid to read from the input
4205  // stream wrapping the line.
4206  if (in.eof () || in.fail ()) {
4207  exMsg << "Unable to get array dimensions line (at all) from line "
4208  << lineNumber << " of input stream. The input stream "
4209  << "claims that it is "
4210  << (in.eof() ? "at end-of-file." : "in a failed state.");
4211  localDimsReadSuccess = 0;
4212  } else {
4213  // Try to get the next line from the input stream.
4214  if (getline (in, line)) {
4215  ++lineNumber; // We did actually read a line.
4216  }
4217  // Is the current line a comment line? Ignore start
4218  // and size; they are only useful for reading the
4219  // actual matrix entries. (We could use them here as
4220  // an optimization, but we've chosen not to.)
4221  size_t start = 0, size = 0;
4222  commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4223  } // whether we failed to read the line at all
4224  } // while the line we just read is a comment line
4225 
4226  //
4227  // Get <numRows> <numCols> from the line we just read.
4228  //
4229  std::istringstream istr (line);
4230 
4231  // Test whether it is even valid to read from the input
4232  // stream wrapping the line.
4233  if (istr.eof () || istr.fail ()) {
4234  exMsg << "Unable to read any data from line " << lineNumber
4235  << " of input; the line should contain the matrix dimensions "
4236  << "\"<numRows> <numCols>\".";
4237  localDimsReadSuccess = 0;
4238  } else { // It's valid to read from the line.
4239  GO theNumRows = 0;
4240  istr >> theNumRows; // Read in the number of rows.
4241  if (istr.fail ()) {
4242  exMsg << "Failed to get number of rows from line "
4243  << lineNumber << " of input; the line should contains the "
4244  << "matrix dimensions \"<numRows> <numCols>\".";
4245  localDimsReadSuccess = 0;
4246  } else { // We successfully read the number of rows
4247  dims[0] = theNumRows; // Save the number of rows
4248  if (istr.eof ()) { // Do we still have data to read?
4249  exMsg << "No more data after number of rows on line "
4250  << lineNumber << " of input; the line should contain the "
4251  << "matrix dimensions \"<numRows> <numCols>\".";
4252  localDimsReadSuccess = 0;
4253  } else { // Still data left to read; read in number of columns.
4254  GO theNumCols = 0;
4255  istr >> theNumCols; // Read in the number of columns
4256  if (istr.fail ()) {
4257  exMsg << "Failed to get number of columns from line "
4258  << lineNumber << " of input; the line should contain "
4259  << "the matrix dimensions \"<numRows> <numCols>\".";
4260  localDimsReadSuccess = 0;
4261  } else { // We successfully read the number of columns
4262  dims[1] = theNumCols; // Save the number of columns
4263  } // if istr.fail ()
4264  } // if istr.eof ()
4265  } // if we read the number of rows
4266  } // if the input stream wrapping the dims line was (in)valid
4267  } // if we successfully read the banner line
4268  } // not in binary format
4269  else { // in binary format
4270  global_size_t numRows, numCols;
4271  in.read(reinterpret_cast<char*>(&numRows), sizeof(numRows));
4272  in.read(reinterpret_cast<char*>(&numCols), sizeof(numCols));
4273  dims[0] = Teuchos::as<GO>(numRows);
4274  dims[1] = Teuchos::as<GO>(numCols);
4275  if ((typeid(ST) == typeid(double)) || Teuchos::ScalarTraits<ST>::isOrdinal) {
4276  dims[2] = 0;
4277  } else if (Teuchos::ScalarTraits<ST>::isComplex) {
4278  dims[2] = 1;
4279  } else {
4280  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
4281  "Unrecognized Matrix Market data type. "
4282  "We should never get here. "
4283  "Please report this bug to the Tpetra developers.");
4284  }
4285  localDimsReadSuccess = true;
4286  localDimsReadSuccess = true;
4287  } // in binary format
4288  } // if (myRank == 0)
4289 
4290  // Broadcast the matrix dimensions, the encoded data type, and
4291  // whether or not Proc 0 succeeded in reading the banner and
4292  // dimensions.
4293  Tuple<GO, 5> bannerDimsReadResult;
4294  if (myRank == 0) {
4295  bannerDimsReadResult[0] = dims[0]; // numRows
4296  bannerDimsReadResult[1] = dims[1]; // numCols
4297  bannerDimsReadResult[2] = dims[2]; // encoded data type
4298  bannerDimsReadResult[3] = localBannerReadSuccess;
4299  bannerDimsReadResult[4] = localDimsReadSuccess;
4300  }
4301  // Broadcast matrix dimensions and the encoded data type from
4302  // Proc 0 to all the MPI processes.
4303  broadcast (*comm, 0, bannerDimsReadResult);
4304 
4305  TEUCHOS_TEST_FOR_EXCEPTION(
4306  bannerDimsReadResult[3] == 0, std::runtime_error,
4307  "Failed to read banner line: " << exMsg.str ());
4308  TEUCHOS_TEST_FOR_EXCEPTION(
4309  bannerDimsReadResult[4] == 0, std::runtime_error,
4310  "Failed to read matrix dimensions line: " << exMsg.str ());
4311  if (myRank != 0) {
4312  dims[0] = bannerDimsReadResult[0];
4313  dims[1] = bannerDimsReadResult[1];
4314  dims[2] = bannerDimsReadResult[2];
4315  }
4316 
4317  // Tpetra objects want the matrix dimensions in these types.
4318  const global_size_t numRows = static_cast<global_size_t> (dims[0]);
4319  const size_t numCols = static_cast<size_t> (dims[1]);
4320 
4321  // Make a "Proc 0 owns everything" Map that we will use to
4322  // read in the multivector entries in the correct order on
4323  // Proc 0. This must be a collective
4324  RCP<const map_type> proc0Map; // "Proc 0 owns everything" Map
4325  if (map.is_null ()) {
4326  // The user didn't supply a Map. Make a contiguous
4327  // distributed Map for them, using the read-in multivector
4328  // dimensions.
4329  map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4330  const size_t localNumRows = (myRank == 0) ? numRows : 0;
4331  // At this point, map exists and has a nonnull node.
4332  proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4333  comm);
4334  }
4335  else { // The user supplied a Map.
4336  proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4337  }
4338 
4339  // Make a multivector X owned entirely by Proc 0.
4340  RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
4341 
4342  //
4343  // On Proc 0, read the Matrix Market data from the input
4344  // stream into the multivector X.
4345  //
4346  int localReadDataSuccess = 1;
4347  if (myRank == 0) {
4348  if (! binary) {
4349  try {
4350  if (debug) {
4351  *err << myRank << ": readDenseImpl: Reading matrix data (dense)"
4352  << endl;
4353  }
4354 
4355  // Make sure that we can get a 1-D view of X.
4356  TEUCHOS_TEST_FOR_EXCEPTION(
4357  ! X->isConstantStride (), std::logic_error,
4358  "Can't get a 1-D view of the entries of the MultiVector X on "
4359  "Process 0, because the stride between the columns of X is not "
4360  "constant. This shouldn't happen because we just created X and "
4361  "haven't filled it in yet. Please report this bug to the Tpetra "
4362  "developers.");
4363 
4364  // Get a writeable 1-D view of the entries of X. Rank 0
4365  // owns all of them. The view will expire at the end of
4366  // scope, so (if necessary) it will be written back to X
4367  // at this time.
4368  ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4369  TEUCHOS_TEST_FOR_EXCEPTION(
4370  as<global_size_t> (X_view.size ()) < numRows * numCols,
4371  std::logic_error,
4372  "The view of X has size " << X_view.size() << " which is not enough to "
4373  "accommodate the expected number of entries numRows*numCols = "
4374  << numRows << "*" << numCols << " = " << numRows*numCols << ". "
4375  "Please report this bug to the Tpetra developers.");
4376  const size_t stride = X->getStride ();
4377 
4378  // The third element of the dimensions Tuple encodes the data
4379  // type reported by the Banner: "real" == 0, "complex" == 1,
4380  // "integer" == 0 (same as "real"), "pattern" == 2. We do not
4381  // allow dense matrices to be pattern matrices, so dims[2] ==
4382  // 0 or 1. We've already checked for this above.
4383  const bool isComplex = (dims[2] == 1);
4384  size_type count = 0, curRow = 0, curCol = 0;
4385 
4386  std::string line;
4387  while (getline (in, line)) {
4388  ++lineNumber;
4389  // Is the current line a comment line? If it's not,
4390  // line.substr(start,size) contains the data.
4391  size_t start = 0, size = 0;
4392  const bool commentLine =
4393  checkCommentLine (line, start, size, lineNumber, tolerant);
4394  if (! commentLine) {
4395  // Make sure we have room in which to put the new matrix
4396  // entry. We check this only after checking for a
4397  // comment line, because there may be one or more
4398  // comment lines at the end of the file. In tolerant
4399  // mode, we simply ignore any extra data.
4400  if (count >= X_view.size()) {
4401  if (tolerant) {
4402  break;
4403  }
4404  else {
4405  TEUCHOS_TEST_FOR_EXCEPTION(
4406  count >= X_view.size(),
4407  std::runtime_error,
4408  "The Matrix Market input stream has more data in it than "
4409  "its metadata reported. Current line number is "
4410  << lineNumber << ".");
4411  }
4412  }
4413 
4414  // mfh 19 Dec 2012: Ignore everything up to the initial
4415  // colon. writeDense() has the option to print out the
4416  // global row index in front of each entry, followed by
4417  // a colon and space.
4418  {
4419  const size_t pos = line.substr (start, size).find (':');
4420  if (pos != std::string::npos) {
4421  start = pos+1;
4422  }
4423  }
4424  std::istringstream istr (line.substr (start, size));
4425  // Does the line contain anything at all? Can we
4426  // safely read from the input stream wrapping the
4427  // line?
4428  if (istr.eof() || istr.fail()) {
4429  // In tolerant mode, simply ignore the line.
4430  if (tolerant) {
4431  break;
4432  }
4433  // We repeat the full test here so the exception
4434  // message is more informative.
4435  TEUCHOS_TEST_FOR_EXCEPTION(
4436  ! tolerant && (istr.eof() || istr.fail()),
4437  std::runtime_error,
4438  "Line " << lineNumber << " of the Matrix Market file is "
4439  "empty, or we cannot read from it for some other reason.");
4440  }
4441  // Current matrix entry to read in.
4442  ST val = STS::zero();
4443  // Real and imaginary parts of the current matrix entry.
4444  // The imaginary part is zero if the matrix is real-valued.
4445  MT real = STM::zero(), imag = STM::zero();
4446 
4447  // isComplex refers to the input stream's data, not to
4448  // the scalar type S. It's OK to read real-valued
4449  // data into a matrix storing complex-valued data; in
4450  // that case, all entries' imaginary parts are zero.
4451  if (isComplex) {
4452  // STS::real() and STS::imag() return a copy of
4453  // their respective components, not a writeable
4454  // reference. Otherwise we could just assign to
4455  // them using the istream extraction operator (>>).
4456  // That's why we have separate magnitude type "real"
4457  // and "imag" variables.
4458 
4459  // Attempt to read the real part of the current entry.
4460  istr >> real;
4461  if (istr.fail()) {
4462  TEUCHOS_TEST_FOR_EXCEPTION(
4463  ! tolerant && istr.eof(), std::runtime_error,
4464  "Failed to get the real part of a complex-valued matrix "
4465  "entry from line " << lineNumber << " of the Matrix Market "
4466  "file.");
4467  // In tolerant mode, just skip bad lines.
4468  if (tolerant) {
4469  break;
4470  }
4471  } else if (istr.eof()) {
4472  TEUCHOS_TEST_FOR_EXCEPTION(
4473  ! tolerant && istr.eof(), std::runtime_error,
4474  "Missing imaginary part of a complex-valued matrix entry "
4475  "on line " << lineNumber << " of the Matrix Market file.");
4476  // In tolerant mode, let any missing imaginary part be 0.
4477  } else {
4478  // Attempt to read the imaginary part of the current
4479  // matrix entry.
4480  istr >> imag;
4481  TEUCHOS_TEST_FOR_EXCEPTION(
4482  ! tolerant && istr.fail(), std::runtime_error,
4483  "Failed to get the imaginary part of a complex-valued "
4484  "matrix entry from line " << lineNumber << " of the "
4485  "Matrix Market file.");
4486  // In tolerant mode, let any missing or corrupted
4487  // imaginary part be 0.
4488  }
4489  } else { // Matrix Market file contains real-valued data.
4490  // Attempt to read the current matrix entry.
4491  istr >> real;
4492  TEUCHOS_TEST_FOR_EXCEPTION(
4493  ! tolerant && istr.fail(), std::runtime_error,
4494  "Failed to get a real-valued matrix entry from line "
4495  << lineNumber << " of the Matrix Market file.");
4496  // In tolerant mode, simply ignore the line if
4497  // we failed to read a matrix entry.
4498  if (istr.fail() && tolerant) {
4499  break;
4500  }
4501  }
4502  // In tolerant mode, we simply let pass through whatever
4503  // data we got.
4504  TEUCHOS_TEST_FOR_EXCEPTION(
4505  ! tolerant && istr.fail(), std::runtime_error,
4506  "Failed to read matrix data from line " << lineNumber
4507  << " of the Matrix Market file.");
4508 
4509  // Assign val = ST(real, imag).
4510  Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4511 
4512  curRow = count % numRows;
4513  curCol = count / numRows;
4514  X_view[curRow + curCol*stride] = val;
4515  ++count;
4516  } // if not a comment line
4517  } // while there are still lines in the file, get the next one
4518 
4519  TEUCHOS_TEST_FOR_EXCEPTION(
4520  ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
4521  std::runtime_error,
4522  "The Matrix Market metadata reports that the dense matrix is "
4523  << numRows << " x " << numCols << ", and thus has "
4524  << numRows*numCols << " total entries, but we only found " << count
4525  << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
4526  } catch (std::exception& e) {
4527  exMsg << e.what ();
4528  localReadDataSuccess = 0;
4529  }
4530  } // not binary file
4531  else {
4532  if (debug) {
4533  *err << myRank << ": readDenseImpl: Reading matrix data (dense)"
4534  << endl;
4535  }
4536 
4537  // Make sure that we can get a 1-D view of X.
4538  TEUCHOS_TEST_FOR_EXCEPTION(
4539  ! X->isConstantStride (), std::logic_error,
4540  "Can't get a 1-D view of the entries of the MultiVector X on "
4541  "Process 0, because the stride between the columns of X is not "
4542  "constant. This shouldn't happen because we just created X and "
4543  "haven't filled it in yet. Please report this bug to the Tpetra "
4544  "developers.");
4545 
4546  // Get a writeable 1-D view of the entries of X. Rank 0
4547  // owns all of them. The view will expire at the end of
4548  // scope, so (if necessary) it will be written back to X
4549  // at this time.
4550  auto X_view = X->getLocalViewHost (Access::OverwriteAll);
4551 
4552  TEUCHOS_TEST_FOR_EXCEPTION(
4553  as<global_size_t> (X_view.extent(0)) < numRows,
4554  std::logic_error,
4555  "The view of X has " << X_view.extent(0) << " rows which is not enough to "
4556  "accommodate the expected number of entries numRows = "
4557  << numRows << ". "
4558  "Please report this bug to the Tpetra developers.");
4559  TEUCHOS_TEST_FOR_EXCEPTION(
4560  as<global_size_t> (X_view.extent(1)) < numCols,
4561  std::logic_error,
4562  "The view of X has " << X_view.extent(1) << " colums which is not enough to "
4563  "accommodate the expected number of entries numRows = "
4564  << numCols << ". "
4565  "Please report this bug to the Tpetra developers.");
4566 
4567  for (size_t curRow = 0; curRow < numRows; ++curRow) {
4568  for (size_t curCol = 0; curCol < numCols; ++curCol) {
4569  if (Teuchos::ScalarTraits<ST>::isOrdinal){
4570  global_size_t val;
4571  in.read(reinterpret_cast<char*>(&val), sizeof(val));
4572  X_view(curRow, curCol) = val;
4573  } else {
4574  double val;
4575  in.read(reinterpret_cast<char*>(&val), sizeof(val));
4576  X_view(curRow, curCol) = val;
4577  }
4578  }
4579  }
4580  } // binary
4581  } // if (myRank == 0)
4582 
4583  if (debug) {
4584  *err << myRank << ": readDenseImpl: done reading data" << endl;
4585  }
4586 
4587  // Synchronize on whether Proc 0 successfully read the data.
4588  int globalReadDataSuccess = localReadDataSuccess;
4589  broadcast (*comm, 0, outArg (globalReadDataSuccess));
4590  TEUCHOS_TEST_FOR_EXCEPTION(
4591  globalReadDataSuccess == 0, std::runtime_error,
4592  "Failed to read the multivector's data: " << exMsg.str ());
4593 
4594  // If there's only one MPI process and the user didn't supply
4595  // a Map (i.e., pMap is null), we're done. Set pMap to the
4596  // Map used to distribute X, and return X.
4597  if (comm->getSize () == 1 && map.is_null ()) {
4598  map = proc0Map;
4599  if (! err.is_null ()) {
4600  err->popTab ();
4601  }
4602  if (debug) {
4603  *err << myRank << ": readDenseImpl: done" << endl;
4604  }
4605  if (! err.is_null ()) {
4606  err->popTab ();
4607  }
4608  return X;
4609  }
4610 
4611  if (debug) {
4612  *err << myRank << ": readDenseImpl: Creating target MV" << endl;
4613  }
4614 
4615  // Make a multivector Y with the distributed map pMap.
4616  RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
4617 
4618  if (debug) {
4619  *err << myRank << ": readDenseImpl: Creating Export" << endl;
4620  }
4621 
4622  // Make an Export object that will export X to Y. First
4623  // argument is the source map, second argument is the target
4624  // map.
4625  Export<LO, GO, NT> exporter (proc0Map, map, err);
4626 
4627  if (debug) {
4628  *err << myRank << ": readDenseImpl: Exporting" << endl;
4629  }
4630  // Export X into Y.
4631  Y->doExport (*X, exporter, INSERT);
4632 
4633  if (! err.is_null ()) {
4634  err->popTab ();
4635  }
4636  if (debug) {
4637  *err << myRank << ": readDenseImpl: done" << endl;
4638  }
4639  if (! err.is_null ()) {
4640  err->popTab ();
4641  }
4642 
4643  // Y is distributed over all process(es) in the communicator.
4644  return Y;
4645  }
4646 
4647 
4648  template<class VectorScalarType>
4649  static Teuchos::RCP<Tpetra::Vector<VectorScalarType,
4652  node_type> >
4653  readVectorImpl (std::istream& in,
4654  const trcp_tcomm_t& comm,
4655  Teuchos::RCP<const map_type>& map,
4656  const Teuchos::RCP<Teuchos::FancyOStream>& err,
4657  const bool tolerant=false,
4658  const bool debug=false)
4659  {
4660  using Teuchos::MatrixMarket::Banner;
4661  using Teuchos::MatrixMarket::checkCommentLine;
4662  using Teuchos::as;
4663  using Teuchos::broadcast;
4664  using Teuchos::outArg;
4665  using Teuchos::RCP;
4666  using Teuchos::Tuple;
4667  using std::endl;
4668  typedef VectorScalarType ST;
4669  typedef local_ordinal_type LO;
4670  typedef global_ordinal_type GO;
4671  typedef node_type NT;
4672  typedef Teuchos::ScalarTraits<ST> STS;
4673  typedef typename STS::magnitudeType MT;
4674  typedef Teuchos::ScalarTraits<MT> STM;
4675  typedef Tpetra::Vector<ST, LO, GO, NT> MV;
4676 
4677  // Rank 0 is the only (MPI) process allowed to read from the
4678  // input stream.
4679  const int myRank = comm->getRank ();
4680 
4681  if (! err.is_null ()) {
4682  err->pushTab ();
4683  }
4684  if (debug) {
4685  *err << myRank << ": readVectorImpl" << endl;
4686  }
4687  if (! err.is_null ()) {
4688  err->pushTab ();
4689  }
4690 
4691  // mfh 17 Feb 2013: It's not strictly necessary that the Comm
4692  // instances be identical and that the Node instances be
4693  // identical. The essential condition is more complicated to
4694  // test and isn't the same for all Node types. Thus, we just
4695  // leave it up to the user.
4696 
4697  // // If map is nonnull, check the precondition that its
4698  // // communicator resp. node equal comm resp. node. Checking
4699  // // now avoids doing a lot of file reading before we detect the
4700  // // violated precondition.
4701  // TEUCHOS_TEST_FOR_EXCEPTION(
4702  // ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
4703  // std::invalid_argument, "If you supply a nonnull Map, the Map's "
4704  // "communicator and node must equal the supplied communicator resp. "
4705  // "node.");
4706 
4707  // Process 0 will read in the matrix dimensions from the file,
4708  // and broadcast them to all ranks in the given communicator.
4709  // There are only 2 dimensions in the matrix, but we use the
4710  // third element of the Tuple to encode the banner's reported
4711  // data type: "real" == 0, "complex" == 1, and "integer" == 0
4712  // (same as "real"). We don't allow pattern matrices (i.e.,
4713  // graphs) since they only make sense for sparse data.
4714  Tuple<GO, 3> dims;
4715  dims[0] = 0;
4716  dims[1] = 0;
4717 
4718  // Current line number in the input stream. Only valid on
4719  // Proc 0. Various calls will modify this depending on the
4720  // number of lines that are read from the input stream.
4721  size_t lineNumber = 1;
4722 
4723  // Capture errors and their messages on Proc 0.
4724  std::ostringstream exMsg;
4725  int localBannerReadSuccess = 1;
4726  int localDimsReadSuccess = 1;
4727 
4728  // Only Proc 0 gets to read matrix data from the input stream.
4729  if (myRank == 0) {
4730  if (debug) {
4731  *err << myRank << ": readVectorImpl: Reading banner line (dense)" << endl;
4732  }
4733 
4734  // The "Banner" tells you whether the input stream
4735  // represents a dense matrix, the symmetry type of the
4736  // matrix, and the type of the data it contains.
4737  RCP<const Banner> pBanner;
4738  try {
4739  pBanner = readBanner (in, lineNumber, tolerant, debug);
4740  } catch (std::exception& e) {
4741  exMsg << e.what ();
4742  localBannerReadSuccess = 0;
4743  }
4744  // Make sure the input stream is the right kind of data.
4745  if (localBannerReadSuccess) {
4746  if (pBanner->matrixType () != "array") {
4747  exMsg << "The Matrix Market file does not contain dense matrix "
4748  "data. Its banner (first) line says that its matrix type is \""
4749  << pBanner->matrixType () << "\", rather that the required "
4750  "\"array\".";
4751  localBannerReadSuccess = 0;
4752  } else if (pBanner->dataType() == "pattern") {
4753  exMsg << "The Matrix Market file's banner (first) "
4754  "line claims that the matrix's data type is \"pattern\". This does "
4755  "not make sense for a dense matrix, yet the file reports the matrix "
4756  "as dense. The only valid data types for a dense matrix are "
4757  "\"real\", \"complex\", and \"integer\".";
4758  localBannerReadSuccess = 0;
4759  } else {
4760  // Encode the data type reported by the Banner as the
4761  // third element of the dimensions Tuple.
4762  dims[2] = encodeDataType (pBanner->dataType ());
4763  }
4764  } // if we successfully read the banner line
4765 
4766  // At this point, we've successfully read the banner line.
4767  // Now read the dimensions line.
4768  if (localBannerReadSuccess) {
4769  if (debug) {
4770  *err << myRank << ": readVectorImpl: Reading dimensions line (dense)" << endl;
4771  }
4772  // Keep reading lines from the input stream until we find
4773  // a non-comment line, or until we run out of lines. The
4774  // latter is an error, since every "array" format Matrix
4775  // Market file must have a dimensions line after the
4776  // banner (even if the matrix has zero rows or columns, or
4777  // zero entries).
4778  std::string line;
4779  bool commentLine = true;
4780 
4781  while (commentLine) {
4782  // Test whether it is even valid to read from the input
4783  // stream wrapping the line.
4784  if (in.eof () || in.fail ()) {
4785  exMsg << "Unable to get array dimensions line (at all) from line "
4786  << lineNumber << " of input stream. The input stream "
4787  << "claims that it is "
4788  << (in.eof() ? "at end-of-file." : "in a failed state.");
4789  localDimsReadSuccess = 0;
4790  } else {
4791  // Try to get the next line from the input stream.
4792  if (getline (in, line)) {
4793  ++lineNumber; // We did actually read a line.
4794  }
4795  // Is the current line a comment line? Ignore start
4796  // and size; they are only useful for reading the
4797  // actual matrix entries. (We could use them here as
4798  // an optimization, but we've chosen not to.)
4799  size_t start = 0, size = 0;
4800  commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4801  } // whether we failed to read the line at all
4802  } // while the line we just read is a comment line
4803 
4804  //
4805  // Get <numRows> <numCols> from the line we just read.
4806  //
4807  std::istringstream istr (line);
4808 
4809  // Test whether it is even valid to read from the input
4810  // stream wrapping the line.
4811  if (istr.eof () || istr.fail ()) {
4812  exMsg << "Unable to read any data from line " << lineNumber
4813  << " of input; the line should contain the matrix dimensions "
4814  << "\"<numRows> <numCols>\".";
4815  localDimsReadSuccess = 0;
4816  } else { // It's valid to read from the line.
4817  GO theNumRows = 0;
4818  istr >> theNumRows; // Read in the number of rows.
4819  if (istr.fail ()) {
4820  exMsg << "Failed to get number of rows from line "
4821  << lineNumber << " of input; the line should contains the "
4822  << "matrix dimensions \"<numRows> <numCols>\".";
4823  localDimsReadSuccess = 0;
4824  } else { // We successfully read the number of rows
4825  dims[0] = theNumRows; // Save the number of rows
4826  if (istr.eof ()) { // Do we still have data to read?
4827  exMsg << "No more data after number of rows on line "
4828  << lineNumber << " of input; the line should contain the "
4829  << "matrix dimensions \"<numRows> <numCols>\".";
4830  localDimsReadSuccess = 0;
4831  } else { // Still data left to read; read in number of columns.
4832  GO theNumCols = 0;
4833  istr >> theNumCols; // Read in the number of columns
4834  if (istr.fail ()) {
4835  exMsg << "Failed to get number of columns from line "
4836  << lineNumber << " of input; the line should contain "
4837  << "the matrix dimensions \"<numRows> <numCols>\".";
4838  localDimsReadSuccess = 0;
4839  } else { // We successfully read the number of columns
4840  dims[1] = theNumCols; // Save the number of columns
4841  } // if istr.fail ()
4842  } // if istr.eof ()
4843  } // if we read the number of rows
4844  } // if the input stream wrapping the dims line was (in)valid
4845  } // if we successfully read the banner line
4846  } // if (myRank == 0)
4847 
4848  // Check if file has a Vector
4849  if (dims[1]!=1) {
4850  exMsg << "File does not contain a 1D Vector";
4851  localDimsReadSuccess = 0;
4852  }
4853 
4854  // Broadcast the matrix dimensions, the encoded data type, and
4855  // whether or not Proc 0 succeeded in reading the banner and
4856  // dimensions.
4857  Tuple<GO, 5> bannerDimsReadResult;
4858  if (myRank == 0) {
4859  bannerDimsReadResult[0] = dims[0]; // numRows
4860  bannerDimsReadResult[1] = dims[1]; // numCols
4861  bannerDimsReadResult[2] = dims[2]; // encoded data type
4862  bannerDimsReadResult[3] = localBannerReadSuccess;
4863  bannerDimsReadResult[4] = localDimsReadSuccess;
4864  }
4865 
4866  // Broadcast matrix dimensions and the encoded data type from
4867  // Proc 0 to all the MPI processes.
4868  broadcast (*comm, 0, bannerDimsReadResult);
4869 
4870  TEUCHOS_TEST_FOR_EXCEPTION(
4871  bannerDimsReadResult[3] == 0, std::runtime_error,
4872  "Failed to read banner line: " << exMsg.str ());
4873  TEUCHOS_TEST_FOR_EXCEPTION(
4874  bannerDimsReadResult[4] == 0, std::runtime_error,
4875  "Failed to read matrix dimensions line: " << exMsg.str ());
4876  if (myRank != 0) {
4877  dims[0] = bannerDimsReadResult[0];
4878  dims[1] = bannerDimsReadResult[1];
4879  dims[2] = bannerDimsReadResult[2];
4880  }
4881 
4882  // Tpetra objects want the matrix dimensions in these types.
4883  const global_size_t numRows = static_cast<global_size_t> (dims[0]);
4884  const size_t numCols = static_cast<size_t> (dims[1]);
4885 
4886  // Make a "Proc 0 owns everything" Map that we will use to
4887  // read in the multivector entries in the correct order on
4888  // Proc 0. This must be a collective
4889  RCP<const map_type> proc0Map; // "Proc 0 owns everything" Map
4890  if (map.is_null ()) {
4891  // The user didn't supply a Map. Make a contiguous
4892  // distributed Map for them, using the read-in multivector
4893  // dimensions.
4894  map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4895  const size_t localNumRows = (myRank == 0) ? numRows : 0;
4896  // At this point, map exists and has a nonnull node.
4897  proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4898  comm);
4899  }
4900  else { // The user supplied a Map.
4901  proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4902  }
4903 
4904  // Make a multivector X owned entirely by Proc 0.
4905  RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
4906 
4907  //
4908  // On Proc 0, read the Matrix Market data from the input
4909  // stream into the multivector X.
4910  //
4911  int localReadDataSuccess = 1;
4912  if (myRank == 0) {
4913  try {
4914  if (debug) {
4915  *err << myRank << ": readVectorImpl: Reading matrix data (dense)"
4916  << endl;
4917  }
4918 
4919  // Make sure that we can get a 1-D view of X.
4920  TEUCHOS_TEST_FOR_EXCEPTION(
4921  ! X->isConstantStride (), std::logic_error,
4922  "Can't get a 1-D view of the entries of the MultiVector X on "
4923  "Process 0, because the stride between the columns of X is not "
4924  "constant. This shouldn't happen because we just created X and "
4925  "haven't filled it in yet. Please report this bug to the Tpetra "
4926  "developers.");
4927 
4928  // Get a writeable 1-D view of the entries of X. Rank 0
4929  // owns all of them. The view will expire at the end of
4930  // scope, so (if necessary) it will be written back to X
4931  // at this time.
4932  Teuchos::ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4933  TEUCHOS_TEST_FOR_EXCEPTION(
4934  as<global_size_t> (X_view.size ()) < numRows * numCols,
4935  std::logic_error,
4936  "The view of X has size " << X_view << " which is not enough to "
4937  "accommodate the expected number of entries numRows*numCols = "
4938  << numRows << "*" << numCols << " = " << numRows*numCols << ". "
4939  "Please report this bug to the Tpetra developers.");
4940  const size_t stride = X->getStride ();
4941 
4942  // The third element of the dimensions Tuple encodes the data
4943  // type reported by the Banner: "real" == 0, "complex" == 1,
4944  // "integer" == 0 (same as "real"), "pattern" == 2. We do not
4945  // allow dense matrices to be pattern matrices, so dims[2] ==
4946  // 0 or 1. We've already checked for this above.
4947  const bool isComplex = (dims[2] == 1);
4948  size_type count = 0, curRow = 0, curCol = 0;
4949 
4950  std::string line;
4951  while (getline (in, line)) {
4952  ++lineNumber;
4953  // Is the current line a comment line? If it's not,
4954  // line.substr(start,size) contains the data.
4955  size_t start = 0, size = 0;
4956  const bool commentLine =
4957  checkCommentLine (line, start, size, lineNumber, tolerant);
4958  if (! commentLine) {
4959  // Make sure we have room in which to put the new matrix
4960  // entry. We check this only after checking for a
4961  // comment line, because there may be one or more
4962  // comment lines at the end of the file. In tolerant
4963  // mode, we simply ignore any extra data.
4964  if (count >= X_view.size()) {
4965  if (tolerant) {
4966  break;
4967  }
4968  else {
4969  TEUCHOS_TEST_FOR_EXCEPTION(
4970  count >= X_view.size(),
4971  std::runtime_error,
4972  "The Matrix Market input stream has more data in it than "
4973  "its metadata reported. Current line number is "
4974  << lineNumber << ".");
4975  }
4976  }
4977 
4978  // mfh 19 Dec 2012: Ignore everything up to the initial
4979  // colon. writeDense() has the option to print out the
4980  // global row index in front of each entry, followed by
4981  // a colon and space.
4982  {
4983  const size_t pos = line.substr (start, size).find (':');
4984  if (pos != std::string::npos) {
4985  start = pos+1;
4986  }
4987  }
4988  std::istringstream istr (line.substr (start, size));
4989  // Does the line contain anything at all? Can we
4990  // safely read from the input stream wrapping the
4991  // line?
4992  if (istr.eof() || istr.fail()) {
4993  // In tolerant mode, simply ignore the line.
4994  if (tolerant) {
4995  break;
4996  }
4997  // We repeat the full test here so the exception
4998  // message is more informative.
4999  TEUCHOS_TEST_FOR_EXCEPTION(
5000  ! tolerant && (istr.eof() || istr.fail()),
5001  std::runtime_error,
5002  "Line " << lineNumber << " of the Matrix Market file is "
5003  "empty, or we cannot read from it for some other reason.");
5004  }
5005  // Current matrix entry to read in.
5006  ST val = STS::zero();
5007  // Real and imaginary parts of the current matrix entry.
5008  // The imaginary part is zero if the matrix is real-valued.
5009  MT real = STM::zero(), imag = STM::zero();
5010 
5011  // isComplex refers to the input stream's data, not to
5012  // the scalar type S. It's OK to read real-valued
5013  // data into a matrix storing complex-valued data; in
5014  // that case, all entries' imaginary parts are zero.
5015  if (isComplex) {
5016  // STS::real() and STS::imag() return a copy of
5017  // their respective components, not a writeable
5018  // reference. Otherwise we could just assign to
5019  // them using the istream extraction operator (>>).
5020  // That's why we have separate magnitude type "real"
5021  // and "imag" variables.
5022 
5023  // Attempt to read the real part of the current entry.
5024  istr >> real;
5025  if (istr.fail()) {
5026  TEUCHOS_TEST_FOR_EXCEPTION(
5027  ! tolerant && istr.eof(), std::runtime_error,
5028  "Failed to get the real part of a complex-valued matrix "
5029  "entry from line " << lineNumber << " of the Matrix Market "
5030  "file.");
5031  // In tolerant mode, just skip bad lines.
5032  if (tolerant) {
5033  break;
5034  }
5035  } else if (istr.eof()) {
5036  TEUCHOS_TEST_FOR_EXCEPTION(
5037  ! tolerant && istr.eof(), std::runtime_error,
5038  "Missing imaginary part of a complex-valued matrix entry "
5039  "on line " << lineNumber << " of the Matrix Market file.");
5040  // In tolerant mode, let any missing imaginary part be 0.
5041  } else {
5042  // Attempt to read the imaginary part of the current
5043  // matrix entry.
5044  istr >> imag;
5045  TEUCHOS_TEST_FOR_EXCEPTION(
5046  ! tolerant && istr.fail(), std::runtime_error,
5047  "Failed to get the imaginary part of a complex-valued "
5048  "matrix entry from line " << lineNumber << " of the "
5049  "Matrix Market file.");
5050  // In tolerant mode, let any missing or corrupted
5051  // imaginary part be 0.
5052  }
5053  } else { // Matrix Market file contains real-valued data.
5054  // Attempt to read the current matrix entry.
5055  istr >> real;
5056  TEUCHOS_TEST_FOR_EXCEPTION(
5057  ! tolerant && istr.fail(), std::runtime_error,
5058  "Failed to get a real-valued matrix entry from line "
5059  << lineNumber << " of the Matrix Market file.");
5060  // In tolerant mode, simply ignore the line if
5061  // we failed to read a matrix entry.
5062  if (istr.fail() && tolerant) {
5063  break;
5064  }
5065  }
5066  // In tolerant mode, we simply let pass through whatever
5067  // data we got.
5068  TEUCHOS_TEST_FOR_EXCEPTION(
5069  ! tolerant && istr.fail(), std::runtime_error,
5070  "Failed to read matrix data from line " << lineNumber
5071  << " of the Matrix Market file.");
5072 
5073  // Assign val = ST(real, imag).
5074  Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
5075 
5076  curRow = count % numRows;
5077  curCol = count / numRows;
5078  X_view[curRow + curCol*stride] = val;
5079  ++count;
5080  } // if not a comment line
5081  } // while there are still lines in the file, get the next one
5082 
5083  TEUCHOS_TEST_FOR_EXCEPTION(
5084  ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
5085  std::runtime_error,
5086  "The Matrix Market metadata reports that the dense matrix is "
5087  << numRows << " x " << numCols << ", and thus has "
5088  << numRows*numCols << " total entries, but we only found " << count
5089  << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
5090  } catch (std::exception& e) {
5091  exMsg << e.what ();
5092  localReadDataSuccess = 0;
5093  }
5094  } // if (myRank == 0)
5095 
5096  if (debug) {
5097  *err << myRank << ": readVectorImpl: done reading data" << endl;
5098  }
5099 
5100  // Synchronize on whether Proc 0 successfully read the data.
5101  int globalReadDataSuccess = localReadDataSuccess;
5102  broadcast (*comm, 0, outArg (globalReadDataSuccess));
5103  TEUCHOS_TEST_FOR_EXCEPTION(
5104  globalReadDataSuccess == 0, std::runtime_error,
5105  "Failed to read the multivector's data: " << exMsg.str ());
5106 
5107  // If there's only one MPI process and the user didn't supply
5108  // a Map (i.e., pMap is null), we're done. Set pMap to the
5109  // Map used to distribute X, and return X.
5110  if (comm->getSize () == 1 && map.is_null ()) {
5111  map = proc0Map;
5112  if (! err.is_null ()) {
5113  err->popTab ();
5114  }
5115  if (debug) {
5116  *err << myRank << ": readVectorImpl: done" << endl;
5117  }
5118  if (! err.is_null ()) {
5119  err->popTab ();
5120  }
5121  return X;
5122  }
5123 
5124  if (debug) {
5125  *err << myRank << ": readVectorImpl: Creating target MV" << endl;
5126  }
5127 
5128  // Make a multivector Y with the distributed map pMap.
5129  RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
5130 
5131  if (debug) {
5132  *err << myRank << ": readVectorImpl: Creating Export" << endl;
5133  }
5134 
5135  // Make an Export object that will export X to Y. First
5136  // argument is the source map, second argument is the target
5137  // map.
5138  Export<LO, GO, NT> exporter (proc0Map, map, err);
5139 
5140  if (debug) {
5141  *err << myRank << ": readVectorImpl: Exporting" << endl;
5142  }
5143  // Export X into Y.
5144  Y->doExport (*X, exporter, INSERT);
5145 
5146  if (! err.is_null ()) {
5147  err->popTab ();
5148  }
5149  if (debug) {
5150  *err << myRank << ": readVectorImpl: done" << endl;
5151  }
5152  if (! err.is_null ()) {
5153  err->popTab ();
5154  }
5155 
5156  // Y is distributed over all process(es) in the communicator.
5157  return Y;
5158  }
5159 
5160  public:
5181  static Teuchos::RCP<const map_type>
5182  readMap (std::istream& in,
5183  const trcp_tcomm_t& comm,
5184  const bool tolerant=false,
5185  const bool debug=false,
5186  const bool binary=false)
5187  {
5188  Teuchos::RCP<Teuchos::FancyOStream> err =
5189  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
5190  return readMap (in, comm, err, tolerant, debug, binary);
5191  }
5192 
5193 
5220  static Teuchos::RCP<const map_type>
5221  readMap (std::istream& in,
5222  const trcp_tcomm_t& comm,
5223  const Teuchos::RCP<Teuchos::FancyOStream>& err,
5224  const bool tolerant=false,
5225  const bool debug=false,
5226  const bool binary=false)
5227  {
5228  using Teuchos::arcp;
5229  using Teuchos::Array;
5230  using Teuchos::ArrayRCP;
5231  using Teuchos::as;
5232  using Teuchos::broadcast;
5233  using Teuchos::Comm;
5234  using Teuchos::CommRequest;
5235  using Teuchos::inOutArg;
5236  using Teuchos::ireceive;
5237  using Teuchos::outArg;
5238  using Teuchos::RCP;
5239  using Teuchos::receive;
5240  using Teuchos::reduceAll;
5241  using Teuchos::REDUCE_MIN;
5242  using Teuchos::isend;
5243  using Teuchos::SerialComm;
5244  using Teuchos::toString;
5245  using Teuchos::wait;
5246  using std::endl;
5247  typedef Tpetra::global_size_t GST;
5248  typedef ptrdiff_t int_type; // Can hold int and GO
5249  typedef local_ordinal_type LO;
5250  typedef global_ordinal_type GO;
5251  typedef node_type NT;
5253 
5254  const int numProcs = comm->getSize ();
5255  const int myRank = comm->getRank ();
5256 
5257  if (err.is_null ()) {
5258  err->pushTab ();
5259  }
5260  if (debug) {
5261  std::ostringstream os;
5262  os << myRank << ": readMap: " << endl;
5263  *err << os.str ();
5264  }
5265  if (err.is_null ()) {
5266  err->pushTab ();
5267  }
5268 
5269  // Tag for receive-size / send-size messages. writeMap used
5270  // tags 1337 and 1338; we count up from there.
5271  const int sizeTag = 1339;
5272  // Tag for receive-data / send-data messages.
5273  const int dataTag = 1340;
5274 
5275  // These are for sends on Process 0, and for receives on all
5276  // other processes. sizeReq is for the {receive,send}-size
5277  // message, and dataReq is for the message containing the
5278  // actual GIDs to belong to the receiving process.
5279  RCP<CommRequest<int> > sizeReq;
5280  RCP<CommRequest<int> > dataReq;
5281 
5282  // Each process will have to receive the number of GIDs to
5283  // expect. Thus, we can post the receives now, and cancel
5284  // them if something should go wrong in the meantime.
5285  ArrayRCP<int_type> numGidsToRecv (1);
5286  numGidsToRecv[0] = 0;
5287  if (myRank != 0) {
5288  sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
5289  }
5290 
5291  int readSuccess = 1;
5292  std::ostringstream exMsg;
5293  RCP<MV> data; // Will only be valid on Proc 0
5294  if (myRank == 0) {
5295  // If we want to reuse readDenseImpl, we have to make a
5296  // communicator that only contains Proc 0. Otherwise,
5297  // readDenseImpl will redistribute the data to all
5298  // processes. While we eventually want that, neither we nor
5299  // readDenseImpl know the correct Map to use at the moment.
5300  // That depends on the second column of the multivector.
5301  RCP<const Comm<int> > proc0Comm (new SerialComm<int> ());
5302  try {
5303  RCP<const map_type> dataMap;
5304  // This is currently the only place where we use the
5305  // 'tolerant' argument. Later, if we want to be clever,
5306  // we could have tolerant mode allow PIDs out of order.
5307  data = readDenseImpl<GO> (in, proc0Comm, dataMap, err, tolerant, debug, binary);
5308  (void) dataMap; // Silence "unused" warnings
5309  if (data.is_null ()) {
5310  readSuccess = 0;
5311  exMsg << "readDenseImpl() returned null." << endl;
5312  }
5313  } catch (std::exception& e) {
5314  readSuccess = 0;
5315  exMsg << e.what () << endl;
5316  }
5317  }
5318 
5319  // Map from PID to all the GIDs for that PID.
5320  // Only populated on Process 0.
5321  std::map<int, Array<GO> > pid2gids;
5322 
5323  // The index base must be the global minimum GID.
5324  // We will compute this on Process 0 and broadcast,
5325  // so that all processes can set up the Map.
5326  int_type globalNumGIDs = 0;
5327 
5328  // The index base must be the global minimum GID.
5329  // We will compute this on Process 0 and broadcast,
5330  // so that all processes can set up the Map.
5331  GO indexBase = 0;
5332 
5333  // Process 0: If the above read of the MultiVector succeeded,
5334  // extract the GIDs and PIDs into pid2gids, and find the
5335  // global min GID.
5336  if (myRank == 0 && readSuccess == 1) {
5337  if (data->getNumVectors () == 2) { // Map format 1.0
5338  ArrayRCP<const GO> GIDs = data->getData (0);
5339  ArrayRCP<const GO> PIDs = data->getData (1); // convert to int
5340  globalNumGIDs = GIDs.size ();
5341 
5342  // Start computing the global min GID, while collecting
5343  // the GIDs for each PID.
5344  if (globalNumGIDs > 0) {
5345  const int pid = static_cast<int> (PIDs[0]);
5346 
5347  if (pid < 0 || pid >= numProcs) {
5348  readSuccess = 0;
5349  exMsg << "Tpetra::MatrixMarket::readMap: "
5350  << "Encountered invalid PID " << pid << "." << endl;
5351  }
5352  else {
5353  const GO gid = GIDs[0];
5354  pid2gids[pid].push_back (gid);
5355  indexBase = gid; // the current min GID
5356  }
5357  }
5358  if (readSuccess == 1) {
5359  // Collect the rest of the GIDs for each PID, and compute
5360  // the global min GID.
5361  for (size_type k = 1; k < globalNumGIDs; ++k) {
5362  const int pid = static_cast<int> (PIDs[k]);
5363  if (pid < 0 || pid >= numProcs) {
5364  readSuccess = 0;
5365  exMsg << "Tpetra::MatrixMarket::readMap: "
5366  << "Encountered invalid PID " << pid << "." << endl;
5367  }
5368  else {
5369  const int_type gid = GIDs[k];
5370  pid2gids[pid].push_back (gid);
5371  if (gid < indexBase) {
5372  indexBase = gid; // the current min GID
5373  }
5374  }
5375  }
5376  }
5377  }
5378  else if (data->getNumVectors () == 1) { // Map format 2.0
5379  if (data->getGlobalLength () % 2 != static_cast<GST> (0)) {
5380  readSuccess = 0;
5381  exMsg << "Tpetra::MatrixMarket::readMap: Input data has the "
5382  "wrong format (for Map format 2.0). The global number of rows "
5383  "in the MultiVector must be even (divisible by 2)." << endl;
5384  }
5385  else {
5386  ArrayRCP<const GO> theData = data->getData (0);
5387  globalNumGIDs = static_cast<GO> (data->getGlobalLength ()) /
5388  static_cast<GO> (2);
5389 
5390  // Start computing the global min GID, while
5391  // collecting the GIDs for each PID.
5392  if (globalNumGIDs > 0) {
5393  const int pid = static_cast<int> (theData[1]);
5394  if (pid < 0 || pid >= numProcs) {
5395  readSuccess = 0;
5396  exMsg << "Tpetra::MatrixMarket::readMap: "
5397  << "Encountered invalid PID " << pid << "." << endl;
5398  }
5399  else {
5400  const GO gid = theData[0];
5401  pid2gids[pid].push_back (gid);
5402  indexBase = gid; // the current min GID
5403  }
5404  }
5405  // Collect the rest of the GIDs for each PID, and
5406  // compute the global min GID.
5407  for (int_type k = 1; k < globalNumGIDs; ++k) {
5408  const int pid = static_cast<int> (theData[2*k + 1]);
5409  if (pid < 0 || pid >= numProcs) {
5410  readSuccess = 0;
5411  exMsg << "Tpetra::MatrixMarket::readMap: "
5412  << "Encountered invalid PID " << pid << "." << endl;
5413  }
5414  else {
5415  const GO gid = theData[2*k];
5416  pid2gids[pid].push_back (gid);
5417  if (gid < indexBase) {
5418  indexBase = gid; // the current min GID
5419  }
5420  }
5421  } // for each GID
5422  } // if the amount of data is correct
5423  }
5424  else {
5425  readSuccess = 0;
5426  exMsg << "Tpetra::MatrixMarket::readMap: Input data must have "
5427  "either 1 column (for the new Map format 2.0) or 2 columns (for "
5428  "the old Map format 1.0).";
5429  }
5430  } // myRank is zero
5431 
5432  // Broadcast the indexBase, the global number of GIDs, and the
5433  // current success status. Use int_type for all of these.
5434  {
5435  int_type readResults[3];
5436  readResults[0] = static_cast<int_type> (indexBase);
5437  readResults[1] = static_cast<int_type> (globalNumGIDs);
5438  readResults[2] = static_cast<int_type> (readSuccess);
5439  broadcast<int, int_type> (*comm, 0, 3, readResults);
5440 
5441  indexBase = static_cast<GO> (readResults[0]);
5442  globalNumGIDs = static_cast<int_type> (readResults[1]);
5443  readSuccess = static_cast<int> (readResults[2]);
5444  }
5445 
5446  // Unwinding the stack will invoke sizeReq's destructor, which
5447  // will cancel the receive-size request on all processes that
5448  // posted it.
5449  TEUCHOS_TEST_FOR_EXCEPTION(
5450  readSuccess != 1, std::runtime_error,
5451  "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
5452  "following exception message: " << exMsg.str ());
5453 
5454  if (myRank == 0) {
5455  // Proc 0: Send each process' number of GIDs to that process.
5456  for (int p = 1; p < numProcs; ++p) {
5457  ArrayRCP<int_type> numGidsToSend (1);
5458 
5459  auto it = pid2gids.find (p);
5460  if (it == pid2gids.end ()) {
5461  numGidsToSend[0] = 0;
5462  } else {
5463  numGidsToSend[0] = it->second.size ();
5464  }
5465  sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
5466  wait<int> (*comm, outArg (sizeReq));
5467  }
5468  }
5469  else {
5470  // Wait on the receive-size message to finish.
5471  wait<int> (*comm, outArg (sizeReq));
5472  }
5473 
5474  // Allocate / get the array for my GIDs.
5475  // Only Process 0 will have its actual GIDs at this point.
5476  ArrayRCP<GO> myGids;
5477  int_type myNumGids = 0;
5478  if (myRank == 0) {
5479  GO* myGidsRaw = NULL;
5480 
5481  typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
5482  if (it != pid2gids.end ()) {
5483  myGidsRaw = it->second.getRawPtr ();
5484  myNumGids = it->second.size ();
5485  // Nonowning ArrayRCP just views the Array.
5486  myGids = arcp<GO> (myGidsRaw, 0, myNumGids, false);
5487  }
5488  }
5489  else { // myRank != 0
5490  myNumGids = numGidsToRecv[0];
5491  myGids = arcp<GO> (myNumGids);
5492  }
5493 
5494  if (myRank != 0) {
5495  // Post receive for data, now that we know how much data we
5496  // will receive. Only post receive if my process actually
5497  // has nonzero GIDs.
5498  if (myNumGids > 0) {
5499  dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
5500  }
5501  }
5502 
5503  for (int p = 1; p < numProcs; ++p) {
5504  if (myRank == 0) {
5505  ArrayRCP<GO> sendGids; // to send to Process p
5506  GO* sendGidsRaw = NULL;
5507  int_type numSendGids = 0;
5508 
5509  typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
5510  if (it != pid2gids.end ()) {
5511  numSendGids = it->second.size ();
5512  sendGidsRaw = it->second.getRawPtr ();
5513  sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids, false);
5514  }
5515  // Only send if that process actually has nonzero GIDs.
5516  if (numSendGids > 0) {
5517  dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
5518  }
5519  wait<int> (*comm, outArg (dataReq));
5520  }
5521  else if (myRank == p) {
5522  // Wait on my receive of GIDs to finish.
5523  wait<int> (*comm, outArg (dataReq));
5524  }
5525  } // for each process rank p in 1, 2, ..., numProcs-1
5526 
5527  if (debug) {
5528  std::ostringstream os;
5529  os << myRank << ": readMap: creating Map" << endl;
5530  *err << os.str ();
5531  }
5532  const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
5533  RCP<const map_type> newMap;
5534 
5535  // Create the Map; test whether the constructor threw. This
5536  // avoids deadlock and makes error reporting more readable.
5537 
5538  int lclSuccess = 1;
5539  int gblSuccess = 0; // output argument
5540  std::ostringstream errStrm;
5541  try {
5542  newMap = rcp (new map_type (INVALID, myGids (), indexBase, comm));
5543  }
5544  catch (std::exception& e) {
5545  lclSuccess = 0;
5546  errStrm << "Process " << comm->getRank () << " threw an exception: "
5547  << e.what () << std::endl;
5548  }
5549  catch (...) {
5550  lclSuccess = 0;
5551  errStrm << "Process " << comm->getRank () << " threw an exception "
5552  "not a subclass of std::exception" << std::endl;
5553  }
5554  Teuchos::reduceAll<int, int> (*comm, Teuchos::REDUCE_MIN,
5555  lclSuccess, Teuchos::outArg (gblSuccess));
5556  if (gblSuccess != 1) {
5557  Tpetra::Details::gathervPrint (std::cerr, errStrm.str (), *comm);
5558  }
5559  TEUCHOS_TEST_FOR_EXCEPTION(gblSuccess != 1, std::runtime_error, "Map constructor failed!");
5560 
5561  if (err.is_null ()) {
5562  err->popTab ();
5563  }
5564  if (debug) {
5565  std::ostringstream os;
5566  os << myRank << ": readMap: done" << endl;
5567  *err << os.str ();
5568  }
5569  if (err.is_null ()) {
5570  err->popTab ();
5571  }
5572  return newMap;
5573  }
5574 
5575 
5576  private:
5577 
5588  static int
5589  encodeDataType (const std::string& dataType)
5590  {
5591  if (dataType == "real" || dataType == "integer") {
5592  return 0;
5593  } else if (dataType == "complex") {
5594  return 1;
5595  } else if (dataType == "pattern") {
5596  return 2;
5597  } else {
5598  // We should never get here, since Banner validates the
5599  // reported data type and ensures it is one of the accepted
5600  // values.
5601  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
5602  "Unrecognized Matrix Market data type \"" << dataType
5603  << "\". We should never get here. "
5604  "Please report this bug to the Tpetra developers.");
5605  }
5606  }
5607 
5608  public:
5609 
5639  static Teuchos::RCP<sparse_matrix_type>
5640  readSparsePerRank (const std::string& filename_prefix,
5641  const std::string& filename_suffix,
5642  const Teuchos::RCP<const map_type>& rowMap,
5643  Teuchos::RCP<const map_type>& colMap,
5644  const Teuchos::RCP<const map_type>& domainMap,
5645  const Teuchos::RCP<const map_type>& rangeMap,
5646  const bool callFillComplete=true,
5647  const bool tolerant=false,
5648  const int ranksToReadAtOnce=8,
5649  const bool debug=false)
5650  {
5651  using ST = scalar_type;
5652  using LO = local_ordinal_type;
5653  using GO = global_ordinal_type;
5654  using STS = typename Teuchos::ScalarTraits<ST>;
5655  using Teuchos::RCP;
5656  using Teuchos::ArrayRCP;
5657  using Teuchos::arcp;
5658  using Teuchos::rcp;
5659 
5660  // Sanity Checks
5661  // Fast checks for invalid input. We can't check other
5662  // attributes of the Maps until we've read in the matrix
5663  // dimensions.
5664  TEUCHOS_TEST_FOR_EXCEPTION(
5665  rowMap.is_null (), std::invalid_argument,
5666  "Row Map must be nonnull.");
5667  Teuchos::RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
5668  TEUCHOS_TEST_FOR_EXCEPTION
5669  (comm.is_null (), std::invalid_argument,
5670  "The input row map's communicator (Teuchos::Comm object) is null.");
5671  TEUCHOS_TEST_FOR_EXCEPTION(
5672  rangeMap.is_null (), std::invalid_argument,
5673  "Range Map must be nonnull.");
5674  TEUCHOS_TEST_FOR_EXCEPTION(
5675  domainMap.is_null (), std::invalid_argument,
5676  "Domain Map must be nonnull.");
5677  TEUCHOS_TEST_FOR_EXCEPTION(
5678  domainMap->getComm().getRawPtr() != comm.getRawPtr(),
5679  std::invalid_argument,
5680  "The specified domain Map's communicator (domainMap->getComm())"
5681  "differs from row Map's communicator");
5682  TEUCHOS_TEST_FOR_EXCEPTION(
5683  rangeMap->getComm().getRawPtr() != comm.getRawPtr(),
5684  std::invalid_argument,
5685  "The specified range Map's communicator (rangeMap->getComm())"
5686  "differs from row Map's communicator");
5687 
5688  // Setup
5689  const int myRank = comm->getRank();
5690  const int numProc = comm->getSize();
5691  std::string filename = filename_prefix + std::to_string(myRank) + filename_suffix;
5692 
5693  // Bounds check the writing limits
5694  int rank_limit = std::min(std::max(ranksToReadAtOnce,1),numProc);
5695 
5696  // Data structures for constructor
5697  ArrayRCP<size_t> numEntriesPerRow;
5698  ArrayRCP<size_t> rowPtr;
5699  ArrayRCP<global_ordinal_type> colInd;
5700  ArrayRCP<scalar_type> values;
5701  std::ostringstream errMsg;
5702 
5704  // Start the reading of the banners to get
5705  // local row / nnz counts and then read the
5706  // data. We'll pack everything into big ol'
5707  // rowptr/colind/values ArrayRCPs
5708  bool success = true;
5709  int bannerIsCorrect = 1, readSuccess = 1;
5710  LO numRows, numCols, numNonzeros;
5711  for(int base_rank = 0; base_rank < numProc; base_rank += rank_limit) {
5712  int stop = std::min(base_rank+rank_limit,numProc);
5713 
5714  // Is my rank in this batch?
5715  if(base_rank <= myRank && myRank < stop) {
5716  // My turn to read
5717  std::ifstream in(filename);
5718  using Teuchos::MatrixMarket::Banner;
5719  size_t lineNumber = 1;
5720  RCP<const Banner> pBanner;
5721  try {
5722  pBanner = readBanner (in, lineNumber, tolerant, debug);
5723  }
5724  catch (std::exception& e) {
5725  errMsg << "Attempt to read the Matrix Market file's Banner line "
5726  "threw an exception: " << e.what();
5727  bannerIsCorrect = 0;
5728  }
5729  if (bannerIsCorrect) {
5730  // Validate the Banner for the case of a sparse matrix.
5731  // We validate on Proc 0, since it reads the Banner.
5732 
5733  // In intolerant mode, the matrix type must be "coordinate".
5734  if (! tolerant && pBanner->matrixType() != "coordinate") {
5735  bannerIsCorrect = 0;
5736  errMsg << "The Matrix Market input file must contain a "
5737  "\"coordinate\"-format sparse matrix in order to create a "
5738  "Tpetra::CrsMatrix object from it, but the file's matrix "
5739  "type is \"" << pBanner->matrixType() << "\" instead.";
5740  }
5741  // In tolerant mode, we allow the matrix type to be
5742  // anything other than "array" (which would mean that
5743  // the file contains a dense matrix).
5744  if (tolerant && pBanner->matrixType() == "array") {
5745  bannerIsCorrect = 0;
5746  errMsg << "Matrix Market file must contain a \"coordinate\"-"
5747  "format sparse matrix in order to create a Tpetra::CrsMatrix "
5748  "object from it, but the file's matrix type is \"array\" "
5749  "instead. That probably means the file contains dense matrix "
5750  "data.";
5751  }
5752  }
5753 
5754  // Unpacked coordinate matrix dimensions
5755  using Teuchos::MatrixMarket::readCoordinateDimensions;
5756  success = readCoordinateDimensions (in, numRows, numCols,
5757  numNonzeros, lineNumber,
5758  tolerant);
5759 
5760  // Sanity checking of headers
5761  TEUCHOS_TEST_FOR_EXCEPTION(numRows != (LO)rowMap->getLocalNumElements(), std::invalid_argument,
5762  "# rows in file does not match rowmap.");
5763  TEUCHOS_TEST_FOR_EXCEPTION(!colMap.is_null() && numCols != (LO)colMap->getLocalNumElements(), std::invalid_argument,
5764  "# rows in file does not match colmap.");
5765 
5766 
5767  // Read the data
5768  typedef Teuchos::MatrixMarket::Raw::Adder<scalar_type,global_ordinal_type> raw_adder_type;
5769  bool tolerant_required = true;
5770  Teuchos::RCP<raw_adder_type> pRaw =
5771  Teuchos::rcp (new raw_adder_type (numRows,numCols,numNonzeros,tolerant_required,debug));
5772  RCP<adder_type> pAdder = Teuchos::rcp (new adder_type (pRaw, pBanner->symmType ()));
5773 
5774  if (debug) {
5775  std::cerr << "-- Reading matrix data" << std::endl;
5776  }
5777 
5778  try {
5779  // Reader for "coordinate" format sparse matrix data.
5780  typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
5781  global_ordinal_type, scalar_type, STS::isComplex> reader_type;
5782  reader_type reader (pAdder);
5783 
5784  // Read the sparse matrix entries.
5785  std::pair<bool, std::vector<size_t> > results = reader.read (in, lineNumber, tolerant_required, debug);
5786 
5787 
5788  readSuccess = results.first ? 1 : 0;
5789  }
5790  catch (std::exception& e) {
5791  readSuccess = 0;
5792  errMsg << e.what();
5793  }
5794 
5796  // Create the CSR Arrays
5797  typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,global_ordinal_type> element_type;
5798 
5799  // Additively merge duplicate matrix entries.
5800  pAdder->getAdder()->merge ();
5801 
5802  // Get a temporary const view of the merged matrix entries.
5803  const std::vector<element_type>& entries = pAdder->getAdder()->getEntries();
5804 
5805  // Number of matrix entries (after merging).
5806  const size_t numEntries = (size_t)entries.size();
5807 
5808  if (debug) {
5809  std::cerr << "----- Proc "<<myRank<<": Matrix has numRows=" << numRows
5810  << " rows and numEntries=" << numEntries
5811  << " entries." << std::endl;
5812  }
5813 
5814 
5815  // Make space for the CSR matrix data. Converting to
5816  // CSR is easier if we fill numEntriesPerRow with zeros
5817  // at first.
5818  numEntriesPerRow = arcp<size_t> (numRows);
5819  std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
5820  rowPtr = arcp<size_t> (numRows+1);
5821  std::fill (rowPtr.begin(), rowPtr.end(), 0);
5822  colInd = arcp<global_ordinal_type> (numEntries);
5823  values = arcp<scalar_type> (numEntries);
5824 
5825  // Convert from array-of-structs coordinate format to CSR
5826  // (compressed sparse row) format.
5827  global_ordinal_type l_prvRow = 0;
5828  size_t curPos = 0;
5829  LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
5830  rowPtr[0] = 0;
5831  LO indexBase = rowMap->getIndexBase();
5832  for (curPos = 0; curPos < numEntries; ++curPos) {
5833  const element_type& curEntry = entries[curPos];
5834  const global_ordinal_type curRow = curEntry.rowIndex() + indexBase;
5835  LO l_curRow = rowMap->getLocalElement(curRow);
5836 
5837 
5838  TEUCHOS_TEST_FOR_EXCEPTION(l_curRow == INVALID,std::logic_error,
5839  "Current global row "<< curRow << " is invalid.");
5840 
5841  TEUCHOS_TEST_FOR_EXCEPTION(l_curRow < l_prvRow, std::logic_error,
5842  "Row indices are out of order, even though they are supposed "
5843  "to be sorted. curRow = " << l_curRow << ", prvRow = "
5844  << l_prvRow << ", at curPos = " << curPos << ". Please report "
5845  "this bug to the Tpetra developers.");
5846  if (l_curRow > l_prvRow) {
5847  for (LO r = l_prvRow+1; r <= l_curRow; ++r) {
5848  rowPtr[r] = curPos;
5849  }
5850  l_prvRow = l_curRow;
5851  }
5852  numEntriesPerRow[l_curRow]++;
5853  colInd[curPos] = curEntry.colIndex() + indexBase;
5854  values[curPos] = curEntry.value();
5855 
5856  }
5857  // rowPtr has one more entry than numEntriesPerRow. The
5858  // last entry of rowPtr is the number of entries in
5859  // colInd and values.
5860  rowPtr[numRows] = numEntries;
5861 
5862  }// end base_rank <= myRank < stop
5863 
5864  // Barrier between batches to keep the filesystem happy
5865  comm->barrier();
5866 
5867  }//end outer rank loop
5868 
5869 
5870  // Call the matrix constructor and fill. This isn't particularly efficient
5871  RCP<sparse_matrix_type> A;
5872  if(colMap.is_null()) {
5873  A=rcp(new sparse_matrix_type(rowMap,numEntriesPerRow()));
5874  for(size_t i=0; i<rowMap->getLocalNumElements(); i++) {
5875  GO g_row = rowMap->getGlobalElement(i);
5876  size_t start = rowPtr[i];
5877  size_t size = rowPtr[i+1] - rowPtr[i];
5878  if(size>0) {
5879  A->insertGlobalValues(g_row,size,&values[start],&colInd[start]);
5880  }
5881  }
5882  }
5883  else {
5884  throw std::runtime_error("Reading with a column map is not yet implemented");
5885  }
5886  RCP<const map_type> myDomainMap = domainMap.is_null() ? rowMap : domainMap;
5887  RCP<const map_type> myRangeMap = rangeMap.is_null() ? rowMap : rangeMap;
5888 
5889  A->fillComplete(myDomainMap,myRangeMap);
5890 
5891  if(!readSuccess)
5892  success = false;
5893  TEUCHOS_TEST_FOR_EXCEPTION(success == false, std::runtime_error,
5894  "Read failed.");
5895 
5896  return A;
5897  }// end readSparsePerRank
5898 
5899 
5900  }; // class Reader
5901 
5930  template<class SparseMatrixType>
5931  class Writer {
5932  public:
5934  typedef SparseMatrixType sparse_matrix_type;
5935  typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
5936 
5938  typedef typename SparseMatrixType::scalar_type scalar_type;
5940  typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
5946  typedef typename SparseMatrixType::global_ordinal_type global_ordinal_type;
5948  typedef typename SparseMatrixType::node_type node_type;
5949 
5951  typedef MultiVector<scalar_type,
5959 
5962 
5964  using trcp_tcomm_t = Teuchos::RCP<const Teuchos::Comm<int>>;
5965 
5997  static void
5998  writeSparseFile (const std::string& filename,
5999  const sparse_matrix_type& matrix,
6000  const std::string& matrixName,
6001  const std::string& matrixDescription,
6002  const bool debug=false)
6003  {
6004  trcp_tcomm_t comm = matrix.getComm ();
6005  TEUCHOS_TEST_FOR_EXCEPTION
6006  (comm.is_null (), std::invalid_argument,
6007  "The input matrix's communicator (Teuchos::Comm object) is null.");
6008  const int myRank = comm->getRank ();
6009 
6010  auto out = Writer::openOutFileOnRankZero(comm, filename, myRank, true);
6011 
6012  writeSparse (out, matrix, matrixName, matrixDescription, debug);
6013  // We can rely on the destructor of the output stream to close
6014  // the file on scope exit, even if writeSparse() throws an
6015  // exception.
6016  }
6017 
6019  static void
6020  writeSparseFile (const std::string& filename,
6021  const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6022  const std::string& matrixName,
6023  const std::string& matrixDescription,
6024  const bool debug=false)
6025  {
6026  TEUCHOS_TEST_FOR_EXCEPTION
6027  (pMatrix.is_null (), std::invalid_argument,
6028  "The input matrix is null.");
6029  writeSparseFile (filename, *pMatrix, matrixName,
6030  matrixDescription, debug);
6031  }
6032 
6052  static void
6053  writeSparseFile (const std::string& filename,
6054  const sparse_matrix_type& matrix,
6055  const bool debug=false)
6056  {
6057  writeSparseFile (filename, matrix, "", "", debug);
6058  }
6059 
6061  static void
6062  writeSparseFile (const std::string& filename,
6063  const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6064  const bool debug=false)
6065  {
6066  writeSparseFile (filename, *pMatrix, "", "", debug);
6067  }
6068 
6099  static void
6100  writeSparse (std::ostream& out,
6101  const sparse_matrix_type& matrix,
6102  const std::string& matrixName,
6103  const std::string& matrixDescription,
6104  const bool debug=false)
6105  {
6106  using Teuchos::ArrayView;
6107  using Teuchos::Comm;
6108  using Teuchos::FancyOStream;
6109  using Teuchos::getFancyOStream;
6110  using Teuchos::null;
6111  using Teuchos::RCP;
6112  using Teuchos::rcpFromRef;
6113  using std::cerr;
6114  using std::endl;
6115  using ST = scalar_type;
6116  using LO = local_ordinal_type;
6117  using GO = global_ordinal_type;
6118  using STS = typename Teuchos::ScalarTraits<ST>;
6119 
6120  // Make the output stream write floating-point numbers in
6121  // scientific notation. It will politely put the output
6122  // stream back to its state on input, when this scope
6123  // terminates.
6124  Teuchos::SetScientific<ST> sci (out);
6125 
6126  // Get the matrix's communicator.
6127  trcp_tcomm_t comm = matrix.getComm ();
6128  TEUCHOS_TEST_FOR_EXCEPTION(
6129  comm.is_null (), std::invalid_argument,
6130  "The input matrix's communicator (Teuchos::Comm object) is null.");
6131  const int myRank = comm->getRank ();
6132 
6133  // Optionally, make a stream for debugging output.
6134  RCP<FancyOStream> err =
6135  debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6136  if (debug) {
6137  std::ostringstream os;
6138  os << myRank << ": writeSparse" << endl;
6139  *err << os.str ();
6140  comm->barrier ();
6141  os << "-- " << myRank << ": past barrier" << endl;
6142  *err << os.str ();
6143  }
6144 
6145  // Whether to print debugging output to stderr.
6146  const bool debugPrint = debug && myRank == 0;
6147 
6148  RCP<const map_type> rowMap = matrix.getRowMap ();
6149  RCP<const map_type> colMap = matrix.getColMap ();
6150  RCP<const map_type> domainMap = matrix.getDomainMap ();
6151  RCP<const map_type> rangeMap = matrix.getRangeMap ();
6152 
6153  const global_size_t numRows = rangeMap->getGlobalNumElements ();
6154  const global_size_t numCols = domainMap->getGlobalNumElements ();
6155 
6156  if (debug && myRank == 0) {
6157  std::ostringstream os;
6158  os << "-- Input sparse matrix is:"
6159  << "---- " << numRows << " x " << numCols << endl
6160  << "---- "
6161  << (matrix.isGloballyIndexed() ? "Globally" : "Locally")
6162  << " indexed." << endl
6163  << "---- Its row map has " << rowMap->getGlobalNumElements ()
6164  << " elements." << endl
6165  << "---- Its col map has " << colMap->getGlobalNumElements ()
6166  << " elements." << endl;
6167  *err << os.str ();
6168  }
6169  // Make the "gather" row map, where Proc 0 owns all rows and
6170  // the other procs own no rows.
6171  const size_t localNumRows = (myRank == 0) ? numRows : 0;
6172  if (debug) {
6173  std::ostringstream os;
6174  os << "-- " << myRank << ": making gatherRowMap" << endl;
6175  *err << os.str ();
6176  }
6177  RCP<const map_type> gatherRowMap =
6178  Details::computeGatherMap (rowMap, err, debug);
6179 
6180  // Since the matrix may in general be non-square, we need to
6181  // make a column map as well. In this case, the column map
6182  // contains all the columns of the original matrix, because we
6183  // are gathering the whole matrix onto Proc 0. We call
6184  // computeGatherMap to preserve the original order of column
6185  // indices over all the processes.
6186  const size_t localNumCols = (myRank == 0) ? numCols : 0;
6187  RCP<const map_type> gatherColMap =
6188  Details::computeGatherMap (colMap, err, debug);
6189 
6190  // Current map is the source map, gather map is the target map.
6191  typedef Import<LO, GO, node_type> import_type;
6192  import_type importer (rowMap, gatherRowMap);
6193 
6194  // Create a new CrsMatrix to hold the result of the import.
6195  // The constructor needs a column map as well as a row map,
6196  // for the case that the matrix is not square.
6197  RCP<sparse_matrix_type> newMatrix =
6198  rcp (new sparse_matrix_type (gatherRowMap, gatherColMap,
6199  static_cast<size_t> (0)));
6200  // Import the sparse matrix onto Proc 0.
6201  newMatrix->doImport (matrix, importer, INSERT);
6202 
6203  // fillComplete() needs the domain and range maps for the case
6204  // that the matrix is not square.
6205  {
6206  RCP<const map_type> gatherDomainMap =
6207  rcp (new map_type (numCols, localNumCols,
6208  domainMap->getIndexBase (),
6209  comm));
6210  RCP<const map_type> gatherRangeMap =
6211  rcp (new map_type (numRows, localNumRows,
6212  rangeMap->getIndexBase (),
6213  comm));
6214  newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
6215  }
6216 
6217  if (debugPrint) {
6218  cerr << "-- Output sparse matrix is:"
6219  << "---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
6220  << " x "
6221  << newMatrix->getDomainMap ()->getGlobalNumElements ()
6222  << " with "
6223  << newMatrix->getGlobalNumEntries () << " entries;" << endl
6224  << "---- "
6225  << (newMatrix->isGloballyIndexed () ? "Globally" : "Locally")
6226  << " indexed." << endl
6227  << "---- Its row map has "
6228  << newMatrix->getRowMap ()->getGlobalNumElements ()
6229  << " elements, with index base "
6230  << newMatrix->getRowMap ()->getIndexBase () << "." << endl
6231  << "---- Its col map has "
6232  << newMatrix->getColMap ()->getGlobalNumElements ()
6233  << " elements, with index base "
6234  << newMatrix->getColMap ()->getIndexBase () << "." << endl
6235  << "---- Element count of output matrix's column Map may differ "
6236  << "from that of the input matrix's column Map, if some columns "
6237  << "of the matrix contain no entries." << endl;
6238  }
6239 
6240  //
6241  // Print the metadata and the matrix entries on Rank 0.
6242  //
6243  if (myRank == 0) {
6244  // Print the Matrix Market banner line. CrsMatrix stores
6245  // data nonsymmetrically ("general"). This implies that
6246  // readSparse() on a symmetrically stored input file,
6247  // followed by writeSparse() on the resulting sparse matrix,
6248  // will result in an output file with a different banner
6249  // line than the original input file.
6250  out << "%%MatrixMarket matrix coordinate "
6251  << (STS::isComplex ? "complex" : "real")
6252  << " general" << endl;
6253 
6254  // Print comments (the matrix name and / or description).
6255  if (matrixName != "") {
6256  printAsComment (out, matrixName);
6257  }
6258  if (matrixDescription != "") {
6259  printAsComment (out, matrixDescription);
6260  }
6261 
6262  // Print the Matrix Market header (# rows, # columns, #
6263  // nonzeros). Use the range resp. domain map for the number
6264  // of rows resp. columns, since Tpetra::CrsMatrix uses the
6265  // column map for the number of columns. That only
6266  // corresponds to the "linear-algebraic" number of columns
6267  // when the column map is uniquely owned (a.k.a. one-to-one),
6268  // which only happens if the matrix is (block) diagonal.
6269  out << newMatrix->getRangeMap ()->getGlobalNumElements () << " "
6270  << newMatrix->getDomainMap ()->getGlobalNumElements () << " "
6271  << newMatrix->getGlobalNumEntries () << endl;
6272 
6273  // The Matrix Market format expects one-based row and column
6274  // indices. We'll convert the indices on output from
6275  // whatever index base they use to one-based indices.
6276  const GO rowIndexBase = gatherRowMap->getIndexBase ();
6277  const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
6278  //
6279  // Print the entries of the matrix.
6280  //
6281  // newMatrix can never be globally indexed, since we called
6282  // fillComplete() on it. We include code for both cases
6283  // (globally or locally indexed) just in case that ever
6284  // changes.
6285  if (newMatrix->isGloballyIndexed()) {
6286  // We know that the "gather" row Map is contiguous, so we
6287  // don't need to get the list of GIDs.
6288  const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6289  const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6290  for (GO globalRowIndex = minAllGlobalIndex;
6291  globalRowIndex <= maxAllGlobalIndex; // inclusive range
6292  ++globalRowIndex) {
6293  typename sparse_matrix_type::global_inds_host_view_type ind;
6294  typename sparse_matrix_type::values_host_view_type val;
6295  newMatrix->getGlobalRowView (globalRowIndex, ind, val);
6296  for (size_t ii = 0; ii < ind.extent(0); ii++) {
6297  const GO globalColIndex = ind(ii);
6298  // Convert row and column indices to 1-based.
6299  // This works because the global index type is signed.
6300  out << (globalRowIndex + 1 - rowIndexBase) << " "
6301  << (globalColIndex + 1 - colIndexBase) << " ";
6302  if (STS::isComplex) {
6303  out << STS::real (val(ii)) << " " << STS::imag (val(ii));
6304  } else {
6305  out << val(ii);
6306  }
6307  out << endl;
6308  } // For each entry in the current row
6309  } // For each row of the "gather" matrix
6310  }
6311  else { // newMatrix is locally indexed
6312  using OTG = Teuchos::OrdinalTraits<GO>;
6313  for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
6314  localRowIndex <= gatherRowMap->getMaxLocalIndex();
6315  ++localRowIndex) {
6316  // Convert from local to global row index.
6317  const GO globalRowIndex =
6318  gatherRowMap->getGlobalElement (localRowIndex);
6319  TEUCHOS_TEST_FOR_EXCEPTION(
6320  globalRowIndex == OTG::invalid(), std::logic_error,
6321  "Failed to convert the supposed local row index "
6322  << localRowIndex << " into a global row index. "
6323  "Please report this bug to the Tpetra developers.");
6324  typename sparse_matrix_type::local_inds_host_view_type ind;
6325  typename sparse_matrix_type::values_host_view_type val;
6326  newMatrix->getLocalRowView (localRowIndex, ind, val);
6327  for (size_t ii = 0; ii < ind.extent(0); ii++) {
6328  // Convert the column index from local to global.
6329  const GO globalColIndex =
6330  newMatrix->getColMap()->getGlobalElement (ind(ii));
6331  TEUCHOS_TEST_FOR_EXCEPTION(
6332  globalColIndex == OTG::invalid(), std::logic_error,
6333  "On local row " << localRowIndex << " of the sparse matrix: "
6334  "Failed to convert the supposed local column index "
6335  << ind(ii) << " into a global column index. Please report "
6336  "this bug to the Tpetra developers.");
6337  // Convert row and column indices to 1-based.
6338  // This works because the global index type is signed.
6339  out << (globalRowIndex + 1 - rowIndexBase) << " "
6340  << (globalColIndex + 1 - colIndexBase) << " ";
6341  if (STS::isComplex) {
6342  out << STS::real (val(ii)) << " " << STS::imag (val(ii));
6343  } else {
6344  out << val(ii);
6345  }
6346  out << endl;
6347  } // For each entry in the current row
6348  } // For each row of the "gather" matrix
6349  } // Whether the "gather" matrix is locally or globally indexed
6350  } // If my process' rank is 0
6351  }
6352 
6354  static void
6355  writeSparse (std::ostream& out,
6356  const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6357  const std::string& matrixName,
6358  const std::string& matrixDescription,
6359  const bool debug=false)
6360  {
6361  TEUCHOS_TEST_FOR_EXCEPTION
6362  (pMatrix.is_null (), std::invalid_argument,
6363  "The input matrix is null.");
6364  writeSparse (out, *pMatrix, matrixName, matrixDescription, debug);
6365  }
6366 
6397  static void
6398  writeSparseGraph (std::ostream& out,
6399  const crs_graph_type& graph,
6400  const std::string& graphName,
6401  const std::string& graphDescription,
6402  const bool debug=false)
6403  {
6404  using Teuchos::ArrayView;
6405  using Teuchos::Comm;
6406  using Teuchos::FancyOStream;
6407  using Teuchos::getFancyOStream;
6408  using Teuchos::null;
6409  using Teuchos::RCP;
6410  using Teuchos::rcpFromRef;
6411  using std::cerr;
6412  using std::endl;
6413  typedef local_ordinal_type LO;
6414  typedef global_ordinal_type GO;
6415 
6416  // Get the graph's communicator. Processes on which the
6417  // graph's Map or communicator is null don't participate in
6418  // this operation. This function shouldn't even be called on
6419  // those processes.
6420  auto rowMap = graph.getRowMap ();
6421  if (rowMap.is_null ()) {
6422  return;
6423  }
6424  auto comm = rowMap->getComm ();
6425  if (comm.is_null ()) {
6426  return;
6427  }
6428  const int myRank = comm->getRank ();
6429 
6430  // Optionally, make a stream for debugging output.
6431  RCP<FancyOStream> err =
6432  debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
6433  if (debug) {
6434  std::ostringstream os;
6435  os << myRank << ": writeSparseGraph" << endl;
6436  *err << os.str ();
6437  comm->barrier ();
6438  os << "-- " << myRank << ": past barrier" << endl;
6439  *err << os.str ();
6440  }
6441 
6442  // Whether to print debugging output to stderr.
6443  const bool debugPrint = debug && myRank == 0;
6444 
6445  // We've already gotten the rowMap above.
6446  auto colMap = graph.getColMap ();
6447  auto domainMap = graph.getDomainMap ();
6448  auto rangeMap = graph.getRangeMap ();
6449 
6450  const global_size_t numRows = rangeMap->getGlobalNumElements ();
6451  const global_size_t numCols = domainMap->getGlobalNumElements ();
6452 
6453  if (debug && myRank == 0) {
6454  std::ostringstream os;
6455  os << "-- Input sparse graph is:"
6456  << "---- " << numRows << " x " << numCols << " with "
6457  << graph.getGlobalNumEntries () << " entries;" << endl
6458  << "---- "
6459  << (graph.isGloballyIndexed () ? "Globally" : "Locally")
6460  << " indexed." << endl
6461  << "---- Its row Map has " << rowMap->getGlobalNumElements ()
6462  << " elements." << endl
6463  << "---- Its col Map has " << colMap->getGlobalNumElements ()
6464  << " elements." << endl;
6465  *err << os.str ();
6466  }
6467  // Make the "gather" row map, where Proc 0 owns all rows and
6468  // the other procs own no rows.
6469  const size_t localNumRows = (myRank == 0) ? numRows : 0;
6470  if (debug) {
6471  std::ostringstream os;
6472  os << "-- " << myRank << ": making gatherRowMap" << endl;
6473  *err << os.str ();
6474  }
6475  auto gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
6476 
6477  // Since the graph may in general be non-square, we need to
6478  // make a column map as well. In this case, the column map
6479  // contains all the columns of the original graph, because we
6480  // are gathering the whole graph onto Proc 0. We call
6481  // computeGatherMap to preserve the original order of column
6482  // indices over all the processes.
6483  const size_t localNumCols = (myRank == 0) ? numCols : 0;
6484  auto gatherColMap = Details::computeGatherMap (colMap, err, debug);
6485 
6486  // Current map is the source map, gather map is the target map.
6487  Import<LO, GO, node_type> importer (rowMap, gatherRowMap);
6488 
6489  // Create a new CrsGraph to hold the result of the import.
6490  // The constructor needs a column map as well as a row map,
6491  // for the case that the graph is not square.
6492  crs_graph_type newGraph (gatherRowMap, gatherColMap,
6493  static_cast<size_t> (0));
6494  // Import the sparse graph onto Proc 0.
6495  newGraph.doImport (graph, importer, INSERT);
6496 
6497  // fillComplete() needs the domain and range maps for the case
6498  // that the graph is not square.
6499  {
6500  RCP<const map_type> gatherDomainMap =
6501  rcp (new map_type (numCols, localNumCols,
6502  domainMap->getIndexBase (),
6503  comm));
6504  RCP<const map_type> gatherRangeMap =
6505  rcp (new map_type (numRows, localNumRows,
6506  rangeMap->getIndexBase (),
6507  comm));
6508  newGraph.fillComplete (gatherDomainMap, gatherRangeMap);
6509  }
6510 
6511  if (debugPrint) {
6512  cerr << "-- Output sparse graph is:"
6513  << "---- " << newGraph.getRangeMap ()->getGlobalNumElements ()
6514  << " x "
6515  << newGraph.getDomainMap ()->getGlobalNumElements ()
6516  << " with "
6517  << newGraph.getGlobalNumEntries () << " entries;" << endl
6518  << "---- "
6519  << (newGraph.isGloballyIndexed () ? "Globally" : "Locally")
6520  << " indexed." << endl
6521  << "---- Its row map has "
6522  << newGraph.getRowMap ()->getGlobalNumElements ()
6523  << " elements, with index base "
6524  << newGraph.getRowMap ()->getIndexBase () << "." << endl
6525  << "---- Its col map has "
6526  << newGraph.getColMap ()->getGlobalNumElements ()
6527  << " elements, with index base "
6528  << newGraph.getColMap ()->getIndexBase () << "." << endl
6529  << "---- Element count of output graph's column Map may differ "
6530  << "from that of the input matrix's column Map, if some columns "
6531  << "of the matrix contain no entries." << endl;
6532  }
6533 
6534  //
6535  // Print the metadata and the graph entries on Process 0 of
6536  // the graph's communicator.
6537  //
6538  if (myRank == 0) {
6539  // Print the Matrix Market banner line. CrsGraph stores
6540  // data nonsymmetrically ("general"). This implies that
6541  // readSparseGraph() on a symmetrically stored input file,
6542  // followed by writeSparseGraph() on the resulting sparse
6543  // graph, will result in an output file with a different
6544  // banner line than the original input file.
6545  out << "%%MatrixMarket matrix coordinate pattern general" << endl;
6546 
6547  // Print comments (the graph name and / or description).
6548  if (graphName != "") {
6549  printAsComment (out, graphName);
6550  }
6551  if (graphDescription != "") {
6552  printAsComment (out, graphDescription);
6553  }
6554 
6555  // Print the Matrix Market header (# rows, # columns, #
6556  // stored entries). Use the range resp. domain map for the
6557  // number of rows resp. columns, since Tpetra::CrsGraph uses
6558  // the column map for the number of columns. That only
6559  // corresponds to the "linear-algebraic" number of columns
6560  // when the column map is uniquely owned
6561  // (a.k.a. one-to-one), which only happens if the graph is
6562  // block diagonal (one block per process).
6563  out << newGraph.getRangeMap ()->getGlobalNumElements () << " "
6564  << newGraph.getDomainMap ()->getGlobalNumElements () << " "
6565  << newGraph.getGlobalNumEntries () << endl;
6566 
6567  // The Matrix Market format expects one-based row and column
6568  // indices. We'll convert the indices on output from
6569  // whatever index base they use to one-based indices.
6570  const GO rowIndexBase = gatherRowMap->getIndexBase ();
6571  const GO colIndexBase = newGraph.getColMap()->getIndexBase ();
6572  //
6573  // Print the entries of the graph.
6574  //
6575  // newGraph can never be globally indexed, since we called
6576  // fillComplete() on it. We include code for both cases
6577  // (globally or locally indexed) just in case that ever
6578  // changes.
6579  if (newGraph.isGloballyIndexed ()) {
6580  // We know that the "gather" row Map is contiguous, so we
6581  // don't need to get the list of GIDs.
6582  const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
6583  const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
6584  for (GO globalRowIndex = minAllGlobalIndex;
6585  globalRowIndex <= maxAllGlobalIndex; // inclusive range
6586  ++globalRowIndex) {
6587  typename crs_graph_type::global_inds_host_view_type ind;
6588  newGraph.getGlobalRowView (globalRowIndex, ind);
6589  for (size_t ii = 0; ii < ind.extent(0); ii++) {
6590  const GO globalColIndex = ind(ii);
6591  // Convert row and column indices to 1-based.
6592  // This works because the global index type is signed.
6593  out << (globalRowIndex + 1 - rowIndexBase) << " "
6594  << (globalColIndex + 1 - colIndexBase) << " ";
6595  out << endl;
6596  } // For each entry in the current row
6597  } // For each row of the "gather" graph
6598  }
6599  else { // newGraph is locally indexed
6600  typedef Teuchos::OrdinalTraits<GO> OTG;
6601  for (LO localRowIndex = gatherRowMap->getMinLocalIndex ();
6602  localRowIndex <= gatherRowMap->getMaxLocalIndex ();
6603  ++localRowIndex) {
6604  // Convert from local to global row index.
6605  const GO globalRowIndex =
6606  gatherRowMap->getGlobalElement (localRowIndex);
6607  TEUCHOS_TEST_FOR_EXCEPTION
6608  (globalRowIndex == OTG::invalid (), std::logic_error, "Failed "
6609  "to convert the supposed local row index " << localRowIndex <<
6610  " into a global row index. Please report this bug to the "
6611  "Tpetra developers.");
6612  typename crs_graph_type::local_inds_host_view_type ind;
6613  newGraph.getLocalRowView (localRowIndex, ind);
6614  for (size_t ii = 0; ii < ind.extent(0); ii++) {
6615  // Convert the column index from local to global.
6616  const GO globalColIndex =
6617  newGraph.getColMap ()->getGlobalElement (ind(ii));
6618  TEUCHOS_TEST_FOR_EXCEPTION(
6619  globalColIndex == OTG::invalid(), std::logic_error,
6620  "On local row " << localRowIndex << " of the sparse graph: "
6621  "Failed to convert the supposed local column index "
6622  << ind(ii) << " into a global column index. Please report "
6623  "this bug to the Tpetra developers.");
6624  // Convert row and column indices to 1-based.
6625  // This works because the global index type is signed.
6626  out << (globalRowIndex + 1 - rowIndexBase) << " "
6627  << (globalColIndex + 1 - colIndexBase) << " ";
6628  out << endl;
6629  } // For each entry in the current row
6630  } // For each row of the "gather" graph
6631  } // Whether the "gather" graph is locally or globally indexed
6632  } // If my process' rank is 0
6633  }
6634 
6640  static void
6641  writeSparseGraph (std::ostream& out,
6642  const crs_graph_type& graph,
6643  const bool debug=false)
6644  {
6645  writeSparseGraph (out, graph, "", "", debug);
6646  }
6647 
6682  static void
6683  writeSparseGraphFile (const std::string& filename,
6684  const crs_graph_type& graph,
6685  const std::string& graphName,
6686  const std::string& graphDescription,
6687  const bool debug=false)
6688  {
6689  auto comm = graph.getComm ();
6690  if (comm.is_null ()) {
6691  // Processes on which the communicator is null shouldn't
6692  // even call this function. The convention is that
6693  // processes on which the object's communicator is null do
6694  // not participate in collective operations involving the
6695  // object.
6696  return;
6697  }
6698  const int myRank = comm->getRank ();
6699 
6700  auto out = Writer::openOutFileOnRankZero(comm, filename, myRank, true);
6701 
6702  writeSparseGraph (out, graph, graphName, graphDescription, debug);
6703  // We can rely on the destructor of the output stream to close
6704  // the file on scope exit, even if writeSparseGraph() throws
6705  // an exception.
6706  }
6707 
6712  static void
6713  writeSparseGraphFile (const std::string& filename,
6714  const crs_graph_type& graph,
6715  const bool debug=false)
6716  {
6717  writeSparseGraphFile (filename, graph, "", "", debug);
6718  }
6719 
6728  static void
6729  writeSparseGraphFile (const std::string& filename,
6730  const Teuchos::RCP<const crs_graph_type>& pGraph,
6731  const std::string& graphName,
6732  const std::string& graphDescription,
6733  const bool debug=false)
6734  {
6735  writeSparseGraphFile (filename, *pGraph, graphName, graphDescription, debug);
6736  }
6737 
6747  static void
6748  writeSparseGraphFile (const std::string& filename,
6749  const Teuchos::RCP<const crs_graph_type>& pGraph,
6750  const bool debug=false)
6751  {
6752  writeSparseGraphFile (filename, *pGraph, "", "", debug);
6753  }
6754 
6777  static void
6778  writeSparse (std::ostream& out,
6779  const sparse_matrix_type& matrix,
6780  const bool debug=false)
6781  {
6782  writeSparse (out, matrix, "", "", debug);
6783  }
6784 
6786  static void
6787  writeSparse (std::ostream& out,
6788  const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
6789  const bool debug=false)
6790  {
6791  writeSparse (out, *pMatrix, "", "", debug);
6792  }
6793 
6822  static void
6823  writeDenseFile (const std::string& filename,
6824  const multivector_type& X,
6825  const std::string& matrixName,
6826  const std::string& matrixDescription,
6827  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6828  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6829  {
6830  trcp_tcomm_t comm = Writer::getComm(X.getMap());
6831  const int myRank = Writer::getRank(comm);
6832 
6833  auto out = Writer::openOutFileOnRankZero(comm, filename, myRank, true);
6834 
6835  writeDense (out, X, matrixName, matrixDescription, err, dbg);
6836  // We can rely on the destructor of the output stream to close
6837  // the file on scope exit, even if writeDense() throws an
6838  // exception.
6839  }
6840 
6846  static void
6847  writeDenseFile (const std::string& filename,
6848  const Teuchos::RCP<const multivector_type>& X,
6849  const std::string& matrixName,
6850  const std::string& matrixDescription,
6851  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6852  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6853  {
6854  TEUCHOS_TEST_FOR_EXCEPTION(
6855  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
6856  "writeDenseFile: The input MultiVector X is null.");
6857  writeDenseFile (filename, *X, matrixName, matrixDescription, err, dbg);
6858  }
6859 
6865  static void
6866  writeDenseFile (const std::string& filename,
6867  const multivector_type& X,
6868  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6869  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6870  {
6871  writeDenseFile (filename, X, "", "", err, dbg);
6872  }
6873 
6879  static void
6880  writeDenseFile (const std::string& filename,
6881  const Teuchos::RCP<const multivector_type>& X,
6882  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6883  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6884  {
6885  TEUCHOS_TEST_FOR_EXCEPTION(
6886  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
6887  "writeDenseFile: The input MultiVector X is null.");
6888  writeDenseFile (filename, *X, err, dbg);
6889  }
6890 
6891 
6922  static void
6923  writeDense (std::ostream& out,
6924  const multivector_type& X,
6925  const std::string& matrixName,
6926  const std::string& matrixDescription,
6927  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6928  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6929  {
6930  using Teuchos::Comm;
6931  using Teuchos::outArg;
6932  using Teuchos::REDUCE_MAX;
6933  using Teuchos::reduceAll;
6934  using std::endl;
6935 
6936  trcp_tcomm_t comm = Writer::getComm(X.getMap());
6937  const int myRank = Writer::getRank(comm);
6938 
6939  // If the caller provides a nonnull debug output stream, we
6940  // print debugging output to it. This is a local thing; we
6941  // don't have to check across processes.
6942  const bool debug = ! dbg.is_null ();
6943  if (debug) {
6944  dbg->pushTab ();
6945  std::ostringstream os;
6946  os << myRank << ": writeDense" << endl;
6947  *dbg << os.str ();
6948  dbg->pushTab ();
6949  }
6950  // Print the Matrix Market header.
6951  writeDenseHeader (out, X, matrixName, matrixDescription, err, dbg);
6952 
6953  // Print each column one at a time. This is a (perhaps)
6954  // temporary fix for Bug 6288.
6955  const size_t numVecs = X.getNumVectors ();
6956  for (size_t j = 0; j < numVecs; ++j) {
6957  writeDenseColumn (out, * (X.getVector (j)), err, dbg);
6958  }
6959 
6960  if (debug) {
6961  dbg->popTab ();
6962  std::ostringstream os;
6963  os << myRank << ": writeDense: Done" << endl;
6964  *dbg << os.str ();
6965  dbg->popTab ();
6966  }
6967  }
6968 
6969  private:
6975  static std::ofstream openOutFileOnRankZero(
6976  const trcp_tcomm_t& comm,
6977  const std::string& filename, const int rank, const bool safe = true,
6978  const std::ios_base::openmode mode = std::ios_base::out
6979  ){
6980  // Placeholder for the output stream.
6981  std::ofstream out;
6982 
6983  // State that will make all ranks throw if the root rank wasn't able to open the stream (using @c int for broadcasting).
6984  int all_should_stop = 0;
6985 
6986  // Try to open the file and update the state.
6987  if(rank == 0) {
6988  out.open(filename, mode);
6989  all_should_stop = !out && safe;
6990  }
6991 
6992  // Broadcast the stream state and throw from all ranks if needed.
6993  if(comm) Teuchos::broadcast(*comm, 0, &all_should_stop);
6994 
6995  TEUCHOS_TEST_FOR_EXCEPTION(
6996  all_should_stop,
6997  std::runtime_error,
6998  "Could not open output file '" + filename + "' on root rank 0."
6999  );
7000 
7001  return out;
7002  }
7003 
7029  static void
7030  writeDenseHeader (std::ostream& out,
7031  const multivector_type& X,
7032  const std::string& matrixName,
7033  const std::string& matrixDescription,
7034  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7035  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7036  {
7037  using Teuchos::Comm;
7038  using Teuchos::outArg;
7039  using Teuchos::RCP;
7040  using Teuchos::REDUCE_MAX;
7041  using Teuchos::reduceAll;
7042  using std::endl;
7043  typedef Teuchos::ScalarTraits<scalar_type> STS;
7044  const char prefix[] = "Tpetra::MatrixMarket::writeDenseHeader: ";
7045 
7046  trcp_tcomm_t comm = Writer::getComm(X.getMap());
7047  const int myRank = Writer::getRank(comm);
7048  int lclErr = 0; // whether this MPI process has seen an error
7049  int gblErr = 0; // whether we know if some MPI process has seen an error
7050 
7051  // If the caller provides a nonnull debug output stream, we
7052  // print debugging output to it. This is a local thing; we
7053  // don't have to check across processes.
7054  const bool debug = ! dbg.is_null ();
7055 
7056  if (debug) {
7057  dbg->pushTab ();
7058  std::ostringstream os;
7059  os << myRank << ": writeDenseHeader" << endl;
7060  *dbg << os.str ();
7061  dbg->pushTab ();
7062  }
7063 
7064  //
7065  // Process 0: Write the MatrixMarket header.
7066  //
7067  if (myRank == 0) {
7068  try {
7069  // Print the Matrix Market header. MultiVector stores data
7070  // nonsymmetrically, hence "general" in the banner line.
7071  // Print first to a temporary string output stream, and then
7072  // write it to the main output stream, so that at least the
7073  // header output has transactional semantics. We can't
7074  // guarantee transactional semantics for the whole output,
7075  // since that would not be memory scalable. (This could be
7076  // done in the file system by using a temporary file; we
7077  // don't do this, but users could.)
7078  std::ostringstream hdr;
7079  {
7080  std::string dataType;
7081  if (STS::isComplex) {
7082  dataType = "complex";
7083  } else if (STS::isOrdinal) {
7084  dataType = "integer";
7085  } else {
7086  dataType = "real";
7087  }
7088  hdr << "%%MatrixMarket matrix array " << dataType << " general"
7089  << endl;
7090  }
7091 
7092  // Print comments (the matrix name and / or description).
7093  if (matrixName != "") {
7094  printAsComment (hdr, matrixName);
7095  }
7096  if (matrixDescription != "") {
7097  printAsComment (hdr, matrixDescription);
7098  }
7099  // Print the Matrix Market dimensions header for dense matrices.
7100  hdr << X.getGlobalLength () << " " << X.getNumVectors () << endl;
7101 
7102  // Write the MatrixMarket header to the output stream.
7103  out << hdr.str ();
7104  } catch (std::exception& e) {
7105  if (! err.is_null ()) {
7106  *err << prefix << "While writing the Matrix Market header, "
7107  "Process 0 threw an exception: " << e.what () << endl;
7108  }
7109  lclErr = 1;
7110  }
7111  } // if I am Process 0
7112 
7113  // Establish global agreement on the error state. It wouldn't
7114  // be good for other processes to keep going, if Process 0
7115  // finds out that it can't write to the given output stream.
7116  reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
7117  TEUCHOS_TEST_FOR_EXCEPTION(
7118  gblErr == 1, std::runtime_error, prefix << "Some error occurred "
7119  "which prevented this method from completing.");
7120 
7121  if (debug) {
7122  dbg->popTab ();
7123  *dbg << myRank << ": writeDenseHeader: Done" << endl;
7124  dbg->popTab ();
7125  }
7126  }
7127 
7145  static void
7146  writeDenseColumn (std::ostream& out,
7147  const multivector_type& X,
7148  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7149  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7150  {
7151  using Teuchos::arcp;
7152  using Teuchos::Array;
7153  using Teuchos::ArrayRCP;
7154  using Teuchos::ArrayView;
7155  using Teuchos::Comm;
7156  using Teuchos::CommRequest;
7157  using Teuchos::ireceive;
7158  using Teuchos::isend;
7159  using Teuchos::outArg;
7160  using Teuchos::REDUCE_MAX;
7161  using Teuchos::reduceAll;
7162  using Teuchos::RCP;
7163  using Teuchos::TypeNameTraits;
7164  using Teuchos::wait;
7165  using std::endl;
7166  typedef Teuchos::ScalarTraits<scalar_type> STS;
7167 
7168  const Comm<int>& comm = * (X.getMap ()->getComm ());
7169  const int myRank = comm.getRank ();
7170  const int numProcs = comm.getSize ();
7171  int lclErr = 0; // whether this MPI process has seen an error
7172  int gblErr = 0; // whether we know if some MPI process has seen an error
7173 
7174  // If the caller provides a nonnull debug output stream, we
7175  // print debugging output to it. This is a local thing; we
7176  // don't have to check across processes.
7177  const bool debug = ! dbg.is_null ();
7178 
7179  if (debug) {
7180  dbg->pushTab ();
7181  std::ostringstream os;
7182  os << myRank << ": writeDenseColumn" << endl;
7183  *dbg << os.str ();
7184  dbg->pushTab ();
7185  }
7186 
7187  // Make the output stream write floating-point numbers in
7188  // scientific notation. It will politely put the output
7189  // stream back to its state on input, when this scope
7190  // terminates.
7191  Teuchos::SetScientific<scalar_type> sci (out);
7192 
7193  const size_t myNumRows = X.getLocalLength ();
7194  const size_t numCols = X.getNumVectors ();
7195  // Use a different tag for the "size" messages than for the
7196  // "data" messages, in order to help us debug any mix-ups.
7197  const int sizeTag = 1337;
7198  const int dataTag = 1338;
7199 
7200  // Process 0 pipelines nonblocking receives with file output.
7201  //
7202  // Constraints:
7203  // - Process 0 can't post a receive for another process'
7204  // actual data, until it posts and waits on the receive
7205  // from that process with the amount of data to receive.
7206  // (We could just post receives with a max data size, but
7207  // I feel uncomfortable about that.)
7208  // - The C++ standard library doesn't allow nonblocking
7209  // output to an std::ostream. (Thus, we have to start a
7210  // receive or send before starting the write, and hope
7211  // that MPI completes it in the background.)
7212  //
7213  // Process 0: Post receive-size receives from Processes 1 and 2.
7214  // Process 1: Post send-size send to Process 0.
7215  // Process 2: Post send-size send to Process 0.
7216  //
7217  // All processes: Pack my entries.
7218  //
7219  // Process 1:
7220  // - Post send-data send to Process 0.
7221  // - Wait on my send-size send to Process 0.
7222  //
7223  // Process 0:
7224  // - Print MatrixMarket header.
7225  // - Print my entries.
7226  // - Wait on receive-size receive from Process 1.
7227  // - Post receive-data receive from Process 1.
7228  //
7229  // For each process p = 1, 2, ... numProcs-1:
7230  // If I am Process 0:
7231  // - Post receive-size receive from Process p + 2
7232  // - Wait on receive-size receive from Process p + 1
7233  // - Post receive-data receive from Process p + 1
7234  // - Wait on receive-data receive from Process p
7235  // - Write data from Process p.
7236  // Else if I am Process p:
7237  // - Wait on my send-data send.
7238  // Else if I am Process p+1:
7239  // - Post send-data send to Process 0.
7240  // - Wait on my send-size send.
7241  // Else if I am Process p+2:
7242  // - Post send-size send to Process 0.
7243  //
7244  // Pipelining has three goals here:
7245  // 1. Overlap communication (the receives) with file I/O
7246  // 2. Give Process 0 a chance to prepost some receives,
7247  // before sends show up, by packing local data before
7248  // posting sends
7249  // 3. Don't post _all_ receives or _all_ sends, because that
7250  // wouldn't be memory scalable. (Just because we can't
7251  // see how much memory MPI consumes, doesn't mean that it
7252  // doesn't consume any!)
7253 
7254  // These are used on every process. sendReqSize[0] holds the
7255  // number of rows on this process, and sendReqBuf holds this
7256  // process' data. Process 0 packs into sendReqBuf, but
7257  // doesn't send; it only uses that for printing. All other
7258  // processes send both of these to Process 0.
7259  RCP<CommRequest<int> > sendReqSize, sendReqData;
7260 
7261  // These are used only on Process 0, for received data. Keep
7262  // 3 of each, and treat the arrays as circular buffers. When
7263  // receiving from Process p, the corresponding array index
7264  // here is p % 3.
7265  Array<ArrayRCP<size_t> > recvSizeBufs (3);
7266  Array<ArrayRCP<scalar_type> > recvDataBufs (3);
7267  Array<RCP<CommRequest<int> > > recvSizeReqs (3);
7268  Array<RCP<CommRequest<int> > > recvDataReqs (3);
7269 
7270  // Buffer for nonblocking send of the "send size."
7271  ArrayRCP<size_t> sendDataSize (1);
7272  sendDataSize[0] = myNumRows;
7273 
7274  if (myRank == 0) {
7275  if (debug) {
7276  std::ostringstream os;
7277  os << myRank << ": Post receive-size receives from "
7278  "Procs 1 and 2: tag = " << sizeTag << endl;
7279  *dbg << os.str ();
7280  }
7281  // Process 0: Post receive-size receives from Processes 1 and 2.
7282  recvSizeBufs[0].resize (1);
7283  // Set these three to an invalid value as a flag. If we
7284  // don't get these messages, then the invalid value will
7285  // remain, so we can test for it.
7286  (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7287  recvSizeBufs[1].resize (1);
7288  (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7289  recvSizeBufs[2].resize (1);
7290  (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
7291  if (numProcs > 1) {
7292  recvSizeReqs[1] =
7293  ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
7294  }
7295  if (numProcs > 2) {
7296  recvSizeReqs[2] =
7297  ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
7298  }
7299  }
7300  else if (myRank == 1 || myRank == 2) {
7301  if (debug) {
7302  std::ostringstream os;
7303  os << myRank << ": Post send-size send: size = "
7304  << sendDataSize[0] << ", tag = " << sizeTag << endl;
7305  *dbg << os.str ();
7306  }
7307  // Prime the pipeline by having Processes 1 and 2 start
7308  // their send-size sends. We don't want _all_ the processes
7309  // to start their send-size sends, because that wouldn't be
7310  // memory scalable.
7311  sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7312  }
7313  else {
7314  if (debug) {
7315  std::ostringstream os;
7316  os << myRank << ": Not posting my send-size send yet" << endl;
7317  *dbg << os.str ();
7318  }
7319  }
7320 
7321  //
7322  // Pack my entries, in column-major order.
7323  //
7324  if (debug) {
7325  std::ostringstream os;
7326  os << myRank << ": Pack my entries" << endl;
7327  *dbg << os.str ();
7328  }
7329  ArrayRCP<scalar_type> sendDataBuf;
7330  try {
7331  sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
7332  X.get1dCopy (sendDataBuf (), myNumRows);
7333  }
7334  catch (std::exception& e) {
7335  lclErr = 1;
7336  if (! err.is_null ()) {
7337  std::ostringstream os;
7338  os << "Process " << myRank << ": Attempt to pack my MultiVector "
7339  "entries threw an exception: " << e.what () << endl;
7340  *err << os.str ();
7341  }
7342  }
7343  if (debug) {
7344  std::ostringstream os;
7345  os << myRank << ": Done packing my entries" << endl;
7346  *dbg << os.str ();
7347  }
7348 
7349  //
7350  // Process 1: post send-data send to Process 0.
7351  //
7352  if (myRank == 1) {
7353  if (debug) {
7354  *dbg << myRank << ": Post send-data send: tag = " << dataTag
7355  << endl;
7356  }
7357  sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7358  }
7359 
7360  //
7361  // Process 0: Write my entries.
7362  //
7363  if (myRank == 0) {
7364  if (debug) {
7365  std::ostringstream os;
7366  os << myRank << ": Write my entries" << endl;
7367  *dbg << os.str ();
7368  }
7369 
7370  // Write Process 0's data to the output stream.
7371  // Matrix Market prints dense matrices in column-major order.
7372  const size_t printNumRows = myNumRows;
7373  ArrayView<const scalar_type> printData = sendDataBuf ();
7374  const size_t printStride = printNumRows;
7375  if (static_cast<size_t> (printData.size ()) < printStride * numCols) {
7376  lclErr = 1;
7377  if (! err.is_null ()) {
7378  std::ostringstream os;
7379  os << "Process " << myRank << ": My MultiVector data's size "
7380  << printData.size () << " does not match my local dimensions "
7381  << printStride << " x " << numCols << "." << endl;
7382  *err << os.str ();
7383  }
7384  }
7385  else {
7386  // Matrix Market dense format wants one number per line.
7387  // It wants each complex number as two real numbers (real
7388  // resp. imaginary parts) with a space between.
7389  for (size_t col = 0; col < numCols; ++col) {
7390  for (size_t row = 0; row < printNumRows; ++row) {
7391  if (STS::isComplex) {
7392  out << STS::real (printData[row + col * printStride]) << " "
7393  << STS::imag (printData[row + col * printStride]) << endl;
7394  } else {
7395  out << printData[row + col * printStride] << endl;
7396  }
7397  }
7398  }
7399  }
7400  }
7401 
7402  if (myRank == 0) {
7403  // Wait on receive-size receive from Process 1.
7404  const int recvRank = 1;
7405  const int circBufInd = recvRank % 3;
7406  if (debug) {
7407  std::ostringstream os;
7408  os << myRank << ": Wait on receive-size receive from Process "
7409  << recvRank << endl;
7410  *dbg << os.str ();
7411  }
7412  if (numProcs > 1) {
7413  wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7414 
7415  // We received the number of rows of data. (The data
7416  // come in two columns.)
7417  size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7418  if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7419  lclErr = 1;
7420  if (! err.is_null ()) {
7421  std::ostringstream os;
7422  os << myRank << ": Result of receive-size receive from Process "
7423  << recvRank << " is Teuchos::OrdinalTraits<size_t>::invalid() "
7424  << "= " << Teuchos::OrdinalTraits<size_t>::invalid () << ". "
7425  "This should never happen, and suggests that the receive never "
7426  "got posted. Please report this bug to the Tpetra developers."
7427  << endl;
7428  *err << os.str ();
7429  }
7430 
7431  // If we're going to continue after error, set the
7432  // number of rows to receive to a reasonable size. This
7433  // may cause MPI_ERR_TRUNCATE if the sending process is
7434  // sending more than 0 rows, but that's better than MPI
7435  // overflowing due to the huge positive value that is
7436  // Teuchos::OrdinalTraits<size_t>::invalid().
7437  recvNumRows = 0;
7438  }
7439 
7440  // Post receive-data receive from Process 1.
7441  recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7442  if (debug) {
7443  std::ostringstream os;
7444  os << myRank << ": Post receive-data receive from Process "
7445  << recvRank << ": tag = " << dataTag << ", buffer size = "
7446  << recvDataBufs[circBufInd].size () << endl;
7447  *dbg << os.str ();
7448  }
7449  if (! recvSizeReqs[circBufInd].is_null ()) {
7450  lclErr = 1;
7451  if (! err.is_null ()) {
7452  std::ostringstream os;
7453  os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
7454  "null, before posting the receive-data receive from Process "
7455  << recvRank << ". This should never happen. Please report "
7456  "this bug to the Tpetra developers." << endl;
7457  *err << os.str ();
7458  }
7459  }
7460  recvDataReqs[circBufInd] =
7461  ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7462  recvRank, dataTag, comm);
7463  } // numProcs > 1
7464  }
7465  else if (myRank == 1) {
7466  // Wait on my send-size send.
7467  if (debug) {
7468  std::ostringstream os;
7469  os << myRank << ": Wait on my send-size send" << endl;
7470  *dbg << os.str ();
7471  }
7472  wait<int> (comm, outArg (sendReqSize));
7473  }
7474 
7475  //
7476  // Pipeline loop
7477  //
7478  for (int p = 1; p < numProcs; ++p) {
7479  if (myRank == 0) {
7480  if (p + 2 < numProcs) {
7481  // Post receive-size receive from Process p + 2.
7482  const int recvRank = p + 2;
7483  const int circBufInd = recvRank % 3;
7484  if (debug) {
7485  std::ostringstream os;
7486  os << myRank << ": Post receive-size receive from Process "
7487  << recvRank << ": tag = " << sizeTag << endl;
7488  *dbg << os.str ();
7489  }
7490  if (! recvSizeReqs[circBufInd].is_null ()) {
7491  lclErr = 1;
7492  if (! err.is_null ()) {
7493  std::ostringstream os;
7494  os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
7495  << "null, for the receive-size receive from Process "
7496  << recvRank << "! This may mean that this process never "
7497  << "finished waiting for the receive from Process "
7498  << (recvRank - 3) << "." << endl;
7499  *err << os.str ();
7500  }
7501  }
7502  recvSizeReqs[circBufInd] =
7503  ireceive<int, size_t> (recvSizeBufs[circBufInd],
7504  recvRank, sizeTag, comm);
7505  }
7506 
7507  if (p + 1 < numProcs) {
7508  const int recvRank = p + 1;
7509  const int circBufInd = recvRank % 3;
7510 
7511  // Wait on receive-size receive from Process p + 1.
7512  if (debug) {
7513  std::ostringstream os;
7514  os << myRank << ": Wait on receive-size receive from Process "
7515  << recvRank << endl;
7516  *dbg << os.str ();
7517  }
7518  wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
7519 
7520  // We received the number of rows of data. (The data
7521  // come in two columns.)
7522  size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
7523  if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7524  lclErr = 1;
7525  if (! err.is_null ()) {
7526  std::ostringstream os;
7527  os << myRank << ": Result of receive-size receive from Process "
7528  << recvRank << " is Teuchos::OrdinalTraits<size_t>::invalid() "
7529  << "= " << Teuchos::OrdinalTraits<size_t>::invalid () << ". "
7530  "This should never happen, and suggests that the receive never "
7531  "got posted. Please report this bug to the Tpetra developers."
7532  << endl;
7533  *err << os.str ();
7534  }
7535  // If we're going to continue after error, set the
7536  // number of rows to receive to a reasonable size.
7537  // This may cause MPI_ERR_TRUNCATE if the sending
7538  // process sends more than 0 rows, but that's better
7539  // than MPI overflowing due to the huge positive value
7540  // Teuchos::OrdinalTraits<size_t>::invalid().
7541  recvNumRows = 0;
7542  }
7543 
7544  // Post receive-data receive from Process p + 1.
7545  recvDataBufs[circBufInd].resize (recvNumRows * numCols);
7546  if (debug) {
7547  std::ostringstream os;
7548  os << myRank << ": Post receive-data receive from Process "
7549  << recvRank << ": tag = " << dataTag << ", buffer size = "
7550  << recvDataBufs[circBufInd].size () << endl;
7551  *dbg << os.str ();
7552  }
7553  if (! recvDataReqs[circBufInd].is_null ()) {
7554  lclErr = 1;
7555  if (! err.is_null ()) {
7556  std::ostringstream os;
7557  os << myRank << ": recvDataReqs[" << circBufInd << "] is not "
7558  << "null, for the receive-data receive from Process "
7559  << recvRank << "! This may mean that this process never "
7560  << "finished waiting for the receive from Process "
7561  << (recvRank - 3) << "." << endl;
7562  *err << os.str ();
7563  }
7564  }
7565  recvDataReqs[circBufInd] =
7566  ireceive<int, scalar_type> (recvDataBufs[circBufInd],
7567  recvRank, dataTag, comm);
7568  }
7569 
7570  // Wait on receive-data receive from Process p.
7571  const int recvRank = p;
7572  const int circBufInd = recvRank % 3;
7573  if (debug) {
7574  std::ostringstream os;
7575  os << myRank << ": Wait on receive-data receive from Process "
7576  << recvRank << endl;
7577  *dbg << os.str ();
7578  }
7579  wait<int> (comm, outArg (recvDataReqs[circBufInd]));
7580 
7581  // Write Process p's data. Number of rows lives in
7582  // recvSizeBufs[circBufInd], and the actual data live in
7583  // recvDataBufs[circBufInd]. Do this after posting receives,
7584  // in order to expose overlap of comm. with file I/O.
7585  if (debug) {
7586  std::ostringstream os;
7587  os << myRank << ": Write entries from Process " << recvRank
7588  << endl;
7589  *dbg << os.str () << endl;
7590  }
7591  size_t printNumRows = (recvSizeBufs[circBufInd])[0];
7592  if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7593  lclErr = 1;
7594  if (! err.is_null ()) {
7595  std::ostringstream os;
7596  os << myRank << ": Result of receive-size receive from Process "
7597  << recvRank << " was Teuchos::OrdinalTraits<size_t>::"
7598  "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
7599  << ". This should never happen, and suggests that its "
7600  "receive-size receive was never posted. "
7601  "Please report this bug to the Tpetra developers." << endl;
7602  *err << os.str ();
7603  }
7604  // If we're going to continue after error, set the
7605  // number of rows to print to a reasonable size.
7606  printNumRows = 0;
7607  }
7608  if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
7609  lclErr = 1;
7610  if (! err.is_null ()) {
7611  std::ostringstream os;
7612  os << myRank << ": Result of receive-size receive from Proc "
7613  << recvRank << " was " << printNumRows << " > 0, but "
7614  "recvDataBufs[" << circBufInd << "] is null. This should "
7615  "never happen. Please report this bug to the Tpetra "
7616  "developers." << endl;
7617  *err << os.str ();
7618  }
7619  // If we're going to continue after error, set the
7620  // number of rows to print to a reasonable size.
7621  printNumRows = 0;
7622  }
7623 
7624  // Write the received data to the output stream.
7625  // Matrix Market prints dense matrices in column-major order.
7626  ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
7627  const size_t printStride = printNumRows;
7628  // Matrix Market dense format wants one number per line.
7629  // It wants each complex number as two real numbers (real
7630  // resp. imaginary parts) with a space between.
7631  for (size_t col = 0; col < numCols; ++col) {
7632  for (size_t row = 0; row < printNumRows; ++row) {
7633  if (STS::isComplex) {
7634  out << STS::real (printData[row + col * printStride]) << " "
7635  << STS::imag (printData[row + col * printStride]) << endl;
7636  } else {
7637  out << printData[row + col * printStride] << endl;
7638  }
7639  }
7640  }
7641  }
7642  else if (myRank == p) { // Process p
7643  // Wait on my send-data send.
7644  if (debug) {
7645  std::ostringstream os;
7646  os << myRank << ": Wait on my send-data send" << endl;
7647  *dbg << os.str ();
7648  }
7649  wait<int> (comm, outArg (sendReqData));
7650  }
7651  else if (myRank == p + 1) { // Process p + 1
7652  // Post send-data send to Process 0.
7653  if (debug) {
7654  std::ostringstream os;
7655  os << myRank << ": Post send-data send: tag = " << dataTag
7656  << endl;
7657  *dbg << os.str ();
7658  }
7659  sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
7660  // Wait on my send-size send.
7661  if (debug) {
7662  std::ostringstream os;
7663  os << myRank << ": Wait on my send-size send" << endl;
7664  *dbg << os.str ();
7665  }
7666  wait<int> (comm, outArg (sendReqSize));
7667  }
7668  else if (myRank == p + 2) { // Process p + 2
7669  // Post send-size send to Process 0.
7670  if (debug) {
7671  std::ostringstream os;
7672  os << myRank << ": Post send-size send: size = "
7673  << sendDataSize[0] << ", tag = " << sizeTag << endl;
7674  *dbg << os.str ();
7675  }
7676  sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
7677  }
7678  }
7679 
7680  // Establish global agreement on the error state.
7681  reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
7682  TEUCHOS_TEST_FOR_EXCEPTION(
7683  gblErr == 1, std::runtime_error, "Tpetra::MatrixMarket::writeDense "
7684  "experienced some kind of error and was unable to complete.");
7685 
7686  if (debug) {
7687  dbg->popTab ();
7688  *dbg << myRank << ": writeDenseColumn: Done" << endl;
7689  dbg->popTab ();
7690  }
7691  }
7692 
7693  public:
7694 
7700  static void
7701  writeDense (std::ostream& out,
7702  const Teuchos::RCP<const multivector_type>& X,
7703  const std::string& matrixName,
7704  const std::string& matrixDescription,
7705  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7706  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7707  {
7708  TEUCHOS_TEST_FOR_EXCEPTION(
7709  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
7710  "writeDense: The input MultiVector X is null.");
7711  writeDense (out, *X, matrixName, matrixDescription, err, dbg);
7712  }
7713 
7719  static void
7720  writeDense (std::ostream& out,
7721  const multivector_type& X,
7722  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7723  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7724  {
7725  writeDense (out, X, "", "", err, dbg);
7726  }
7727 
7733  static void
7734  writeDense (std::ostream& out,
7735  const Teuchos::RCP<const multivector_type>& X,
7736  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
7737  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
7738  {
7739  TEUCHOS_TEST_FOR_EXCEPTION(
7740  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
7741  "writeDense: The input MultiVector X is null.");
7742  writeDense (out, *X, "", "", err, dbg);
7743  }
7744 
7764  static void
7765  writeMap (std::ostream& out, const map_type& map, const bool debug=false)
7766  {
7767  Teuchos::RCP<Teuchos::FancyOStream> err =
7768  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
7769  writeMap (out, map, err, debug);
7770  }
7771 
7780  static void
7781  writeMap (std::ostream& out,
7782  const map_type& map,
7783  const Teuchos::RCP<Teuchos::FancyOStream>& err,
7784  const bool debug=false)
7785  {
7786  using Teuchos::Array;
7787  using Teuchos::ArrayRCP;
7788  using Teuchos::ArrayView;
7789  using Teuchos::Comm;
7790  using Teuchos::CommRequest;
7791  using Teuchos::ireceive;
7792  using Teuchos::isend;
7793  using Teuchos::RCP;
7794  using Teuchos::TypeNameTraits;
7795  using Teuchos::wait;
7796  using std::endl;
7797  typedef global_ordinal_type GO;
7798  typedef int pid_type;
7799 
7800  // Treat the Map as a 1-column "multivector." This differs
7801  // from the previous two-column format, in which column 0 held
7802  // the GIDs, and column 1 held the corresponding PIDs. It
7803  // differs because printing that format requires knowing the
7804  // entire first column -- that is, all the GIDs -- in advance.
7805  // Sending messages from each process one at a time saves
7806  // memory, but it means that Process 0 doesn't ever have all
7807  // the GIDs at once.
7808  //
7809  // We pack the entries as ptrdiff_t, since this should be the
7810  // biggest signed built-in integer type that can hold any GO
7811  // or pid_type (= int) quantity without overflow. Test this
7812  // assumption at run time.
7813  typedef ptrdiff_t int_type;
7814  TEUCHOS_TEST_FOR_EXCEPTION(
7815  sizeof (GO) > sizeof (int_type), std::logic_error,
7816  "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
7817  << " is too big for ptrdiff_t. sizeof(GO) = " << sizeof (GO)
7818  << " > sizeof(ptrdiff_t) = " << sizeof (ptrdiff_t) << ".");
7819  TEUCHOS_TEST_FOR_EXCEPTION(
7820  sizeof (pid_type) > sizeof (int_type), std::logic_error,
7821  "The (MPI) process rank type pid_type=" <<
7822  TypeNameTraits<pid_type>::name () << " is too big for ptrdiff_t. "
7823  "sizeof(pid_type) = " << sizeof (pid_type) << " > sizeof(ptrdiff_t)"
7824  " = " << sizeof (ptrdiff_t) << ".");
7825 
7826  const Comm<int>& comm = * (map.getComm ());
7827  const int myRank = comm.getRank ();
7828  const int numProcs = comm.getSize ();
7829 
7830  if (! err.is_null ()) {
7831  err->pushTab ();
7832  }
7833  if (debug) {
7834  std::ostringstream os;
7835  os << myRank << ": writeMap" << endl;
7836  *err << os.str ();
7837  }
7838  if (! err.is_null ()) {
7839  err->pushTab ();
7840  }
7841 
7842  const size_t myNumRows = map.getLocalNumElements ();
7843  // Use a different tag for the "size" messages than for the
7844  // "data" messages, in order to help us debug any mix-ups.
7845  const int sizeTag = 1337;
7846  const int dataTag = 1338;
7847 
7848  // Process 0 pipelines nonblocking receives with file output.
7849  //
7850  // Constraints: