Tpetra parallel linear algebra  Version of the Day
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 
56 #include "Tpetra_ConfigDefs.hpp"
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_SymmetrizingAdder.hpp"
63 #include "Teuchos_MatrixMarket_assignScalar.hpp"
64 #include "Teuchos_MatrixMarket_Banner.hpp"
65 #include "Teuchos_MatrixMarket_CoordDataReader.hpp"
66 #include "Teuchos_MatrixMarket_SetScientific.hpp"
67 
68 #include <algorithm>
69 #include <fstream>
70 #include <iostream>
71 #include <iterator>
72 #include <vector>
73 #include <stdexcept>
74 
75 namespace Tpetra {
105  namespace MatrixMarket {
161  template<class SparseMatrixType>
162  class Reader {
163  public:
165  typedef SparseMatrixType sparse_matrix_type;
166  typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
167 
170  typedef typename SparseMatrixType::scalar_type scalar_type;
173  typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
181  typedef typename SparseMatrixType::global_ordinal_type
184  typedef typename SparseMatrixType::node_type node_type;
185 
187  typedef MultiVector<scalar_type,
190  node_type> multivector_type;
191 
193  typedef Vector<scalar_type,
196  node_type> vector_type;
197 
198  typedef Teuchos::Comm<int> comm_type;
200 
201  // DEPRECATED typedefs for backwards compatibility.
202  typedef Teuchos::RCP<const comm_type> comm_ptr;
203  typedef Teuchos::RCP<const map_type> map_ptr;
204  typedef Teuchos::RCP<node_type> node_ptr;
205 
206  private:
212  typedef Teuchos::ArrayRCP<int>::size_type size_type;
213 
225  static Teuchos::RCP<const map_type>
226  makeRangeMap (const Teuchos::RCP<const comm_type>& pComm,
227  const Teuchos::RCP<node_type>& pNode,
228  const global_ordinal_type numRows)
229  {
230  // Return a conventional, uniformly partitioned, contiguous map.
231  if (pNode.is_null ()) {
232  return rcp (new map_type (static_cast<global_size_t> (numRows),
233  static_cast<global_ordinal_type> (0),
234  pComm, GloballyDistributed));
235  }
236  else {
237  return rcp (new map_type (static_cast<global_size_t> (numRows),
238  static_cast<global_ordinal_type> (0),
239  pComm, GloballyDistributed, pNode));
240  }
241  }
242 
271  static RCP<const map_type>
272  makeRowMap (const RCP<const map_type>& pRowMap,
273  const RCP<const comm_type>& pComm,
274  const RCP<node_type>& pNode,
275  const global_ordinal_type numRows)
276  {
277  // If the caller didn't provide a map, return a conventional,
278  // uniformly partitioned, contiguous map.
279  if (pRowMap.is_null ()) {
280  if (pNode.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  return rcp (new map_type (static_cast<global_size_t> (numRows),
287  static_cast<global_ordinal_type> (0),
288  pComm, GloballyDistributed, pNode));
289  }
290  }
291  else {
292  TEUCHOS_TEST_FOR_EXCEPTION
293  (! pRowMap->isDistributed () && pComm->getSize () > 1,
294  std::invalid_argument, "The specified row map is not distributed, "
295  "but the given communicator includes more than one process (in "
296  "fact, there are " << pComm->getSize () << " processes).");
297  TEUCHOS_TEST_FOR_EXCEPTION
298  (pRowMap->getComm () != pComm, std::invalid_argument,
299  "The specified row Map's communicator (pRowMap->getComm()) "
300  "differs from the given separately supplied communicator pComm.");
301  TEUCHOS_TEST_FOR_EXCEPTION
302  (pRowMap->getNode () != pNode, std::invalid_argument,
303  "The specified row Map's node (pRowMap->getNode()) differs from "
304  "the given separately supplied node pNode.");
305  return pRowMap;
306  }
307  }
308 
323  static Teuchos::RCP<const map_type>
324  makeDomainMap (const Teuchos::RCP<const map_type>& pRangeMap,
325  const global_ordinal_type numRows,
326  const global_ordinal_type numCols)
327  {
328  // Abbreviations so that the map creation call isn't too long.
329  typedef local_ordinal_type LO;
330  typedef global_ordinal_type GO;
331  typedef node_type NT;
332 
333  if (numRows == numCols) {
334  return pRangeMap;
335  } else {
336  return createUniformContigMapWithNode<LO,GO,NT> (numCols,
337  pRangeMap->getComm (),
338  pRangeMap->getNode ());
339  }
340  }
341 
414  static void
415  distribute (ArrayRCP<size_t>& myNumEntriesPerRow,
416  ArrayRCP<size_t>& myRowPtr,
417  ArrayRCP<global_ordinal_type>& myColInd,
418  ArrayRCP<scalar_type>& myValues,
419  const RCP<const map_type>& pRowMap,
420  ArrayRCP<size_t>& numEntriesPerRow,
421  ArrayRCP<size_t>& rowPtr,
422  ArrayRCP<global_ordinal_type>& colInd,
423  ArrayRCP<scalar_type>& values,
424  const bool debug=false)
425  {
426  using Teuchos::as;
427  using Teuchos::Comm;
428  using Teuchos::CommRequest;
429  using Teuchos::RCP;
430  using Teuchos::receive;
431  using Teuchos::send;
432  using std::cerr;
433  using std::endl;
434 
435  const bool extraDebug = false;
436  RCP<const comm_type> pComm = pRowMap->getComm ();
437  const int numProcs = pComm->getSize ();
438  const int myRank = pComm->getRank ();
439  const int rootRank = 0;
440 
441  // Type abbreviations to make the code more concise.
442  typedef global_ordinal_type GO;
443 
444  // List of the global indices of my rows. They may or may
445  // not be contiguous, and the row map need not be one-to-one.
446  ArrayView<const GO> myRows = pRowMap->getNodeElementList();
447  const size_type myNumRows = myRows.size();
448  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<size_t>(myNumRows) !=
449  pRowMap->getNodeNumElements(),
450  std::logic_error,
451  "pRowMap->getNodeElementList().size() = "
452  << myNumRows
453  << " != pRowMap->getNodeNumElements() = "
454  << pRowMap->getNodeNumElements() << ". "
455  "Please report this bug to the Tpetra developers.");
456  TEUCHOS_TEST_FOR_EXCEPTION(myRank == 0 && numEntriesPerRow.size() < myNumRows,
457  std::logic_error,
458  "On Proc 0: numEntriesPerRow.size() = "
459  << numEntriesPerRow.size()
460  << " != pRowMap->getNodeElementList().size() = "
461  << myNumRows << ". Please report this bug to the "
462  "Tpetra developers.");
463 
464  // Space for my proc's number of entries per row. Will be
465  // filled in below.
466  myNumEntriesPerRow = arcp<size_t> (myNumRows);
467 
468  if (myRank != rootRank) {
469  // Tell the root how many rows we have. If we're sending
470  // none, then we don't have anything else to send, nor does
471  // the root have to receive anything else.
472  send (*pComm, myNumRows, rootRank);
473  if (myNumRows != 0) {
474  // Now send my rows' global indices. Hopefully the cast
475  // to int doesn't overflow. This is unlikely, since it
476  // should fit in a LO, even though it is a GO.
477  send (*pComm, static_cast<int> (myNumRows),
478  myRows.getRawPtr(), rootRank);
479 
480  // I (this proc) don't care if my global row indices are
481  // contiguous, though the root proc does (since otherwise
482  // it needs to pack noncontiguous data into contiguous
483  // storage before sending). That's why we don't check
484  // for contiguousness here.
485 
486  // Ask the root process for my part of the array of the
487  // number of entries per row.
488  receive (*pComm, rootRank,
489  static_cast<int> (myNumRows),
490  myNumEntriesPerRow.getRawPtr());
491 
492  // Use the resulting array to figure out how many column
493  // indices and values I should ask from the root process.
494  const local_ordinal_type myNumEntries =
495  std::accumulate (myNumEntriesPerRow.begin(),
496  myNumEntriesPerRow.end(), 0);
497 
498  // Make space for my entries of the sparse matrix. Note
499  // that they don't have to be sorted by row index.
500  // Iterating through all my rows requires computing a
501  // running sum over myNumEntriesPerRow.
502  myColInd = arcp<GO> (myNumEntries);
503  myValues = arcp<scalar_type> (myNumEntries);
504  if (myNumEntries > 0) {
505  // Ask for that many column indices and values, if
506  // there are any.
507  receive (*pComm, rootRank,
508  static_cast<int> (myNumEntries),
509  myColInd.getRawPtr());
510  receive (*pComm, rootRank,
511  static_cast<int> (myNumEntries),
512  myValues.getRawPtr());
513  }
514  } // If I own at least one row
515  } // If I am not the root processor
516  else { // I _am_ the root processor
517  if (debug) {
518  cerr << "-- Proc 0: Copying my data from global arrays" << endl;
519  }
520  // Proc 0 still needs to (allocate, if not done already)
521  // and fill its part of the matrix (my*).
522  for (size_type k = 0; k < myNumRows; ++k) {
523  const GO myCurRow = myRows[k];
524  const local_ordinal_type numEntriesInThisRow = numEntriesPerRow[myCurRow];
525  myNumEntriesPerRow[k] = numEntriesInThisRow;
526  }
527  if (extraDebug && debug) {
528  cerr << "Proc " << pRowMap->getComm ()->getRank ()
529  << ": myNumEntriesPerRow[0.." << (myNumRows-1) << "] = [";
530  for (size_type k = 0; k < myNumRows; ++k) {
531  cerr << myNumEntriesPerRow[k];
532  if (k < myNumRows-1) {
533  cerr << " ";
534  }
535  }
536  cerr << "]" << endl;
537  }
538  // The total number of matrix entries that my proc owns.
539  const local_ordinal_type myNumEntries =
540  std::accumulate (myNumEntriesPerRow.begin(),
541  myNumEntriesPerRow.end(), 0);
542  if (debug) {
543  cerr << "-- Proc 0: I own " << myNumRows << " rows and "
544  << myNumEntries << " entries" << endl;
545  }
546  myColInd = arcp<GO> (myNumEntries);
547  myValues = arcp<scalar_type> (myNumEntries);
548 
549  // Copy Proc 0's part of the matrix into the my* arrays.
550  // It's important that myCurPos be updated _before_ k,
551  // otherwise myCurPos will get the wrong number of entries
552  // per row (it should be for the row in the just-completed
553  // iteration, not for the next iteration's row).
554  local_ordinal_type myCurPos = 0;
555  for (size_type k = 0; k < myNumRows;
556  myCurPos += myNumEntriesPerRow[k], ++k) {
557  const local_ordinal_type curNumEntries = myNumEntriesPerRow[k];
558  const GO myRow = myRows[k];
559  const size_t curPos = rowPtr[myRow];
560  // Only copy if there are entries to copy, in order not
561  // to construct empty ranges for the ArrayRCP views.
562  if (curNumEntries > 0) {
563  ArrayView<GO> colIndView = colInd (curPos, curNumEntries);
564  ArrayView<GO> myColIndView = myColInd (myCurPos, curNumEntries);
565  std::copy (colIndView.begin(), colIndView.end(),
566  myColIndView.begin());
567 
568  ArrayView<scalar_type> valuesView =
569  values (curPos, curNumEntries);
570  ArrayView<scalar_type> myValuesView =
571  myValues (myCurPos, curNumEntries);
572  std::copy (valuesView.begin(), valuesView.end(),
573  myValuesView.begin());
574  }
575  }
576 
577  // Proc 0 processes each other proc p in turn.
578  for (int p = 1; p < numProcs; ++p) {
579  if (debug) {
580  cerr << "-- Proc 0: Processing proc " << p << endl;
581  }
582 
583  size_type theirNumRows = 0;
584  // Ask Proc p how many rows it has. If it doesn't
585  // have any, we can move on to the next proc. This
586  // has to be a standard receive so that we can avoid
587  // the degenerate case of sending zero data.
588  receive (*pComm, p, &theirNumRows);
589  if (debug) {
590  cerr << "-- Proc 0: Proc " << p << " owns "
591  << theirNumRows << " rows" << endl;
592  }
593  if (theirNumRows != 0) {
594  // Ask Proc p which rows it owns. The resulting global
595  // row indices are not guaranteed to be contiguous or
596  // sorted. Global row indices are themselves indices
597  // into the numEntriesPerRow array.
598  ArrayRCP<GO> theirRows = arcp<GO> (theirNumRows);
599  receive (*pComm, p, as<int> (theirNumRows),
600  theirRows.getRawPtr ());
601  // Extra test to make sure that the rows we received
602  // are all sensible. This is a good idea since we are
603  // going to use the global row indices we've received
604  // to index into the numEntriesPerRow array. Better to
605  // catch any bugs here and print a sensible error
606  // message, rather than segfault and print a cryptic
607  // error message.
608  {
609  const global_size_t numRows = pRowMap->getGlobalNumElements ();
610  const GO indexBase = pRowMap->getIndexBase ();
611  bool theirRowsValid = true;
612  for (size_type k = 0; k < theirNumRows; ++k) {
613  if (theirRows[k] < indexBase ||
614  as<global_size_t> (theirRows[k] - indexBase) >= numRows) {
615  theirRowsValid = false;
616  }
617  }
618  if (! theirRowsValid) {
619  TEUCHOS_TEST_FOR_EXCEPTION(
620  ! theirRowsValid, std::logic_error,
621  "Proc " << p << " has at least one invalid row index. "
622  "Here are all of them: " <<
623  Teuchos::toString (theirRows ()) << ". Valid row index "
624  "range (zero-based): [0, " << (numRows - 1) << "].");
625  }
626  }
627 
628  // Perhaps we could save a little work if we check
629  // whether Proc p's row indices are contiguous. That
630  // would make lookups in the global data arrays
631  // faster. For now, we just implement the general
632  // case and don't prematurely optimize. (Remember
633  // that you're making Proc 0 read the whole file, so
634  // you've already lost scalability.)
635 
636  // Compute the number of entries in each of Proc p's
637  // rows. (Proc p will compute its row pointer array
638  // on its own, after it gets the data from Proc 0.)
639  ArrayRCP<size_t> theirNumEntriesPerRow;
640  theirNumEntriesPerRow = arcp<size_t> (theirNumRows);
641  for (size_type k = 0; k < theirNumRows; ++k) {
642  theirNumEntriesPerRow[k] = numEntriesPerRow[theirRows[k]];
643  }
644 
645  // Tell Proc p the number of entries in each of its
646  // rows. Hopefully the cast to int doesn't overflow.
647  // This is unlikely, since it should fit in a LO,
648  // even though it is a GO.
649  send (*pComm, static_cast<int> (theirNumRows),
650  theirNumEntriesPerRow.getRawPtr(), p);
651 
652  // Figure out how many entries Proc p owns.
653  const local_ordinal_type theirNumEntries =
654  std::accumulate (theirNumEntriesPerRow.begin(),
655  theirNumEntriesPerRow.end(), 0);
656 
657  if (debug) {
658  cerr << "-- Proc 0: Proc " << p << " owns "
659  << theirNumEntries << " entries" << endl;
660  }
661 
662  // If there are no entries to send, then we're done
663  // with Proc p.
664  if (theirNumEntries == 0) {
665  continue;
666  }
667 
668  // Construct (views of) proc p's column indices and
669  // values. Later, we might like to optimize for the
670  // (common) contiguous case, for which we don't need to
671  // copy data into separate "their*" arrays (we can just
672  // use contiguous views of the global arrays).
673  ArrayRCP<GO> theirColInd (theirNumEntries);
674  ArrayRCP<scalar_type> theirValues (theirNumEntries);
675  // Copy Proc p's part of the matrix into the their*
676  // arrays. It's important that theirCurPos be updated
677  // _before_ k, otherwise theirCurPos will get the wrong
678  // number of entries per row (it should be for the row
679  // in the just-completed iteration, not for the next
680  // iteration's row).
681  local_ordinal_type theirCurPos = 0;
682  for (size_type k = 0; k < theirNumRows;
683  theirCurPos += theirNumEntriesPerRow[k], k++) {
684  const local_ordinal_type curNumEntries = theirNumEntriesPerRow[k];
685  const GO theirRow = theirRows[k];
686  const local_ordinal_type curPos = rowPtr[theirRow];
687 
688  // Only copy if there are entries to copy, in order
689  // not to construct empty ranges for the ArrayRCP
690  // views.
691  if (curNumEntries > 0) {
692  ArrayView<GO> colIndView =
693  colInd (curPos, curNumEntries);
694  ArrayView<GO> theirColIndView =
695  theirColInd (theirCurPos, curNumEntries);
696  std::copy (colIndView.begin(), colIndView.end(),
697  theirColIndView.begin());
698 
699  ArrayView<scalar_type> valuesView =
700  values (curPos, curNumEntries);
701  ArrayView<scalar_type> theirValuesView =
702  theirValues (theirCurPos, curNumEntries);
703  std::copy (valuesView.begin(), valuesView.end(),
704  theirValuesView.begin());
705  }
706  }
707  // Send Proc p its column indices and values.
708  // Hopefully the cast to int doesn't overflow. This
709  // is unlikely, since it should fit in a LO, even
710  // though it is a GO.
711  send (*pComm, static_cast<int> (theirNumEntries),
712  theirColInd.getRawPtr(), p);
713  send (*pComm, static_cast<int> (theirNumEntries),
714  theirValues.getRawPtr(), p);
715 
716  if (debug) {
717  cerr << "-- Proc 0: Finished with proc " << p << endl;
718  }
719  } // If proc p owns at least one row
720  } // For each proc p not the root proc 0
721  } // If I'm (not) the root proc 0
722 
723  // Invalidate the input data to save space, since we don't
724  // need it anymore.
725  numEntriesPerRow = null;
726  rowPtr = null;
727  colInd = null;
728  values = null;
729 
730  if (debug && myRank == 0) {
731  cerr << "-- Proc 0: About to fill in myRowPtr" << endl;
732  }
733 
734  // Allocate and fill in myRowPtr (the row pointer array for
735  // my rank's rows). We delay this until the end because we
736  // don't need it to compute anything else in distribute().
737  // Each proc can do this work for itself, since it only needs
738  // myNumEntriesPerRow to do so.
739  myRowPtr = arcp<size_t> (myNumRows+1);
740  myRowPtr[0] = 0;
741  for (size_type k = 1; k < myNumRows+1; ++k) {
742  myRowPtr[k] = myRowPtr[k-1] + myNumEntriesPerRow[k-1];
743  }
744  if (extraDebug && debug) {
745  cerr << "Proc " << Teuchos::rank (*(pRowMap->getComm()))
746  << ": myRowPtr[0.." << myNumRows << "] = [";
747  for (size_type k = 0; k < myNumRows+1; ++k) {
748  cerr << myRowPtr[k];
749  if (k < myNumRows) {
750  cerr << " ";
751  }
752  }
753  cerr << "]" << endl << endl;
754  }
755 
756  if (debug && myRank == 0) {
757  cerr << "-- Proc 0: Done with distribute" << endl;
758  }
759  }
760 
774  static Teuchos::RCP<sparse_matrix_type>
775  makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
776  Teuchos::ArrayRCP<size_t>& myRowPtr,
777  Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
778  Teuchos::ArrayRCP<scalar_type>& myValues,
779  const Teuchos::RCP<const map_type>& pRowMap,
780  const Teuchos::RCP<const map_type>& pRangeMap,
781  const Teuchos::RCP<const map_type>& pDomainMap,
782  const bool callFillComplete = true)
783  {
784  using Teuchos::ArrayView;
785  using Teuchos::null;
786  using Teuchos::RCP;
787  using Teuchos::rcp;
788  using std::cerr;
789  using std::endl;
790  // Typedef to make certain type declarations shorter.
791  typedef global_ordinal_type GO;
792 
793  // The row pointer array always has at least one entry, even
794  // if the matrix has zero rows. myNumEntriesPerRow, myColInd,
795  // and myValues would all be empty arrays in that degenerate
796  // case, but the row and domain maps would still be nonnull
797  // (though they would be trivial maps).
798  TEUCHOS_TEST_FOR_EXCEPTION(myRowPtr.is_null(), std::logic_error,
799  "makeMatrix: myRowPtr array is null. "
800  "Please report this bug to the Tpetra developers.");
801  TEUCHOS_TEST_FOR_EXCEPTION(pDomainMap.is_null(), std::logic_error,
802  "makeMatrix: domain map is null. "
803  "Please report this bug to the Tpetra developers.");
804  TEUCHOS_TEST_FOR_EXCEPTION(pRangeMap.is_null(), std::logic_error,
805  "makeMatrix: range map is null. "
806  "Please report this bug to the Tpetra developers.");
807  TEUCHOS_TEST_FOR_EXCEPTION(pRowMap.is_null(), std::logic_error,
808  "makeMatrix: row map is null. "
809  "Please report this bug to the Tpetra developers.");
810 
811  // Construct the CrsMatrix, using the row map, with the
812  // constructor specifying the number of nonzeros for each row.
813  // Create with DynamicProfile, so that the fillComplete() can
814  // do first-touch reallocation (a NUMA (Non-Uniform Memory
815  // Access) optimization on multicore CPUs).
816  RCP<sparse_matrix_type> A =
817  rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
818  DynamicProfile));
819 
820  // List of the global indices of my rows.
821  // They may or may not be contiguous.
822  ArrayView<const GO> myRows = pRowMap->getNodeElementList ();
823  const size_type myNumRows = myRows.size ();
824 
825  // Add this processor's matrix entries to the CrsMatrix.
826  const GO indexBase = pRowMap->getIndexBase ();
827  for (size_type i = 0; i < myNumRows; ++i) {
828  const size_type myCurPos = myRowPtr[i];
829  const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
830  ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
831  ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
832 
833  // Modify the column indices in place to have the right index base.
834  for (size_type k = 0; k < curNumEntries; ++k) {
835  curColInd[k] += indexBase;
836  }
837  // Avoid constructing empty views of ArrayRCP objects.
838  if (curNumEntries > 0) {
839  A->insertGlobalValues (myRows[i], curColInd, curValues);
840  }
841  }
842  // We've entered in all our matrix entries, so we can delete
843  // the original data. This will save memory when we call
844  // fillComplete(), so that we never keep more than two copies
845  // of the matrix's data in memory at once.
846  myNumEntriesPerRow = null;
847  myRowPtr = null;
848  myColInd = null;
849  myValues = null;
850 
851  if (callFillComplete) {
852  A->fillComplete (pDomainMap, pRangeMap);
853  }
854  return A;
855  }
856 
862  static Teuchos::RCP<sparse_matrix_type>
863  makeMatrix (ArrayRCP<size_t>& myNumEntriesPerRow,
864  ArrayRCP<size_t>& myRowPtr,
865  ArrayRCP<global_ordinal_type>& myColInd,
866  ArrayRCP<scalar_type>& myValues,
867  const Teuchos::RCP<const map_type>& pRowMap,
868  const Teuchos::RCP<const map_type>& pRangeMap,
869  const Teuchos::RCP<const map_type>& pDomainMap,
870  const RCP<Teuchos::ParameterList>& constructorParams,
871  const RCP<Teuchos::ParameterList>& fillCompleteParams)
872  {
873  using Teuchos::RCP;
874  using Teuchos::rcp;
875  using std::cerr;
876  using std::endl;
877  // Typedef to make certain type declarations shorter.
878  typedef global_ordinal_type GO;
879 
880  // The row pointer array always has at least one entry, even
881  // if the matrix has zero rows. myNumEntriesPerRow, myColInd,
882  // and myValues would all be empty arrays in that degenerate
883  // case, but the row and domain maps would still be nonnull
884  // (though they would be trivial maps).
885  TEUCHOS_TEST_FOR_EXCEPTION(
886  myRowPtr.is_null(), std::logic_error,
887  "makeMatrix: myRowPtr array is null. "
888  "Please report this bug to the Tpetra developers.");
889  TEUCHOS_TEST_FOR_EXCEPTION(
890  pDomainMap.is_null(), std::logic_error,
891  "makeMatrix: domain map is null. "
892  "Please report this bug to the Tpetra developers.");
893  TEUCHOS_TEST_FOR_EXCEPTION(
894  pRangeMap.is_null(), std::logic_error,
895  "makeMatrix: range map is null. "
896  "Please report this bug to the Tpetra developers.");
897  TEUCHOS_TEST_FOR_EXCEPTION(
898  pRowMap.is_null(), std::logic_error,
899  "makeMatrix: row map is null. "
900  "Please report this bug to the Tpetra developers.");
901 
902  // Construct the CrsMatrix, using the row map, with the
903  // constructor specifying the number of nonzeros for each row.
904  // Create with DynamicProfile, so that the fillComplete() can
905  // do first-touch reallocation (a NUMA (Non-Uniform Memory
906  // Access) optimization on multicore CPUs).
907  RCP<sparse_matrix_type> A =
908  rcp (new sparse_matrix_type (pRowMap, myNumEntriesPerRow,
909  DynamicProfile, constructorParams));
910 
911  // List of the global indices of my rows.
912  // They may or may not be contiguous.
913  ArrayView<const GO> myRows = pRowMap->getNodeElementList();
914  const size_type myNumRows = myRows.size();
915 
916  // Add this processor's matrix entries to the CrsMatrix.
917  const GO indexBase = pRowMap->getIndexBase ();
918  for (size_type i = 0; i < myNumRows; ++i) {
919  const size_type myCurPos = myRowPtr[i];
920  const local_ordinal_type curNumEntries = myNumEntriesPerRow[i];
921  ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
922  ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
923 
924  // Modify the column indices in place to have the right index base.
925  for (size_type k = 0; k < curNumEntries; ++k) {
926  curColInd[k] += indexBase;
927  }
928  if (curNumEntries > 0) {
929  A->insertGlobalValues (myRows[i], curColInd, curValues);
930  }
931  }
932  // We've entered in all our matrix entries, so we can delete
933  // the original data. This will save memory when we call
934  // fillComplete(), so that we never keep more than two copies
935  // of the matrix's data in memory at once.
936  myNumEntriesPerRow = null;
937  myRowPtr = null;
938  myColInd = null;
939  myValues = null;
940 
941  A->fillComplete (pDomainMap, pRangeMap, fillCompleteParams);
942  return A;
943  }
944 
949  static Teuchos::RCP<sparse_matrix_type>
950  makeMatrix (Teuchos::ArrayRCP<size_t>& myNumEntriesPerRow,
951  Teuchos::ArrayRCP<size_t>& myRowPtr,
952  Teuchos::ArrayRCP<global_ordinal_type>& myColInd,
953  Teuchos::ArrayRCP<scalar_type>& myValues,
954  const Teuchos::RCP<const map_type>& rowMap,
955  Teuchos::RCP<const map_type>& colMap,
956  const Teuchos::RCP<const map_type>& domainMap,
957  const Teuchos::RCP<const map_type>& rangeMap,
958  const bool callFillComplete = true)
959  {
960  using Teuchos::ArrayView;
961  using Teuchos::as;
962  using Teuchos::null;
963  using Teuchos::RCP;
964  using Teuchos::rcp;
965  typedef global_ordinal_type GO;
966  typedef typename ArrayView<const GO>::size_type size_type;
967 
968  // Construct the CrsMatrix.
969  //
970  // Create with DynamicProfile, so that the fillComplete() can
971  // do first-touch reallocation.
972  RCP<sparse_matrix_type> A; // the matrix to return.
973  if (colMap.is_null ()) { // the user didn't provide a column Map
974  A = rcp (new sparse_matrix_type (rowMap, myNumEntriesPerRow, DynamicProfile));
975  } else { // the user provided a column Map
976  A = rcp (new sparse_matrix_type (rowMap, colMap, myNumEntriesPerRow, DynamicProfile));
977  }
978 
979  // List of the global indices of my rows.
980  // They may or may not be contiguous.
981  ArrayView<const GO> myRows = rowMap->getNodeElementList ();
982  const size_type myNumRows = myRows.size ();
983 
984  // Add this process' matrix entries to the CrsMatrix.
985  const GO indexBase = rowMap->getIndexBase ();
986  for (size_type i = 0; i < myNumRows; ++i) {
987  const size_type myCurPos = myRowPtr[i];
988  const size_type curNumEntries = as<size_type> (myNumEntriesPerRow[i]);
989  ArrayView<GO> curColInd = myColInd.view (myCurPos, curNumEntries);
990  ArrayView<scalar_type> curValues = myValues.view (myCurPos, curNumEntries);
991 
992  // Modify the column indices in place to have the right index base.
993  for (size_type k = 0; k < curNumEntries; ++k) {
994  curColInd[k] += indexBase;
995  }
996  if (curNumEntries > 0) {
997  A->insertGlobalValues (myRows[i], curColInd, curValues);
998  }
999  }
1000  // We've entered in all our matrix entries, so we can delete
1001  // the original data. This will save memory when we call
1002  // fillComplete(), so that we never keep more than two copies
1003  // of the matrix's data in memory at once.
1004  myNumEntriesPerRow = null;
1005  myRowPtr = null;
1006  myColInd = null;
1007  myValues = null;
1008 
1009  if (callFillComplete) {
1010  A->fillComplete (domainMap, rangeMap);
1011  if (colMap.is_null ()) {
1012  colMap = A->getColMap ();
1013  }
1014  }
1015  return A;
1016  }
1017 
1018  private:
1019 
1036  static RCP<const Teuchos::MatrixMarket::Banner>
1037  readBanner (std::istream& in,
1038  size_t& lineNumber,
1039  const bool tolerant=false,
1040  const bool debug=false)
1041  {
1042  using Teuchos::MatrixMarket::Banner;
1043  using std::cerr;
1044  using std::endl;
1045  typedef Teuchos::ScalarTraits<scalar_type> STS;
1046 
1047  RCP<Banner> pBanner; // On output, if successful: the read-in Banner.
1048  std::string line; // If read from stream successful: the Banner line
1049 
1050  // Try to read a line from the input stream.
1051  const bool readFailed = ! getline(in, line);
1052  TEUCHOS_TEST_FOR_EXCEPTION(readFailed, std::invalid_argument,
1053  "Failed to get Matrix Market banner line from input.");
1054 
1055  // We read a line from the input stream.
1056  lineNumber++;
1057 
1058  // Assume that the line we found is the Banner line.
1059  try {
1060  pBanner = rcp (new Banner (line, tolerant));
1061  } catch (std::exception& e) {
1062  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
1063  "Matrix Market banner line contains syntax error(s): "
1064  << e.what());
1065  }
1066  TEUCHOS_TEST_FOR_EXCEPTION(pBanner->objectType() != "matrix",
1067  std::invalid_argument, "The Matrix Market file does not contain "
1068  "matrix data. Its Banner (first) line says that its object type is \""
1069  << pBanner->matrixType() << "\", rather than the required \"matrix\".");
1070 
1071  // Validate the data type of the matrix, with respect to the
1072  // Scalar type of the CrsMatrix entries.
1073  TEUCHOS_TEST_FOR_EXCEPTION(
1074  ! STS::isComplex && pBanner->dataType() == "complex",
1075  std::invalid_argument,
1076  "The Matrix Market file contains complex-valued data, but you are "
1077  "trying to read it into a matrix containing entries of the real-"
1078  "valued Scalar type \""
1079  << Teuchos::TypeNameTraits<scalar_type>::name() << "\".");
1080  TEUCHOS_TEST_FOR_EXCEPTION(
1081  pBanner->dataType() != "real" &&
1082  pBanner->dataType() != "complex" &&
1083  pBanner->dataType() != "integer",
1084  std::invalid_argument,
1085  "When reading Matrix Market data into a Tpetra::CrsMatrix, the "
1086  "Matrix Market file may not contain a \"pattern\" matrix. A "
1087  "pattern matrix is really just a graph with no weights. It "
1088  "should be stored in a CrsGraph, not a CrsMatrix.");
1089 
1090  return pBanner;
1091  }
1092 
1115  static Tuple<global_ordinal_type, 3>
1116  readCoordDims (std::istream& in,
1117  size_t& lineNumber,
1118  const Teuchos::RCP<const Teuchos::MatrixMarket::Banner>& pBanner,
1119  const Teuchos::RCP<const comm_type>& pComm,
1120  const bool tolerant = false,
1121  const bool debug = false)
1122  {
1123  using Teuchos::MatrixMarket::readCoordinateDimensions;
1124  using Teuchos::Tuple;
1125 
1126  // Packed coordinate matrix dimensions (numRows, numCols,
1127  // numNonzeros); computed on Rank 0 and broadcasted to all MPI
1128  // ranks.
1129  Tuple<global_ordinal_type, 3> dims;
1130 
1131  // Read in the coordinate matrix dimensions from the input
1132  // stream. "success" tells us whether reading in the
1133  // coordinate matrix dimensions succeeded ("Guilty unless
1134  // proven innocent").
1135  bool success = false;
1136  if (pComm->getRank() == 0) {
1137  TEUCHOS_TEST_FOR_EXCEPTION(pBanner->matrixType() != "coordinate",
1138  std::invalid_argument, "The Tpetra::CrsMatrix Matrix Market reader "
1139  "only accepts \"coordinate\" (sparse) matrix data.");
1140  // Unpacked coordinate matrix dimensions
1141  global_ordinal_type numRows, numCols, numNonzeros;
1142  // Only MPI Rank 0 reads from the input stream
1143  success = readCoordinateDimensions (in, numRows, numCols,
1144  numNonzeros, lineNumber,
1145  tolerant);
1146  // Pack up the data into a Tuple so we can send them with
1147  // one broadcast instead of three.
1148  dims[0] = numRows;
1149  dims[1] = numCols;
1150  dims[2] = numNonzeros;
1151  }
1152  // Only Rank 0 did the reading, so it decides success.
1153  //
1154  // FIXME (mfh 02 Feb 2011) Teuchos::broadcast doesn't know how
1155  // to send bools. For now, we convert to/from int instead,
1156  // using the usual "true is 1, false is 0" encoding.
1157  {
1158  int the_success = success ? 1 : 0; // only matters on MPI Rank 0
1159  Teuchos::broadcast (*pComm, 0, &the_success);
1160  success = (the_success == 1);
1161  }
1162  if (success) {
1163  // Broadcast (numRows, numCols, numNonzeros) from Rank 0
1164  // to all the other MPI ranks.
1165  Teuchos::broadcast (*pComm, 0, dims);
1166  }
1167  else {
1168  // Perhaps in tolerant mode, we could set all the
1169  // dimensions to zero for now, and deduce correct
1170  // dimensions by reading all of the file's entries and
1171  // computing the max(row index) and max(column index).
1172  // However, for now we just error out in that case.
1173  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
1174  "Error reading Matrix Market sparse matrix: failed to read "
1175  "coordinate matrix dimensions.");
1176  }
1177  return dims;
1178  }
1179 
1190  typedef Teuchos::MatrixMarket::SymmetrizingAdder<Teuchos::MatrixMarket::Raw::Adder<scalar_type, global_ordinal_type> > adder_type;
1191 
1217  static RCP<adder_type>
1218  makeAdder (const Teuchos::RCP<const 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,
1226  global_ordinal_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 
1238  public:
1262  static Teuchos::RCP<sparse_matrix_type>
1263  readSparseFile (const std::string& filename,
1264  const RCP<const Comm<int> >& pComm,
1265  const bool callFillComplete=true,
1266  const bool tolerant=false,
1267  const bool debug=false)
1268  {
1269  return readSparseFile (filename, pComm, Teuchos::null, callFillComplete, tolerant, debug);
1270  }
1271 
1273  static Teuchos::RCP<sparse_matrix_type>
1274  readSparseFile (const std::string& filename,
1275  const RCP<const Comm<int> >& pComm,
1276  const RCP<node_type>& pNode,
1277  const bool callFillComplete=true,
1278  const bool tolerant=false,
1279  const bool debug=false)
1280  {
1281  const int myRank = pComm->getRank ();
1282  std::ifstream in;
1283 
1284  // Only open the file on Rank 0.
1285  if (myRank == 0) {
1286  in.open (filename.c_str ());
1287  }
1288  // FIXME (mfh 16 Jun 2015) Do a broadcast to make sure that
1289  // opening the file succeeded, before continuing. That will
1290  // avoid hangs if the read doesn't work. On the other hand,
1291  // readSparse could do that too, by checking the status of the
1292  // std::ostream.
1293 
1294  return readSparse (in, pComm, pNode, callFillComplete, tolerant, debug);
1295  // We can rely on the destructor of the input stream to close
1296  // the file on scope exit, even if readSparse() throws an
1297  // exception.
1298  }
1299 
1328  static Teuchos::RCP<sparse_matrix_type>
1329  readSparseFile (const std::string& filename,
1330  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1331  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1332  const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1333  const bool tolerant=false,
1334  const bool debug=false)
1335  {
1336  return readSparseFile (filename, pComm, Teuchos::null,
1337  constructorParams, fillCompleteParams,
1338  tolerant, debug);
1339  }
1340 
1342  static Teuchos::RCP<sparse_matrix_type>
1343  readSparseFile (const std::string& filename,
1344  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1345  const Teuchos::RCP<node_type>& pNode,
1346  const Teuchos::RCP<Teuchos::ParameterList>& constructorParams,
1347  const Teuchos::RCP<Teuchos::ParameterList>& fillCompleteParams,
1348  const bool tolerant=false,
1349  const bool debug=false)
1350  {
1351  std::ifstream in;
1352  if (pComm->getRank () == 0) { // only open on Process 0
1353  in.open (filename.c_str ());
1354  }
1355  return readSparse (in, pComm, pNode, constructorParams,
1356  fillCompleteParams, tolerant, debug);
1357  }
1358 
1396  static Teuchos::RCP<sparse_matrix_type>
1397  readSparseFile (const std::string& filename,
1398  const RCP<const map_type>& rowMap,
1399  RCP<const map_type>& colMap,
1400  const RCP<const map_type>& domainMap,
1401  const RCP<const map_type>& rangeMap,
1402  const bool callFillComplete=true,
1403  const bool tolerant=false,
1404  const bool debug=false)
1405  {
1406  using Teuchos::broadcast;
1407  using Teuchos::outArg;
1408  TEUCHOS_TEST_FOR_EXCEPTION(
1409  rowMap.is_null (), std::invalid_argument,
1410  "Row Map must be nonnull.");
1411 
1412  RCP<const Comm<int> > comm = rowMap->getComm ();
1413  const int myRank = comm->getRank ();
1414 
1415  // Only open the file on Process 0. Test carefully to make
1416  // sure that the file opened successfully (and broadcast that
1417  // result to all processes to prevent a hang on exception
1418  // throw), since it's a common mistake to misspell a filename.
1419  std::ifstream in;
1420  int opened = 0;
1421  if (myRank == 0) {
1422  try {
1423  in.open (filename.c_str ());
1424  opened = 1;
1425  }
1426  catch (...) {
1427  opened = 0;
1428  }
1429  }
1430  broadcast<int, int> (*comm, 0, outArg (opened));
1431  TEUCHOS_TEST_FOR_EXCEPTION(
1432  opened == 0, std::runtime_error,
1433  "readSparseFile: Failed to open file \"" << filename << "\" on "
1434  "Process 0.");
1435  return readSparse (in, rowMap, colMap, domainMap, rangeMap,
1436  callFillComplete, tolerant, debug);
1437  }
1438 
1464  static Teuchos::RCP<sparse_matrix_type>
1465  readSparse (std::istream& in,
1466  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1467  const bool callFillComplete=true,
1468  const bool tolerant=false,
1469  const bool debug=false)
1470  {
1471  return readSparse (in, pComm, Teuchos::null, callFillComplete, tolerant, debug);
1472  }
1473 
1475  static Teuchos::RCP<sparse_matrix_type>
1476  readSparse (std::istream& in,
1477  const Teuchos::RCP<const Teuchos::Comm<int> >& pComm,
1478  const Teuchos::RCP<node_type>& pNode,
1479  const bool callFillComplete=true,
1480  const bool tolerant=false,
1481  const bool debug=false)
1482  {
1483  using Teuchos::MatrixMarket::Banner;
1484  using Teuchos::arcp;
1485  using Teuchos::ArrayRCP;
1486  using Teuchos::broadcast;
1487  using Teuchos::null;
1488  using Teuchos::ptr;
1489  using Teuchos::RCP;
1490  using Teuchos::REDUCE_MAX;
1491  using Teuchos::reduceAll;
1492  using Teuchos::Tuple;
1493  using std::cerr;
1494  using std::endl;
1495  typedef Teuchos::ScalarTraits<scalar_type> STS;
1496 
1497  const bool extraDebug = false;
1498  const int myRank = pComm->getRank ();
1499  const int rootRank = 0;
1500 
1501  // Current line number in the input stream. Various calls
1502  // will modify this depending on the number of lines that are
1503  // read from the input stream. Only Rank 0 modifies this.
1504  size_t lineNumber = 1;
1505 
1506  if (debug && myRank == rootRank) {
1507  cerr << "Matrix Market reader: readSparse:" << endl
1508  << "-- Reading banner line" << endl;
1509  }
1510 
1511  // The "Banner" tells you whether the input stream represents
1512  // a sparse matrix, the symmetry type of the matrix, and the
1513  // type of the data it contains.
1514  //
1515  // pBanner will only be nonnull on MPI Rank 0. It will be
1516  // null on all other MPI processes.
1517  RCP<const Banner> pBanner;
1518  {
1519  // We read and validate the Banner on Proc 0, but broadcast
1520  // the validation result to all processes.
1521  // Teuchos::broadcast doesn't currently work with bool, so
1522  // we use int (true -> 1, false -> 0).
1523  int bannerIsCorrect = 1;
1524  std::ostringstream errMsg;
1525 
1526  if (myRank == rootRank) {
1527  // Read the Banner line from the input stream.
1528  try {
1529  pBanner = readBanner (in, lineNumber, tolerant, debug);
1530  }
1531  catch (std::exception& e) {
1532  errMsg << "Attempt to read the Matrix Market file's Banner line "
1533  "threw an exception: " << e.what();
1534  bannerIsCorrect = 0;
1535  }
1536 
1537  if (bannerIsCorrect) {
1538  // Validate the Banner for the case of a sparse matrix.
1539  // We validate on Proc 0, since it reads the Banner.
1540 
1541  // In intolerant mode, the matrix type must be "coordinate".
1542  if (! tolerant && pBanner->matrixType() != "coordinate") {
1543  bannerIsCorrect = 0;
1544  errMsg << "The Matrix Market input file must contain a "
1545  "\"coordinate\"-format sparse matrix in order to create a "
1546  "Tpetra::CrsMatrix object from it, but the file's matrix "
1547  "type is \"" << pBanner->matrixType() << "\" instead.";
1548  }
1549  // In tolerant mode, we allow the matrix type to be
1550  // anything other than "array" (which would mean that
1551  // the file contains a dense matrix).
1552  if (tolerant && pBanner->matrixType() == "array") {
1553  bannerIsCorrect = 0;
1554  errMsg << "Matrix Market file must contain a \"coordinate\"-"
1555  "format sparse matrix in order to create a Tpetra::CrsMatrix "
1556  "object from it, but the file's matrix type is \"array\" "
1557  "instead. That probably means the file contains dense matrix "
1558  "data.";
1559  }
1560  }
1561  } // Proc 0: Done reading the Banner, hopefully successfully.
1562 
1563  // Broadcast from Proc 0 whether the Banner was read correctly.
1564  broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
1565 
1566  // If the Banner is invalid, all processes throw an
1567  // exception. Only Proc 0 gets the exception message, but
1568  // that's OK, since the main point is to "stop the world"
1569  // (rather than throw an exception on one process and leave
1570  // the others hanging).
1571  TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
1572  std::invalid_argument, errMsg.str ());
1573  } // Done reading the Banner line and broadcasting success.
1574  if (debug && myRank == rootRank) {
1575  cerr << "-- Reading dimensions line" << endl;
1576  }
1577 
1578  // Read the matrix dimensions from the Matrix Market metadata.
1579  // dims = (numRows, numCols, numEntries). Proc 0 does the
1580  // reading, but it broadcasts the results to all MPI
1581  // processes. Thus, readCoordDims() is a collective
1582  // operation. It does a collective check for correctness too.
1583  Tuple<global_ordinal_type, 3> dims =
1584  readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
1585 
1586  if (debug && myRank == rootRank) {
1587  cerr << "-- Making Adder for collecting matrix data" << endl;
1588  }
1589 
1590  // "Adder" object for collecting all the sparse matrix entries
1591  // from the input stream. This is only nonnull on Proc 0.
1592  RCP<adder_type> pAdder =
1593  makeAdder (pComm, pBanner, dims, tolerant, debug);
1594 
1595  if (debug && myRank == rootRank) {
1596  cerr << "-- Reading matrix data" << endl;
1597  }
1598  //
1599  // Read the matrix entries from the input stream on Proc 0.
1600  //
1601  {
1602  // We use readSuccess to broadcast the results of the read
1603  // (succeeded or not) to all MPI processes. Since
1604  // Teuchos::broadcast doesn't currently know how to send
1605  // bools, we convert to int (true -> 1, false -> 0).
1606  int readSuccess = 1;
1607  std::ostringstream errMsg; // Exception message (only valid on Proc 0)
1608  if (myRank == rootRank) {
1609  try {
1610  // Reader for "coordinate" format sparse matrix data.
1611  typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
1612  global_ordinal_type, scalar_type, STS::isComplex> reader_type;
1613  reader_type reader (pAdder);
1614 
1615  // Read the sparse matrix entries.
1616  std::pair<bool, std::vector<size_t> > results =
1617  reader.read (in, lineNumber, tolerant, debug);
1618  readSuccess = results.first ? 1 : 0;
1619  }
1620  catch (std::exception& e) {
1621  readSuccess = 0;
1622  errMsg << e.what();
1623  }
1624  }
1625  broadcast (*pComm, rootRank, ptr (&readSuccess));
1626 
1627  // It would be nice to add a "verbose" flag, so that in
1628  // tolerant mode, we could log any bad line number(s) on
1629  // Proc 0. For now, we just throw if the read fails to
1630  // succeed.
1631  //
1632  // Question: If we're in tolerant mode, and if the read did
1633  // not succeed, should we attempt to call fillComplete()?
1634  TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
1635  "Failed to read the Matrix Market sparse matrix file: "
1636  << errMsg.str());
1637  } // Done reading the matrix entries (stored on Proc 0 for now)
1638 
1639  if (debug && myRank == rootRank) {
1640  cerr << "-- Successfully read the Matrix Market data" << endl;
1641  }
1642 
1643  // In tolerant mode, we need to rebroadcast the matrix
1644  // dimensions, since they may be different after reading the
1645  // actual matrix data. We only need to broadcast the number
1646  // of rows and columns. Only Rank 0 needs to know the actual
1647  // global number of entries, since (a) we need to merge
1648  // duplicates on Rank 0 first anyway, and (b) when we
1649  // distribute the entries, each rank other than Rank 0 will
1650  // only need to know how many entries it owns, not the total
1651  // number of entries.
1652  if (tolerant) {
1653  if (debug && myRank == rootRank) {
1654  cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
1655  << endl
1656  << "----- Dimensions before: "
1657  << dims[0] << " x " << dims[1]
1658  << endl;
1659  }
1660  // Packed coordinate matrix dimensions (numRows, numCols).
1661  Tuple<global_ordinal_type, 2> updatedDims;
1662  if (myRank == rootRank) {
1663  // If one or more bottom rows of the matrix contain no
1664  // entries, then the Adder will report that the number
1665  // of rows is less than that specified in the
1666  // metadata. We allow this case, and favor the
1667  // metadata so that the zero row(s) will be included.
1668  updatedDims[0] =
1669  std::max (dims[0], pAdder->getAdder()->numRows());
1670  updatedDims[1] = pAdder->getAdder()->numCols();
1671  }
1672  broadcast (*pComm, rootRank, updatedDims);
1673  dims[0] = updatedDims[0];
1674  dims[1] = updatedDims[1];
1675  if (debug && myRank == rootRank) {
1676  cerr << "----- Dimensions after: " << dims[0] << " x "
1677  << dims[1] << endl;
1678  }
1679  }
1680  else {
1681  // In strict mode, we require that the matrix's metadata and
1682  // its actual data agree, at least somewhat. In particular,
1683  // the number of rows must agree, since otherwise we cannot
1684  // distribute the matrix correctly.
1685 
1686  // Teuchos::broadcast() doesn't know how to broadcast bools,
1687  // so we use an int with the standard 1 == true, 0 == false
1688  // encoding.
1689  int dimsMatch = 1;
1690  if (myRank == rootRank) {
1691  // If one or more bottom rows of the matrix contain no
1692  // entries, then the Adder will report that the number of
1693  // rows is less than that specified in the metadata. We
1694  // allow this case, and favor the metadata, but do not
1695  // allow the Adder to think there are more rows in the
1696  // matrix than the metadata says.
1697  if (dims[0] < pAdder->getAdder ()->numRows ()) {
1698  dimsMatch = 0;
1699  }
1700  }
1701  broadcast (*pComm, 0, ptr (&dimsMatch));
1702  if (dimsMatch == 0) {
1703  // We're in an error state anyway, so we might as well
1704  // work a little harder to print an informative error
1705  // message.
1706  //
1707  // Broadcast the Adder's idea of the matrix dimensions
1708  // from Proc 0 to all processes.
1709  Tuple<global_ordinal_type, 2> addersDims;
1710  if (myRank == rootRank) {
1711  addersDims[0] = pAdder->getAdder()->numRows();
1712  addersDims[1] = pAdder->getAdder()->numCols();
1713  }
1714  broadcast (*pComm, 0, addersDims);
1715  TEUCHOS_TEST_FOR_EXCEPTION(
1716  dimsMatch == 0, std::runtime_error,
1717  "The matrix metadata says that the matrix is " << dims[0] << " x "
1718  << dims[1] << ", but the actual data says that the matrix is "
1719  << addersDims[0] << " x " << addersDims[1] << ". That means the "
1720  "data includes more rows than reported in the metadata. This "
1721  "is not allowed when parsing in strict mode. Parse the matrix in "
1722  "tolerant mode to ignore the metadata when it disagrees with the "
1723  "data.");
1724  }
1725  } // Matrix dimensions (# rows, # cols, # entries) agree.
1726 
1727  if (debug && myRank == rootRank) {
1728  cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
1729  }
1730 
1731  // Now that we've read in all the matrix entries from the
1732  // input stream into the adder on Proc 0, post-process them
1733  // into CSR format (still on Proc 0). This will facilitate
1734  // distributing them to all the processors.
1735  //
1736  // These arrays represent the global matrix data as a CSR
1737  // matrix (with numEntriesPerRow as redundant but convenient
1738  // metadata, since it's computable from rowPtr and vice
1739  // versa). They are valid only on Proc 0.
1740  ArrayRCP<size_t> numEntriesPerRow;
1741  ArrayRCP<size_t> rowPtr;
1742  ArrayRCP<global_ordinal_type> colInd;
1743  ArrayRCP<scalar_type> values;
1744 
1745  // Proc 0 first merges duplicate entries, and then converts
1746  // the coordinate-format matrix data to CSR.
1747  {
1748  int mergeAndConvertSucceeded = 1;
1749  std::ostringstream errMsg;
1750 
1751  if (myRank == rootRank) {
1752  try {
1753  typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
1754  global_ordinal_type> element_type;
1755 
1756  // Number of rows in the matrix. If we are in tolerant
1757  // mode, we've already synchronized dims with the actual
1758  // matrix data. If in strict mode, we should use dims
1759  // (as read from the file's metadata) rather than the
1760  // matrix data to determine the dimensions. (The matrix
1761  // data will claim fewer rows than the metadata, if one
1762  // or more rows have no entries stored in the file.)
1763  const size_type numRows = dims[0];
1764 
1765  // Additively merge duplicate matrix entries.
1766  pAdder->getAdder()->merge ();
1767 
1768  // Get a temporary const view of the merged matrix entries.
1769  const std::vector<element_type>& entries =
1770  pAdder->getAdder()->getEntries();
1771 
1772  // Number of matrix entries (after merging).
1773  const size_t numEntries = (size_t)entries.size();
1774 
1775  if (debug) {
1776  cerr << "----- Proc 0: Matrix has numRows=" << numRows
1777  << " rows and numEntries=" << numEntries
1778  << " entries." << endl;
1779  }
1780 
1781  // Make space for the CSR matrix data. Converting to
1782  // CSR is easier if we fill numEntriesPerRow with zeros
1783  // at first.
1784  numEntriesPerRow = arcp<size_t> (numRows);
1785  std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
1786  rowPtr = arcp<size_t> (numRows+1);
1787  std::fill (rowPtr.begin(), rowPtr.end(), 0);
1788  colInd = arcp<global_ordinal_type> (numEntries);
1789  values = arcp<scalar_type> (numEntries);
1790 
1791  // Convert from array-of-structs coordinate format to CSR
1792  // (compressed sparse row) format.
1793  global_ordinal_type prvRow = 0;
1794  size_t curPos = 0;
1795  rowPtr[0] = 0;
1796  for (curPos = 0; curPos < numEntries; ++curPos) {
1797  const element_type& curEntry = entries[curPos];
1798  const global_ordinal_type curRow = curEntry.rowIndex();
1799  TEUCHOS_TEST_FOR_EXCEPTION(
1800  curRow < prvRow, std::logic_error,
1801  "Row indices are out of order, even though they are supposed "
1802  "to be sorted. curRow = " << curRow << ", prvRow = "
1803  << prvRow << ", at curPos = " << curPos << ". Please report "
1804  "this bug to the Tpetra developers.");
1805  if (curRow > prvRow) {
1806  for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
1807  rowPtr[r] = curPos;
1808  }
1809  prvRow = curRow;
1810  }
1811  numEntriesPerRow[curRow]++;
1812  colInd[curPos] = curEntry.colIndex();
1813  values[curPos] = curEntry.value();
1814  }
1815  // rowPtr has one more entry than numEntriesPerRow. The
1816  // last entry of rowPtr is the number of entries in
1817  // colInd and values.
1818  rowPtr[numRows] = numEntries;
1819  } // Finished conversion to CSR format
1820  catch (std::exception& e) {
1821  mergeAndConvertSucceeded = 0;
1822  errMsg << "Failed to merge sparse matrix entries and convert to "
1823  "CSR format: " << e.what();
1824  }
1825 
1826  if (debug && mergeAndConvertSucceeded) {
1827  // Number of rows in the matrix.
1828  const size_type numRows = dims[0];
1829  const size_type maxToDisplay = 100;
1830 
1831  cerr << "----- Proc 0: numEntriesPerRow[0.."
1832  << (numEntriesPerRow.size()-1) << "] ";
1833  if (numRows > maxToDisplay) {
1834  cerr << "(only showing first and last few entries) ";
1835  }
1836  cerr << "= [";
1837  if (numRows > 0) {
1838  if (numRows > maxToDisplay) {
1839  for (size_type k = 0; k < 2; ++k) {
1840  cerr << numEntriesPerRow[k] << " ";
1841  }
1842  cerr << "... ";
1843  for (size_type k = numRows-2; k < numRows-1; ++k) {
1844  cerr << numEntriesPerRow[k] << " ";
1845  }
1846  }
1847  else {
1848  for (size_type k = 0; k < numRows-1; ++k) {
1849  cerr << numEntriesPerRow[k] << " ";
1850  }
1851  }
1852  cerr << numEntriesPerRow[numRows-1];
1853  } // numRows > 0
1854  cerr << "]" << endl;
1855 
1856  cerr << "----- Proc 0: rowPtr ";
1857  if (numRows > maxToDisplay) {
1858  cerr << "(only showing first and last few entries) ";
1859  }
1860  cerr << "= [";
1861  if (numRows > maxToDisplay) {
1862  for (size_type k = 0; k < 2; ++k) {
1863  cerr << rowPtr[k] << " ";
1864  }
1865  cerr << "... ";
1866  for (size_type k = numRows-2; k < numRows; ++k) {
1867  cerr << rowPtr[k] << " ";
1868  }
1869  }
1870  else {
1871  for (size_type k = 0; k < numRows; ++k) {
1872  cerr << rowPtr[k] << " ";
1873  }
1874  }
1875  cerr << rowPtr[numRows] << "]" << endl;
1876  }
1877  } // if myRank == rootRank
1878  } // Done converting sparse matrix data to CSR format
1879 
1880  // Now we're done with the Adder, so we can release the
1881  // reference ("free" it) to save space. This only actually
1882  // does anything on Rank 0, since pAdder is null on all the
1883  // other MPI processes.
1884  pAdder = null;
1885 
1886  if (debug && myRank == rootRank) {
1887  cerr << "-- Making range, domain, and row maps" << endl;
1888  }
1889 
1890  // Make the maps that describe the matrix's range and domain,
1891  // and the distribution of its rows. Creating a Map is a
1892  // collective operation, so we don't have to do a broadcast of
1893  // a success Boolean.
1894  RCP<const map_type> pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
1895  RCP<const map_type> pDomainMap =
1896  makeDomainMap (pRangeMap, dims[0], dims[1]);
1897  RCP<const map_type> pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
1898 
1899  if (debug && myRank == rootRank) {
1900  cerr << "-- Distributing the matrix data" << endl;
1901  }
1902 
1903  // Distribute the matrix data. Each processor has to add the
1904  // rows that it owns. If you try to make Proc 0 call
1905  // insertGlobalValues() for _all_ the rows, not just those it
1906  // owns, then fillComplete() will compute the number of
1907  // columns incorrectly. That's why Proc 0 has to distribute
1908  // the matrix data and why we make all the processors (not
1909  // just Proc 0) call insertGlobalValues() on their own data.
1910  //
1911  // These arrays represent each processor's part of the matrix
1912  // data, in "CSR" format (sort of, since the row indices might
1913  // not be contiguous).
1914  ArrayRCP<size_t> myNumEntriesPerRow;
1915  ArrayRCP<size_t> myRowPtr;
1916  ArrayRCP<global_ordinal_type> myColInd;
1917  ArrayRCP<scalar_type> myValues;
1918  // Distribute the matrix data. This is a collective operation.
1919  distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
1920  numEntriesPerRow, rowPtr, colInd, values, debug);
1921 
1922  if (debug && myRank == rootRank) {
1923  cerr << "-- Inserting matrix entries on each processor";
1924  if (callFillComplete) {
1925  cerr << " and calling fillComplete()";
1926  }
1927  cerr << endl;
1928  }
1929  // Each processor inserts its part of the matrix data, and
1930  // then they all call fillComplete(). This method invalidates
1931  // the my* distributed matrix data before calling
1932  // fillComplete(), in order to save space. In general, we
1933  // never store more than two copies of the matrix's entries in
1934  // memory at once, which is no worse than what Tpetra
1935  // promises.
1936  RCP<sparse_matrix_type> pMatrix =
1937  makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
1938  pRowMap, pRangeMap, pDomainMap, callFillComplete);
1939  // Only use a reduce-all in debug mode to check if pMatrix is
1940  // null. Otherwise, just throw an exception. We never expect
1941  // a null pointer here, so we can save a communication.
1942  if (debug) {
1943  int localIsNull = pMatrix.is_null () ? 1 : 0;
1944  int globalIsNull = 0;
1945  reduceAll (*pComm, REDUCE_MAX, localIsNull, ptr (&globalIsNull));
1946  TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
1947  "Reader::makeMatrix() returned a null pointer on at least one "
1948  "process. Please report this bug to the Tpetra developers.");
1949  }
1950  else {
1951  TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
1952  "Reader::makeMatrix() returned a null pointer. "
1953  "Please report this bug to the Tpetra developers.");
1954  }
1955 
1956  // We can't get the dimensions of the matrix until after
1957  // fillComplete() is called. Thus, we can't do the sanity
1958  // check (dimensions read from the Matrix Market data,
1959  // vs. dimensions reported by the CrsMatrix) unless the user
1960  // asked makeMatrix() to call fillComplete().
1961  //
1962  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
1963  // what one might think it does, so you have to ask the range
1964  // resp. domain map for the number of rows resp. columns.
1965  if (callFillComplete) {
1966  const int numProcs = pComm->getSize ();
1967 
1968  if (extraDebug && debug) {
1969  const global_size_t globalNumRows =
1970  pRangeMap->getGlobalNumElements ();
1971  const global_size_t globalNumCols =
1972  pDomainMap->getGlobalNumElements ();
1973  if (myRank == rootRank) {
1974  cerr << "-- Matrix is "
1975  << globalNumRows << " x " << globalNumCols
1976  << " with " << pMatrix->getGlobalNumEntries()
1977  << " entries, and index base "
1978  << pMatrix->getIndexBase() << "." << endl;
1979  }
1980  pComm->barrier ();
1981  for (int p = 0; p < numProcs; ++p) {
1982  if (myRank == p) {
1983  cerr << "-- Proc " << p << " owns "
1984  << pMatrix->getNodeNumCols() << " columns, and "
1985  << pMatrix->getNodeNumEntries() << " entries." << endl;
1986  }
1987  pComm->barrier ();
1988  }
1989  } // if (extraDebug && debug)
1990  } // if (callFillComplete)
1991 
1992  if (debug && myRank == rootRank) {
1993  cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
1994  << endl;
1995  }
1996  return pMatrix;
1997  }
1998 
1999 
2028  static Teuchos::RCP<sparse_matrix_type>
2029  readSparse (std::istream& in,
2030  const RCP<const Comm<int> >& pComm,
2031  const RCP<Teuchos::ParameterList>& constructorParams,
2032  const RCP<Teuchos::ParameterList>& fillCompleteParams,
2033  const bool tolerant=false,
2034  const bool debug=false)
2035  {
2036  return readSparse (in, pComm, Teuchos::null, constructorParams,
2037  fillCompleteParams, tolerant, debug);
2038  }
2039 
2041  static Teuchos::RCP<sparse_matrix_type>
2042  readSparse (std::istream& in,
2043  const RCP<const Comm<int> >& pComm,
2044  const RCP<node_type>& pNode,
2045  const RCP<Teuchos::ParameterList>& constructorParams,
2046  const RCP<Teuchos::ParameterList>& fillCompleteParams,
2047  const bool tolerant=false,
2048  const bool debug=false)
2049  {
2050  using Teuchos::MatrixMarket::Banner;
2051  using Teuchos::broadcast;
2052  using Teuchos::ptr;
2053  using Teuchos::reduceAll;
2054  using std::cerr;
2055  using std::endl;
2056  typedef Teuchos::ScalarTraits<scalar_type> STS;
2057 
2058  const bool extraDebug = false;
2059  const int myRank = pComm->getRank ();
2060  const int rootRank = 0;
2061 
2062  // Current line number in the input stream. Various calls
2063  // will modify this depending on the number of lines that are
2064  // read from the input stream. Only Rank 0 modifies this.
2065  size_t lineNumber = 1;
2066 
2067  if (debug && myRank == rootRank) {
2068  cerr << "Matrix Market reader: readSparse:" << endl
2069  << "-- Reading banner line" << endl;
2070  }
2071 
2072  // The "Banner" tells you whether the input stream represents
2073  // a sparse matrix, the symmetry type of the matrix, and the
2074  // type of the data it contains.
2075  //
2076  // pBanner will only be nonnull on MPI Rank 0. It will be
2077  // null on all other MPI processes.
2078  RCP<const Banner> pBanner;
2079  {
2080  // We read and validate the Banner on Proc 0, but broadcast
2081  // the validation result to all processes.
2082  // Teuchos::broadcast doesn't currently work with bool, so
2083  // we use int (true -> 1, false -> 0).
2084  int bannerIsCorrect = 1;
2085  std::ostringstream errMsg;
2086 
2087  if (myRank == rootRank) {
2088  // Read the Banner line from the input stream.
2089  try {
2090  pBanner = readBanner (in, lineNumber, tolerant, debug);
2091  }
2092  catch (std::exception& e) {
2093  errMsg << "Attempt to read the Matrix Market file's Banner line "
2094  "threw an exception: " << e.what();
2095  bannerIsCorrect = 0;
2096  }
2097 
2098  if (bannerIsCorrect) {
2099  // Validate the Banner for the case of a sparse matrix.
2100  // We validate on Proc 0, since it reads the Banner.
2101 
2102  // In intolerant mode, the matrix type must be "coordinate".
2103  if (! tolerant && pBanner->matrixType() != "coordinate") {
2104  bannerIsCorrect = 0;
2105  errMsg << "The Matrix Market input file must contain a "
2106  "\"coordinate\"-format sparse matrix in order to create a "
2107  "Tpetra::CrsMatrix object from it, but the file's matrix "
2108  "type is \"" << pBanner->matrixType() << "\" instead.";
2109  }
2110  // In tolerant mode, we allow the matrix type to be
2111  // anything other than "array" (which would mean that
2112  // the file contains a dense matrix).
2113  if (tolerant && pBanner->matrixType() == "array") {
2114  bannerIsCorrect = 0;
2115  errMsg << "Matrix Market file must contain a \"coordinate\"-"
2116  "format sparse matrix in order to create a Tpetra::CrsMatrix "
2117  "object from it, but the file's matrix type is \"array\" "
2118  "instead. That probably means the file contains dense matrix "
2119  "data.";
2120  }
2121  }
2122  } // Proc 0: Done reading the Banner, hopefully successfully.
2123 
2124  // Broadcast from Proc 0 whether the Banner was read correctly.
2125  broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2126 
2127  // If the Banner is invalid, all processes throw an
2128  // exception. Only Proc 0 gets the exception message, but
2129  // that's OK, since the main point is to "stop the world"
2130  // (rather than throw an exception on one process and leave
2131  // the others hanging).
2132  TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2133  std::invalid_argument, errMsg.str ());
2134  } // Done reading the Banner line and broadcasting success.
2135  if (debug && myRank == rootRank) {
2136  cerr << "-- Reading dimensions line" << endl;
2137  }
2138 
2139  // Read the matrix dimensions from the Matrix Market metadata.
2140  // dims = (numRows, numCols, numEntries). Proc 0 does the
2141  // reading, but it broadcasts the results to all MPI
2142  // processes. Thus, readCoordDims() is a collective
2143  // operation. It does a collective check for correctness too.
2144  Tuple<global_ordinal_type, 3> dims =
2145  readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2146 
2147  if (debug && myRank == rootRank) {
2148  cerr << "-- Making Adder for collecting matrix data" << endl;
2149  }
2150 
2151  // "Adder" object for collecting all the sparse matrix entries
2152  // from the input stream. This is only nonnull on Proc 0.
2153  RCP<adder_type> pAdder =
2154  makeAdder (pComm, pBanner, dims, tolerant, debug);
2155 
2156  if (debug && myRank == rootRank) {
2157  cerr << "-- Reading matrix data" << endl;
2158  }
2159  //
2160  // Read the matrix entries from the input stream on Proc 0.
2161  //
2162  {
2163  // We use readSuccess to broadcast the results of the read
2164  // (succeeded or not) to all MPI processes. Since
2165  // Teuchos::broadcast doesn't currently know how to send
2166  // bools, we convert to int (true -> 1, false -> 0).
2167  int readSuccess = 1;
2168  std::ostringstream errMsg; // Exception message (only valid on Proc 0)
2169  if (myRank == rootRank) {
2170  try {
2171  // Reader for "coordinate" format sparse matrix data.
2172  typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2173  global_ordinal_type, scalar_type, STS::isComplex> reader_type;
2174  reader_type reader (pAdder);
2175 
2176  // Read the sparse matrix entries.
2177  std::pair<bool, std::vector<size_t> > results =
2178  reader.read (in, lineNumber, tolerant, debug);
2179  readSuccess = results.first ? 1 : 0;
2180  }
2181  catch (std::exception& e) {
2182  readSuccess = 0;
2183  errMsg << e.what();
2184  }
2185  }
2186  broadcast (*pComm, rootRank, ptr (&readSuccess));
2187 
2188  // It would be nice to add a "verbose" flag, so that in
2189  // tolerant mode, we could log any bad line number(s) on
2190  // Proc 0. For now, we just throw if the read fails to
2191  // succeed.
2192  //
2193  // Question: If we're in tolerant mode, and if the read did
2194  // not succeed, should we attempt to call fillComplete()?
2195  TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2196  "Failed to read the Matrix Market sparse matrix file: "
2197  << errMsg.str());
2198  } // Done reading the matrix entries (stored on Proc 0 for now)
2199 
2200  if (debug && myRank == rootRank) {
2201  cerr << "-- Successfully read the Matrix Market data" << endl;
2202  }
2203 
2204  // In tolerant mode, we need to rebroadcast the matrix
2205  // dimensions, since they may be different after reading the
2206  // actual matrix data. We only need to broadcast the number
2207  // of rows and columns. Only Rank 0 needs to know the actual
2208  // global number of entries, since (a) we need to merge
2209  // duplicates on Rank 0 first anyway, and (b) when we
2210  // distribute the entries, each rank other than Rank 0 will
2211  // only need to know how many entries it owns, not the total
2212  // number of entries.
2213  if (tolerant) {
2214  if (debug && myRank == rootRank) {
2215  cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
2216  << endl
2217  << "----- Dimensions before: "
2218  << dims[0] << " x " << dims[1]
2219  << endl;
2220  }
2221  // Packed coordinate matrix dimensions (numRows, numCols).
2222  Tuple<global_ordinal_type, 2> updatedDims;
2223  if (myRank == rootRank) {
2224  // If one or more bottom rows of the matrix contain no
2225  // entries, then the Adder will report that the number
2226  // of rows is less than that specified in the
2227  // metadata. We allow this case, and favor the
2228  // metadata so that the zero row(s) will be included.
2229  updatedDims[0] =
2230  std::max (dims[0], pAdder->getAdder()->numRows());
2231  updatedDims[1] = pAdder->getAdder()->numCols();
2232  }
2233  broadcast (*pComm, rootRank, updatedDims);
2234  dims[0] = updatedDims[0];
2235  dims[1] = updatedDims[1];
2236  if (debug && myRank == rootRank) {
2237  cerr << "----- Dimensions after: " << dims[0] << " x "
2238  << dims[1] << endl;
2239  }
2240  }
2241  else {
2242  // In strict mode, we require that the matrix's metadata and
2243  // its actual data agree, at least somewhat. In particular,
2244  // the number of rows must agree, since otherwise we cannot
2245  // distribute the matrix correctly.
2246 
2247  // Teuchos::broadcast() doesn't know how to broadcast bools,
2248  // so we use an int with the standard 1 == true, 0 == false
2249  // encoding.
2250  int dimsMatch = 1;
2251  if (myRank == rootRank) {
2252  // If one or more bottom rows of the matrix contain no
2253  // entries, then the Adder will report that the number of
2254  // rows is less than that specified in the metadata. We
2255  // allow this case, and favor the metadata, but do not
2256  // allow the Adder to think there are more rows in the
2257  // matrix than the metadata says.
2258  if (dims[0] < pAdder->getAdder ()->numRows ()) {
2259  dimsMatch = 0;
2260  }
2261  }
2262  broadcast (*pComm, 0, ptr (&dimsMatch));
2263  if (dimsMatch == 0) {
2264  // We're in an error state anyway, so we might as well
2265  // work a little harder to print an informative error
2266  // message.
2267  //
2268  // Broadcast the Adder's idea of the matrix dimensions
2269  // from Proc 0 to all processes.
2270  Tuple<global_ordinal_type, 2> addersDims;
2271  if (myRank == rootRank) {
2272  addersDims[0] = pAdder->getAdder()->numRows();
2273  addersDims[1] = pAdder->getAdder()->numCols();
2274  }
2275  broadcast (*pComm, 0, addersDims);
2276  TEUCHOS_TEST_FOR_EXCEPTION(
2277  dimsMatch == 0, std::runtime_error,
2278  "The matrix metadata says that the matrix is " << dims[0] << " x "
2279  << dims[1] << ", but the actual data says that the matrix is "
2280  << addersDims[0] << " x " << addersDims[1] << ". That means the "
2281  "data includes more rows than reported in the metadata. This "
2282  "is not allowed when parsing in strict mode. Parse the matrix in "
2283  "tolerant mode to ignore the metadata when it disagrees with the "
2284  "data.");
2285  }
2286  } // Matrix dimensions (# rows, # cols, # entries) agree.
2287 
2288  if (debug && myRank == rootRank) {
2289  cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
2290  }
2291 
2292  // Now that we've read in all the matrix entries from the
2293  // input stream into the adder on Proc 0, post-process them
2294  // into CSR format (still on Proc 0). This will facilitate
2295  // distributing them to all the processors.
2296  //
2297  // These arrays represent the global matrix data as a CSR
2298  // matrix (with numEntriesPerRow as redundant but convenient
2299  // metadata, since it's computable from rowPtr and vice
2300  // versa). They are valid only on Proc 0.
2301  ArrayRCP<size_t> numEntriesPerRow;
2302  ArrayRCP<size_t> rowPtr;
2303  ArrayRCP<global_ordinal_type> colInd;
2304  ArrayRCP<scalar_type> values;
2305 
2306  // Proc 0 first merges duplicate entries, and then converts
2307  // the coordinate-format matrix data to CSR.
2308  {
2309  int mergeAndConvertSucceeded = 1;
2310  std::ostringstream errMsg;
2311 
2312  if (myRank == rootRank) {
2313  try {
2314  typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
2315  global_ordinal_type> element_type;
2316 
2317  // Number of rows in the matrix. If we are in tolerant
2318  // mode, we've already synchronized dims with the actual
2319  // matrix data. If in strict mode, we should use dims
2320  // (as read from the file's metadata) rather than the
2321  // matrix data to determine the dimensions. (The matrix
2322  // data will claim fewer rows than the metadata, if one
2323  // or more rows have no entries stored in the file.)
2324  const size_type numRows = dims[0];
2325 
2326  // Additively merge duplicate matrix entries.
2327  pAdder->getAdder()->merge ();
2328 
2329  // Get a temporary const view of the merged matrix entries.
2330  const std::vector<element_type>& entries =
2331  pAdder->getAdder()->getEntries();
2332 
2333  // Number of matrix entries (after merging).
2334  const size_t numEntries = (size_t)entries.size();
2335 
2336  if (debug) {
2337  cerr << "----- Proc 0: Matrix has numRows=" << numRows
2338  << " rows and numEntries=" << numEntries
2339  << " entries." << endl;
2340  }
2341 
2342  // Make space for the CSR matrix data. Converting to
2343  // CSR is easier if we fill numEntriesPerRow with zeros
2344  // at first.
2345  numEntriesPerRow = arcp<size_t> (numRows);
2346  std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2347  rowPtr = arcp<size_t> (numRows+1);
2348  std::fill (rowPtr.begin(), rowPtr.end(), 0);
2349  colInd = arcp<global_ordinal_type> (numEntries);
2350  values = arcp<scalar_type> (numEntries);
2351 
2352  // Convert from array-of-structs coordinate format to CSR
2353  // (compressed sparse row) format.
2354  global_ordinal_type prvRow = 0;
2355  size_t curPos = 0;
2356  rowPtr[0] = 0;
2357  for (curPos = 0; curPos < numEntries; ++curPos) {
2358  const element_type& curEntry = entries[curPos];
2359  const global_ordinal_type curRow = curEntry.rowIndex();
2360  TEUCHOS_TEST_FOR_EXCEPTION(
2361  curRow < prvRow, std::logic_error,
2362  "Row indices are out of order, even though they are supposed "
2363  "to be sorted. curRow = " << curRow << ", prvRow = "
2364  << prvRow << ", at curPos = " << curPos << ". Please report "
2365  "this bug to the Tpetra developers.");
2366  if (curRow > prvRow) {
2367  for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
2368  rowPtr[r] = curPos;
2369  }
2370  prvRow = curRow;
2371  }
2372  numEntriesPerRow[curRow]++;
2373  colInd[curPos] = curEntry.colIndex();
2374  values[curPos] = curEntry.value();
2375  }
2376  // rowPtr has one more entry than numEntriesPerRow. The
2377  // last entry of rowPtr is the number of entries in
2378  // colInd and values.
2379  rowPtr[numRows] = numEntries;
2380  } // Finished conversion to CSR format
2381  catch (std::exception& e) {
2382  mergeAndConvertSucceeded = 0;
2383  errMsg << "Failed to merge sparse matrix entries and convert to "
2384  "CSR format: " << e.what();
2385  }
2386 
2387  if (debug && mergeAndConvertSucceeded) {
2388  // Number of rows in the matrix.
2389  const size_type numRows = dims[0];
2390  const size_type maxToDisplay = 100;
2391 
2392  cerr << "----- Proc 0: numEntriesPerRow[0.."
2393  << (numEntriesPerRow.size()-1) << "] ";
2394  if (numRows > maxToDisplay) {
2395  cerr << "(only showing first and last few entries) ";
2396  }
2397  cerr << "= [";
2398  if (numRows > 0) {
2399  if (numRows > maxToDisplay) {
2400  for (size_type k = 0; k < 2; ++k) {
2401  cerr << numEntriesPerRow[k] << " ";
2402  }
2403  cerr << "... ";
2404  for (size_type k = numRows-2; k < numRows-1; ++k) {
2405  cerr << numEntriesPerRow[k] << " ";
2406  }
2407  }
2408  else {
2409  for (size_type k = 0; k < numRows-1; ++k) {
2410  cerr << numEntriesPerRow[k] << " ";
2411  }
2412  }
2413  cerr << numEntriesPerRow[numRows-1];
2414  } // numRows > 0
2415  cerr << "]" << endl;
2416 
2417  cerr << "----- Proc 0: rowPtr ";
2418  if (numRows > maxToDisplay) {
2419  cerr << "(only showing first and last few entries) ";
2420  }
2421  cerr << "= [";
2422  if (numRows > maxToDisplay) {
2423  for (size_type k = 0; k < 2; ++k) {
2424  cerr << rowPtr[k] << " ";
2425  }
2426  cerr << "... ";
2427  for (size_type k = numRows-2; k < numRows; ++k) {
2428  cerr << rowPtr[k] << " ";
2429  }
2430  }
2431  else {
2432  for (size_type k = 0; k < numRows; ++k) {
2433  cerr << rowPtr[k] << " ";
2434  }
2435  }
2436  cerr << rowPtr[numRows] << "]" << endl;
2437  }
2438  } // if myRank == rootRank
2439  } // Done converting sparse matrix data to CSR format
2440 
2441  // Now we're done with the Adder, so we can release the
2442  // reference ("free" it) to save space. This only actually
2443  // does anything on Rank 0, since pAdder is null on all the
2444  // other MPI processes.
2445  pAdder = null;
2446 
2447  if (debug && myRank == rootRank) {
2448  cerr << "-- Making range, domain, and row maps" << endl;
2449  }
2450 
2451  // Make the maps that describe the matrix's range and domain,
2452  // and the distribution of its rows. Creating a Map is a
2453  // collective operation, so we don't have to do a broadcast of
2454  // a success Boolean.
2455  RCP<const map_type> pRangeMap = makeRangeMap (pComm, pNode, dims[0]);
2456  RCP<const map_type> pDomainMap =
2457  makeDomainMap (pRangeMap, dims[0], dims[1]);
2458  RCP<const map_type> pRowMap = makeRowMap (null, pComm, pNode, dims[0]);
2459 
2460  if (debug && myRank == rootRank) {
2461  cerr << "-- Distributing the matrix data" << endl;
2462  }
2463 
2464  // Distribute the matrix data. Each processor has to add the
2465  // rows that it owns. If you try to make Proc 0 call
2466  // insertGlobalValues() for _all_ the rows, not just those it
2467  // owns, then fillComplete() will compute the number of
2468  // columns incorrectly. That's why Proc 0 has to distribute
2469  // the matrix data and why we make all the processors (not
2470  // just Proc 0) call insertGlobalValues() on their own data.
2471  //
2472  // These arrays represent each processor's part of the matrix
2473  // data, in "CSR" format (sort of, since the row indices might
2474  // not be contiguous).
2475  ArrayRCP<size_t> myNumEntriesPerRow;
2476  ArrayRCP<size_t> myRowPtr;
2477  ArrayRCP<global_ordinal_type> myColInd;
2478  ArrayRCP<scalar_type> myValues;
2479  // Distribute the matrix data. This is a collective operation.
2480  distribute (myNumEntriesPerRow, myRowPtr, myColInd, myValues, pRowMap,
2481  numEntriesPerRow, rowPtr, colInd, values, debug);
2482 
2483  if (debug && myRank == rootRank) {
2484  cerr << "-- Inserting matrix entries on each process "
2485  "and calling fillComplete()" << endl;
2486  }
2487  // Each processor inserts its part of the matrix data, and
2488  // then they all call fillComplete(). This method invalidates
2489  // the my* distributed matrix data before calling
2490  // fillComplete(), in order to save space. In general, we
2491  // never store more than two copies of the matrix's entries in
2492  // memory at once, which is no worse than what Tpetra
2493  // promises.
2494  Teuchos::RCP<sparse_matrix_type> pMatrix =
2495  makeMatrix (myNumEntriesPerRow, myRowPtr, myColInd, myValues,
2496  pRowMap, pRangeMap, pDomainMap, constructorParams,
2497  fillCompleteParams);
2498  // Only use a reduce-all in debug mode to check if pMatrix is
2499  // null. Otherwise, just throw an exception. We never expect
2500  // a null pointer here, so we can save a communication.
2501  if (debug) {
2502  int localIsNull = pMatrix.is_null () ? 1 : 0;
2503  int globalIsNull = 0;
2504  reduceAll (*pComm, Teuchos::REDUCE_MAX, localIsNull, ptr (&globalIsNull));
2505  TEUCHOS_TEST_FOR_EXCEPTION(globalIsNull != 0, std::logic_error,
2506  "Reader::makeMatrix() returned a null pointer on at least one "
2507  "process. Please report this bug to the Tpetra developers.");
2508  }
2509  else {
2510  TEUCHOS_TEST_FOR_EXCEPTION(pMatrix.is_null(), std::logic_error,
2511  "Reader::makeMatrix() returned a null pointer. "
2512  "Please report this bug to the Tpetra developers.");
2513  }
2514 
2515  // Sanity check for dimensions (read from the Matrix Market
2516  // data, vs. dimensions reported by the CrsMatrix).
2517  //
2518  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
2519  // what one might think it does, so you have to ask the range
2520  // resp. domain map for the number of rows resp. columns.
2521  if (extraDebug && debug) {
2522  const int numProcs = pComm->getSize ();
2523  const global_size_t globalNumRows =
2524  pRangeMap->getGlobalNumElements();
2525  const global_size_t globalNumCols =
2526  pDomainMap->getGlobalNumElements();
2527  if (myRank == rootRank) {
2528  cerr << "-- Matrix is "
2529  << globalNumRows << " x " << globalNumCols
2530  << " with " << pMatrix->getGlobalNumEntries()
2531  << " entries, and index base "
2532  << pMatrix->getIndexBase() << "." << endl;
2533  }
2534  pComm->barrier ();
2535  for (int p = 0; p < numProcs; ++p) {
2536  if (myRank == p) {
2537  cerr << "-- Proc " << p << " owns "
2538  << pMatrix->getNodeNumCols() << " columns, and "
2539  << pMatrix->getNodeNumEntries() << " entries." << endl;
2540  }
2541  pComm->barrier ();
2542  }
2543  } // if (extraDebug && debug)
2544 
2545  if (debug && myRank == rootRank) {
2546  cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
2547  << endl;
2548  }
2549  return pMatrix;
2550  }
2551 
2592  static Teuchos::RCP<sparse_matrix_type>
2593  readSparse (std::istream& in,
2594  const Teuchos::RCP<const map_type>& rowMap,
2595  Teuchos::RCP<const map_type>& colMap,
2596  const Teuchos::RCP<const map_type>& domainMap,
2597  const Teuchos::RCP<const map_type>& rangeMap,
2598  const bool callFillComplete=true,
2599  const bool tolerant=false,
2600  const bool debug=false)
2601  {
2602  using Teuchos::MatrixMarket::Banner;
2603  using Teuchos::as;
2604  using Teuchos::broadcast;
2605  using Teuchos::Comm;
2606  using Teuchos::ptr;
2607  using Teuchos::RCP;
2608  using Teuchos::reduceAll;
2609  using Teuchos::Tuple;
2610  using std::cerr;
2611  using std::endl;
2612  typedef Teuchos::ScalarTraits<scalar_type> STS;
2613 
2614  RCP<const Comm<int> > pComm = rowMap->getComm ();
2615  const int myRank = pComm->getRank ();
2616  const int rootRank = 0;
2617  const bool extraDebug = false;
2618 
2619  // Fast checks for invalid input. We can't check other
2620  // attributes of the Maps until we've read in the matrix
2621  // dimensions.
2622  TEUCHOS_TEST_FOR_EXCEPTION(
2623  rowMap.is_null (), std::invalid_argument,
2624  "Row Map must be nonnull.");
2625  TEUCHOS_TEST_FOR_EXCEPTION(
2626  rangeMap.is_null (), std::invalid_argument,
2627  "Range Map must be nonnull.");
2628  TEUCHOS_TEST_FOR_EXCEPTION(
2629  domainMap.is_null (), std::invalid_argument,
2630  "Domain Map must be nonnull.");
2631  TEUCHOS_TEST_FOR_EXCEPTION(
2632  rowMap->getComm().getRawPtr() != pComm.getRawPtr(),
2633  std::invalid_argument,
2634  "The specified row Map's communicator (rowMap->getComm())"
2635  "differs from the given separately supplied communicator pComm.");
2636  TEUCHOS_TEST_FOR_EXCEPTION(
2637  domainMap->getComm().getRawPtr() != pComm.getRawPtr(),
2638  std::invalid_argument,
2639  "The specified domain Map's communicator (domainMap->getComm())"
2640  "differs from the given separately supplied communicator pComm.");
2641  TEUCHOS_TEST_FOR_EXCEPTION(
2642  rangeMap->getComm().getRawPtr() != pComm.getRawPtr(),
2643  std::invalid_argument,
2644  "The specified range Map's communicator (rangeMap->getComm())"
2645  "differs from the given separately supplied communicator pComm.");
2646 
2647  // Current line number in the input stream. Various calls
2648  // will modify this depending on the number of lines that are
2649  // read from the input stream. Only Rank 0 modifies this.
2650  size_t lineNumber = 1;
2651 
2652  if (debug && myRank == rootRank) {
2653  cerr << "Matrix Market reader: readSparse:" << endl
2654  << "-- Reading banner line" << endl;
2655  }
2656 
2657  // The "Banner" tells you whether the input stream represents
2658  // a sparse matrix, the symmetry type of the matrix, and the
2659  // type of the data it contains.
2660  //
2661  // pBanner will only be nonnull on MPI Rank 0. It will be
2662  // null on all other MPI processes.
2663  RCP<const Banner> pBanner;
2664  {
2665  // We read and validate the Banner on Proc 0, but broadcast
2666  // the validation result to all processes.
2667  // Teuchos::broadcast doesn't currently work with bool, so
2668  // we use int (true -> 1, false -> 0).
2669  int bannerIsCorrect = 1;
2670  std::ostringstream errMsg;
2671 
2672  if (myRank == rootRank) {
2673  // Read the Banner line from the input stream.
2674  try {
2675  pBanner = readBanner (in, lineNumber, tolerant, debug);
2676  }
2677  catch (std::exception& e) {
2678  errMsg << "Attempt to read the Matrix Market file's Banner line "
2679  "threw an exception: " << e.what();
2680  bannerIsCorrect = 0;
2681  }
2682 
2683  if (bannerIsCorrect) {
2684  // Validate the Banner for the case of a sparse matrix.
2685  // We validate on Proc 0, since it reads the Banner.
2686 
2687  // In intolerant mode, the matrix type must be "coordinate".
2688  if (! tolerant && pBanner->matrixType() != "coordinate") {
2689  bannerIsCorrect = 0;
2690  errMsg << "The Matrix Market input file must contain a "
2691  "\"coordinate\"-format sparse matrix in order to create a "
2692  "Tpetra::CrsMatrix object from it, but the file's matrix "
2693  "type is \"" << pBanner->matrixType() << "\" instead.";
2694  }
2695  // In tolerant mode, we allow the matrix type to be
2696  // anything other than "array" (which would mean that
2697  // the file contains a dense matrix).
2698  if (tolerant && pBanner->matrixType() == "array") {
2699  bannerIsCorrect = 0;
2700  errMsg << "Matrix Market file must contain a \"coordinate\"-"
2701  "format sparse matrix in order to create a Tpetra::CrsMatrix "
2702  "object from it, but the file's matrix type is \"array\" "
2703  "instead. That probably means the file contains dense matrix "
2704  "data.";
2705  }
2706  }
2707  } // Proc 0: Done reading the Banner, hopefully successfully.
2708 
2709  // Broadcast from Proc 0 whether the Banner was read correctly.
2710  broadcast (*pComm, rootRank, ptr (&bannerIsCorrect));
2711 
2712  // If the Banner is invalid, all processes throw an
2713  // exception. Only Proc 0 gets the exception message, but
2714  // that's OK, since the main point is to "stop the world"
2715  // (rather than throw an exception on one process and leave
2716  // the others hanging).
2717  TEUCHOS_TEST_FOR_EXCEPTION(bannerIsCorrect == 0,
2718  std::invalid_argument, errMsg.str ());
2719  } // Done reading the Banner line and broadcasting success.
2720  if (debug && myRank == rootRank) {
2721  cerr << "-- Reading dimensions line" << endl;
2722  }
2723 
2724  // Read the matrix dimensions from the Matrix Market metadata.
2725  // dims = (numRows, numCols, numEntries). Proc 0 does the
2726  // reading, but it broadcasts the results to all MPI
2727  // processes. Thus, readCoordDims() is a collective
2728  // operation. It does a collective check for correctness too.
2729  Tuple<global_ordinal_type, 3> dims =
2730  readCoordDims (in, lineNumber, pBanner, pComm, tolerant, debug);
2731 
2732  if (debug && myRank == rootRank) {
2733  cerr << "-- Making Adder for collecting matrix data" << endl;
2734  }
2735 
2736  // "Adder" object for collecting all the sparse matrix entries
2737  // from the input stream. This is only nonnull on Proc 0.
2738  // The Adder internally converts the one-based indices (native
2739  // Matrix Market format) into zero-based indices.
2740  RCP<adder_type> pAdder =
2741  makeAdder (pComm, pBanner, dims, tolerant, debug);
2742 
2743  if (debug && myRank == rootRank) {
2744  cerr << "-- Reading matrix data" << endl;
2745  }
2746  //
2747  // Read the matrix entries from the input stream on Proc 0.
2748  //
2749  {
2750  // We use readSuccess to broadcast the results of the read
2751  // (succeeded or not) to all MPI processes. Since
2752  // Teuchos::broadcast doesn't currently know how to send
2753  // bools, we convert to int (true -> 1, false -> 0).
2754  int readSuccess = 1;
2755  std::ostringstream errMsg; // Exception message (only valid on Proc 0)
2756  if (myRank == rootRank) {
2757  try {
2758  // Reader for "coordinate" format sparse matrix data.
2759  typedef Teuchos::MatrixMarket::CoordDataReader<adder_type,
2760  global_ordinal_type, scalar_type, STS::isComplex> reader_type;
2761  reader_type reader (pAdder);
2762 
2763  // Read the sparse matrix entries.
2764  std::pair<bool, std::vector<size_t> > results =
2765  reader.read (in, lineNumber, tolerant, debug);
2766  readSuccess = results.first ? 1 : 0;
2767  }
2768  catch (std::exception& e) {
2769  readSuccess = 0;
2770  errMsg << e.what();
2771  }
2772  }
2773  broadcast (*pComm, rootRank, ptr (&readSuccess));
2774 
2775  // It would be nice to add a "verbose" flag, so that in
2776  // tolerant mode, we could log any bad line number(s) on
2777  // Proc 0. For now, we just throw if the read fails to
2778  // succeed.
2779  //
2780  // Question: If we're in tolerant mode, and if the read did
2781  // not succeed, should we attempt to call fillComplete()?
2782  TEUCHOS_TEST_FOR_EXCEPTION(readSuccess == 0, std::runtime_error,
2783  "Failed to read the Matrix Market sparse matrix file: "
2784  << errMsg.str());
2785  } // Done reading the matrix entries (stored on Proc 0 for now)
2786 
2787  if (debug && myRank == rootRank) {
2788  cerr << "-- Successfully read the Matrix Market data" << endl;
2789  }
2790 
2791  // In tolerant mode, we need to rebroadcast the matrix
2792  // dimensions, since they may be different after reading the
2793  // actual matrix data. We only need to broadcast the number
2794  // of rows and columns. Only Rank 0 needs to know the actual
2795  // global number of entries, since (a) we need to merge
2796  // duplicates on Rank 0 first anyway, and (b) when we
2797  // distribute the entries, each rank other than Rank 0 will
2798  // only need to know how many entries it owns, not the total
2799  // number of entries.
2800  if (tolerant) {
2801  if (debug && myRank == rootRank) {
2802  cerr << "-- Tolerant mode: rebroadcasting matrix dimensions"
2803  << endl
2804  << "----- Dimensions before: "
2805  << dims[0] << " x " << dims[1]
2806  << endl;
2807  }
2808  // Packed coordinate matrix dimensions (numRows, numCols).
2809  Tuple<global_ordinal_type, 2> updatedDims;
2810  if (myRank == rootRank) {
2811  // If one or more bottom rows of the matrix contain no
2812  // entries, then the Adder will report that the number
2813  // of rows is less than that specified in the
2814  // metadata. We allow this case, and favor the
2815  // metadata so that the zero row(s) will be included.
2816  updatedDims[0] =
2817  std::max (dims[0], pAdder->getAdder()->numRows());
2818  updatedDims[1] = pAdder->getAdder()->numCols();
2819  }
2820  broadcast (*pComm, rootRank, updatedDims);
2821  dims[0] = updatedDims[0];
2822  dims[1] = updatedDims[1];
2823  if (debug && myRank == rootRank) {
2824  cerr << "----- Dimensions after: " << dims[0] << " x "
2825  << dims[1] << endl;
2826  }
2827  }
2828  else {
2829  // In strict mode, we require that the matrix's metadata and
2830  // its actual data agree, at least somewhat. In particular,
2831  // the number of rows must agree, since otherwise we cannot
2832  // distribute the matrix correctly.
2833 
2834  // Teuchos::broadcast() doesn't know how to broadcast bools,
2835  // so we use an int with the standard 1 == true, 0 == false
2836  // encoding.
2837  int dimsMatch = 1;
2838  if (myRank == rootRank) {
2839  // If one or more bottom rows of the matrix contain no
2840  // entries, then the Adder will report that the number of
2841  // rows is less than that specified in the metadata. We
2842  // allow this case, and favor the metadata, but do not
2843  // allow the Adder to think there are more rows in the
2844  // matrix than the metadata says.
2845  if (dims[0] < pAdder->getAdder ()->numRows ()) {
2846  dimsMatch = 0;
2847  }
2848  }
2849  broadcast (*pComm, 0, ptr (&dimsMatch));
2850  if (dimsMatch == 0) {
2851  // We're in an error state anyway, so we might as well
2852  // work a little harder to print an informative error
2853  // message.
2854  //
2855  // Broadcast the Adder's idea of the matrix dimensions
2856  // from Proc 0 to all processes.
2857  Tuple<global_ordinal_type, 2> addersDims;
2858  if (myRank == rootRank) {
2859  addersDims[0] = pAdder->getAdder()->numRows();
2860  addersDims[1] = pAdder->getAdder()->numCols();
2861  }
2862  broadcast (*pComm, 0, addersDims);
2863  TEUCHOS_TEST_FOR_EXCEPTION(
2864  dimsMatch == 0, std::runtime_error,
2865  "The matrix metadata says that the matrix is " << dims[0] << " x "
2866  << dims[1] << ", but the actual data says that the matrix is "
2867  << addersDims[0] << " x " << addersDims[1] << ". That means the "
2868  "data includes more rows than reported in the metadata. This "
2869  "is not allowed when parsing in strict mode. Parse the matrix in "
2870  "tolerant mode to ignore the metadata when it disagrees with the "
2871  "data.");
2872  }
2873  } // Matrix dimensions (# rows, # cols, # entries) agree.
2874 
2875  if (debug && myRank == rootRank) {
2876  cerr << "-- Converting matrix data into CSR format on Proc 0" << endl;
2877  }
2878 
2879  // Now that we've read in all the matrix entries from the
2880  // input stream into the adder on Proc 0, post-process them
2881  // into CSR format (still on Proc 0). This will facilitate
2882  // distributing them to all the processors.
2883  //
2884  // These arrays represent the global matrix data as a CSR
2885  // matrix (with numEntriesPerRow as redundant but convenient
2886  // metadata, since it's computable from rowPtr and vice
2887  // versa). They are valid only on Proc 0.
2888  ArrayRCP<size_t> numEntriesPerRow;
2889  ArrayRCP<size_t> rowPtr;
2890  ArrayRCP<global_ordinal_type> colInd;
2891  ArrayRCP<scalar_type> values;
2892 
2893  // Proc 0 first merges duplicate entries, and then converts
2894  // the coordinate-format matrix data to CSR.
2895  {
2896  int mergeAndConvertSucceeded = 1;
2897  std::ostringstream errMsg;
2898 
2899  if (myRank == rootRank) {
2900  try {
2901  typedef Teuchos::MatrixMarket::Raw::Element<scalar_type,
2902  global_ordinal_type> element_type;
2903 
2904  // Number of rows in the matrix. If we are in tolerant
2905  // mode, we've already synchronized dims with the actual
2906  // matrix data. If in strict mode, we should use dims
2907  // (as read from the file's metadata) rather than the
2908  // matrix data to determine the dimensions. (The matrix
2909  // data will claim fewer rows than the metadata, if one
2910  // or more rows have no entries stored in the file.)
2911  const size_type numRows = dims[0];
2912 
2913  // Additively merge duplicate matrix entries.
2914  pAdder->getAdder()->merge ();
2915 
2916  // Get a temporary const view of the merged matrix entries.
2917  const std::vector<element_type>& entries =
2918  pAdder->getAdder()->getEntries();
2919 
2920  // Number of matrix entries (after merging).
2921  const size_t numEntries = (size_t)entries.size();
2922 
2923  if (debug) {
2924  cerr << "----- Proc 0: Matrix has numRows=" << numRows
2925  << " rows and numEntries=" << numEntries
2926  << " entries." << endl;
2927  }
2928 
2929  // Make space for the CSR matrix data. Converting to
2930  // CSR is easier if we fill numEntriesPerRow with zeros
2931  // at first.
2932  numEntriesPerRow = arcp<size_t> (numRows);
2933  std::fill (numEntriesPerRow.begin(), numEntriesPerRow.end(), 0);
2934  rowPtr = arcp<size_t> (numRows+1);
2935  std::fill (rowPtr.begin(), rowPtr.end(), 0);
2936  colInd = arcp<global_ordinal_type> (numEntries);
2937  values = arcp<scalar_type> (numEntries);
2938 
2939  // Convert from array-of-structs coordinate format to CSR
2940  // (compressed sparse row) format.
2941  global_ordinal_type prvRow = 0;
2942  size_t curPos = 0;
2943  rowPtr[0] = 0;
2944  for (curPos = 0; curPos < numEntries; ++curPos) {
2945  const element_type& curEntry = entries[curPos];
2946  const global_ordinal_type curRow = curEntry.rowIndex();
2947  TEUCHOS_TEST_FOR_EXCEPTION(
2948  curRow < prvRow, std::logic_error,
2949  "Row indices are out of order, even though they are supposed "
2950  "to be sorted. curRow = " << curRow << ", prvRow = "
2951  << prvRow << ", at curPos = " << curPos << ". Please report "
2952  "this bug to the Tpetra developers.");
2953  if (curRow > prvRow) {
2954  for (global_ordinal_type r = prvRow+1; r <= curRow; ++r) {
2955  rowPtr[r] = curPos;
2956  }
2957  prvRow = curRow;
2958  }
2959  numEntriesPerRow[curRow]++;
2960  colInd[curPos] = curEntry.colIndex();
2961  values[curPos] = curEntry.value();
2962  }
2963  // rowPtr has one more entry than numEntriesPerRow. The
2964  // last entry of rowPtr is the number of entries in
2965  // colInd and values.
2966  rowPtr[numRows] = numEntries;
2967  } // Finished conversion to CSR format
2968  catch (std::exception& e) {
2969  mergeAndConvertSucceeded = 0;
2970  errMsg << "Failed to merge sparse matrix entries and convert to "
2971  "CSR format: " << e.what();
2972  }
2973 
2974  if (debug && mergeAndConvertSucceeded) {
2975  // Number of rows in the matrix.
2976  const size_type numRows = dims[0];
2977  const size_type maxToDisplay = 100;
2978 
2979  cerr << "----- Proc 0: numEntriesPerRow[0.."
2980  << (numEntriesPerRow.size()-1) << "] ";
2981  if (numRows > maxToDisplay) {
2982  cerr << "(only showing first and last few entries) ";
2983  }
2984  cerr << "= [";
2985  if (numRows > 0) {
2986  if (numRows > maxToDisplay) {
2987  for (size_type k = 0; k < 2; ++k) {
2988  cerr << numEntriesPerRow[k] << " ";
2989  }
2990  cerr << "... ";
2991  for (size_type k = numRows-2; k < numRows-1; ++k) {
2992  cerr << numEntriesPerRow[k] << " ";
2993  }
2994  }
2995  else {
2996  for (size_type k = 0; k < numRows-1; ++k) {
2997  cerr << numEntriesPerRow[k] << " ";
2998  }
2999  }
3000  cerr << numEntriesPerRow[numRows-1];
3001  } // numRows > 0
3002  cerr << "]" << endl;
3003 
3004  cerr << "----- Proc 0: rowPtr ";
3005  if (numRows > maxToDisplay) {
3006  cerr << "(only showing first and last few entries) ";
3007  }
3008  cerr << "= [";
3009  if (numRows > maxToDisplay) {
3010  for (size_type k = 0; k < 2; ++k) {
3011  cerr << rowPtr[k] << " ";
3012  }
3013  cerr << "... ";
3014  for (size_type k = numRows-2; k < numRows; ++k) {
3015  cerr << rowPtr[k] << " ";
3016  }
3017  }
3018  else {
3019  for (size_type k = 0; k < numRows; ++k) {
3020  cerr << rowPtr[k] << " ";
3021  }
3022  }
3023  cerr << rowPtr[numRows] << "]" << endl;
3024 
3025  cerr << "----- Proc 0: colInd = [";
3026  for (size_t k = 0; k < rowPtr[numRows]; ++k) {
3027  cerr << colInd[k] << " ";
3028  }
3029  cerr << "]" << endl;
3030  }
3031  } // if myRank == rootRank
3032  } // Done converting sparse matrix data to CSR format
3033 
3034  // Now we're done with the Adder, so we can release the
3035  // reference ("free" it) to save space. This only actually
3036  // does anything on Rank 0, since pAdder is null on all the
3037  // other MPI processes.
3038  pAdder = null;
3039 
3040  // Verify details of the Maps. Don't count the global number
3041  // of entries in the row Map, since that number doesn't
3042  // correctly count overlap.
3043  if (debug && myRank == rootRank) {
3044  cerr << "-- Verifying Maps" << endl;
3045  }
3046  TEUCHOS_TEST_FOR_EXCEPTION(
3047  as<global_size_t> (dims[0]) != rangeMap->getGlobalNumElements(),
3048  std::invalid_argument,
3049  "The range Map has " << rangeMap->getGlobalNumElements ()
3050  << " entries, but the matrix has a global number of rows " << dims[0]
3051  << ".");
3052  TEUCHOS_TEST_FOR_EXCEPTION(
3053  as<global_size_t> (dims[1]) != domainMap->getGlobalNumElements (),
3054  std::invalid_argument,
3055  "The domain Map has " << domainMap->getGlobalNumElements ()
3056  << " entries, but the matrix has a global number of columns "
3057  << dims[1] << ".");
3058 
3059  // Create a row Map which is entirely owned on Proc 0.
3060  RCP<Teuchos::FancyOStream> err = debug ?
3061  Teuchos::getFancyOStream (Teuchos::rcpFromRef (cerr)) : null;
3062 
3063  RCP<const map_type> gatherRowMap = Details::computeGatherMap (rowMap, err, debug);
3064  ArrayView<const global_ordinal_type> myRows =
3065  gatherRowMap->getNodeElementList ();
3066  const size_type myNumRows = myRows.size ();
3067  const global_ordinal_type indexBase = gatherRowMap->getIndexBase ();
3068 
3069  ArrayRCP<size_t> gatherNumEntriesPerRow = arcp<size_t>(myNumRows);
3070  for (size_type i_ = 0; i_ < myNumRows; i_++) {
3071  gatherNumEntriesPerRow[i_] = numEntriesPerRow[myRows[i_]-indexBase];
3072  }
3073 
3074  // Create a matrix using this Map, and fill in on Proc 0. We
3075  // know how many entries there are in each row, so we can use
3076  // static profile.
3077  RCP<sparse_matrix_type> A_proc0 =
3078  rcp (new sparse_matrix_type (gatherRowMap, gatherNumEntriesPerRow,
3080  if (myRank == rootRank) {
3081  if (debug) {
3082  cerr << "-- Proc 0: Filling gather matrix" << endl;
3083  }
3084  if (debug) {
3085  cerr << "---- Rows: " << Teuchos::toString (myRows) << endl;
3086  }
3087 
3088  // Add Proc 0's matrix entries to the CrsMatrix.
3089  for (size_type i_ = 0; i_ < myNumRows; ++i_) {
3090  size_type i = myRows[i_] - indexBase;
3091 
3092  const size_type curPos = as<size_type> (rowPtr[i]);
3093  const local_ordinal_type curNumEntries = numEntriesPerRow[i];
3094  ArrayView<global_ordinal_type> curColInd =
3095  colInd.view (curPos, curNumEntries);
3096  ArrayView<scalar_type> curValues =
3097  values.view (curPos, curNumEntries);
3098 
3099  // Modify the column indices in place to have the right index base.
3100  for (size_type k = 0; k < curNumEntries; ++k) {
3101  curColInd[k] += indexBase;
3102  }
3103  if (debug) {
3104  cerr << "------ Columns: " << Teuchos::toString (curColInd) << endl;
3105  cerr << "------ Values: " << Teuchos::toString (curValues) << endl;
3106  }
3107  // Avoid constructing empty views of ArrayRCP objects.
3108  if (curNumEntries > 0) {
3109  A_proc0->insertGlobalValues (myRows[i_], curColInd, curValues);
3110  }
3111  }
3112  // Now we can save space by deallocating numEntriesPerRow,
3113  // rowPtr, colInd, and values, since we've already put those
3114  // data in the matrix.
3115  numEntriesPerRow = null;
3116  rowPtr = null;
3117  colInd = null;
3118  values = null;
3119  } // if myRank == rootRank
3120 
3121  RCP<sparse_matrix_type> A;
3122  if (colMap.is_null ()) {
3123  A = rcp (new sparse_matrix_type (rowMap, 0));
3124  } else {
3125  A = rcp (new sparse_matrix_type (rowMap, colMap, 0));
3126  }
3128  export_type exp (gatherRowMap, rowMap);
3129  A->doExport (*A_proc0, exp, INSERT);
3130 
3131  if (callFillComplete) {
3132  A->fillComplete (domainMap, rangeMap);
3133  }
3134 
3135  // We can't get the dimensions of the matrix until after
3136  // fillComplete() is called. Thus, we can't do the sanity
3137  // check (dimensions read from the Matrix Market data,
3138  // vs. dimensions reported by the CrsMatrix) unless the user
3139  // asked us to call fillComplete().
3140  //
3141  // Note that pMatrix->getGlobalNum{Rows,Cols}() does _not_ do
3142  // what one might think it does, so you have to ask the range
3143  // resp. domain map for the number of rows resp. columns.
3144  if (callFillComplete) {
3145  const int numProcs = pComm->getSize ();
3146 
3147  if (extraDebug && debug) {
3148  const global_size_t globalNumRows = rangeMap->getGlobalNumElements ();
3149  const global_size_t globalNumCols = domainMap->getGlobalNumElements ();
3150  if (myRank == rootRank) {
3151  cerr << "-- Matrix is "
3152  << globalNumRows << " x " << globalNumCols
3153  << " with " << A->getGlobalNumEntries()
3154  << " entries, and index base "
3155  << A->getIndexBase() << "." << endl;
3156  }
3157  pComm->barrier ();
3158  for (int p = 0; p < numProcs; ++p) {
3159  if (myRank == p) {
3160  cerr << "-- Proc " << p << " owns "
3161  << A->getNodeNumCols() << " columns, and "
3162  << A->getNodeNumEntries() << " entries." << endl;
3163  }
3164  pComm->barrier ();
3165  }
3166  } // if (extraDebug && debug)
3167  } // if (callFillComplete)
3168 
3169  if (debug && myRank == rootRank) {
3170  cerr << "-- Done creating the CrsMatrix from the Matrix Market data"
3171  << endl;
3172  }
3173  return A;
3174  }
3175 
3205  static RCP<multivector_type>
3206  readDenseFile (const std::string& filename,
3207  const RCP<const comm_type>& comm,
3208  RCP<const map_type>& map,
3209  const bool tolerant=false,
3210  const bool debug=false)
3211  {
3212  std::ifstream in;
3213  if (comm->getRank () == 0) { // Only open the file on Proc 0.
3214  in.open (filename.c_str ()); // Destructor closes safely
3215  }
3216  return readDense (in, comm, map, tolerant, debug);
3217  }
3218 
3220  static RCP<multivector_type>
3221  readDenseFile (const std::string& filename,
3222  const RCP<const comm_type>& comm,
3223  const RCP<node_type>& node,
3224  RCP<const map_type>& map,
3225  const bool tolerant=false,
3226  const bool debug=false)
3227  {
3228  std::ifstream in;
3229  if (comm->getRank () == 0) { // Only open the file on Proc 0.
3230  in.open (filename.c_str ()); // Destructor closes safely
3231  }
3232  return readDense (in, comm, node, map, tolerant, debug);
3233  }
3234 
3264  static RCP<vector_type>
3265  readVectorFile (const std::string& filename,
3266  const RCP<const comm_type>& comm,
3267  RCP<const map_type>& map,
3268  const bool tolerant=false,
3269  const bool debug=false)
3270  {
3271  std::ifstream in;
3272  if (comm->getRank () == 0) { // Only open the file on Proc 0.
3273  in.open (filename.c_str ()); // Destructor closes safely
3274  }
3275  return readVector (in, comm, map, tolerant, debug);
3276  }
3277 
3280  static RCP<vector_type>
3281  readVectorFile (const std::string& filename,
3282  const RCP<const comm_type>& comm,
3283  const RCP<node_type>& node,
3284  RCP<const map_type>& map,
3285  const bool tolerant=false,
3286  const bool debug=false)
3287  {
3288  std::ifstream in;
3289  if (comm->getRank () == 0) { // Only open the file on Proc 0.
3290  in.open (filename.c_str ()); // Destructor closes safely
3291  }
3292  return readVector (in, comm, node, map, tolerant, debug);
3293  }
3294 
3362  static RCP<multivector_type>
3363  readDense (std::istream& in,
3364  const RCP<const comm_type>& comm,
3365  RCP<const map_type>& map,
3366  const bool tolerant=false,
3367  const bool debug=false)
3368  {
3369  return readDense (in, comm, Teuchos::null, map, tolerant, debug);
3370  }
3371 
3373  static RCP<multivector_type>
3374  readDense (std::istream& in,
3375  const RCP<const comm_type>& comm,
3376  const RCP<node_type>& node,
3377  RCP<const map_type>& map,
3378  const bool tolerant=false,
3379  const bool debug=false)
3380  {
3381  Teuchos::RCP<Teuchos::FancyOStream> err =
3382  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
3383  return readDenseImpl<scalar_type> (in, comm, node, map,
3384  err, tolerant, debug);
3385  }
3386 
3388  static RCP<vector_type>
3389  readVector (std::istream& in,
3390  const RCP<const comm_type>& comm,
3391  RCP<const map_type>& map,
3392  const bool tolerant=false,
3393  const bool debug=false)
3394  {
3395  Teuchos::RCP<Teuchos::FancyOStream> err =
3396  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
3397  return readVectorImpl<scalar_type> (in, comm, Teuchos::null, map,
3398  err, tolerant, debug);
3399  }
3400 
3402  static RCP<vector_type>
3403  readVector (std::istream& in,
3404  const RCP<const comm_type>& comm,
3405  const RCP<node_type>& node,
3406  RCP<const map_type>& map,
3407  const bool tolerant=false,
3408  const bool debug=false)
3409  {
3410  Teuchos::RCP<Teuchos::FancyOStream> err =
3411  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
3412  return readVectorImpl<scalar_type> (in, comm, node, map,
3413  err, tolerant, debug);
3414  }
3415 
3437  static RCP<const map_type>
3438  readMapFile (const std::string& filename,
3439  const RCP<const comm_type>& comm,
3440  const RCP<node_type>& node,
3441  const bool tolerant=false,
3442  const bool debug=false)
3443  {
3444  using Teuchos::inOutArg;
3445  using Teuchos::broadcast;
3446  std::ifstream in;
3447 
3448  int success = 1;
3449  if (comm->getRank () == 0) { // Only open the file on Proc 0.
3450  in.open (filename.c_str ()); // Destructor closes safely
3451  if (! in) {
3452  success = 0;
3453  }
3454  }
3455  broadcast<int, int> (*comm, 0, inOutArg (success));
3456  TEUCHOS_TEST_FOR_EXCEPTION(
3457  success == 0, std::runtime_error,
3458  "Tpetra::MatrixMarket::Reader::readMapFile: "
3459  "Failed to read file \"" << filename << "\" on Process 0.");
3460  return readMap (in, comm, node, tolerant, debug);
3461  }
3462 
3463  private:
3464  template<class MultiVectorScalarType>
3465  static RCP<Tpetra::MultiVector<MultiVectorScalarType,
3468  node_type> >
3469  readDenseImpl (std::istream& in,
3470  const RCP<const comm_type>& comm,
3471  const RCP<node_type>& node,
3472  RCP<const map_type>& map,
3473  const Teuchos::RCP<Teuchos::FancyOStream>& err,
3474  const bool tolerant=false,
3475  const bool debug=false)
3476  {
3477  using Teuchos::MatrixMarket::Banner;
3478  using Teuchos::MatrixMarket::checkCommentLine;
3479  using Teuchos::as;
3480  using Teuchos::broadcast;
3481  using Teuchos::outArg;
3482  using std::endl;
3483  typedef MultiVectorScalarType ST;
3484  typedef local_ordinal_type LO;
3485  typedef global_ordinal_type GO;
3486  typedef node_type NT;
3487  typedef Teuchos::ScalarTraits<ST> STS;
3488  typedef typename STS::magnitudeType MT;
3489  typedef Teuchos::ScalarTraits<MT> STM;
3491 
3492  // Rank 0 is the only (MPI) process allowed to read from the
3493  // input stream.
3494  const int myRank = comm->getRank ();
3495 
3496  if (! err.is_null ()) {
3497  err->pushTab ();
3498  }
3499  if (debug) {
3500  *err << myRank << ": readDenseImpl" << endl;
3501  }
3502  if (! err.is_null ()) {
3503  err->pushTab ();
3504  }
3505 
3506  // mfh 17 Feb 2013: It's not strictly necessary that the Comm
3507  // instances be identical and that the Node instances be
3508  // identical. The essential condition is more complicated to
3509  // test and isn't the same for all Node types. Thus, we just
3510  // leave it up to the user.
3511 
3512  // // If map is nonnull, check the precondition that its
3513  // // communicator resp. node equal comm resp. node. Checking
3514  // // now avoids doing a lot of file reading before we detect the
3515  // // violated precondition.
3516  // TEUCHOS_TEST_FOR_EXCEPTION(
3517  // ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
3518  // std::invalid_argument, "If you supply a nonnull Map, the Map's "
3519  // "communicator and node must equal the supplied communicator resp. "
3520  // "node.");
3521 
3522  // Process 0 will read in the matrix dimensions from the file,
3523  // and broadcast them to all ranks in the given communicator.
3524  // There are only 2 dimensions in the matrix, but we use the
3525  // third element of the Tuple to encode the banner's reported
3526  // data type: "real" == 0, "complex" == 1, and "integer" == 0
3527  // (same as "real"). We don't allow pattern matrices (i.e.,
3528  // graphs) since they only make sense for sparse data.
3529  Tuple<GO, 3> dims;
3530  dims[0] = 0;
3531  dims[1] = 0;
3532 
3533  // Current line number in the input stream. Only valid on
3534  // Proc 0. Various calls will modify this depending on the
3535  // number of lines that are read from the input stream.
3536  size_t lineNumber = 1;
3537 
3538  // Capture errors and their messages on Proc 0.
3539  std::ostringstream exMsg;
3540  int localBannerReadSuccess = 1;
3541  int localDimsReadSuccess = 1;
3542 
3543  // Only Proc 0 gets to read matrix data from the input stream.
3544  if (myRank == 0) {
3545  if (debug) {
3546  *err << myRank << ": readDenseImpl: Reading banner line (dense)" << endl;
3547  }
3548 
3549  // The "Banner" tells you whether the input stream
3550  // represents a dense matrix, the symmetry type of the
3551  // matrix, and the type of the data it contains.
3552  RCP<const Banner> pBanner;
3553  try {
3554  pBanner = readBanner (in, lineNumber, tolerant, debug);
3555  } catch (std::exception& e) {
3556  exMsg << e.what ();
3557  localBannerReadSuccess = 0;
3558  }
3559  // Make sure the input stream is the right kind of data.
3560  if (localBannerReadSuccess) {
3561  if (pBanner->matrixType () != "array") {
3562  exMsg << "The Matrix Market file does not contain dense matrix "
3563  "data. Its banner (first) line says that its matrix type is \""
3564  << pBanner->matrixType () << "\", rather that the required "
3565  "\"array\".";
3566  localBannerReadSuccess = 0;
3567  } else if (pBanner->dataType() == "pattern") {
3568  exMsg << "The Matrix Market file's banner (first) "
3569  "line claims that the matrix's data type is \"pattern\". This does "
3570  "not make sense for a dense matrix, yet the file reports the matrix "
3571  "as dense. The only valid data types for a dense matrix are "
3572  "\"real\", \"complex\", and \"integer\".";
3573  localBannerReadSuccess = 0;
3574  } else {
3575  // Encode the data type reported by the Banner as the
3576  // third element of the dimensions Tuple.
3577  dims[2] = encodeDataType (pBanner->dataType ());
3578  }
3579  } // if we successfully read the banner line
3580 
3581  // At this point, we've successfully read the banner line.
3582  // Now read the dimensions line.
3583  if (localBannerReadSuccess) {
3584  if (debug) {
3585  *err << myRank << ": readDenseImpl: Reading dimensions line (dense)" << endl;
3586  }
3587  // Keep reading lines from the input stream until we find
3588  // a non-comment line, or until we run out of lines. The
3589  // latter is an error, since every "array" format Matrix
3590  // Market file must have a dimensions line after the
3591  // banner (even if the matrix has zero rows or columns, or
3592  // zero entries).
3593  std::string line;
3594  bool commentLine = true;
3595 
3596  while (commentLine) {
3597  // Test whether it is even valid to read from the input
3598  // stream wrapping the line.
3599  if (in.eof () || in.fail ()) {
3600  exMsg << "Unable to get array dimensions line (at all) from line "
3601  << lineNumber << " of input stream. The input stream "
3602  << "claims that it is "
3603  << (in.eof() ? "at end-of-file." : "in a failed state.");
3604  localDimsReadSuccess = 0;
3605  } else {
3606  // Try to get the next line from the input stream.
3607  if (getline (in, line)) {
3608  ++lineNumber; // We did actually read a line.
3609  }
3610  // Is the current line a comment line? Ignore start
3611  // and size; they are only useful for reading the
3612  // actual matrix entries. (We could use them here as
3613  // an optimization, but we've chosen not to.)
3614  size_t start = 0, size = 0;
3615  commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
3616  } // whether we failed to read the line at all
3617  } // while the line we just read is a comment line
3618 
3619  //
3620  // Get <numRows> <numCols> from the line we just read.
3621  //
3622  std::istringstream istr (line);
3623 
3624  // Test whether it is even valid to read from the input
3625  // stream wrapping the line.
3626  if (istr.eof () || istr.fail ()) {
3627  exMsg << "Unable to read any data from line " << lineNumber
3628  << " of input; the line should contain the matrix dimensions "
3629  << "\"<numRows> <numCols>\".";
3630  localDimsReadSuccess = 0;
3631  } else { // It's valid to read from the line.
3632  GO theNumRows = 0;
3633  istr >> theNumRows; // Read in the number of rows.
3634  if (istr.fail ()) {
3635  exMsg << "Failed to get number of rows from line "
3636  << lineNumber << " of input; the line should contains the "
3637  << "matrix dimensions \"<numRows> <numCols>\".";
3638  localDimsReadSuccess = 0;
3639  } else { // We successfully read the number of rows
3640  dims[0] = theNumRows; // Save the number of rows
3641  if (istr.eof ()) { // Do we still have data to read?
3642  exMsg << "No more data after number of rows on line "
3643  << lineNumber << " of input; the line should contain the "
3644  << "matrix dimensions \"<numRows> <numCols>\".";
3645  localDimsReadSuccess = 0;
3646  } else { // Still data left to read; read in number of columns.
3647  GO theNumCols = 0;
3648  istr >> theNumCols; // Read in the number of columns
3649  if (istr.fail ()) {
3650  exMsg << "Failed to get number of columns from line "
3651  << lineNumber << " of input; the line should contain "
3652  << "the matrix dimensions \"<numRows> <numCols>\".";
3653  localDimsReadSuccess = 0;
3654  } else { // We successfully read the number of columns
3655  dims[1] = theNumCols; // Save the number of columns
3656  } // if istr.fail ()
3657  } // if istr.eof ()
3658  } // if we read the number of rows
3659  } // if the input stream wrapping the dims line was (in)valid
3660  } // if we successfully read the banner line
3661  } // if (myRank == 0)
3662 
3663  // Broadcast the matrix dimensions, the encoded data type, and
3664  // whether or not Proc 0 succeeded in reading the banner and
3665  // dimensions.
3666  Tuple<GO, 5> bannerDimsReadResult;
3667  if (myRank == 0) {
3668  bannerDimsReadResult[0] = dims[0]; // numRows
3669  bannerDimsReadResult[1] = dims[1]; // numCols
3670  bannerDimsReadResult[2] = dims[2]; // encoded data type
3671  bannerDimsReadResult[3] = localBannerReadSuccess;
3672  bannerDimsReadResult[4] = localDimsReadSuccess;
3673  }
3674  // Broadcast matrix dimensions and the encoded data type from
3675  // Proc 0 to all the MPI processes.
3676  broadcast (*comm, 0, bannerDimsReadResult);
3677 
3678  TEUCHOS_TEST_FOR_EXCEPTION(
3679  bannerDimsReadResult[3] == 0, std::runtime_error,
3680  "Failed to read banner line: " << exMsg.str ());
3681  TEUCHOS_TEST_FOR_EXCEPTION(
3682  bannerDimsReadResult[4] == 0, std::runtime_error,
3683  "Failed to read matrix dimensions line: " << exMsg.str ());
3684  if (myRank != 0) {
3685  dims[0] = bannerDimsReadResult[0];
3686  dims[1] = bannerDimsReadResult[1];
3687  dims[2] = bannerDimsReadResult[2];
3688  }
3689 
3690  // Tpetra objects want the matrix dimensions in these types.
3691  const global_size_t numRows = static_cast<global_size_t> (dims[0]);
3692  const size_t numCols = static_cast<size_t> (dims[1]);
3693 
3694  // Make a "Proc 0 owns everything" Map that we will use to
3695  // read in the multivector entries in the correct order on
3696  // Proc 0. This must be a collective
3697  RCP<const map_type> proc0Map; // "Proc 0 owns everything" Map
3698  if (map.is_null ()) {
3699  // The user didn't supply a Map. Make a contiguous
3700  // distributed Map for them, using the read-in multivector
3701  // dimensions.
3702  if (node.is_null ()) {
3703  map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
3704  } else {
3705  map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
3706  }
3707  const size_t localNumRows = (myRank == 0) ? numRows : 0;
3708  // At this point, map exists and has a nonnull node.
3709  proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
3710  comm, map->getNode ());
3711  }
3712  else { // The user supplied a Map.
3713  proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
3714  }
3715 
3716  // Make a multivector X owned entirely by Proc 0.
3717  RCP<MV> X = createMultiVector<ST, LO, GO, NT> (proc0Map, numCols);
3718 
3719  //
3720  // On Proc 0, read the Matrix Market data from the input
3721  // stream into the multivector X.
3722  //
3723  int localReadDataSuccess = 1;
3724  if (myRank == 0) {
3725  try {
3726  if (debug) {
3727  *err << myRank << ": readDenseImpl: Reading matrix data (dense)"
3728  << endl;
3729  }
3730 
3731  // Make sure that we can get a 1-D view of X.
3732  TEUCHOS_TEST_FOR_EXCEPTION(
3733  ! X->isConstantStride (), std::logic_error,
3734  "Can't get a 1-D view of the entries of the MultiVector X on "
3735  "Process 0, because the stride between the columns of X is not "
3736  "constant. This shouldn't happen because we just created X and "
3737  "haven't filled it in yet. Please report this bug to the Tpetra "
3738  "developers.");
3739 
3740  // Get a writeable 1-D view of the entries of X. Rank 0
3741  // owns all of them. The view will expire at the end of
3742  // scope, so (if necessary) it will be written back to X
3743  // at this time.
3744  ArrayRCP<ST> X_view = X->get1dViewNonConst ();
3745  TEUCHOS_TEST_FOR_EXCEPTION(
3746  as<global_size_t> (X_view.size ()) < numRows * numCols,
3747  std::logic_error,
3748  "The view of X has size " << X_view << " which is not enough to "
3749  "accommodate the expected number of entries numRows*numCols = "
3750  << numRows << "*" << numCols << " = " << numRows*numCols << ". "
3751  "Please report this bug to the Tpetra developers.");
3752  const size_t stride = X->getStride ();
3753 
3754  // The third element of the dimensions Tuple encodes the data
3755  // type reported by the Banner: "real" == 0, "complex" == 1,
3756  // "integer" == 0 (same as "real"), "pattern" == 2. We do not
3757  // allow dense matrices to be pattern matrices, so dims[2] ==
3758  // 0 or 1. We've already checked for this above.
3759  const bool isComplex = (dims[2] == 1);
3760  size_type count = 0, curRow = 0, curCol = 0;
3761 
3762  std::string line;
3763  while (getline (in, line)) {
3764  ++lineNumber;
3765  // Is the current line a comment line? If it's not,
3766  // line.substr(start,size) contains the data.
3767  size_t start = 0, size = 0;
3768  const bool commentLine =
3769  checkCommentLine (line, start, size, lineNumber, tolerant);
3770  if (! commentLine) {
3771  // Make sure we have room in which to put the new matrix
3772  // entry. We check this only after checking for a
3773  // comment line, because there may be one or more
3774  // comment lines at the end of the file. In tolerant
3775  // mode, we simply ignore any extra data.
3776  if (count >= X_view.size()) {
3777  if (tolerant) {
3778  break;
3779  }
3780  else {
3781  TEUCHOS_TEST_FOR_EXCEPTION(
3782  count >= X_view.size(),
3783  std::runtime_error,
3784  "The Matrix Market input stream has more data in it than "
3785  "its metadata reported. Current line number is "
3786  << lineNumber << ".");
3787  }
3788  }
3789 
3790  // mfh 19 Dec 2012: Ignore everything up to the initial
3791  // colon. writeDense() has the option to print out the
3792  // global row index in front of each entry, followed by
3793  // a colon and space.
3794  {
3795  const size_t pos = line.substr (start, size).find (':');
3796  if (pos != std::string::npos) {
3797  start = pos+1;
3798  }
3799  }
3800  std::istringstream istr (line.substr (start, size));
3801  // Does the line contain anything at all? Can we
3802  // safely read from the input stream wrapping the
3803  // line?
3804  if (istr.eof() || istr.fail()) {
3805  // In tolerant mode, simply ignore the line.
3806  if (tolerant) {
3807  break;
3808  }
3809  // We repeat the full test here so the exception
3810  // message is more informative.
3811  TEUCHOS_TEST_FOR_EXCEPTION(
3812  ! tolerant && (istr.eof() || istr.fail()),
3813  std::runtime_error,
3814  "Line " << lineNumber << " of the Matrix Market file is "
3815  "empty, or we cannot read from it for some other reason.");
3816  }
3817  // Current matrix entry to read in.
3818  ST val = STS::zero();
3819  // Real and imaginary parts of the current matrix entry.
3820  // The imaginary part is zero if the matrix is real-valued.
3821  MT real = STM::zero(), imag = STM::zero();
3822 
3823  // isComplex refers to the input stream's data, not to
3824  // the scalar type S. It's OK to read real-valued
3825  // data into a matrix storing complex-valued data; in
3826  // that case, all entries' imaginary parts are zero.
3827  if (isComplex) {
3828  // STS::real() and STS::imag() return a copy of
3829  // their respective components, not a writeable
3830  // reference. Otherwise we could just assign to
3831  // them using the istream extraction operator (>>).
3832  // That's why we have separate magnitude type "real"
3833  // and "imag" variables.
3834 
3835  // Attempt to read the real part of the current entry.
3836  istr >> real;
3837  if (istr.fail()) {
3838  TEUCHOS_TEST_FOR_EXCEPTION(
3839  ! tolerant && istr.eof(), std::runtime_error,
3840  "Failed to get the real part of a complex-valued matrix "
3841  "entry from line " << lineNumber << " of the Matrix Market "
3842  "file.");
3843  // In tolerant mode, just skip bad lines.
3844  if (tolerant) {
3845  break;
3846  }
3847  } else if (istr.eof()) {
3848  TEUCHOS_TEST_FOR_EXCEPTION(
3849  ! tolerant && istr.eof(), std::runtime_error,
3850  "Missing imaginary part of a complex-valued matrix entry "
3851  "on line " << lineNumber << " of the Matrix Market file.");
3852  // In tolerant mode, let any missing imaginary part be 0.
3853  } else {
3854  // Attempt to read the imaginary part of the current
3855  // matrix entry.
3856  istr >> imag;
3857  TEUCHOS_TEST_FOR_EXCEPTION(
3858  ! tolerant && istr.fail(), std::runtime_error,
3859  "Failed to get the imaginary part of a complex-valued "
3860  "matrix entry from line " << lineNumber << " of the "
3861  "Matrix Market file.");
3862  // In tolerant mode, let any missing or corrupted
3863  // imaginary part be 0.
3864  }
3865  } else { // Matrix Market file contains real-valued data.
3866  // Attempt to read the current matrix entry.
3867  istr >> real;
3868  TEUCHOS_TEST_FOR_EXCEPTION(
3869  ! tolerant && istr.fail(), std::runtime_error,
3870  "Failed to get a real-valued matrix entry from line "
3871  << lineNumber << " of the Matrix Market file.");
3872  // In tolerant mode, simply ignore the line if
3873  // we failed to read a matrix entry.
3874  if (istr.fail() && tolerant) {
3875  break;
3876  }
3877  }
3878  // In tolerant mode, we simply let pass through whatever
3879  // data we got.
3880  TEUCHOS_TEST_FOR_EXCEPTION(
3881  ! tolerant && istr.fail(), std::runtime_error,
3882  "Failed to read matrix data from line " << lineNumber
3883  << " of the Matrix Market file.");
3884 
3885  // Assign val = ST(real, imag).
3886  Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
3887 
3888  curRow = count % numRows;
3889  curCol = count / numRows;
3890  X_view[curRow + curCol*stride] = val;
3891  ++count;
3892  } // if not a comment line
3893  } // while there are still lines in the file, get the next one
3894 
3895  TEUCHOS_TEST_FOR_EXCEPTION(
3896  ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
3897  std::runtime_error,
3898  "The Matrix Market metadata reports that the dense matrix is "
3899  << numRows << " x " << numCols << ", and thus has "
3900  << numRows*numCols << " total entries, but we only found " << count
3901  << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
3902  } catch (std::exception& e) {
3903  exMsg << e.what ();
3904  localReadDataSuccess = 0;
3905  }
3906  } // if (myRank == 0)
3907 
3908  if (debug) {
3909  *err << myRank << ": readDenseImpl: done reading data" << endl;
3910  }
3911 
3912  // Synchronize on whether Proc 0 successfully read the data.
3913  int globalReadDataSuccess = localReadDataSuccess;
3914  broadcast (*comm, 0, outArg (globalReadDataSuccess));
3915  TEUCHOS_TEST_FOR_EXCEPTION(
3916  globalReadDataSuccess == 0, std::runtime_error,
3917  "Failed to read the multivector's data: " << exMsg.str ());
3918 
3919  // If there's only one MPI process and the user didn't supply
3920  // a Map (i.e., pMap is null), we're done. Set pMap to the
3921  // Map used to distribute X, and return X.
3922  if (comm->getSize () == 1 && map.is_null ()) {
3923  map = proc0Map;
3924  if (! err.is_null ()) {
3925  err->popTab ();
3926  }
3927  if (debug) {
3928  *err << myRank << ": readDenseImpl: done" << endl;
3929  }
3930  if (! err.is_null ()) {
3931  err->popTab ();
3932  }
3933  return X;
3934  }
3935 
3936  if (debug) {
3937  *err << myRank << ": readDenseImpl: Creating target MV" << endl;
3938  }
3939 
3940  // Make a multivector Y with the distributed map pMap.
3941  RCP<MV> Y = createMultiVector<ST, LO, GO, NT> (map, numCols);
3942 
3943  if (debug) {
3944  *err << myRank << ": readDenseImpl: Creating Export" << endl;
3945  }
3946 
3947  // Make an Export object that will export X to Y. First
3948  // argument is the source map, second argument is the target
3949  // map.
3950  Export<LO, GO, NT> exporter (proc0Map, map, err);
3951 
3952  if (debug) {
3953  *err << myRank << ": readDenseImpl: Exporting" << endl;
3954  }
3955  // Export X into Y.
3956  Y->doExport (*X, exporter, INSERT);
3957 
3958  if (! err.is_null ()) {
3959  err->popTab ();
3960  }
3961  if (debug) {
3962  *err << myRank << ": readDenseImpl: done" << endl;
3963  }
3964  if (! err.is_null ()) {
3965  err->popTab ();
3966  }
3967 
3968  // Y is distributed over all process(es) in the communicator.
3969  return Y;
3970  }
3971 
3972 
3973  template<class VectorScalarType>
3974  static RCP<Tpetra::Vector<VectorScalarType,
3977  node_type> >
3978  readVectorImpl (std::istream& in,
3979  const RCP<const comm_type>& comm,
3980  const RCP<node_type>& node, // allowed to be null
3981  RCP<const map_type>& map,
3982  const Teuchos::RCP<Teuchos::FancyOStream>& err,
3983  const bool tolerant=false,
3984  const bool debug=false)
3985  {
3986  using Teuchos::MatrixMarket::Banner;
3987  using Teuchos::MatrixMarket::checkCommentLine;
3988  using Teuchos::as;
3989  using Teuchos::broadcast;
3990  using Teuchos::outArg;
3991  using std::endl;
3992  typedef VectorScalarType ST;
3993  typedef local_ordinal_type LO;
3994  typedef global_ordinal_type GO;
3995  typedef node_type NT;
3996  typedef Teuchos::ScalarTraits<ST> STS;
3997  typedef typename STS::magnitudeType MT;
3998  typedef Teuchos::ScalarTraits<MT> STM;
3999  typedef Tpetra::Vector<ST, LO, GO, NT> MV;
4000 
4001  // Rank 0 is the only (MPI) process allowed to read from the
4002  // input stream.
4003  const int myRank = comm->getRank ();
4004 
4005  if (! err.is_null ()) {
4006  err->pushTab ();
4007  }
4008  if (debug) {
4009  *err << myRank << ": readVectorImpl" << endl;
4010  }
4011  if (! err.is_null ()) {
4012  err->pushTab ();
4013  }
4014 
4015  // mfh 17 Feb 2013: It's not strictly necessary that the Comm
4016  // instances be identical and that the Node instances be
4017  // identical. The essential condition is more complicated to
4018  // test and isn't the same for all Node types. Thus, we just
4019  // leave it up to the user.
4020 
4021  // // If map is nonnull, check the precondition that its
4022  // // communicator resp. node equal comm resp. node. Checking
4023  // // now avoids doing a lot of file reading before we detect the
4024  // // violated precondition.
4025  // TEUCHOS_TEST_FOR_EXCEPTION(
4026  // ! map.is_null() && (map->getComm() != comm || map->getNode () != node,
4027  // std::invalid_argument, "If you supply a nonnull Map, the Map's "
4028  // "communicator and node must equal the supplied communicator resp. "
4029  // "node.");
4030 
4031  // Process 0 will read in the matrix dimensions from the file,
4032  // and broadcast them to all ranks in the given communicator.
4033  // There are only 2 dimensions in the matrix, but we use the
4034  // third element of the Tuple to encode the banner's reported
4035  // data type: "real" == 0, "complex" == 1, and "integer" == 0
4036  // (same as "real"). We don't allow pattern matrices (i.e.,
4037  // graphs) since they only make sense for sparse data.
4038  Tuple<GO, 3> dims;
4039  dims[0] = 0;
4040  dims[1] = 0;
4041 
4042  // Current line number in the input stream. Only valid on
4043  // Proc 0. Various calls will modify this depending on the
4044  // number of lines that are read from the input stream.
4045  size_t lineNumber = 1;
4046 
4047  // Capture errors and their messages on Proc 0.
4048  std::ostringstream exMsg;
4049  int localBannerReadSuccess = 1;
4050  int localDimsReadSuccess = 1;
4051 
4052  // Only Proc 0 gets to read matrix data from the input stream.
4053  if (myRank == 0) {
4054  if (debug) {
4055  *err << myRank << ": readVectorImpl: Reading banner line (dense)" << endl;
4056  }
4057 
4058  // The "Banner" tells you whether the input stream
4059  // represents a dense matrix, the symmetry type of the
4060  // matrix, and the type of the data it contains.
4061  RCP<const Banner> pBanner;
4062  try {
4063  pBanner = readBanner (in, lineNumber, tolerant, debug);
4064  } catch (std::exception& e) {
4065  exMsg << e.what ();
4066  localBannerReadSuccess = 0;
4067  }
4068  // Make sure the input stream is the right kind of data.
4069  if (localBannerReadSuccess) {
4070  if (pBanner->matrixType () != "array") {
4071  exMsg << "The Matrix Market file does not contain dense matrix "
4072  "data. Its banner (first) line says that its matrix type is \""
4073  << pBanner->matrixType () << "\", rather that the required "
4074  "\"array\".";
4075  localBannerReadSuccess = 0;
4076  } else if (pBanner->dataType() == "pattern") {
4077  exMsg << "The Matrix Market file's banner (first) "
4078  "line claims that the matrix's data type is \"pattern\". This does "
4079  "not make sense for a dense matrix, yet the file reports the matrix "
4080  "as dense. The only valid data types for a dense matrix are "
4081  "\"real\", \"complex\", and \"integer\".";
4082  localBannerReadSuccess = 0;
4083  } else {
4084  // Encode the data type reported by the Banner as the
4085  // third element of the dimensions Tuple.
4086  dims[2] = encodeDataType (pBanner->dataType ());
4087  }
4088  } // if we successfully read the banner line
4089 
4090  // At this point, we've successfully read the banner line.
4091  // Now read the dimensions line.
4092  if (localBannerReadSuccess) {
4093  if (debug) {
4094  *err << myRank << ": readVectorImpl: Reading dimensions line (dense)" << endl;
4095  }
4096  // Keep reading lines from the input stream until we find
4097  // a non-comment line, or until we run out of lines. The
4098  // latter is an error, since every "array" format Matrix
4099  // Market file must have a dimensions line after the
4100  // banner (even if the matrix has zero rows or columns, or
4101  // zero entries).
4102  std::string line;
4103  bool commentLine = true;
4104 
4105  while (commentLine) {
4106  // Test whether it is even valid to read from the input
4107  // stream wrapping the line.
4108  if (in.eof () || in.fail ()) {
4109  exMsg << "Unable to get array dimensions line (at all) from line "
4110  << lineNumber << " of input stream. The input stream "
4111  << "claims that it is "
4112  << (in.eof() ? "at end-of-file." : "in a failed state.");
4113  localDimsReadSuccess = 0;
4114  } else {
4115  // Try to get the next line from the input stream.
4116  if (getline (in, line)) {
4117  ++lineNumber; // We did actually read a line.
4118  }
4119  // Is the current line a comment line? Ignore start
4120  // and size; they are only useful for reading the
4121  // actual matrix entries. (We could use them here as
4122  // an optimization, but we've chosen not to.)
4123  size_t start = 0, size = 0;
4124  commentLine = checkCommentLine (line, start, size, lineNumber, tolerant);
4125  } // whether we failed to read the line at all
4126  } // while the line we just read is a comment line
4127 
4128  //
4129  // Get <numRows> <numCols> from the line we just read.
4130  //
4131  std::istringstream istr (line);
4132 
4133  // Test whether it is even valid to read from the input
4134  // stream wrapping the line.
4135  if (istr.eof () || istr.fail ()) {
4136  exMsg << "Unable to read any data from line " << lineNumber
4137  << " of input; the line should contain the matrix dimensions "
4138  << "\"<numRows> <numCols>\".";
4139  localDimsReadSuccess = 0;
4140  } else { // It's valid to read from the line.
4141  GO theNumRows = 0;
4142  istr >> theNumRows; // Read in the number of rows.
4143  if (istr.fail ()) {
4144  exMsg << "Failed to get number of rows from line "
4145  << lineNumber << " of input; the line should contains the "
4146  << "matrix dimensions \"<numRows> <numCols>\".";
4147  localDimsReadSuccess = 0;
4148  } else { // We successfully read the number of rows
4149  dims[0] = theNumRows; // Save the number of rows
4150  if (istr.eof ()) { // Do we still have data to read?
4151  exMsg << "No more data after number of rows on line "
4152  << lineNumber << " of input; the line should contain the "
4153  << "matrix dimensions \"<numRows> <numCols>\".";
4154  localDimsReadSuccess = 0;
4155  } else { // Still data left to read; read in number of columns.
4156  GO theNumCols = 0;
4157  istr >> theNumCols; // Read in the number of columns
4158  if (istr.fail ()) {
4159  exMsg << "Failed to get number of columns from line "
4160  << lineNumber << " of input; the line should contain "
4161  << "the matrix dimensions \"<numRows> <numCols>\".";
4162  localDimsReadSuccess = 0;
4163  } else { // We successfully read the number of columns
4164  dims[1] = theNumCols; // Save the number of columns
4165  } // if istr.fail ()
4166  } // if istr.eof ()
4167  } // if we read the number of rows
4168  } // if the input stream wrapping the dims line was (in)valid
4169  } // if we successfully read the banner line
4170  } // if (myRank == 0)
4171 
4172  // Check if file has a Vector
4173  if (dims[1]!=1) {
4174  exMsg << "File does not contain a 1D Vector";
4175  localDimsReadSuccess = 0;
4176  }
4177 
4178  // Broadcast the matrix dimensions, the encoded data type, and
4179  // whether or not Proc 0 succeeded in reading the banner and
4180  // dimensions.
4181  Tuple<GO, 5> bannerDimsReadResult;
4182  if (myRank == 0) {
4183  bannerDimsReadResult[0] = dims[0]; // numRows
4184  bannerDimsReadResult[1] = dims[1]; // numCols
4185  bannerDimsReadResult[2] = dims[2]; // encoded data type
4186  bannerDimsReadResult[3] = localBannerReadSuccess;
4187  bannerDimsReadResult[4] = localDimsReadSuccess;
4188  }
4189 
4190  // Broadcast matrix dimensions and the encoded data type from
4191  // Proc 0 to all the MPI processes.
4192  broadcast (*comm, 0, bannerDimsReadResult);
4193 
4194  TEUCHOS_TEST_FOR_EXCEPTION(
4195  bannerDimsReadResult[3] == 0, std::runtime_error,
4196  "Failed to read banner line: " << exMsg.str ());
4197  TEUCHOS_TEST_FOR_EXCEPTION(
4198  bannerDimsReadResult[4] == 0, std::runtime_error,
4199  "Failed to read matrix dimensions line: " << exMsg.str ());
4200  if (myRank != 0) {
4201  dims[0] = bannerDimsReadResult[0];
4202  dims[1] = bannerDimsReadResult[1];
4203  dims[2] = bannerDimsReadResult[2];
4204  }
4205 
4206  // Tpetra objects want the matrix dimensions in these types.
4207  const global_size_t numRows = static_cast<global_size_t> (dims[0]);
4208  const size_t numCols = static_cast<size_t> (dims[1]);
4209 
4210  // Make a "Proc 0 owns everything" Map that we will use to
4211  // read in the multivector entries in the correct order on
4212  // Proc 0. This must be a collective
4213  RCP<const map_type> proc0Map; // "Proc 0 owns everything" Map
4214  if (map.is_null ()) {
4215  // The user didn't supply a Map. Make a contiguous
4216  // distributed Map for them, using the read-in multivector
4217  // dimensions.
4218  if (node.is_null ()) {
4219  map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm);
4220  } else {
4221  map = createUniformContigMapWithNode<LO, GO, NT> (numRows, comm, node);
4222  }
4223  const size_t localNumRows = (myRank == 0) ? numRows : 0;
4224  // At this point, map exists and has a nonnull node.
4225  proc0Map = createContigMapWithNode<LO, GO, NT> (numRows, localNumRows,
4226  comm, map->getNode ());
4227  }
4228  else { // The user supplied a Map.
4229  proc0Map = Details::computeGatherMap<map_type> (map, err, debug);
4230  }
4231 
4232  // Make a multivector X owned entirely by Proc 0.
4233  RCP<MV> X = createVector<ST, LO, GO, NT> (proc0Map);
4234 
4235  //
4236  // On Proc 0, read the Matrix Market data from the input
4237  // stream into the multivector X.
4238  //
4239  int localReadDataSuccess = 1;
4240  if (myRank == 0) {
4241  try {
4242  if (debug) {
4243  *err << myRank << ": readVectorImpl: Reading matrix data (dense)"
4244  << endl;
4245  }
4246 
4247  // Make sure that we can get a 1-D view of X.
4248  TEUCHOS_TEST_FOR_EXCEPTION(
4249  ! X->isConstantStride (), std::logic_error,
4250  "Can't get a 1-D view of the entries of the MultiVector X on "
4251  "Process 0, because the stride between the columns of X is not "
4252  "constant. This shouldn't happen because we just created X and "
4253  "haven't filled it in yet. Please report this bug to the Tpetra "
4254  "developers.");
4255 
4256  // Get a writeable 1-D view of the entries of X. Rank 0
4257  // owns all of them. The view will expire at the end of
4258  // scope, so (if necessary) it will be written back to X
4259  // at this time.
4260  ArrayRCP<ST> X_view = X->get1dViewNonConst ();
4261  TEUCHOS_TEST_FOR_EXCEPTION(
4262  as<global_size_t> (X_view.size ()) < numRows * numCols,
4263  std::logic_error,
4264  "The view of X has size " << X_view << " which is not enough to "
4265  "accommodate the expected number of entries numRows*numCols = "
4266  << numRows << "*" << numCols << " = " << numRows*numCols << ". "
4267  "Please report this bug to the Tpetra developers.");
4268  const size_t stride = X->getStride ();
4269 
4270  // The third element of the dimensions Tuple encodes the data
4271  // type reported by the Banner: "real" == 0, "complex" == 1,
4272  // "integer" == 0 (same as "real"), "pattern" == 2. We do not
4273  // allow dense matrices to be pattern matrices, so dims[2] ==
4274  // 0 or 1. We've already checked for this above.
4275  const bool isComplex = (dims[2] == 1);
4276  size_type count = 0, curRow = 0, curCol = 0;
4277 
4278  std::string line;
4279  while (getline (in, line)) {
4280  ++lineNumber;
4281  // Is the current line a comment line? If it's not,
4282  // line.substr(start,size) contains the data.
4283  size_t start = 0, size = 0;
4284  const bool commentLine =
4285  checkCommentLine (line, start, size, lineNumber, tolerant);
4286  if (! commentLine) {
4287  // Make sure we have room in which to put the new matrix
4288  // entry. We check this only after checking for a
4289  // comment line, because there may be one or more
4290  // comment lines at the end of the file. In tolerant
4291  // mode, we simply ignore any extra data.
4292  if (count >= X_view.size()) {
4293  if (tolerant) {
4294  break;
4295  }
4296  else {
4297  TEUCHOS_TEST_FOR_EXCEPTION(
4298  count >= X_view.size(),
4299  std::runtime_error,
4300  "The Matrix Market input stream has more data in it than "
4301  "its metadata reported. Current line number is "
4302  << lineNumber << ".");
4303  }
4304  }
4305 
4306  // mfh 19 Dec 2012: Ignore everything up to the initial
4307  // colon. writeDense() has the option to print out the
4308  // global row index in front of each entry, followed by
4309  // a colon and space.
4310  {
4311  const size_t pos = line.substr (start, size).find (':');
4312  if (pos != std::string::npos) {
4313  start = pos+1;
4314  }
4315  }
4316  std::istringstream istr (line.substr (start, size));
4317  // Does the line contain anything at all? Can we
4318  // safely read from the input stream wrapping the
4319  // line?
4320  if (istr.eof() || istr.fail()) {
4321  // In tolerant mode, simply ignore the line.
4322  if (tolerant) {
4323  break;
4324  }
4325  // We repeat the full test here so the exception
4326  // message is more informative.
4327  TEUCHOS_TEST_FOR_EXCEPTION(
4328  ! tolerant && (istr.eof() || istr.fail()),
4329  std::runtime_error,
4330  "Line " << lineNumber << " of the Matrix Market file is "
4331  "empty, or we cannot read from it for some other reason.");
4332  }
4333  // Current matrix entry to read in.
4334  ST val = STS::zero();
4335  // Real and imaginary parts of the current matrix entry.
4336  // The imaginary part is zero if the matrix is real-valued.
4337  MT real = STM::zero(), imag = STM::zero();
4338 
4339  // isComplex refers to the input stream's data, not to
4340  // the scalar type S. It's OK to read real-valued
4341  // data into a matrix storing complex-valued data; in
4342  // that case, all entries' imaginary parts are zero.
4343  if (isComplex) {
4344  // STS::real() and STS::imag() return a copy of
4345  // their respective components, not a writeable
4346  // reference. Otherwise we could just assign to
4347  // them using the istream extraction operator (>>).
4348  // That's why we have separate magnitude type "real"
4349  // and "imag" variables.
4350 
4351  // Attempt to read the real part of the current entry.
4352  istr >> real;
4353  if (istr.fail()) {
4354  TEUCHOS_TEST_FOR_EXCEPTION(
4355  ! tolerant && istr.eof(), std::runtime_error,
4356  "Failed to get the real part of a complex-valued matrix "
4357  "entry from line " << lineNumber << " of the Matrix Market "
4358  "file.");
4359  // In tolerant mode, just skip bad lines.
4360  if (tolerant) {
4361  break;
4362  }
4363  } else if (istr.eof()) {
4364  TEUCHOS_TEST_FOR_EXCEPTION(
4365  ! tolerant && istr.eof(), std::runtime_error,
4366  "Missing imaginary part of a complex-valued matrix entry "
4367  "on line " << lineNumber << " of the Matrix Market file.");
4368  // In tolerant mode, let any missing imaginary part be 0.
4369  } else {
4370  // Attempt to read the imaginary part of the current
4371  // matrix entry.
4372  istr >> imag;
4373  TEUCHOS_TEST_FOR_EXCEPTION(
4374  ! tolerant && istr.fail(), std::runtime_error,
4375  "Failed to get the imaginary part of a complex-valued "
4376  "matrix entry from line " << lineNumber << " of the "
4377  "Matrix Market file.");
4378  // In tolerant mode, let any missing or corrupted
4379  // imaginary part be 0.
4380  }
4381  } else { // Matrix Market file contains real-valued data.
4382  // Attempt to read the current matrix entry.
4383  istr >> real;
4384  TEUCHOS_TEST_FOR_EXCEPTION(
4385  ! tolerant && istr.fail(), std::runtime_error,
4386  "Failed to get a real-valued matrix entry from line "
4387  << lineNumber << " of the Matrix Market file.");
4388  // In tolerant mode, simply ignore the line if
4389  // we failed to read a matrix entry.
4390  if (istr.fail() && tolerant) {
4391  break;
4392  }
4393  }
4394  // In tolerant mode, we simply let pass through whatever
4395  // data we got.
4396  TEUCHOS_TEST_FOR_EXCEPTION(
4397  ! tolerant && istr.fail(), std::runtime_error,
4398  "Failed to read matrix data from line " << lineNumber
4399  << " of the Matrix Market file.");
4400 
4401  // Assign val = ST(real, imag).
4402  Teuchos::MatrixMarket::details::assignScalar<ST> (val, real, imag);
4403 
4404  curRow = count % numRows;
4405  curCol = count / numRows;
4406  X_view[curRow + curCol*stride] = val;
4407  ++count;
4408  } // if not a comment line
4409  } // while there are still lines in the file, get the next one
4410 
4411  TEUCHOS_TEST_FOR_EXCEPTION(
4412  ! tolerant && static_cast<global_size_t> (count) < numRows * numCols,
4413  std::runtime_error,
4414  "The Matrix Market metadata reports that the dense matrix is "
4415  << numRows << " x " << numCols << ", and thus has "
4416  << numRows*numCols << " total entries, but we only found " << count
4417  << " entr" << (count == 1 ? "y" : "ies") << " in the file.");
4418  } catch (std::exception& e) {
4419  exMsg << e.what ();
4420  localReadDataSuccess = 0;
4421  }
4422  } // if (myRank == 0)
4423 
4424  if (debug) {
4425  *err << myRank << ": readVectorImpl: done reading data" << endl;
4426  }
4427 
4428  // Synchronize on whether Proc 0 successfully read the data.
4429  int globalReadDataSuccess = localReadDataSuccess;
4430  broadcast (*comm, 0, outArg (globalReadDataSuccess));
4431  TEUCHOS_TEST_FOR_EXCEPTION(
4432  globalReadDataSuccess == 0, std::runtime_error,
4433  "Failed to read the multivector's data: " << exMsg.str ());
4434 
4435  // If there's only one MPI process and the user didn't supply
4436  // a Map (i.e., pMap is null), we're done. Set pMap to the
4437  // Map used to distribute X, and return X.
4438  if (comm->getSize () == 1 && map.is_null ()) {
4439  map = proc0Map;
4440  if (! err.is_null ()) {
4441  err->popTab ();
4442  }
4443  if (debug) {
4444  *err << myRank << ": readVectorImpl: done" << endl;
4445  }
4446  if (! err.is_null ()) {
4447  err->popTab ();
4448  }
4449  return X;
4450  }
4451 
4452  if (debug) {
4453  *err << myRank << ": readVectorImpl: Creating target MV" << endl;
4454  }
4455 
4456  // Make a multivector Y with the distributed map pMap.
4457  RCP<MV> Y = createVector<ST, LO, GO, NT> (map);
4458 
4459  if (debug) {
4460  *err << myRank << ": readVectorImpl: Creating Export" << endl;
4461  }
4462 
4463  // Make an Export object that will export X to Y. First
4464  // argument is the source map, second argument is the target
4465  // map.
4466  Export<LO, GO, NT> exporter (proc0Map, map, err);
4467 
4468  if (debug) {
4469  *err << myRank << ": readVectorImpl: Exporting" << endl;
4470  }
4471  // Export X into Y.
4472  Y->doExport (*X, exporter, INSERT);
4473 
4474  if (! err.is_null ()) {
4475  err->popTab ();
4476  }
4477  if (debug) {
4478  *err << myRank << ": readVectorImpl: done" << endl;
4479  }
4480  if (! err.is_null ()) {
4481  err->popTab ();
4482  }
4483 
4484  // Y is distributed over all process(es) in the communicator.
4485  return Y;
4486  }
4487 
4488  public:
4509  static RCP<const map_type>
4510  readMap (std::istream& in,
4511  const RCP<const comm_type>& comm,
4512  const RCP<node_type>& node,
4513  const bool tolerant=false,
4514  const bool debug=false)
4515  {
4516  Teuchos::RCP<Teuchos::FancyOStream> err =
4517  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
4518  return readMap (in, comm, node, err, tolerant, debug);
4519  }
4520 
4547  static RCP<const map_type>
4548  readMap (std::istream& in,
4549  const RCP<const comm_type>& comm,
4550  const RCP<node_type>& node,
4551  const Teuchos::RCP<Teuchos::FancyOStream>& err,
4552  const bool tolerant=false,
4553  const bool debug=false)
4554  {
4555  using Teuchos::arcp;
4556  using Teuchos::Array;
4557  using Teuchos::ArrayRCP;
4558  using Teuchos::as;
4559  using Teuchos::broadcast;
4560  using Teuchos::Comm;
4561  using Teuchos::CommRequest;
4562  using Teuchos::inOutArg;
4563  using Teuchos::ireceive;
4564  using Teuchos::outArg;
4565  using Teuchos::receive;
4566  using Teuchos::reduceAll;
4567  using Teuchos::REDUCE_MIN;
4568  using Teuchos::isend;
4569  using Teuchos::SerialComm;
4570  using Teuchos::toString;
4571  using Teuchos::wait;
4572  using std::endl;
4573  typedef Tpetra::global_size_t GST;
4574  typedef ptrdiff_t int_type; // Can hold int and GO
4575  typedef local_ordinal_type LO;
4576  typedef global_ordinal_type GO;
4577  typedef node_type NT;
4579 
4580  const int numProcs = comm->getSize ();
4581  const int myRank = comm->getRank ();
4582 
4583  if (err.is_null ()) {
4584  err->pushTab ();
4585  }
4586  if (debug) {
4587  std::ostringstream os;
4588  os << myRank << ": readMap: " << endl;
4589  *err << os.str ();
4590  }
4591  if (err.is_null ()) {
4592  err->pushTab ();
4593  }
4594 
4595  // Tag for receive-size / send-size messages. writeMap used
4596  // tags 1337 and 1338; we count up from there.
4597  const int sizeTag = 1339;
4598  // Tag for receive-data / send-data messages.
4599  const int dataTag = 1340;
4600 
4601  // These are for sends on Process 0, and for receives on all
4602  // other processes. sizeReq is for the {receive,send}-size
4603  // message, and dataReq is for the message containing the
4604  // actual GIDs to belong to the receiving process.
4605  RCP<CommRequest<int> > sizeReq;
4606  RCP<CommRequest<int> > dataReq;
4607 
4608  // Each process will have to receive the number of GIDs to
4609  // expect. Thus, we can post the receives now, and cancel
4610  // them if something should go wrong in the meantime.
4611  ArrayRCP<int_type> numGidsToRecv (1);
4612  numGidsToRecv[0] = 0;
4613  if (myRank != 0) {
4614  sizeReq = ireceive<int, int_type> (numGidsToRecv, 0, sizeTag, *comm);
4615  }
4616 
4617  int readSuccess = 1;
4618  std::ostringstream exMsg;
4619  RCP<MV> data; // Will only be valid on Proc 0
4620  if (myRank == 0) {
4621  // If we want to reuse readDenseImpl, we have to make a
4622  // communicator that only contains Proc 0. Otherwise,
4623  // readDenseImpl will redistribute the data to all
4624  // processes. While we eventually want that, neither we nor
4625  // readDenseImpl know the correct Map to use at the moment.
4626  // That depends on the second column of the multivector.
4627  RCP<const Comm<int> > proc0Comm (new SerialComm<int> ());
4628  try {
4629  RCP<const map_type> dataMap;
4630  // This is currently the only place where we use the
4631  // 'tolerant' argument. Later, if we want to be clever,
4632  // we could have tolerant mode allow PIDs out of order.
4633  data = readDenseImpl<GO> (in, proc0Comm, node, dataMap,
4634  err, tolerant, debug);
4635  (void) dataMap; // Silence "unused" warnings
4636  if (data.is_null ()) {
4637  readSuccess = 0;
4638  exMsg << "readDenseImpl() returned null." << endl;
4639  }
4640  } catch (std::exception& e) {
4641  readSuccess = 0;
4642  exMsg << e.what () << endl;
4643  }
4644  }
4645 
4646  // Map from PID to all the GIDs for that PID.
4647  // Only populated on Process 0.
4648  std::map<int, Array<GO> > pid2gids;
4649 
4650  // The index base must be the global minimum GID.
4651  // We will compute this on Process 0 and broadcast,
4652  // so that all processes can set up the Map.
4653  int_type globalNumGIDs = 0;
4654 
4655  // The index base must be the global minimum GID.
4656  // We will compute this on Process 0 and broadcast,
4657  // so that all processes can set up the Map.
4658  GO indexBase = 0;
4659 
4660  // Process 0: If the above read of the MultiVector succeeded,
4661  // extract the GIDs and PIDs into pid2gids, and find the
4662  // global min GID.
4663  if (myRank == 0 && readSuccess == 1) {
4664  if (data->getNumVectors () == 2) { // Map format 1.0
4665  ArrayRCP<const GO> GIDs = data->getData (0);
4666  ArrayRCP<const GO> PIDs = data->getData (1); // convert to int
4667  globalNumGIDs = GIDs.size ();
4668 
4669  // Start computing the global min GID, while collecting
4670  // the GIDs for each PID.
4671  if (globalNumGIDs > 0) {
4672  const int pid = static_cast<int> (PIDs[0]);
4673 
4674  if (pid < 0 || pid >= numProcs) {
4675  readSuccess = 0;
4676  exMsg << "Tpetra::MatrixMarket::readMap: "
4677  << "Encountered invalid PID " << pid << "." << endl;
4678  }
4679  else {
4680  const GO gid = GIDs[0];
4681  pid2gids[pid].push_back (gid);
4682  indexBase = gid; // the current min GID
4683  }
4684  }
4685  if (readSuccess == 1) {
4686  // Collect the rest of the GIDs for each PID, and compute
4687  // the global min GID.
4688  for (size_type k = 1; k < globalNumGIDs; ++k) {
4689  const int pid = static_cast<int> (PIDs[k]);
4690  if (pid < 0 || pid >= numProcs) {
4691  readSuccess = 0;
4692  exMsg << "Tpetra::MatrixMarket::readMap: "
4693  << "Encountered invalid PID " << pid << "." << endl;
4694  }
4695  else {
4696  const int_type gid = GIDs[k];
4697  pid2gids[pid].push_back (gid);
4698  if (gid < indexBase) {
4699  indexBase = gid; // the current min GID
4700  }
4701  }
4702  }
4703  }
4704  }
4705  else if (data->getNumVectors () == 1) { // Map format 2.0
4706  if (data->getGlobalLength () % 2 != static_cast<GST> (0)) {
4707  readSuccess = 0;
4708  exMsg << "Tpetra::MatrixMarket::readMap: Input data has the "
4709  "wrong format (for Map format 2.0). The global number of rows "
4710  "in the MultiVector must be even (divisible by 2)." << endl;
4711  }
4712  else {
4713  ArrayRCP<const GO> theData = data->getData (0);
4714  globalNumGIDs = static_cast<GO> (data->getGlobalLength ()) /
4715  static_cast<GO> (2);
4716 
4717  // Start computing the global min GID, while
4718  // collecting the GIDs for each PID.
4719  if (globalNumGIDs > 0) {
4720  const int pid = static_cast<int> (theData[1]);
4721  if (pid < 0 || pid >= numProcs) {
4722  readSuccess = 0;
4723  exMsg << "Tpetra::MatrixMarket::readMap: "
4724  << "Encountered invalid PID " << pid << "." << endl;
4725  }
4726  else {
4727  const GO gid = theData[0];
4728  pid2gids[pid].push_back (gid);
4729  indexBase = gid; // the current min GID
4730  }
4731  }
4732  // Collect the rest of the GIDs for each PID, and
4733  // compute the global min GID.
4734  for (int_type k = 1; k < globalNumGIDs; ++k) {
4735  const int pid = static_cast<int> (theData[2*k + 1]);
4736  if (pid < 0 || pid >= numProcs) {
4737  readSuccess = 0;
4738  exMsg << "Tpetra::MatrixMarket::readMap: "
4739  << "Encountered invalid PID " << pid << "." << endl;
4740  }
4741  else {
4742  const GO gid = theData[2*k];
4743  pid2gids[pid].push_back (gid);
4744  if (gid < indexBase) {
4745  indexBase = gid; // the current min GID
4746  }
4747  }
4748  } // for each GID
4749  } // if the amount of data is correct
4750  }
4751  else {
4752  readSuccess = 0;
4753  exMsg << "Tpetra::MatrixMarket::readMap: Input data must have "
4754  "either 1 column (for the new Map format 2.0) or 2 columns (for "
4755  "the old Map format 1.0).";
4756  }
4757  } // myRank is zero
4758 
4759  // Broadcast the indexBase, the global number of GIDs, and the
4760  // current success status. Use int_type for all of these.
4761  {
4762  int_type readResults[3];
4763  readResults[0] = static_cast<int_type> (indexBase);
4764  readResults[1] = static_cast<int_type> (globalNumGIDs);
4765  readResults[2] = static_cast<int_type> (readSuccess);
4766  broadcast<int, int_type> (*comm, 0, 3, readResults);
4767 
4768  indexBase = static_cast<GO> (readResults[0]);
4769  globalNumGIDs = static_cast<int_type> (readResults[1]);
4770  readSuccess = static_cast<int> (readResults[2]);
4771  }
4772 
4773  // Unwinding the stack will invoke sizeReq's destructor, which
4774  // will cancel the receive-size request on all processes that
4775  // posted it.
4776  TEUCHOS_TEST_FOR_EXCEPTION(
4777  readSuccess != 1, std::runtime_error,
4778  "Tpetra::MatrixMarket::readMap: Reading the Map failed with the "
4779  "following exception message: " << exMsg.str ());
4780 
4781  if (myRank == 0) {
4782  // Proc 0: Send each process' number of GIDs to that process.
4783  for (int p = 1; p < numProcs; ++p) {
4784  ArrayRCP<int_type> numGidsToSend (1);
4785 
4786  typename std::map<int, Array<GO> >::const_iterator it = pid2gids.find (p);
4787  if (it == pid2gids.end ()) {
4788  numGidsToSend[0] = 0;
4789  } else {
4790  numGidsToSend[0] = it->second.size ();
4791  }
4792  sizeReq = isend<int, int_type> (numGidsToSend, p, sizeTag, *comm);
4793  wait<int> (*comm, outArg (sizeReq));
4794  }
4795  }
4796  else {
4797  // Wait on the receive-size message to finish.
4798  wait<int> (*comm, outArg (sizeReq));
4799  }
4800 
4801  // Allocate / get the array for my GIDs.
4802  // Only Process 0 will have its actual GIDs at this point.
4803  ArrayRCP<GO> myGids;
4804  int_type myNumGids = 0;
4805  if (myRank == 0) {
4806  GO* myGidsRaw = NULL;
4807 
4808  typename std::map<int, Array<GO> >::iterator it = pid2gids.find (0);
4809  if (it != pid2gids.end ()) {
4810  myGidsRaw = it->second.getRawPtr ();
4811  myNumGids = it->second.size ();
4812  // Nonowning ArrayRCP just views the Array.
4813  myGids = arcp<GO> (myGidsRaw, 0, myNumGids, false);
4814  }
4815  }
4816  else { // myRank != 0
4817  myNumGids = numGidsToRecv[0];
4818  myGids = arcp<GO> (myNumGids);
4819  }
4820 
4821  if (myRank != 0) {
4822  // Post receive for data, now that we know how much data we
4823  // will receive. Only post receive if my process actually
4824  // has nonzero GIDs.
4825  if (myNumGids > 0) {
4826  dataReq = ireceive<int, GO> (myGids, 0, dataTag, *comm);
4827  }
4828  }
4829 
4830  for (int p = 1; p < numProcs; ++p) {
4831  if (myRank == 0) {
4832  ArrayRCP<GO> sendGids; // to send to Process p
4833  GO* sendGidsRaw = NULL;
4834  int_type numSendGids = 0;
4835 
4836  typename std::map<int, Array<GO> >::iterator it = pid2gids.find (p);
4837  if (it != pid2gids.end ()) {
4838  numSendGids = it->second.size ();
4839  sendGidsRaw = it->second.getRawPtr ();
4840  sendGids = arcp<GO> (sendGidsRaw, 0, numSendGids, false);
4841  }
4842  // Only send if that process actually has nonzero GIDs.
4843  if (numSendGids > 0) {
4844  dataReq = isend<int, GO> (sendGids, p, dataTag, *comm);
4845  }
4846  wait<int> (*comm, outArg (dataReq));
4847  }
4848  else if (myRank == p) {
4849  // Wait on my receive of GIDs to finish.
4850  wait<int> (*comm, outArg (dataReq));
4851  }
4852  } // for each process rank p in 1, 2, ..., numProcs-1
4853 
4854  if (debug) {
4855  std::ostringstream os;
4856  os << myRank << ": readMap: creating Map" << endl;
4857  *err << os.str ();
4858  }
4859  const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();
4860  RCP<const map_type> newMap =
4861  rcp (new map_type (INVALID, myGids (), indexBase, comm, node));
4862 
4863  if (err.is_null ()) {
4864  err->popTab ();
4865  }
4866  if (debug) {
4867  std::ostringstream os;
4868  os << myRank << ": readMap: done" << endl;
4869  *err << os.str ();
4870  }
4871  if (err.is_null ()) {
4872  err->popTab ();
4873  }
4874  return newMap;
4875  }
4876 
4877  private:
4878 
4889  static int
4890  encodeDataType (const std::string& dataType)
4891  {
4892  if (dataType == "real" || dataType == "integer") {
4893  return 0;
4894  } else if (dataType == "complex") {
4895  return 1;
4896  } else if (dataType == "pattern") {
4897  return 2;
4898  } else {
4899  // We should never get here, since Banner validates the
4900  // reported data type and ensures it is one of the accepted
4901  // values.
4902  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
4903  "Unrecognized Matrix Market data type \"" << dataType
4904  << "\". We should never get here. "
4905  "Please report this bug to the Tpetra developers.");
4906  }
4907  }
4908  };
4909 
4938  template<class SparseMatrixType>
4939  class Writer {
4940  public:
4942  typedef SparseMatrixType sparse_matrix_type;
4943  typedef Teuchos::RCP<sparse_matrix_type> sparse_matrix_ptr;
4944 
4946  typedef typename SparseMatrixType::scalar_type scalar_type;
4948  typedef typename SparseMatrixType::local_ordinal_type local_ordinal_type;
4954  typedef typename SparseMatrixType::global_ordinal_type global_ordinal_type;
4956  typedef typename SparseMatrixType::node_type node_type;
4957 
4959  typedef MultiVector<scalar_type,
4962  node_type> multivector_type;
4965 
4968 
5000  static void
5001  writeSparseFile (const std::string& filename,
5002  const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5003  const std::string& matrixName,
5004  const std::string& matrixDescription,
5005  const bool debug=false)
5006  {
5007  TEUCHOS_TEST_FOR_EXCEPTION(
5008  pMatrix.is_null (), std::invalid_argument,
5009  "The input matrix is null.");
5010  Teuchos::RCP<const Teuchos::Comm<int> > comm = pMatrix->getComm ();
5011  TEUCHOS_TEST_FOR_EXCEPTION(
5012  comm.is_null (), std::invalid_argument,
5013  "The input matrix's communicator (Teuchos::Comm object) is null.");
5014  const int myRank = comm->getRank ();
5015  std::ofstream out;
5016 
5017  // Only open the file on Rank 0.
5018  if (myRank == 0) {
5019  out.open (filename.c_str ());
5020  }
5021  writeSparse (out, pMatrix, matrixName, matrixDescription, debug);
5022  // We can rely on the destructor of the output stream to close
5023  // the file on scope exit, even if writeSparse() throws an
5024  // exception.
5025  }
5026 
5047  static void
5048  writeSparseFile (const std::string& filename,
5049  const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5050  const bool debug=false)
5051  {
5052  writeSparseFile (filename, pMatrix, "", "", debug);
5053  }
5054 
5085  static void
5086  writeSparse (std::ostream& out,
5087  const Teuchos::RCP<const sparse_matrix_type>& pMatrix,
5088  const std::string& matrixName,
5089  const std::string& matrixDescription,
5090  const bool debug=false)
5091  {
5092  using Teuchos::Comm;
5093  using Teuchos::FancyOStream;
5094  using Teuchos::getFancyOStream;
5095  using Teuchos::null;
5096  using Teuchos::RCP;
5097  using Teuchos::rcpFromRef;
5098  using std::cerr;
5099  using std::endl;
5100  typedef scalar_type ST;
5101  typedef local_ordinal_type LO;
5102  typedef global_ordinal_type GO;
5103  typedef typename Teuchos::ScalarTraits<ST> STS;
5104  typedef typename ArrayView<const LO>::const_iterator lo_iter;
5105  typedef typename ArrayView<const GO>::const_iterator go_iter;
5106  typedef typename ArrayView<const ST>::const_iterator st_iter;
5107 
5108  TEUCHOS_TEST_FOR_EXCEPTION(
5109  pMatrix.is_null (), std::invalid_argument,
5110  "The input matrix is null.");
5111 
5112  // Make the output stream write floating-point numbers in
5113  // scientific notation. It will politely put the output
5114  // stream back to its state on input, when this scope
5115  // terminates.
5116  Teuchos::MatrixMarket::details::SetScientific<ST> sci (out);
5117 
5118  // Get the matrix's communicator.
5119  RCP<const Comm<int> > comm = pMatrix->getComm ();
5120  TEUCHOS_TEST_FOR_EXCEPTION(
5121  comm.is_null (), std::invalid_argument,
5122  "The input matrix's communicator (Teuchos::Comm object) is null.");
5123  const int myRank = comm->getRank ();
5124 
5125  // Optionally, make a stream for debugging output.
5126  RCP<FancyOStream> err =
5127  debug ? getFancyOStream (rcpFromRef (std::cerr)) : null;
5128  if (debug) {
5129  std::ostringstream os;
5130  os << myRank << ": writeSparse" << endl;
5131  *err << os.str ();
5132  comm->barrier ();
5133  os << "-- " << myRank << ": past barrier" << endl;
5134  *err << os.str ();
5135  }
5136 
5137  // Whether to print debugging output to stderr.
5138  const bool debugPrint = debug && myRank == 0;
5139 
5140  RCP<const map_type> rowMap = pMatrix->getRowMap ();
5141  RCP<const map_type> colMap = pMatrix->getColMap ();
5142  RCP<const map_type> domainMap = pMatrix->getDomainMap ();
5143  RCP<const map_type> rangeMap = pMatrix->getRangeMap ();
5144 
5145  const global_size_t numRows = rangeMap->getGlobalNumElements ();
5146  const global_size_t numCols = domainMap->getGlobalNumElements ();
5147 
5148  if (debug && myRank == 0) {
5149  std::ostringstream os;
5150  os << "-- Input sparse matrix is:"
5151  << "---- " << numRows << " x " << numCols << " with "
5152  << pMatrix->getGlobalNumEntries() << " entries;" << endl
5153  << "---- "
5154  << (pMatrix->isGloballyIndexed() ? "Globally" : "Locally")
5155  << " indexed." << endl
5156  << "---- Its row map has " << rowMap->getGlobalNumElements ()
5157  << " elements." << endl
5158  << "---- Its col map has " << colMap->getGlobalNumElements ()
5159  << " elements." << endl;
5160  *err << os.str ();
5161  }
5162  // Make the "gather" row map, where Proc 0 owns all rows and
5163  // the other procs own no rows.
5164  const size_t localNumRows = (myRank == 0) ? numRows : 0;
5165  RCP<node_type> node = rowMap->getNode();
5166  if (debug) {
5167  std::ostringstream os;
5168  os << "-- " << myRank << ": making gatherRowMap" << endl;
5169  *err << os.str ();
5170  }
5171  RCP<const map_type> gatherRowMap =
5172  Details::computeGatherMap (rowMap, err, debug);
5173 
5174  // Since the matrix may in general be non-square, we need to
5175  // make a column map as well. In this case, the column map
5176  // contains all the columns of the original matrix, because we
5177  // are gathering the whole matrix onto Proc 0. We call
5178  // computeGatherMap to preserve the original order of column
5179  // indices over all the processes.
5180  const size_t localNumCols = (myRank == 0) ? numCols : 0;
5181  RCP<const map_type> gatherColMap =
5182  Details::computeGatherMap (colMap, err, debug);
5183 
5184  // Current map is the source map, gather map is the target map.
5185  typedef Import<LO, GO, node_type> import_type;
5186  import_type importer (rowMap, gatherRowMap);
5187 
5188  // Create a new CrsMatrix to hold the result of the import.
5189  // The constructor needs a column map as well as a row map,
5190  // for the case that the matrix is not square.
5191  RCP<sparse_matrix_type> newMatrix =
5192  rcp (new sparse_matrix_type (gatherRowMap, gatherColMap,
5193  static_cast<size_t> (0)));
5194  // Import the sparse matrix onto Proc 0.
5195  newMatrix->doImport (*pMatrix, importer, INSERT);
5196 
5197  // fillComplete() needs the domain and range maps for the case
5198  // that the matrix is not square.
5199  {
5200  RCP<const map_type> gatherDomainMap =
5201  rcp (new map_type (numCols, localNumCols,
5202  domainMap->getIndexBase (),
5203  comm, node));
5204  RCP<const map_type> gatherRangeMap =
5205  rcp (new map_type (numRows, localNumRows,
5206  rangeMap->getIndexBase (),
5207  comm, node));
5208  newMatrix->fillComplete (gatherDomainMap, gatherRangeMap);
5209  }
5210 
5211  if (debugPrint) {
5212  cerr << "-- Output sparse matrix is:"
5213  << "---- " << newMatrix->getRangeMap ()->getGlobalNumElements ()
5214  << " x "
5215  << newMatrix->getDomainMap ()->getGlobalNumElements ()
5216  << " with "
5217  << newMatrix->getGlobalNumEntries () << " entries;" << endl
5218  << "---- "
5219  << (newMatrix->isGloballyIndexed () ? "Globally" : "Locally")
5220  << " indexed." << endl
5221  << "---- Its row map has "
5222  << newMatrix->getRowMap ()->getGlobalNumElements ()
5223  << " elements, with index base "
5224  << newMatrix->getRowMap ()->getIndexBase () << "." << endl
5225  << "---- Its col map has "
5226  << newMatrix->getColMap ()->getGlobalNumElements ()
5227  << " elements, with index base "
5228  << newMatrix->getColMap ()->getIndexBase () << "." << endl
5229  << "---- Element count of output matrix's column Map may differ "
5230  << "from that of the input matrix's column Map, if some columns "
5231  << "of the matrix contain no entries." << endl;
5232  }
5233 
5234  //
5235  // Print the metadata and the matrix entries on Rank 0.
5236  //
5237  if (myRank == 0) {
5238  // Print the Matrix Market banner line. CrsMatrix stores
5239  // data nonsymmetrically ("general"). This implies that
5240  // readSparse() on a symmetrically stored input file,
5241  // followed by writeSparse() on the resulting sparse matrix,
5242  // will result in an output file with a different banner
5243  // line than the original input file.
5244  out << "%%MatrixMarket matrix coordinate "
5245  << (STS::isComplex ? "complex" : "real")
5246  << " general" << endl;
5247 
5248  // Print comments (the matrix name and / or description).
5249  if (matrixName != "") {
5250  printAsComment (out, matrixName);
5251  }
5252  if (matrixDescription != "") {
5253  printAsComment (out, matrixDescription);
5254  }
5255 
5256  // Print the Matrix Market header (# rows, # columns, #
5257  // nonzeros). Use the range resp. domain map for the number
5258  // of rows resp. columns, since Tpetra::CrsMatrix uses the
5259  // column map for the number of columns. That only
5260  // corresponds to the "linear-algebraic" number of columns
5261  // when the column map is uniquely owned (a.k.a. one-to-one),
5262  // which only happens if the matrix is (block) diagonal.
5263  out << newMatrix->getRangeMap ()->getGlobalNumElements () << " "
5264  << newMatrix->getDomainMap ()->getGlobalNumElements () << " "
5265  << newMatrix->getGlobalNumEntries () << endl;
5266 
5267  // The Matrix Market format expects one-based row and column
5268  // indices. We'll convert the indices on output from
5269  // whatever index base they use to one-based indices.
5270  const GO rowIndexBase = gatherRowMap->getIndexBase ();
5271  const GO colIndexBase = newMatrix->getColMap()->getIndexBase ();
5272  //
5273  // Print the entries of the matrix.
5274  //
5275  // newMatrix can never be globally indexed, since we called
5276  // fillComplete() on it. We include code for both cases
5277  // (globally or locally indexed) just in case that ever
5278  // changes.
5279  if (newMatrix->isGloballyIndexed()) {
5280  // We know that the "gather" row Map is contiguous, so we
5281  // don't need to get the list of GIDs.
5282  const GO minAllGlobalIndex = gatherRowMap->getMinAllGlobalIndex ();
5283  const GO maxAllGlobalIndex = gatherRowMap->getMaxAllGlobalIndex ();
5284  for (GO globalRowIndex = minAllGlobalIndex;
5285  globalRowIndex <= maxAllGlobalIndex; // inclusive range
5286  ++globalRowIndex) {
5287  ArrayView<const GO> ind;
5288  ArrayView<const ST> val;
5289  newMatrix->getGlobalRowView (globalRowIndex, ind, val);
5290  go_iter indIter = ind.begin ();
5291  st_iter valIter = val.begin ();
5292  for (; indIter != ind.end() && valIter != val.end();
5293  ++indIter, ++valIter) {
5294  const GO globalColIndex = *indIter;
5295  // Convert row and column indices to 1-based.
5296  // This works because the global index type is signed.
5297  out << (globalRowIndex + 1 - rowIndexBase) << " "
5298  << (globalColIndex + 1 - colIndexBase) << " ";
5299  if (STS::isComplex) {
5300  out << STS::real (*valIter) << " " << STS::imag (*valIter);
5301  } else {
5302  out << *valIter;
5303  }
5304  out << endl;
5305  } // For each entry in the current row
5306  } // For each row of the "gather" matrix
5307  } else { // newMatrix is locally indexed
5308  typedef OrdinalTraits<GO> OTG;
5309  for (LO localRowIndex = gatherRowMap->getMinLocalIndex();
5310  localRowIndex <= gatherRowMap->getMaxLocalIndex();
5311  ++localRowIndex) {
5312  // Convert from local to global row index.
5313  const GO globalRowIndex =
5314  gatherRowMap->getGlobalElement (localRowIndex);
5315  TEUCHOS_TEST_FOR_EXCEPTION(
5316  globalRowIndex == OTG::invalid(), std::logic_error,
5317  "Failed to convert the supposed local row index "
5318  << localRowIndex << " into a global row index. "
5319  "Please report this bug to the Tpetra developers.");
5320  ArrayView<const LO> ind;
5321  ArrayView<const ST> val;
5322  newMatrix->getLocalRowView (localRowIndex, ind, val);
5323  lo_iter indIter = ind.begin ();
5324  st_iter valIter = val.begin ();
5325  for (; indIter != ind.end() && valIter != val.end();
5326  ++indIter, ++valIter) {
5327  // Convert the column index from local to global.
5328  const GO globalColIndex =
5329  newMatrix->getColMap()->getGlobalElement (*indIter);
5330  TEUCHOS_TEST_FOR_EXCEPTION(
5331  globalColIndex == OTG::invalid(), std::logic_error,
5332  "On local row " << localRowIndex << " of the sparse matrix: "
5333  "Failed to convert the supposed local column index "
5334  << *indIter << " into a global column index. Please report "
5335  "this bug to the Tpetra developers.");
5336  // Convert row and column indices to 1-based.
5337  // This works because the global index type is signed.
5338  out << (globalRowIndex + 1 - rowIndexBase) << " "
5339  << (globalColIndex + 1 - colIndexBase) << " ";
5340  if (STS::isComplex) {
5341  out << STS::real (*valIter) << " " << STS::imag (*valIter);
5342  } else {
5343  out << *valIter;
5344  }
5345  out << endl;
5346  } // For each entry in the current row
5347  } // For each row of the "gather" matrix
5348  } // Whether the "gather" matrix is locally or globally indexed
5349  } // If my process' rank is 0
5350  }
5351 
5374  static void
5375  writeSparse (std::ostream& out,
5376  const RCP<const sparse_matrix_type>& pMatrix,
5377  const bool debug=false)
5378  {
5379  writeSparse (out, pMatrix, "", "", debug);
5380  }
5381 
5410  static void
5411  writeDenseFile (const std::string& filename,
5412  const multivector_type& X,
5413  const std::string& matrixName,
5414  const std::string& matrixDescription,
5415  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5416  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5417  {
5418  const int myRank = X.getMap ().is_null () ? 0 :
5419  (X.getMap ()->getComm ().is_null () ? 0 :
5420  X.getMap ()->getComm ()->getRank ());
5421  std::ofstream out;
5422 
5423  if (myRank == 0) { // Only open the file on Process 0.
5424  out.open (filename.c_str());
5425  }
5426 
5427  writeDense (out, X, matrixName, matrixDescription, err, dbg);
5428  // We can rely on the destructor of the output stream to close
5429  // the file on scope exit, even if writeDense() throws an
5430  // exception.
5431  }
5432 
5438  static void
5439  writeDenseFile (const std::string& filename,
5440  const Teuchos::RCP<const multivector_type>& X,
5441  const std::string& matrixName,
5442  const std::string& matrixDescription,
5443  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5444  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5445  {
5446  TEUCHOS_TEST_FOR_EXCEPTION(
5447  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
5448  "writeDenseFile: The input MultiVector X is null.");
5449  writeDenseFile (filename, *X, matrixName, matrixDescription, err, dbg);
5450  }
5451 
5457  static void
5458  writeDenseFile (const std::string& filename,
5459  const multivector_type& X,
5460  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5461  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5462  {
5463  writeDenseFile (filename, X, "", "", err, dbg);
5464  }
5465 
5471  static void
5472  writeDenseFile (const std::string& filename,
5473  const RCP<const multivector_type>& X,
5474  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5475  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5476  {
5477  TEUCHOS_TEST_FOR_EXCEPTION(
5478  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
5479  "writeDenseFile: The input MultiVector X is null.");
5480  writeDenseFile (filename, *X, err, dbg);
5481  }
5482 
5483 
5514  static void
5515  writeDense (std::ostream& out,
5516  const multivector_type& X,
5517  const std::string& matrixName,
5518  const std::string& matrixDescription,
5519  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5520  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5521  {
5522  using Teuchos::Comm;
5523  using Teuchos::outArg;
5524  using Teuchos::REDUCE_MAX;
5525  using Teuchos::reduceAll;
5526  using Teuchos::RCP;
5527  using std::endl;
5528 
5529  RCP<const Comm<int> > comm = X.getMap ().is_null () ?
5530  Teuchos::null : X.getMap ()->getComm ();
5531  const int myRank = comm.is_null () ? 0 : comm->getRank ();
5532 
5533  // If the caller provides a nonnull debug output stream, we
5534  // print debugging output to it. This is a local thing; we
5535  // don't have to check across processes.
5536  const bool debug = ! dbg.is_null ();
5537  if (debug) {
5538  dbg->pushTab ();
5539  std::ostringstream os;
5540  os << myRank << ": writeDense" << endl;
5541  *dbg << os.str ();
5542  dbg->pushTab ();
5543  }
5544  // Print the Matrix Market header.
5545  writeDenseHeader (out, X, matrixName, matrixDescription, err, dbg);
5546 
5547  // Print each column one at a time. This is a (perhaps)
5548  // temporary fix for Bug 6288.
5549  const size_t numVecs = X.getNumVectors ();
5550  for (size_t j = 0; j < numVecs; ++j) {
5551  writeDenseColumn (out, * (X.getVector (j)), err, dbg);
5552  }
5553 
5554  if (debug) {
5555  dbg->popTab ();
5556  std::ostringstream os;
5557  os << myRank << ": writeDense: Done" << endl;
5558  *dbg << os.str ();
5559  dbg->popTab ();
5560  }
5561  }
5562 
5563  private:
5564 
5590  static void
5591  writeDenseHeader (std::ostream& out,
5592  const multivector_type& X,
5593  const std::string& matrixName,
5594  const std::string& matrixDescription,
5595  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5596  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5597  {
5598  using Teuchos::Comm;
5599  using Teuchos::outArg;
5600  using Teuchos::RCP;
5601  using Teuchos::REDUCE_MAX;
5602  using Teuchos::reduceAll;
5603  using std::endl;
5604  typedef typename multivector_type::scalar_type scalar_type;
5605  typedef Teuchos::ScalarTraits<scalar_type> STS;
5606  const char prefix[] = "Tpetra::MatrixMarket::writeDenseHeader: ";
5607 
5608  RCP<const Comm<int> > comm = X.getMap ().is_null () ?
5609  Teuchos::null : X.getMap ()->getComm ();
5610  const int myRank = comm.is_null () ? 0 : comm->getRank ();
5611  int lclErr = 0; // whether this MPI process has seen an error
5612  int gblErr = 0; // whether we know if some MPI process has seen an error
5613 
5614  // If the caller provides a nonnull debug output stream, we
5615  // print debugging output to it. This is a local thing; we
5616  // don't have to check across processes.
5617  const bool debug = ! dbg.is_null ();
5618 
5619  if (debug) {
5620  dbg->pushTab ();
5621  std::ostringstream os;
5622  os << myRank << ": writeDenseHeader" << endl;
5623  *dbg << os.str ();
5624  dbg->pushTab ();
5625  }
5626 
5627  //
5628  // Process 0: Write the MatrixMarket header.
5629  //
5630  if (myRank == 0) {
5631  try {
5632  // Print the Matrix Market header. MultiVector stores data
5633  // nonsymmetrically, hence "general" in the banner line.
5634  // Print first to a temporary string output stream, and then
5635  // write it to the main output stream, so that at least the
5636  // header output has transactional semantics. We can't
5637  // guarantee transactional semantics for the whole output,
5638  // since that would not be memory scalable. (This could be
5639  // done in the file system by using a temporary file; we
5640  // don't do this, but users could.)
5641  std::ostringstream hdr;
5642  {
5643  std::string dataType;
5644  if (STS::isComplex) {
5645  dataType = "complex";
5646  } else if (STS::isOrdinal) {
5647  dataType = "integer";
5648  } else {
5649  dataType = "real";
5650  }
5651  hdr << "%%MatrixMarket matrix array " << dataType << " general"
5652  << endl;
5653  }
5654 
5655  // Print comments (the matrix name and / or description).
5656  if (matrixName != "") {
5657  printAsComment (hdr, matrixName);
5658  }
5659  if (matrixDescription != "") {
5660  printAsComment (hdr, matrixDescription);
5661  }
5662  // Print the Matrix Market dimensions header for dense matrices.
5663  hdr << X.getGlobalLength () << " " << X.getNumVectors () << endl;
5664 
5665  // Write the MatrixMarket header to the output stream.
5666  out << hdr.str ();
5667  } catch (std::exception& e) {
5668  if (! err.is_null ()) {
5669  *err << prefix << "While writing the Matrix Market header, "
5670  "Process 0 threw an exception: " << e.what () << endl;
5671  }
5672  lclErr = 1;
5673  }
5674  } // if I am Process 0
5675 
5676  // Establish global agreement on the error state. It wouldn't
5677  // be good for other processes to keep going, if Process 0
5678  // finds out that it can't write to the given output stream.
5679  reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
5680  TEUCHOS_TEST_FOR_EXCEPTION(
5681  gblErr == 1, std::runtime_error, prefix << "Some error occurred "
5682  "which prevented this method from completing.");
5683 
5684  if (debug) {
5685  dbg->popTab ();
5686  *dbg << myRank << ": writeDenseHeader: Done" << endl;
5687  dbg->popTab ();
5688  }
5689  }
5690 
5708  static void
5709  writeDenseColumn (std::ostream& out,
5710  const multivector_type& X,
5711  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
5712  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
5713  {
5714  using Teuchos::arcp;
5715  using Teuchos::ArrayRCP;
5716  using Teuchos::ArrayView;
5717  using Teuchos::CommRequest;
5718  using Teuchos::ireceive;
5719  using Teuchos::isend;
5720  using Teuchos::outArg;
5721  using Teuchos::REDUCE_MAX;
5722  using Teuchos::reduceAll;
5723  using Teuchos::RCP;
5724  using Teuchos::TypeNameTraits;
5725  using Teuchos::wait;
5726  using std::endl;
5727  typedef typename multivector_type::scalar_type scalar_type;
5728  typedef Teuchos::ScalarTraits<scalar_type> STS;
5729 
5730  const Comm<int>& comm = * (X.getMap ()->getComm ());
5731  const int myRank = comm.getRank ();
5732  const int numProcs = comm.getSize ();
5733  int lclErr = 0; // whether this MPI process has seen an error
5734  int gblErr = 0; // whether we know if some MPI process has seen an error
5735 
5736  // If the caller provides a nonnull debug output stream, we
5737  // print debugging output to it. This is a local thing; we
5738  // don't have to check across processes.
5739  const bool debug = ! dbg.is_null ();
5740 
5741  if (debug) {
5742  dbg->pushTab ();
5743  std::ostringstream os;
5744  os << myRank << ": writeDenseColumn" << endl;
5745  *dbg << os.str ();
5746  dbg->pushTab ();
5747  }
5748 
5749  // Make the output stream write floating-point numbers in
5750  // scientific notation. It will politely put the output
5751  // stream back to its state on input, when this scope
5752  // terminates.
5753  Teuchos::MatrixMarket::details::SetScientific<scalar_type> sci (out);
5754 
5755  const size_t myNumRows = X.getLocalLength ();
5756  const size_t numCols = X.getNumVectors ();
5757  // Use a different tag for the "size" messages than for the
5758  // "data" messages, in order to help us debug any mix-ups.
5759  const int sizeTag = 1337;
5760  const int dataTag = 1338;
5761 
5762  // Process 0 pipelines nonblocking receives with file output.
5763  //
5764  // Constraints:
5765  // - Process 0 can't post a receive for another process'
5766  // actual data, until it posts and waits on the receive
5767  // from that process with the amount of data to receive.
5768  // (We could just post receives with a max data size, but
5769  // I feel uncomfortable about that.)
5770  // - The C++ standard library doesn't allow nonblocking
5771  // output to an std::ostream. (Thus, we have to start a
5772  // receive or send before starting the write, and hope
5773  // that MPI completes it in the background.)
5774  //
5775  // Process 0: Post receive-size receives from Processes 1 and 2.
5776  // Process 1: Post send-size send to Process 0.
5777  // Process 2: Post send-size send to Process 0.
5778  //
5779  // All processes: Pack my entries.
5780  //
5781  // Process 1:
5782  // - Post send-data send to Process 0.
5783  // - Wait on my send-size send to Process 0.
5784  //
5785  // Process 0:
5786  // - Print MatrixMarket header.
5787  // - Print my entries.
5788  // - Wait on receive-size receive from Process 1.
5789  // - Post receive-data receive from Process 1.
5790  //
5791  // For each process p = 1, 2, ... numProcs-1:
5792  // If I am Process 0:
5793  // - Post receive-size receive from Process p + 2
5794  // - Wait on receive-size receive from Process p + 1
5795  // - Post receive-data receive from Process p + 1
5796  // - Wait on receive-data receive from Process p
5797  // - Write data from Process p.
5798  // Else if I am Process p:
5799  // - Wait on my send-data send.
5800  // Else if I am Process p+1:
5801  // - Post send-data send to Process 0.
5802  // - Wait on my send-size send.
5803  // Else if I am Process p+2:
5804  // - Post send-size send to Process 0.
5805  //
5806  // Pipelining has three goals here:
5807  // 1. Overlap communication (the receives) with file I/O
5808  // 2. Give Process 0 a chance to prepost some receives,
5809  // before sends show up, by packing local data before
5810  // posting sends
5811  // 3. Don't post _all_ receives or _all_ sends, because that
5812  // wouldn't be memory scalable. (Just because we can't
5813  // see how much memory MPI consumes, doesn't mean that it
5814  // doesn't consume any!)
5815 
5816  // These are used on every process. sendReqSize[0] holds the
5817  // number of rows on this process, and sendReqBuf holds this
5818  // process' data. Process 0 packs into sendReqBuf, but
5819  // doesn't send; it only uses that for printing. All other
5820  // processes send both of these to Process 0.
5821  RCP<CommRequest<int> > sendReqSize, sendReqData;
5822 
5823  // These are used only on Process 0, for received data. Keep
5824  // 3 of each, and treat the arrays as circular buffers. When
5825  // receiving from Process p, the corresponding array index
5826  // here is p % 3.
5827  Array<ArrayRCP<size_t> > recvSizeBufs (3);
5828  Array<ArrayRCP<scalar_type> > recvDataBufs (3);
5829  Array<RCP<CommRequest<int> > > recvSizeReqs (3);
5830  Array<RCP<CommRequest<int> > > recvDataReqs (3);
5831 
5832  // Buffer for nonblocking send of the "send size."
5833  ArrayRCP<size_t> sendDataSize (1);
5834  sendDataSize[0] = myNumRows;
5835 
5836  if (myRank == 0) {
5837  if (debug) {
5838  std::ostringstream os;
5839  os << myRank << ": Post receive-size receives from "
5840  "Procs 1 and 2: tag = " << sizeTag << endl;
5841  *dbg << os.str ();
5842  }
5843  // Process 0: Post receive-size receives from Processes 1 and 2.
5844  recvSizeBufs[0].resize (1);
5845  // Set these three to an invalid value as a flag. If we
5846  // don't get these messages, then the invalid value will
5847  // remain, so we can test for it.
5848  (recvSizeBufs[0])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
5849  recvSizeBufs[1].resize (1);
5850  (recvSizeBufs[1])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
5851  recvSizeBufs[2].resize (1);
5852  (recvSizeBufs[2])[0] = Teuchos::OrdinalTraits<size_t>::invalid ();
5853  if (numProcs > 1) {
5854  recvSizeReqs[1] =
5855  ireceive<int, size_t> (recvSizeBufs[1], 1, sizeTag, comm);
5856  }
5857  if (numProcs > 2) {
5858  recvSizeReqs[2] =
5859  ireceive<int, size_t> (recvSizeBufs[2], 2, sizeTag, comm);
5860  }
5861  }
5862  else if (myRank == 1 || myRank == 2) {
5863  if (debug) {
5864  std::ostringstream os;
5865  os << myRank << ": Post send-size send: size = "
5866  << sendDataSize[0] << ", tag = " << sizeTag << endl;
5867  *dbg << os.str ();
5868  }
5869  // Prime the pipeline by having Processes 1 and 2 start
5870  // their send-size sends. We don't want _all_ the processes
5871  // to start their send-size sends, because that wouldn't be
5872  // memory scalable.
5873  sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
5874  }
5875  else {
5876  if (debug) {
5877  std::ostringstream os;
5878  os << myRank << ": Not posting my send-size send yet" << endl;
5879  *dbg << os.str ();
5880  }
5881  }
5882 
5883  //
5884  // Pack my entries, in column-major order.
5885  //
5886  if (debug) {
5887  std::ostringstream os;
5888  os << myRank << ": Pack my entries" << endl;
5889  *dbg << os.str ();
5890  }
5891  ArrayRCP<scalar_type> sendDataBuf;
5892  try {
5893  sendDataBuf = arcp<scalar_type> (myNumRows * numCols);
5894  X.get1dCopy (sendDataBuf (), myNumRows);
5895  }
5896  catch (std::exception& e) {
5897  lclErr = 1;
5898  if (! err.is_null ()) {
5899  std::ostringstream os;
5900  os << "Process " << myRank << ": Attempt to pack my MultiVector "
5901  "entries threw an exception: " << e.what () << endl;
5902  *err << os.str ();
5903  }
5904  }
5905  if (debug) {
5906  std::ostringstream os;
5907  os << myRank << ": Done packing my entries" << endl;
5908  *dbg << os.str ();
5909  }
5910 
5911  //
5912  // Process 1: post send-data send to Process 0.
5913  //
5914  if (myRank == 1) {
5915  if (debug) {
5916  *dbg << myRank << ": Post send-data send: tag = " << dataTag
5917  << endl;
5918  }
5919  sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
5920  }
5921 
5922  //
5923  // Process 0: Write my entries.
5924  //
5925  if (myRank == 0) {
5926  if (debug) {
5927  std::ostringstream os;
5928  os << myRank << ": Write my entries" << endl;
5929  *dbg << os.str ();
5930  }
5931 
5932  // Write Process 0's data to the output stream.
5933  // Matrix Market prints dense matrices in column-major order.
5934  const size_t printNumRows = myNumRows;
5935  ArrayView<const scalar_type> printData = sendDataBuf ();
5936  const size_t printStride = printNumRows;
5937  if (static_cast<size_t> (printData.size ()) < printStride * numCols) {
5938  lclErr = 1;
5939  if (! err.is_null ()) {
5940  std::ostringstream os;
5941  os << "Process " << myRank << ": My MultiVector data's size "
5942  << printData.size () << " does not match my local dimensions "
5943  << printStride << " x " << numCols << "." << endl;
5944  *err << os.str ();
5945  }
5946  }
5947  else {
5948  // Matrix Market dense format wants one number per line.
5949  // It wants each complex number as two real numbers (real
5950  // resp. imaginary parts) with a space between.
5951  for (size_t col = 0; col < numCols; ++col) {
5952  for (size_t row = 0; row < printNumRows; ++row) {
5953  if (STS::isComplex) {
5954  out << STS::real (printData[row + col * printStride]) << " "
5955  << STS::imag (printData[row + col * printStride]) << endl;
5956  } else {
5957  out << printData[row + col * printStride] << endl;
5958  }
5959  }
5960  }
5961  }
5962  }
5963 
5964  if (myRank == 0) {
5965  // Wait on receive-size receive from Process 1.
5966  const int recvRank = 1;
5967  const int circBufInd = recvRank % 3;
5968  if (debug) {
5969  std::ostringstream os;
5970  os << myRank << ": Wait on receive-size receive from Process "
5971  << recvRank << endl;
5972  *dbg << os.str ();
5973  }
5974  if (numProcs > 1) {
5975  wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
5976 
5977  // We received the number of rows of data. (The data
5978  // come in two columns.)
5979  size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
5980  if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
5981  lclErr = 1;
5982  if (! err.is_null ()) {
5983  std::ostringstream os;
5984  os << myRank << ": Result of receive-size receive from Process "
5985  << recvRank << " is Teuchos::OrdinalTraits<size_t>::invalid() "
5986  << "= " << Teuchos::OrdinalTraits<size_t>::invalid () << ". "
5987  "This should never happen, and suggests that the receive never "
5988  "got posted. Please report this bug to the Tpetra developers."
5989  << endl;
5990  *err << os.str ();
5991  }
5992 
5993  // If we're going to continue after error, set the
5994  // number of rows to receive to a reasonable size. This
5995  // may cause MPI_ERR_TRUNCATE if the sending process is
5996  // sending more than 0 rows, but that's better than MPI
5997  // overflowing due to the huge positive value that is
5998  // Teuchos::OrdinalTraits<size_t>::invalid().
5999  recvNumRows = 0;
6000  }
6001 
6002  // Post receive-data receive from Process 1.
6003  recvDataBufs[circBufInd].resize (recvNumRows * numCols);
6004  if (debug) {
6005  std::ostringstream os;
6006  os << myRank << ": Post receive-data receive from Process "
6007  << recvRank << ": tag = " << dataTag << ", buffer size = "
6008  << recvDataBufs[circBufInd].size () << endl;
6009  *dbg << os.str ();
6010  }
6011  if (! recvSizeReqs[circBufInd].is_null ()) {
6012  lclErr = 1;
6013  if (! err.is_null ()) {
6014  std::ostringstream os;
6015  os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
6016  "null, before posting the receive-data receive from Process "
6017  << recvRank << ". This should never happen. Please report "
6018  "this bug to the Tpetra developers." << endl;
6019  *err << os.str ();
6020  }
6021  }
6022  recvDataReqs[circBufInd] =
6023  ireceive<int, scalar_type> (recvDataBufs[circBufInd],
6024  recvRank, dataTag, comm);
6025  } // numProcs > 1
6026  }
6027  else if (myRank == 1) {
6028  // Wait on my send-size send.
6029  if (debug) {
6030  std::ostringstream os;
6031  os << myRank << ": Wait on my send-size send" << endl;
6032  *dbg << os.str ();
6033  }
6034  wait<int> (comm, outArg (sendReqSize));
6035  }
6036 
6037  //
6038  // Pipeline loop
6039  //
6040  for (int p = 1; p < numProcs; ++p) {
6041  if (myRank == 0) {
6042  if (p + 2 < numProcs) {
6043  // Post receive-size receive from Process p + 2.
6044  const int recvRank = p + 2;
6045  const int circBufInd = recvRank % 3;
6046  if (debug) {
6047  std::ostringstream os;
6048  os << myRank << ": Post receive-size receive from Process "
6049  << recvRank << ": tag = " << sizeTag << endl;
6050  *dbg << os.str ();
6051  }
6052  if (! recvSizeReqs[circBufInd].is_null ()) {
6053  lclErr = 1;
6054  if (! err.is_null ()) {
6055  std::ostringstream os;
6056  os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
6057  << "null, for the receive-size receive from Process "
6058  << recvRank << "! This may mean that this process never "
6059  << "finished waiting for the receive from Process "
6060  << (recvRank - 3) << "." << endl;
6061  *err << os.str ();
6062  }
6063  }
6064  recvSizeReqs[circBufInd] =
6065  ireceive<int, size_t> (recvSizeBufs[circBufInd],
6066  recvRank, sizeTag, comm);
6067  }
6068 
6069  if (p + 1 < numProcs) {
6070  const int recvRank = p + 1;
6071  const int circBufInd = recvRank % 3;
6072 
6073  // Wait on receive-size receive from Process p + 1.
6074  if (debug) {
6075  std::ostringstream os;
6076  os << myRank << ": Wait on receive-size receive from Process "
6077  << recvRank << endl;
6078  *dbg << os.str ();
6079  }
6080  wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
6081 
6082  // We received the number of rows of data. (The data
6083  // come in two columns.)
6084  size_t recvNumRows = (recvSizeBufs[circBufInd])[0];
6085  if (recvNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
6086  lclErr = 1;
6087  if (! err.is_null ()) {
6088  std::ostringstream os;
6089  os << myRank << ": Result of receive-size receive from Process "
6090  << recvRank << " is Teuchos::OrdinalTraits<size_t>::invalid() "
6091  << "= " << Teuchos::OrdinalTraits<size_t>::invalid () << ". "
6092  "This should never happen, and suggests that the receive never "
6093  "got posted. Please report this bug to the Tpetra developers."
6094  << endl;
6095  *err << os.str ();
6096  }
6097  // If we're going to continue after error, set the
6098  // number of rows to receive to a reasonable size.
6099  // This may cause MPI_ERR_TRUNCATE if the sending
6100  // process sends more than 0 rows, but that's better
6101  // than MPI overflowing due to the huge positive value
6102  // Teuchos::OrdinalTraits<size_t>::invalid().
6103  recvNumRows = 0;
6104  }
6105 
6106  // Post receive-data receive from Process p + 1.
6107  recvDataBufs[circBufInd].resize (recvNumRows * numCols);
6108  if (debug) {
6109  std::ostringstream os;
6110  os << myRank << ": Post receive-data receive from Process "
6111  << recvRank << ": tag = " << dataTag << ", buffer size = "
6112  << recvDataBufs[circBufInd].size () << endl;
6113  *dbg << os.str ();
6114  }
6115  if (! recvDataReqs[circBufInd].is_null ()) {
6116  lclErr = 1;
6117  if (! err.is_null ()) {
6118  std::ostringstream os;
6119  os << myRank << ": recvDataReqs[" << circBufInd << "] is not "
6120  << "null, for the receive-data receive from Process "
6121  << recvRank << "! This may mean that this process never "
6122  << "finished waiting for the receive from Process "
6123  << (recvRank - 3) << "." << endl;
6124  *err << os.str ();
6125  }
6126  }
6127  recvDataReqs[circBufInd] =
6128  ireceive<int, scalar_type> (recvDataBufs[circBufInd],
6129  recvRank, dataTag, comm);
6130  }
6131 
6132  // Wait on receive-data receive from Process p.
6133  const int recvRank = p;
6134  const int circBufInd = recvRank % 3;
6135  if (debug) {
6136  std::ostringstream os;
6137  os << myRank << ": Wait on receive-data receive from Process "
6138  << recvRank << endl;
6139  *dbg << os.str ();
6140  }
6141  wait<int> (comm, outArg (recvDataReqs[circBufInd]));
6142 
6143  // Write Process p's data. Number of rows lives in
6144  // recvSizeBufs[circBufInd], and the actual data live in
6145  // recvDataBufs[circBufInd]. Do this after posting receives,
6146  // in order to expose overlap of comm. with file I/O.
6147  if (debug) {
6148  std::ostringstream os;
6149  os << myRank << ": Write entries from Process " << recvRank
6150  << endl;
6151  *dbg << os.str () << endl;
6152  }
6153  size_t printNumRows = (recvSizeBufs[circBufInd])[0];
6154  if (printNumRows == Teuchos::OrdinalTraits<size_t>::invalid ()) {
6155  lclErr = 1;
6156  if (! err.is_null ()) {
6157  std::ostringstream os;
6158  os << myRank << ": Result of receive-size receive from Process "
6159  << recvRank << " was Teuchos::OrdinalTraits<size_t>::"
6160  "invalid() = " << Teuchos::OrdinalTraits<size_t>::invalid ()
6161  << ". This should never happen, and suggests that its "
6162  "receive-size receive was never posted. "
6163  "Please report this bug to the Tpetra developers." << endl;
6164  *err << os.str ();
6165  }
6166  // If we're going to continue after error, set the
6167  // number of rows to print to a reasonable size.
6168  printNumRows = 0;
6169  }
6170  if (printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
6171  lclErr = 1;
6172  if (! err.is_null ()) {
6173  std::ostringstream os;
6174  os << myRank << ": Result of receive-size receive from Proc "
6175  << recvRank << " was " << printNumRows << " > 0, but "
6176  "recvDataBufs[" << circBufInd << "] is null. This should "
6177  "never happen. Please report this bug to the Tpetra "
6178  "developers." << endl;
6179  *err << os.str ();
6180  }
6181  // If we're going to continue after error, set the
6182  // number of rows to print to a reasonable size.
6183  printNumRows = 0;
6184  }
6185 
6186  // Write the received data to the output stream.
6187  // Matrix Market prints dense matrices in column-major order.
6188  ArrayView<const scalar_type> printData = (recvDataBufs[circBufInd]) ();
6189  const size_t printStride = printNumRows;
6190  // Matrix Market dense format wants one number per line.
6191  // It wants each complex number as two real numbers (real
6192  // resp. imaginary parts) with a space between.
6193  for (size_t col = 0; col < numCols; ++col) {
6194  for (size_t row = 0; row < printNumRows; ++row) {
6195  if (STS::isComplex) {
6196  out << STS::real (printData[row + col * printStride]) << " "
6197  << STS::imag (printData[row + col * printStride]) << endl;
6198  } else {
6199  out << printData[row + col * printStride] << endl;
6200  }
6201  }
6202  }
6203  }
6204  else if (myRank == p) { // Process p
6205  // Wait on my send-data send.
6206  if (debug) {
6207  std::ostringstream os;
6208  os << myRank << ": Wait on my send-data send" << endl;
6209  *dbg << os.str ();
6210  }
6211  wait<int> (comm, outArg (sendReqData));
6212  }
6213  else if (myRank == p + 1) { // Process p + 1
6214  // Post send-data send to Process 0.
6215  if (debug) {
6216  std::ostringstream os;
6217  os << myRank << ": Post send-data send: tag = " << dataTag
6218  << endl;
6219  *dbg << os.str ();
6220  }
6221  sendReqData = isend<int, scalar_type> (sendDataBuf, 0, dataTag, comm);
6222  // Wait on my send-size send.
6223  if (debug) {
6224  std::ostringstream os;
6225  os << myRank << ": Wait on my send-size send" << endl;
6226  *dbg << os.str ();
6227  }
6228  wait<int> (comm, outArg (sendReqSize));
6229  }
6230  else if (myRank == p + 2) { // Process p + 2
6231  // Post send-size send to Process 0.
6232  if (debug) {
6233  std::ostringstream os;
6234  os << myRank << ": Post send-size send: size = "
6235  << sendDataSize[0] << ", tag = " << sizeTag << endl;
6236  *dbg << os.str ();
6237  }
6238  sendReqSize = isend<int, size_t> (sendDataSize, 0, sizeTag, comm);
6239  }
6240  }
6241 
6242  // Establish global agreement on the error state.
6243  reduceAll<int, int> (comm, REDUCE_MAX, lclErr, outArg (gblErr));
6244  TEUCHOS_TEST_FOR_EXCEPTION(
6245  gblErr == 1, std::runtime_error, "Tpetra::MatrixMarket::writeDense "
6246  "experienced some kind of error and was unable to complete.");
6247 
6248  if (debug) {
6249  dbg->popTab ();
6250  *dbg << myRank << ": writeDenseColumn: Done" << endl;
6251  dbg->popTab ();
6252  }
6253  }
6254 
6255  public:
6256 
6262  static void
6263  writeDense (std::ostream& out,
6264  const RCP<const multivector_type>& X,
6265  const std::string& matrixName,
6266  const std::string& matrixDescription,
6267  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6268  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6269  {
6270  TEUCHOS_TEST_FOR_EXCEPTION(
6271  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
6272  "writeDense: The input MultiVector X is null.");
6273  writeDense (out, *X, matrixName, matrixDescription, err, dbg);
6274  }
6275 
6281  static void
6282  writeDense (std::ostream& out,
6283  const multivector_type& X,
6284  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6285  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6286  {
6287  writeDense (out, X, "", "", err, dbg);
6288  }
6289 
6295  static void
6296  writeDense (std::ostream& out,
6297  const RCP<const multivector_type>& X,
6298  const Teuchos::RCP<Teuchos::FancyOStream>& err = Teuchos::null,
6299  const Teuchos::RCP<Teuchos::FancyOStream>& dbg = Teuchos::null)
6300  {
6301  TEUCHOS_TEST_FOR_EXCEPTION(
6302  X.is_null (), std::invalid_argument, "Tpetra::MatrixMarket::"
6303  "writeDense: The input MultiVector X is null.");
6304  writeDense (out, *X, "", "", err, dbg);
6305  }
6306 
6326  static void
6327  writeMap (std::ostream& out, const map_type& map, const bool debug=false)
6328  {
6329  Teuchos::RCP<Teuchos::FancyOStream> err =
6330  Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
6331  writeMap (out, map, err, debug);
6332  }
6333 
6342  static void
6343  writeMap (std::ostream& out,
6344  const map_type& map,
6345  const Teuchos::RCP<Teuchos::FancyOStream>& err,
6346  const bool debug=false)
6347  {
6348  using Teuchos::ArrayRCP;
6349  using Teuchos::ArrayView;
6350  using Teuchos::CommRequest;
6351  using Teuchos::ireceive;
6352  using Teuchos::isend;
6353  using Teuchos::RCP;
6354  using Teuchos::TypeNameTraits;
6355  using Teuchos::wait;
6356  using std::endl;
6357  typedef global_ordinal_type GO;
6358  typedef int pid_type;
6359 
6360  // Treat the Map as a 1-column "multivector." This differs
6361  // from the previous two-column format, in which column 0 held
6362  // the GIDs, and column 1 held the corresponding PIDs. It
6363  // differs because printing that format requires knowing the
6364  // entire first column -- that is, all the GIDs -- in advance.
6365  // Sending messages from each process one at a time saves
6366  // memory, but it means that Process 0 doesn't ever have all
6367  // the GIDs at once.
6368  //
6369  // We pack the entries as ptrdiff_t, since this should be the
6370  // biggest signed built-in integer type that can hold any GO
6371  // or pid_type (= int) quantity without overflow. Test this
6372  // assumption at run time.
6373  typedef ptrdiff_t int_type;
6374  TEUCHOS_TEST_FOR_EXCEPTION(
6375  sizeof (GO) > sizeof (int_type), std::logic_error,
6376  "The global ordinal type GO=" << TypeNameTraits<GO>::name ()
6377  << " is too big for ptrdiff_t. sizeof(GO) = " << sizeof (GO)
6378  << " > sizeof(ptrdiff_t) = " << sizeof (ptrdiff_t) << ".");
6379  TEUCHOS_TEST_FOR_EXCEPTION(
6380  sizeof (pid_type) > sizeof (int_type), std::logic_error,
6381  "The (MPI) process rank type pid_type=" <<
6382  TypeNameTraits<pid_type>::name () << " is too big for ptrdiff_t. "
6383  "sizeof(pid_type) = " << sizeof (pid_type) << " > sizeof(ptrdiff_t)"
6384  " = " << sizeof (ptrdiff_t) << ".");
6385 
6386  const Comm<int>& comm = * (map.getComm ());
6387  const int myRank = comm.getRank ();
6388  const int numProcs = comm.getSize ();
6389 
6390  if (! err.is_null ()) {
6391  err->pushTab ();
6392  }
6393  if (debug) {
6394  std::ostringstream os;
6395  os << myRank << ": writeMap" << endl;
6396  *err << os.str ();
6397  }
6398  if (! err.is_null ()) {
6399  err->pushTab ();
6400  }
6401 
6402  const size_t myNumRows = map.getNodeNumElements ();
6403  // Use a different tag for the "size" messages than for the
6404  // "data" messages, in order to help us debug any mix-ups.
6405  const int sizeTag = 1337;
6406  const int dataTag = 1338;
6407 
6408  // Process 0 pipelines nonblocking receives with file output.
6409  //
6410  // Constraints:
6411  // - Process 0 can't post a receive for another process'
6412  // actual data, until it posts and waits on the receive
6413  // from that process with the amount of data to receive.
6414  // (We could just post receives with a max data size, but
6415  // I feel uncomfortable about that.)
6416  // - The C++ standard library doesn't allow nonblocking
6417  // output to an std::ostream.
6418  //
6419  // Process 0: Post receive-size receives from Processes 1 and 2.
6420  // Process 1: Post send-size send to Process 0.
6421  // Process 2: Post send-size send to Process 0.
6422  //
6423  // All processes: Pack my GIDs and PIDs.
6424  //
6425  // Process 1:
6426  // - Post send-data send to Process 0.
6427  // - Wait on my send-size send to Process 0.
6428  //
6429  // Process 0:
6430  // - Print MatrixMarket header.
6431  // - Print my GIDs and PIDs.
6432  // - Wait on receive-size receive from Process 1.
6433  // - Post receive-data receive from Process 1.
6434  //
6435  // For each process p = 1, 2, ... numProcs-1:
6436  // If I am Process 0:
6437  // - Post receive-size receive from Process p + 2
6438  // - Wait on receive-size receive from Process p + 1
6439  // - Post receive-data receive from Process p + 1
6440  // - Wait on receive-data receive from Process p
6441  // - Write data from Process p.
6442  // Else if I am Process p:
6443  // - Wait on my send-data send.
6444  // Else if I am Process p+1:
6445  // - Post send-data send to Process 0.
6446  // - Wait on my send-size send.
6447  // Else if I am Process p+2:
6448  // - Post send-size send to Process 0.
6449  //
6450  // Pipelining has three goals here:
6451  // 1. Overlap communication (the receives) with file I/O
6452  // 2. Give Process 0 a chance to prepost some receives,
6453  // before sends show up, by packing local data before
6454  // posting sends
6455  // 3. Don't post _all_ receives or _all_ sends, because that
6456  // wouldn't be memory scalable. (Just because we can't
6457  // see how much memory MPI consumes, doesn't mean that it
6458  // doesn't consume any!)
6459 
6460  // These are used on every process. sendReqSize[0] holds the
6461  // number of rows on this process, and sendReqBuf holds this
6462  // process' data. Process 0 packs into sendReqBuf, but
6463  // doesn't send; it only uses that for printing. All other
6464  // processes send both of these to Process 0.
6465  RCP<CommRequest<int> > sendReqSize, sendReqData;
6466 
6467  // These are used only on Process 0, for received data. Keep
6468  // 3 of each, and treat the arrays as circular buffers. When
6469  // receiving from Process p, the corresponding array index
6470  // here is p % 3.
6471  Array<ArrayRCP<int_type> > recvSizeBufs (3);
6472  Array<ArrayRCP<int_type> > recvDataBufs (3);
6473  Array<RCP<CommRequest<int> > > recvSizeReqs (3);
6474  Array<RCP<CommRequest<int> > > recvDataReqs (3);
6475 
6476  // Buffer for nonblocking send of the "send size."
6477  ArrayRCP<int_type> sendDataSize (1);
6478  sendDataSize[0] = myNumRows;
6479 
6480  if (myRank == 0) {
6481  if (debug) {
6482  std::ostringstream os;
6483  os << myRank << ": Post receive-size receives from "
6484  "Procs 1 and 2: tag = " << sizeTag << endl;
6485  *err << os.str ();
6486  }
6487  // Process 0: Post receive-size receives from Processes 1 and 2.
6488  recvSizeBufs[0].resize (1);
6489  (recvSizeBufs[0])[0] = -1; // error flag
6490  recvSizeBufs[1].resize (1);
6491  (recvSizeBufs[1])[0] = -1; // error flag
6492  recvSizeBufs[2].resize (1);
6493  (recvSizeBufs[2])[0] = -1; // error flag
6494  if (numProcs > 1) {
6495  recvSizeReqs[1] =
6496  ireceive<int, int_type> (recvSizeBufs[1], 1, sizeTag, comm);
6497  }
6498  if (numProcs > 2) {
6499  recvSizeReqs[2] =
6500  ireceive<int, int_type> (recvSizeBufs[2], 2, sizeTag, comm);
6501  }
6502  }
6503  else if (myRank == 1 || myRank == 2) {
6504  if (debug) {
6505  std::ostringstream os;
6506  os << myRank << ": Post send-size send: size = "
6507  << sendDataSize[0] << ", tag = " << sizeTag << endl;
6508  *err << os.str ();
6509  }
6510  // Prime the pipeline by having Processes 1 and 2 start
6511  // their send-size sends. We don't want _all_ the processes
6512  // to start their send-size sends, because that wouldn't be
6513  // memory scalable.
6514  sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
6515  }
6516  else {
6517  if (debug) {
6518  std::ostringstream os;
6519  os << myRank << ": Not posting my send-size send yet" << endl;
6520  *err << os.str ();
6521  }
6522  }
6523 
6524  //
6525  // Pack my GIDs and PIDs. Each (GID,PID) pair gets packed
6526  // consecutively, for better locality.
6527  //
6528 
6529  if (debug) {
6530  std::ostringstream os;
6531  os << myRank << ": Pack my GIDs and PIDs" << endl;
6532  *err << os.str ();
6533  }
6534 
6535  ArrayRCP<int_type> sendDataBuf (myNumRows * 2);
6536 
6537  if (map.isContiguous ()) {
6538  const int_type myMinGblIdx =
6539  static_cast<int_type> (map.getMinGlobalIndex ());
6540  for (size_t k = 0; k < myNumRows; ++k) {
6541  const int_type gid = myMinGblIdx + static_cast<int_type> (k);
6542  const int_type pid = static_cast<int_type> (myRank);
6543  sendDataBuf[2*k] = gid;
6544  sendDataBuf[2*k+1] = pid;
6545  }
6546  }
6547  else {
6548  ArrayView<const GO> myGblInds = map.getNodeElementList ();
6549  for (size_t k = 0; k < myNumRows; ++k) {
6550  const int_type gid = static_cast<int_type> (myGblInds[k]);
6551  const int_type pid = static_cast<int_type> (myRank);
6552  sendDataBuf[2*k] = gid;
6553  sendDataBuf[2*k+1] = pid;
6554  }
6555  }
6556 
6557  if (debug) {
6558  std::ostringstream os;
6559  os << myRank << ": Done packing my GIDs and PIDs" << endl;
6560  *err << os.str ();
6561  }
6562 
6563  if (myRank == 1) {
6564  // Process 1: post send-data send to Process 0.
6565  if (debug) {
6566  *err << myRank << ": Post send-data send: tag = " << dataTag
6567  << endl;
6568  }
6569  sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
6570  }
6571 
6572  if (myRank == 0) {
6573  if (debug) {
6574  *err << myRank << ": Write MatrixMarket header" << endl;
6575  }
6576 
6577  // Process 0: Write the MatrixMarket header.
6578  // Description section explains each column.
6579  std::ostringstream hdr;
6580 
6581  // Print the Matrix Market header. MultiVector stores data
6582  // nonsymmetrically, hence "general" in the banner line.
6583  hdr << "%%MatrixMarket matrix array integer general" << endl
6584  << "% Format: Version 2.0" << endl
6585  << "%" << endl
6586  << "% This file encodes a Tpetra::Map." << endl
6587  << "% It is stored as a dense vector, with twice as many " << endl
6588  << "% entries as the global number of GIDs (global indices)." << endl
6589  << "% (GID, PID) pairs are stored contiguously, where the PID " << endl
6590  << "% is the rank of the process owning that GID." << endl
6591  << (2 * map.getGlobalNumElements ()) << " " << 1 << endl;
6592  out << hdr.str ();
6593 
6594  if (debug) {
6595  std::ostringstream os;
6596  os << myRank << ": Write my GIDs and PIDs" << endl;
6597  *err << os.str ();
6598  }
6599 
6600  // Write Process 0's data to the output stream.
6601  // Matrix Market prints dense matrices in column-major order.
6602  const int_type printNumRows = myNumRows;
6603  ArrayView<const int_type> printData = sendDataBuf ();
6604  for (int_type k = 0; k < printNumRows; ++k) {
6605  const int_type gid = printData[2*k];
6606  const int_type pid = printData[2*k+1];
6607  out << gid << endl << pid << endl;
6608  }
6609  }
6610 
6611  if (myRank == 0) {
6612  // Wait on receive-size receive from Process 1.
6613  const int recvRank = 1;
6614  const int circBufInd = recvRank % 3;
6615  if (debug) {
6616  std::ostringstream os;
6617  os << myRank << ": Wait on receive-size receive from Process "
6618  << recvRank << endl;
6619  *err << os.str ();
6620  }
6621  if (numProcs > 1) {
6622  wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
6623 
6624  // We received the number of rows of data. (The data
6625  // come in two columns.)
6626  const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
6627  if (debug && recvNumRows == -1) {
6628  std::ostringstream os;
6629  os << myRank << ": Result of receive-size receive from Process "
6630  << recvRank << " is -1. This should never happen, and "
6631  "suggests that the receive never got posted. Please report "
6632  "this bug to the Tpetra developers." << endl;
6633  *err << os.str ();
6634  }
6635 
6636  // Post receive-data receive from Process 1.
6637  recvDataBufs[circBufInd].resize (recvNumRows * 2);
6638  if (debug) {
6639  std::ostringstream os;
6640  os << myRank << ": Post receive-data receive from Process "
6641  << recvRank << ": tag = " << dataTag << ", buffer size = "
6642  << recvDataBufs[circBufInd].size () << endl;
6643  *err << os.str ();
6644  }
6645  if (! recvSizeReqs[circBufInd].is_null ()) {
6646  std::ostringstream os;
6647  os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
6648  "null, before posting the receive-data receive from Process "
6649  << recvRank << ". This should never happen. Please report "
6650  "this bug to the Tpetra developers." << endl;
6651  *err << os.str ();
6652  }
6653  recvDataReqs[circBufInd] =
6654  ireceive<int, int_type> (recvDataBufs[circBufInd],
6655  recvRank, dataTag, comm);
6656  } // numProcs > 1
6657  }
6658  else if (myRank == 1) {
6659  // Wait on my send-size send.
6660  if (debug) {
6661  std::ostringstream os;
6662  os << myRank << ": Wait on my send-size send" << endl;
6663  *err << os.str ();
6664  }
6665  wait<int> (comm, outArg (sendReqSize));
6666  }
6667 
6668  //
6669  // Pipeline loop
6670  //
6671  for (int p = 1; p < numProcs; ++p) {
6672  if (myRank == 0) {
6673  if (p + 2 < numProcs) {
6674  // Post receive-size receive from Process p + 2.
6675  const int recvRank = p + 2;
6676  const int circBufInd = recvRank % 3;
6677  if (debug) {
6678  std::ostringstream os;
6679  os << myRank << ": Post receive-size receive from Process "
6680  << recvRank << ": tag = " << sizeTag << endl;
6681  *err << os.str ();
6682  }
6683  if (! recvSizeReqs[circBufInd].is_null ()) {
6684  std::ostringstream os;
6685  os << myRank << ": recvSizeReqs[" << circBufInd << "] is not "
6686  << "null, for the receive-size receive from Process "
6687  << recvRank << "! This may mean that this process never "
6688  << "finished waiting for the receive from Process "
6689  << (recvRank - 3) << "." << endl;
6690  *err << os.str ();
6691  }
6692  recvSizeReqs[circBufInd] =
6693  ireceive<int, int_type> (recvSizeBufs[circBufInd],
6694  recvRank, sizeTag, comm);
6695  }
6696 
6697  if (p + 1 < numProcs) {
6698  const int recvRank = p + 1;
6699  const int circBufInd = recvRank % 3;
6700 
6701  // Wait on receive-size receive from Process p + 1.
6702  if (debug) {
6703  std::ostringstream os;
6704  os << myRank << ": Wait on receive-size receive from Process "
6705  << recvRank << endl;
6706  *err << os.str ();
6707  }
6708  wait<int> (comm, outArg (recvSizeReqs[circBufInd]));
6709 
6710  // We received the number of rows of data. (The data
6711  // come in two columns.)
6712  const int_type recvNumRows = (recvSizeBufs[circBufInd])[0];
6713  if (debug && recvNumRows == -1) {
6714  std::ostringstream os;
6715  os << myRank << ": Result of receive-size receive from Process "
6716  << recvRank << " is -1. This should never happen, and "
6717  "suggests that the receive never got posted. Please report "
6718  "this bug to the Tpetra developers." << endl;
6719  *err << os.str ();
6720  }
6721 
6722  // Post receive-data receive from Process p + 1.
6723  recvDataBufs[circBufInd].resize (recvNumRows * 2);
6724  if (debug) {
6725  std::ostringstream os;
6726  os << myRank << ": Post receive-data receive from Process "
6727  << recvRank << ": tag = " << dataTag << ", buffer size = "
6728  << recvDataBufs[circBufInd].size () << endl;
6729  *err << os.str ();
6730  }
6731  if (! recvDataReqs[circBufInd].is_null ()) {
6732  std::ostringstream os;
6733  os << myRank << ": recvDataReqs[" << circBufInd << "] is not "
6734  << "null, for the receive-data receive from Process "
6735  << recvRank << "! This may mean that this process never "
6736  << "finished waiting for the receive from Process "
6737  << (recvRank - 3) << "." << endl;
6738  *err << os.str ();
6739  }
6740  recvDataReqs[circBufInd] =
6741  ireceive<int, int_type> (recvDataBufs[circBufInd],
6742  recvRank, dataTag, comm);
6743  }
6744 
6745  // Wait on receive-data receive from Process p.
6746  const int recvRank = p;
6747  const int circBufInd = recvRank % 3;
6748  if (debug) {
6749  std::ostringstream os;
6750  os << myRank << ": Wait on receive-data receive from Process "
6751  << recvRank << endl;
6752  *err << os.str ();
6753  }
6754  wait<int> (comm, outArg (recvDataReqs[circBufInd]));
6755 
6756  // Write Process p's data. Number of rows lives in
6757  // recvSizeBufs[circBufInd], and the actual data live in
6758  // recvDataBufs[circBufInd]. Do this after posting receives,
6759  // in order to expose overlap of comm. with file I/O.
6760  if (debug) {
6761  std::ostringstream os;
6762  os << myRank << ": Write GIDs and PIDs from Process "
6763  << recvRank << endl;
6764  *err << os.str () << endl;
6765  }
6766  const int_type printNumRows = (recvSizeBufs[circBufInd])[0];
6767  if (debug && printNumRows == -1) {
6768  std::ostringstream os;
6769  os << myRank << ": Result of receive-size receive from Process "
6770  << recvRank << " was -1. This should never happen, and "
6771  "suggests that its receive-size receive was never posted. "
6772  "Please report this bug to the Tpetra developers." << endl;
6773  *err << os.str ();
6774  }
6775  if (debug && printNumRows > 0 && recvDataBufs[circBufInd].is_null ()) {
6776  std::ostringstream os;
6777  os << myRank << ": Result of receive-size receive from Proc "
6778  << recvRank << " was " << printNumRows << " > 0, but "
6779  "recvDataBufs[" << circBufInd << "] is null. This should "
6780  "never happen. Please report this bug to the Tpetra "
6781  "developers." << endl;
6782  *err << os.str ();
6783  }
6784  ArrayView<const int_type> printData = (recvDataBufs[circBufInd]) ();
6785  for (int_type k = 0; k < printNumRows; ++k) {
6786  const int_type gid = printData[2*k];
6787  const int_type pid = printData[2*k+1];
6788  out << gid << endl << pid << endl;
6789  }
6790  }
6791  else if (myRank == p) { // Process p
6792  // Wait on my send-data send.
6793  if (debug) {
6794  std::ostringstream os;
6795  os << myRank << ": Wait on my send-data send" << endl;
6796  *err << os.str ();
6797  }
6798  wait<int> (comm, outArg (sendReqData));
6799  }
6800  else if (myRank == p + 1) { // Process p + 1
6801  // Post send-data send to Process 0.
6802  if (debug) {
6803  std::ostringstream os;
6804  os << myRank << ": Post send-data send: tag = " << dataTag
6805  << endl;
6806  *err << os.str ();
6807  }
6808  sendReqData = isend<int, int_type> (sendDataBuf, 0, dataTag, comm);
6809  // Wait on my send-size send.
6810  if (debug) {
6811  std::ostringstream os;
6812  os << myRank << ": Wait on my send-size send" << endl;
6813  *err << os.str ();
6814  }
6815  wait<int> (comm, outArg (sendReqSize));
6816  }
6817  else if (myRank == p + 2) { // Process p + 2
6818  // Post send-size send to Process 0.
6819  if (debug) {
6820  std::ostringstream os;
6821  os << myRank << ": Post send-size send: size = "
6822  << sendDataSize[0] << ", tag = " << sizeTag << endl;
6823  *err << os.str ();
6824  }
6825  sendReqSize = isend<int, int_type> (sendDataSize, 0, sizeTag, comm);
6826  }
6827  }
6828 
6829  if (! err.is_null ()) {
6830  err->popTab ();
6831  }
6832  if (debug) {
6833  *err << myRank << ": writeMap: Done" << endl;
6834  }
6835  if (! err.is_null ()) {
6836  err->popTab ();
6837  }
6838  }
6839 
6841  static void
6842  writeMapFile (const std::string& filename,
6843  const map_type& map)
6844  {
6845  const int myRank = map.getComm ()->getRank ();
6846  std::ofstream out;
6847  if (myRank == 0) { // Only open the file on Proc 0.
6848  out.open (filename.c_str());
6849  }
6850  writeMap (out, map);
6851  // We can rely on the destructor of the output stream to close
6852  // the file on scope exit, even if writeDense() throws an
6853  // exception.
6854  }
6855 
6856  private:
6880  static void
6881  printAsComment (std::ostream& out, const std::string& str)
6882  {
6883  using std::endl;
6884  std::istringstream inpstream (str);
6885  std::string line;
6886 
6887  while (getline (inpstream, line)) {
6888  if (! line.empty()) {
6889  // Note that getline() doesn't store '\n', so we have to
6890  // append the endline ourselves.
6891  if (line[0] == '%') { // Line starts with a comment character.
6892  out << line << endl;
6893  }
6894  else { // Line doesn't start with a comment character.
6895  out << "%% " << line << endl;
6896  }
6897  }
6898  }
6899  }
6900 
6901  public:
6902 
6921  static void
6922  writeOperator(const std::string& fileName, operator_type const &A) {
6923  Teuchos::ParameterList pl;
6924  writeOperator(fileName, A, pl);
6925  }
6926 
6947  static void
6948  writeOperator (std::ostream& out, const operator_type& A) {
6949  Teuchos::ParameterList pl;
6950  writeOperator (out, A, pl);
6951  }
6952 
6989  static void
6990  writeOperator (const std::string& fileName,
6991  const operator_type& A,
6992  const Teuchos::ParameterList& params)
6993  {
6994  std::ofstream out;
6995  std::string tmpFile = "__TMP__" + fileName;
6996  const int myRank = A.getDomainMap()->getComm()->getRank();
6997  bool precisionChanged=false;
6998  int oldPrecision;
6999  // The number of nonzero entries in a Tpetra::Operator is
7000  // unknown until probing is completed. In order to write a
7001  // MatrixMarket header, we write the matrix to a temporary
7002  // file.
7003  //
7004  // FIXME (mfh 23 May 2015) IT WASN'T MY IDEA TO WRITE TO A
7005  // TEMPORARY FILE.
7006  if (myRank==0) {
7007  if (std::ifstream(tmpFile))
7008  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
7009  "writeOperator: temporary file " << tmpFile << " already exists");
7010  out.open(tmpFile.c_str());
7011  if (params.isParameter("precision")) {
7012  oldPrecision = out.precision(params.get<int>("precision"));
7013  precisionChanged=true;
7014  }
7015  }
7016 
7017  const std::string header = writeOperatorImpl(out, A, params);
7018 
7019  if (myRank==0) {
7020  if (precisionChanged)
7021  out.precision(oldPrecision);
7022  out.close();
7023  out.open(fileName.c_str(), std::ios::binary);
7024  bool printMatrixMarketHeader = true;
7025  if (params.isParameter("print MatrixMarket header"))
7026  printMatrixMarketHeader = params.get<bool>("print MatrixMarket header");
7027  if (printMatrixMarketHeader && myRank == 0) {
7028  // Write header to final file.
7029  out << header;
7030  }
7031  // Append matrix from temporary to final file.
7032  std::ifstream src(tmpFile, std::ios_base::binary);
7033  out << src.rdbuf();
7034  src.close();
7035  // Delete the temporary file.
7036  remove(tmpFile.c_str());
7037  }
7038  }
7039 
7078  static void
7079  writeOperator (std::ostream& out,
7080  const operator_type& A,
7081  const Teuchos::ParameterList& params)
7082  {
7083  const int myRank = A.getDomainMap ()->getComm ()->getRank ();
7084 
7085  // The number of nonzero entries in a Tpetra::Operator is
7086  // unknown until probing is completed. In order to write a
7087  // MatrixMarket header, we write the matrix to a temporary
7088  // output stream.
7089  //
7090  // NOTE (mfh 23 May 2015): Writing to a temporary output
7091  // stream may double the memory usage, depending on whether
7092  // 'out' is a file stream or an in-memory output stream (e.g.,
7093  // std::ostringstream). It might be wise to use a temporary
7094  // file instead. However, please look carefully at POSIX
7095  // functions for safe creation of temporary files. Don't just
7096  // prepend "__TMP__" to the filename and hope for the best.
7097  // Furthermore, it should be valid to call the std::ostream
7098  // overload of this method even when Process 0 does not have
7099  // access to a file system.
7100  std::ostringstream tmpOut;
7101  if (myRank == 0) {
7102  if (params.isParameter ("precision") && params.isType<int> ("precision")) {
7103  (void) tmpOut.precision (params.get<int> ("precision"));
7104  }
7105  }
7106 
7107  const std::string header = writeOperatorImpl (tmpOut, A, params);
7108 
7109  if (myRank == 0) {
7110  bool printMatrixMarketHeader = true;
7111  if (params.isParameter ("print MatrixMarket header") &&
7112  params.isType<bool> ("print MatrixMarket header")) {
7113  printMatrixMarketHeader = params.get<bool> ("print MatrixMarket header");
7114  }
7115  if (printMatrixMarketHeader && myRank == 0) {
7116  out << header; // write header to final output stream
7117  }
7118  // Append matrix from temporary output stream to final output stream.
7119  //
7120  // NOTE (mfh 23 May 2015) This might use a lot of memory.
7121  // However, we should not use temporary files in this
7122  // method. Since it does not access the file system (unlike
7123  // the overload that takes a file name), it should not
7124  // require the file system at all.
7125  //
7126  // If memory usage becomes a problem, one thing we could do
7127  // is write the entries of the Operator one column (or a few
7128  // columns) at a time. The Matrix Market sparse format does
7129  // not impose an order on its entries, so it would be OK to
7130  // write them in that order.
7131  out << tmpOut.str ();
7132  }
7133  }
7134 
7135  private:
7136 
7144  static std::string
7145  writeOperatorImpl (std::ostream& os,
7146  const operator_type& A,
7147  const Teuchos::ParameterList& params)
7148  {
7149  using Teuchos::RCP;
7150  using Teuchos::rcp;
7151  using Teuchos::ArrayRCP;
7152  using Teuchos::Array;
7153 
7154  typedef local_ordinal_type LO;
7155  typedef global_ordinal_type GO;
7156  typedef scalar_type Scalar;
7157  typedef Teuchos::OrdinalTraits<LO> TLOT;
7158  typedef Teuchos::OrdinalTraits<GO> TGOT;
7159  typedef Tpetra::Import<LO, GO, node_type> import_type;
7160  typedef Tpetra::MultiVector<GO, LO, GO, node_type> mv_type_go;
7161 
7162  const map_type& domainMap = *(A.getDomainMap());
7163  RCP<const map_type> rangeMap = A.getRangeMap();
7164  RCP<const Teuchos::Comm<int> > comm = rangeMap->getComm();
7165  const int myRank = comm->getRank();
7166  const size_t numProcs = comm->getSize();
7167 
7168  size_t numMVs = 10;
7169  if (params.isParameter("probing size"))
7170  numMVs = params.get<int>("probing size");
7171 
7172  GO globalNnz = 0;
7173  GO minColGid = domainMap.getMinAllGlobalIndex();
7174  GO maxColGid = domainMap.getMaxAllGlobalIndex();
7175  // Rather than replicating the domainMap on all processors, we instead
7176  // iterate from the min GID to the max GID. If the map is gappy,
7177  // there will be invalid GIDs, i.e., GIDs no one has. This will require
7178  // unnecessary matvecs against potentially zero vectors.
7179  GO numGlobElts = maxColGid - minColGid + TGOT::one();
7180  GO numChunks = numGlobElts / numMVs;
7181  GO rem = numGlobElts % numMVs;
7182  GO indexBase = rangeMap->getIndexBase();
7183 
7184  int offsetToUseInPrinting = 1 - indexBase; // default is 1-based indexing
7185  if (params.isParameter("zero-based indexing")) {
7186  if (params.get<bool>("zero-based indexing") == true)
7187  offsetToUseInPrinting = -indexBase; // If 0-based, use as-is. If 1-based, subtract 1.
7188  }
7189 
7190  // Create map that replicates the range map on pid 0 and is empty for all other pids
7191  size_t numLocalRangeEntries = rangeMap->getNodeNumElements();
7192 
7193  // Create contiguous source map
7194  RCP<const map_type> allGidsMap = rcp(new map_type(TGOT::invalid(), numLocalRangeEntries,
7195  indexBase, comm));
7196  // Create vector based on above map. Populate it with GIDs corresponding to this pid's GIDs in rangeMap.
7197  mv_type_go allGids(allGidsMap,1);
7198  Teuchos::ArrayRCP<GO> allGidsData = allGids.getDataNonConst(0);
7199 
7200  for (size_t i=0; i<numLocalRangeEntries; i++)
7201  allGidsData[i] = rangeMap->getGlobalElement(i);
7202  allGidsData = Teuchos::null;
7203 
7204  // Create target map that is nontrivial only on pid 0
7205  GO numTargetMapEntries=TGOT::zero();
7206  Teuchos::Array<GO> importGidList;
7207  if (myRank==0) {
7208  numTargetMapEntries = rangeMap->getGlobalNumElements();
7209  importGidList.reserve(numTargetMapEntries);
7210  for (GO j=0; j<numTargetMapEntries; ++j) importGidList.push_back(j + indexBase);
7211  } else {
7212  importGidList.reserve(numTargetMapEntries);
7213  }
7214  RCP<map_type> importGidMap = rcp(new map_type(TGOT::invalid(), importGidList(), indexBase, comm));
7215 
7216  // Import all rangeMap GIDs to pid 0
7217  import_type gidImporter(allGidsMap, importGidMap);
7218  mv_type_go importedGids(importGidMap, 1);
7219  importedGids.doImport(allGids, gidImporter, INSERT);
7220 
7221  // The following import map will be non-trivial only on pid 0.
7222  ArrayRCP<const GO> importedGidsData = importedGids.getData(0);
7223  RCP<const map_type> importMap = rcp(new map_type(TGOT::invalid(), importedGidsData(), indexBase, comm) );
7224 
7225  // Importer from original range map to pid 0
7226  import_type importer(rangeMap, importMap);
7227  // Target vector on pid 0
7228  RCP<mv_type> colsOnPid0 = rcp(new mv_type(importMap,numMVs));
7229 
7230  RCP<mv_type> ei = rcp(new mv_type(A.getDomainMap(),numMVs)); //probing vector
7231  RCP<mv_type> colsA = rcp(new mv_type(A.getRangeMap(),numMVs)); //columns of A revealed by probing
7232 
7233  Array<GO> globalColsArray, localColsArray;
7234  globalColsArray.reserve(numMVs);
7235  localColsArray.reserve(numMVs);
7236 
7237  ArrayRCP<ArrayRCP<Scalar> > eiData(numMVs);
7238  for (size_t i=0; i<numMVs; ++i)
7239  eiData[i] = ei->getDataNonConst(i);
7240 
7241  // //////////////////////////////////////
7242  // Discover A by chunks
7243  // //////////////////////////////////////
7244  for (GO k=0; k<numChunks; ++k) {
7245  for (size_t j=0; j<numMVs; ++j ) {
7246  //GO curGlobalCol = maxColGid - numMVs + j + TGOT::one();
7247  GO curGlobalCol = minColGid + k*numMVs + j;
7248  globalColsArray.push_back(curGlobalCol);
7249  //TODO extract the g2l map outside of this loop loop
7250  LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
7251  if (curLocalCol != TLOT::invalid()) {
7252  eiData[j][curLocalCol] = TGOT::one();
7253  localColsArray.push_back(curLocalCol);
7254  }
7255  }
7256  //TODO Do the views eiData need to be released prior to the matvec?
7257 
7258  // probe
7259  A.apply(*ei,*colsA);
7260 
7261  colsOnPid0->doImport(*colsA,importer,INSERT);
7262 
7263  if (myRank==0)
7264  globalNnz += writeColumns(os,*colsOnPid0, numMVs, importedGidsData(),
7265  globalColsArray, offsetToUseInPrinting);
7266 
7267  //zero out the ei's
7268  for (size_t j=0; j<numMVs; ++j ) {
7269  for (int i=0; i<localColsArray.size(); ++i)
7270  eiData[j][localColsArray[i]] = TGOT::zero();
7271  }
7272  globalColsArray.clear();
7273  localColsArray.clear();
7274 
7275  }
7276 
7277  // //////////////////////////////////////
7278  // Handle leftover part of A
7279  // //////////////////////////////////////
7280  if (rem > 0) {
7281  for (int j=0; j<rem; ++j ) {
7282  GO curGlobalCol = maxColGid - rem + j + TGOT::one();
7283  globalColsArray.push_back(curGlobalCol);
7284  //TODO extract the g2l map outside of this loop loop
7285  LO curLocalCol = domainMap.getLocalElement(curGlobalCol);
7286  if (curLocalCol != TLOT::invalid()) {
7287  eiData[j][curLocalCol] = TGOT::one();
7288  localColsArray.push_back(curLocalCol);
7289  }
7290  }
7291  //TODO Do the views eiData need to be released prior to the matvec?
7292 
7293  // probe
7294  A.apply(*ei,*colsA);
7295 
7296  colsOnPid0->doImport(*colsA,importer,INSERT);
7297  if (myRank==0)
7298  globalNnz += writeColumns(os,*colsOnPid0, rem, importedGidsData(),
7299  globalColsArray, offsetToUseInPrinting);
7300 
7301  //zero out the ei's
7302  for (int j=0; j<rem; ++j ) {
7303  for (int i=0; i<localColsArray.size(); ++i)
7304  eiData[j][localColsArray[i]] = TGOT::zero();
7305  }
7306  globalColsArray.clear();
7307  localColsArray.clear();
7308 
7309  }
7310 
7311  // Return the Matrix Market header. It includes the header
7312  // line (that starts with "%%"), some comments, and the triple
7313  // of matrix dimensions and number of nonzero entries. We
7314  // don't actually print this here, because we don't know the
7315  // number of nonzero entries until after probing.
7316  std::ostringstream oss;
7317  if (myRank == 0) {
7318  oss << "%%MatrixMarket matrix coordinate ";
7319  if (Teuchos::ScalarTraits<typename operator_type::scalar_type>::isComplex) {
7320  oss << "complex";
7321  } else {
7322  oss << "real";
7323  }
7324  oss << " general" << std::endl;
7325  oss << "% Tpetra::Operator" << std::endl;
7326  std::time_t now = std::time(NULL);
7327  oss << "% time stamp: " << ctime(&now);
7328  oss << "% written from " << numProcs << " processes" << std::endl;
7329  size_t numRows = rangeMap->getGlobalNumElements();
7330  size_t numCols = domainMap.getGlobalNumElements();
7331  oss << numRows << " " << numCols << " " << globalNnz << std::endl;
7332  }
7333 
7334  return oss.str ();
7335  }
7336 
7337  static global_ordinal_type
7338  writeColumns(std::ostream& os, mv_type const &colsA, size_t const &numCols,
7339  Teuchos::ArrayView<const global_ordinal_type> const &rowGids,
7340  Teuchos::Array<global_ordinal_type> const &colsArray,
7341  global_ordinal_type const & indexBase) {
7342 
7343  typedef global_ordinal_type GO;
7344  typedef scalar_type Scalar;
7345  typedef Teuchos::ScalarTraits<Scalar> STS;
7346 
7347  GO nnz=0;
7348  const Scalar zero = STS::zero();
7349  const size_t numRows = colsA.getGlobalLength();
7350  for (size_t j=0; j<numCols; ++j) {
7351  ArrayRCP<const Scalar> const curCol = colsA.getData(j);
7352  const GO J = colsArray[j];
7353  for (size_t i=0; i<numRows; ++i) {
7354  const Scalar val = curCol[i];
7355  if (val!=zero) {
7356  os << rowGids[i]+indexBase << " " << J+indexBase << " " << val << std::endl;
7357  ++nnz;
7358  }
7359  }
7360  }
7361 
7362  return nnz;
7363 
7364  }
7365 
7366  public:
7367 
7368  }; // class Writer
7369 
7370  } // namespace MatrixMarket
7371 } // namespace Tpetra
7372 
7373 #endif // __MatrixMarket_Tpetra_hpp
size_t getLocalLength() const
Local number of rows on the calling process.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
static RCP< multivector_type > readDense(std::istream &in, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market input stream.
void get1dCopy(const Teuchos::ArrayView< Scalar > &A, const size_t LDA) const
Fill the given array with a copy of this multivector&#39;s local values.
static void writeOperator(std::ostream &out, const operator_type &A)
Write a Tpetra::Operator to an output stream.
GlobalOrdinal getMinGlobalIndex() const
The minimum global index owned by the calling process.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
SparseMatrixType::node_type node_type
The fourth template parameter of CrsMatrix and MultiVector.
static void writeDense(std::ostream &out, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
bool isContiguous() const
True if this Map is distributed contiguously, else false.
static RCP< multivector_type > readDense(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Variant of readDense (see above) that takes a Node.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const RCP< const Comm< int > > &pComm, const RCP< node_type > &pNode, const RCP< Teuchos::ParameterList > &constructorParams, const RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of the above readSparse() method that takes a Kokkos Node.
static void writeDenseFile(const std::string &filename, const RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Variant of readSparseFile above that takes a Node object.
SparseMatrixType::local_ordinal_type local_ordinal_type
SparseMatrixType::global_ordinal_type global_ordinal_type
One or more distributed dense vectors.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream.
SparseMatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the sparse matrix.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< Teuchos::ParameterList > &constructorParams, const Teuchos::RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
static void writeDense(std::ostream &out, const RCP< const multivector_type > &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and or description.
virtual void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const =0
Computes the operator-multivector application.
static void writeDenseFile(const std::string &filename, const Teuchos::RCP< const multivector_type > &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
SparseMatrixType::scalar_type scalar_type
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Accessors for the Teuchos::Comm and Kokkos Node objects.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
Specialization of Tpetra::MultiVector that matches SparseMatrixType.
static RCP< const map_type > readMapFile(const std::string &filename, const RCP< const comm_type > &comm, const RCP< node_type > &node, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given Matrix Market file.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
GlobalOrdinal getMinAllGlobalIndex() const
The minimum global index over all processes in the communicator.
static void writeMap(std::ostream &out, const map_type &map, const bool debug=false)
Print the Map to the given output stream.
static void writeDenseFile(const std::string &filename, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > getVector(const size_t j) const
Return a Vector which is a const view of column j.
SparseMatrixType::global_ordinal_type global_ordinal_type
Type of indices as read from the Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const map_type > &rowMap, Teuchos::RCP< const map_type > &colMap, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream, with provided Maps. ...
GlobalOrdinal getMaxAllGlobalIndex() const
The maximum global index over all processes in the communicator.
size_t global_size_t
Global size_t object.
SparseMatrixType::node_type node_type
The Kokkos Node type; fourth template parameter of Tpetra::CrsMatrix.
static void writeOperator(const std::string &fileName, const operator_type &A, const Teuchos::ParameterList &params)
Write a Tpetra::Operator to a file, with options.
Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Specialization of Tpetra::Map that matches SparseMatrixType.
Insert new values that don&#39;t currently exist.
SparseMatrixType::scalar_type scalar_type
Type of the entries of the sparse matrix.
Scalar scalar_type
This class&#39; first template parameter; the type of each entry in the MultiVector.
static void writeMapFile(const std::string &filename, const map_type &map)
Write the Map to the given file.
Abstract interface for operators (e.g., matrices and preconditioners).
static RCP< vector_type > readVector(std::istream &in, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
static RCP< multivector_type > readDenseFile(const std::string &filename, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Variant of readDenseMatrix (see above) that takes a Node.
From a distributed map build a map with all GIDs on the root node.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
static RCP< multivector_type > readDenseFile(const std::string &filename, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read dense matrix (as a MultiVector) from the given Matrix Market file.
static void writeMap(std::ostream &out, const map_type &map, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool debug=false)
Print the Map to the given output stream out.
Teuchos::ArrayRCP< const Scalar > getData(size_t j) const
Const view of the local values in a particular vector of this multivector.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const RCP< const Comm< int > > &pComm, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const RCP< const Comm< int > > &pComm, const RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparseFile that takes a Node object.
static void writeDense(std::ostream &out, const multivector_type &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static void writeSparse(std::ostream &out, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const std::string &matrixName, const std::string &matrixDescription, const bool debug=false)
Print the sparse matrix in Matrix Market format, with comments.
static void writeSparseFile(const std::string &filename, const Teuchos::RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
global_size_t getGlobalLength() const
Global number of rows in the multivector.
MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > multivector_type
The MultiVector specialization associated with SparseMatrixType.
LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const
The local index corresponding to the given global index.
Teuchos::ArrayView< const GlobalOrdinal > getNodeElementList() const
Return a view of the global indices owned by this process.
Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > vector_type
The Vector specialization associated with SparseMatrixType.
size_t getNodeNumElements() const
The number of elements belonging to the calling process.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const RCP< const Comm< int > > &pComm, const RCP< Teuchos::ParameterList > &constructorParams, const RCP< Teuchos::ParameterList > &fillCompleteParams, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market input stream.
Describes a parallel distribution of objects over processes.
static void writeDense(std::ostream &out, const RCP< const multivector_type > &X, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with no matrix name or description.
static RCP< vector_type > readVector(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read Vector from the given Matrix Market input stream, with a supplied Node.
static Teuchos::RCP< sparse_matrix_type > readSparse(std::istream &in, const Teuchos::RCP< const Teuchos::Comm< int > > &pComm, const Teuchos::RCP< node_type > &pNode, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Variant of readSparse() above that takes a Node object.
A distributed dense vector.
static RCP< const map_type > readMap(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream.
size_t getNumVectors() const
Number of columns in the multivector.
SparseMatrixType sparse_matrix_type
This class&#39; template parameter; a specialization of CrsMatrix.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
SparseMatrixType sparse_matrix_type
Template parameter of this class; specialization of CrsMatrix.
static void writeSparse(std::ostream &out, const RCP< const sparse_matrix_type > &pMatrix, const bool debug=false)
Print the sparse matrix in Matrix Market format.
static void writeDenseFile(const std::string &filename, const multivector_type &X, const std::string &matrixName, const std::string &matrixDescription, const Teuchos::RCP< Teuchos::FancyOStream > &err=Teuchos::null, const Teuchos::RCP< Teuchos::FancyOStream > &dbg=Teuchos::null)
Print the multivector in Matrix Market format, with matrix name and description.
static RCP< vector_type > readVectorFile(const std::string &filename, const RCP< const comm_type > &comm, const RCP< node_type > &node, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Like readVectorFile() (see above), but with a supplied Node object.
static RCP< vector_type > readVectorFile(const std::string &filename, const RCP< const comm_type > &comm, RCP< const map_type > &map, const bool tolerant=false, const bool debug=false)
Read a Vector from the given Matrix Market file.
static void writeOperator(std::ostream &out, const operator_type &A, const Teuchos::ParameterList &params)
Write a Tpetra::Operator to an output stream, with options.
Matrix Market file reader for CrsMatrix and MultiVector.
static RCP< const map_type > readMap(std::istream &in, const RCP< const comm_type > &comm, const RCP< node_type > &node, const Teuchos::RCP< Teuchos::FancyOStream > &err, const bool tolerant=false, const bool debug=false)
Read Map (as a MultiVector) from the given input stream, with optional debugging output stream...
Matrix Market file writer for CrsMatrix and MultiVector.
static Teuchos::RCP< sparse_matrix_type > readSparseFile(const std::string &filename, const RCP< const map_type > &rowMap, RCP< const map_type > &colMap, const RCP< const map_type > &domainMap, const RCP< const map_type > &rangeMap, const bool callFillComplete=true, const bool tolerant=false, const bool debug=false)
Read sparse matrix from the given Matrix Market file, with provided Maps.
Matrix Market file readers and writers for sparse and dense matrices (as CrsMatrix resp...
global_size_t getGlobalNumElements() const
The number of elements in this Map.
static void writeOperator(const std::string &fileName, operator_type const &A)
Write a Tpetra::Operator to a file.