Tpetra parallel linear algebra  Version of the Day
Tpetra_MatrixIO_def.hpp
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Tpetra: Templated Linear Algebra Services Package
6 // Copyright (2008) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 */
43 
44 #ifndef TPETRA_MATRIX_IO_DEF
45 #define TPETRA_MATRIX_IO_DEF
46 
47 #include "Tpetra_CrsMatrix.hpp"
48 #include "Tpetra_MatrixIO.hpp"
49 #include <iostream>
50 
51 namespace Tpetra {
52 namespace Utils {
53 
54 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
55 void
56 generateMatrix (const Teuchos::RCP<Teuchos::ParameterList> &plist,
57  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
58  const Teuchos::RCP<Node> &node,
60 {
61  using Teuchos::as;
62  typedef Teuchos::ScalarTraits<Scalar> ST;
64 
65  TEUCHOS_TEST_FOR_EXCEPTION( plist == Teuchos::null, std::runtime_error,
66  "Tpetra::Utils::generateMatrix(): ParameterList is null.");
67  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::isParameterType<std::string>(*plist,"mat_type") == false, std::runtime_error,
68  "Tpetra::Utils::generateMatrix(): ParameterList did not contain string parameter ""mat_type"".");
69  std::string mat_type = plist->get<std::string>("mat_type");
70  if (mat_type == "Lap3D") {
71  // 3D Laplacian, grid is a cube with dimension gridSize x gridSize x gridSize
72  const GlobalOrdinal gridSize = as<GlobalOrdinal>(plist->get<int>("gridSize",100));
73  const GlobalOrdinal gS2 = gridSize*gridSize;
74  const GlobalOrdinal numRows = gS2*gridSize;
75  Teuchos::RCP<map_type> rowMap;
76  if (node.is_null ()) {
77  rowMap = Teuchos::rcp (new map_type (static_cast<global_size_t> (numRows),
78  static_cast<GlobalOrdinal> (0),
79  comm, GloballyDistributed));
80  }
81  else {
82  rowMap = Teuchos::rcp (new map_type (static_cast<global_size_t> (numRows),
83  static_cast<GlobalOrdinal> (0),
84  comm, GloballyDistributed, node));
85  }
87  // fill matrix, one row at a time
88  Teuchos::Array<GlobalOrdinal> neighbors;
89  Teuchos::Array<Scalar> values(7, -ST::one());
90  values[0] = (Scalar)6;
91  for (GlobalOrdinal r = rowMap->getMinGlobalIndex(); r <= rowMap->getMaxGlobalIndex(); ++r) {
92  neighbors.clear();
93  neighbors.push_back(r); // add diagonal
94  GlobalOrdinal ixy, iz, ix, iy; // (x,y,z) coords and index in xy plane
95  ixy = r%gS2;
96  iz = (r - ixy)/gS2;
97  ix = ixy%gridSize;
98  iy = (ixy - ix)/gridSize;
99  //
100  if ( ix != 0 ) neighbors.push_back( r-1 );
101  if ( ix != gridSize-1 ) neighbors.push_back( r+1 );
102  if ( iy != 0 ) neighbors.push_back( r-gridSize );
103  if ( iy != gridSize-1 ) neighbors.push_back( r+gridSize );
104  if ( iz != 0 ) neighbors.push_back( r-gS2 );
105  if ( iz != gridSize-1 ) neighbors.push_back( r+gS2 );
106  A->insertGlobalValues( r, neighbors(), values(0,neighbors.size()) );
107  }
108  A->fillComplete();
109  }
110  else {
111  TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,
112  "Tpetra::Utils::generateMatrix(): ParameterList specified unsupported ""mat_type"".");
113  }
114 }
115 
116 
117 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118 void
119 readHBMatrix (const std::string &filename,
120  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
121  const Teuchos::RCP<Node> &node,
123  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap,
124  const Teuchos::RCP<Teuchos::ParameterList> &params)
125 {
127 
128  const int myRank = comm->getRank();
129  int numRows,numCols,numNZ;
130  Teuchos::ArrayRCP<Scalar> svals;
131  Teuchos::ArrayRCP<GlobalOrdinal> colinds;
132  Teuchos::ArrayRCP<int> rowptrs;
133  Teuchos::ArrayRCP<size_t> nnzPerRow;
134  int fail = 0;
135  if (myRank == 0) {
136  bool isSymmetric=false;
137  Teuchos::ArrayRCP<double> dvals;
138  Teuchos::ArrayRCP<int> colptrs, rowinds;
139  std::string type;
140  Tpetra::Utils::readHBMatDouble(filename,numRows,numCols,numNZ,type,colptrs,rowinds,dvals);
141  TEUCHOS_TEST_FOR_EXCEPT(type.size() != 3);
142  if (type[0] != 'R' && type[0] != 'r') {
143  // only real matrices right now
144  fail = 1;
145  }
146  if (fail == 0 && numNZ > 0) {
147  if (type[1] == 'S' || type[1] == 's') {
148  isSymmetric = true;
149  }
150  else {
151  isSymmetric = false;
152  }
153  }
154  if (fail == 0 && numNZ > 0) {
155  // find num non-zero per row
156  nnzPerRow = Teuchos::arcp<size_t>(numRows);
157  std::fill(nnzPerRow.begin(), nnzPerRow.end(), 0);
158  for (Teuchos::ArrayRCP<int>::const_iterator ri=rowinds.begin(); ri != rowinds.end(); ++ri) {
159  // count each row index towards its row
160  ++nnzPerRow[*ri-1];
161  }
162  if (isSymmetric) {
163  // count each column toward the corresponding row as well
164  for (int c=0; c < numCols; ++c) {
165  // the diagonal was already counted; neglect it, if it exists
166  for (int i=colptrs[c]-1; i != colptrs[c+1]-1; ++i) {
167  if (rowinds[i] != c+1) {
168  ++nnzPerRow[c];
169  ++numNZ;
170  }
171  }
172  }
173  }
174  // allocate/set new matrix data
175  svals = Teuchos::arcp<Scalar>(numNZ);
176  colinds = Teuchos::arcp<GlobalOrdinal>(numNZ);
177  rowptrs = Teuchos::arcp<int>(numRows+1);
178  rowptrs[0] = 0;
179 #ifdef HAVE_TPETRA_DEBUG
180  Teuchos::ArrayRCP<size_t> nnzPerRow_debug(nnzPerRow.size());
181  std::copy(nnzPerRow.begin(), nnzPerRow.end(), nnzPerRow_debug.begin());
182 #endif
183  for (int j=1; j <= numRows; ++j) {
184  rowptrs[j] = rowptrs[j-1] + nnzPerRow[j-1];
185  nnzPerRow[j-1] = 0;
186  }
187  // translate from column-oriented to row-oriented
188  for (int col=0; col<numCols; ++col) {
189  for (int i=colptrs[col]-1; i != colptrs[col+1]-1; ++i) {
190  const int row = rowinds[i]-1;
191  // add entry to (row,col), with value dvals[i]
192  const size_t entry = rowptrs[row] + nnzPerRow[row];
193  svals[entry] = Teuchos::as<Scalar>(dvals[i]);
194  colinds[entry] = Teuchos::as<GlobalOrdinal>(col);
195  ++nnzPerRow[row];
196  if (isSymmetric && row != col) {
197  // add entry to (col,row), with value dvals[i]
198  const size_t symentry = rowptrs[col] + nnzPerRow[col];
199  svals[symentry] = Teuchos::as<Scalar>(dvals[i]);
200  colinds[symentry] = Teuchos::as<GlobalOrdinal>(row);
201  ++nnzPerRow[col];
202  }
203  }
204  }
205 #ifdef HAVE_TPETRA_DEBUG
206  {
207  bool isequal = true;
208  typename Teuchos::ArrayRCP<size_t>::const_iterator it1, it2;
209  for (it1 = nnzPerRow.begin(), it2 = nnzPerRow_debug.begin(); it1 != nnzPerRow.end(); ++it1, ++it2) {
210  if (*it1 != *it2) {
211  isequal = false;
212  break;
213  }
214  }
215  TEUCHOS_TEST_FOR_EXCEPTION(!isequal || nnzPerRow.size() != nnzPerRow_debug.size(), std::logic_error,
216  "Tpetra::Utils::readHBMatrix(): Logic error.");
217  }
218 #endif
219  }
220  // std::cout << "Matrix " << filename << " of type " << type << ": " << numRows << " by " << numCols << ", " << numNZ << " nonzeros" << std::endl;
221  }
222  // check for read errors
223  broadcast(*comm,0,&fail);
224  TEUCHOS_TEST_FOR_EXCEPTION(fail == 1, std::runtime_error, "Tpetra::Utils::readHBMatrix() can only read Real matrices.");
225  // distribute global matrix info
226  broadcast(*comm,0,&numRows);
227  broadcast(*comm,0,&numCols);
228  // create map with uniform partitioning
229  if (rowMap == Teuchos::null) {
230  if (node.is_null ()) {
231  rowMap = Teuchos::rcp (new map_type (static_cast<global_size_t> (numRows),
232  static_cast<GlobalOrdinal> (0),
233  comm, GloballyDistributed));
234  }
235  else {
236  rowMap = Teuchos::rcp (new map_type (static_cast<global_size_t> (numRows),
237  static_cast<GlobalOrdinal> (0),
238  comm, GloballyDistributed, node));
239  }
240  }
241  else {
242  TEUCHOS_TEST_FOR_EXCEPTION( rowMap->getGlobalNumElements() != (global_size_t)numRows, std::runtime_error,
243  "Tpetra::Utils::readHBMatrix(): specified map has incorrect number of elements.");
244  TEUCHOS_TEST_FOR_EXCEPTION( rowMap->isDistributed() == false && comm->getSize() > 1, std::runtime_error,
245  "Tpetra::Utils::readHBMatrix(): specified map is not distributed.");
246  }
247  Teuchos::ArrayRCP<size_t> myNNZ;
248  if (rowMap->getNodeNumElements()) {
249  myNNZ = Teuchos::arcp<size_t>(rowMap->getNodeNumElements());
250  }
251  if (myRank == 0) {
252  LocalOrdinal numRowsAlreadyDistributed = rowMap->getNodeNumElements();
253  std::copy(nnzPerRow.begin(), nnzPerRow.begin()+numRowsAlreadyDistributed,myNNZ);
254  for (int p=1; p < Teuchos::size(*comm); ++p) {
255  size_t numRowsForP;
256  Teuchos::receive(*comm,p,&numRowsForP);
257  if (numRowsForP) {
258  Teuchos::send<int,size_t>(*comm,numRowsForP,nnzPerRow(numRowsAlreadyDistributed,numRowsForP).getRawPtr(),p);
259  numRowsAlreadyDistributed += numRowsForP;
260  }
261  }
262  }
263  else {
264  const size_t numMyRows = rowMap->getNodeNumElements();
265  Teuchos::send(*comm,numMyRows,0);
266  if (numMyRows) {
267  Teuchos::receive<int,size_t>(*comm,0,numMyRows,myNNZ(0,numMyRows).getRawPtr());
268  }
269  }
270  nnzPerRow = Teuchos::null;
271  // create column map
272  Teuchos::RCP<const map_type> domMap;
273  if (numRows == numCols) {
274  domMap = rowMap;
275  }
276  else {
277  if (node.is_null ()) {
278  domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm);
279  }
280  else {
281  domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm,node);
282  }
283  }
285  // free this locally, A will keep it allocated as long as it is needed by A (up until allocation of nonzeros)
286  myNNZ = Teuchos::null;
287  if (myRank == 0 && numNZ > 0) {
288  for (int r=0; r < numRows; ++r) {
289  const LocalOrdinal nnz = rowptrs[r+1] - rowptrs[r];
290  if (nnz > 0) {
291  Teuchos::ArrayView<const GlobalOrdinal> inds = colinds(rowptrs[r],nnz);
292  Teuchos::ArrayView<const Scalar> vals = svals( rowptrs[r],nnz);
293  A->insertGlobalValues(r, inds, vals);
294  }
295  }
296  }
297  // don't need these anymore
298  colinds = Teuchos::null;
299  svals = Teuchos::null;
300  rowptrs = Teuchos::null;
301  A->fillComplete(domMap,rowMap,params);
302 }
303 } // namespace Utils
304 } // namespace Tpetra
305 
306 //
307 // Explicit instantiation macro
308 //
309 // Must be expanded from within the Tpetra::Utils namespace!
310 //
311 
312 #define TPETRA_MATRIXIO_INSTANT(SCALAR,LO,GO,NODE) \
313  template void \
314  readHBMatrix< SCALAR, LO, GO, NODE > (const std::string&, \
315  const Teuchos::RCP<const Teuchos::Comm<int> > &, \
316  const Teuchos::RCP< NODE > &, \
317  Teuchos::RCP<CrsMatrix< SCALAR, LO, GO, NODE > >&, \
318  Teuchos::RCP<const Tpetra::Map< LO, GO, NODE> >, \
319  const Teuchos::RCP<Teuchos::ParameterList>& ); \
320  \
321  template void \
322  generateMatrix< SCALAR, LO, GO, NODE> (const Teuchos::RCP<Teuchos::ParameterList>&, \
323  const Teuchos::RCP<const Teuchos::Comm<int> > &, \
324  const Teuchos::RCP< NODE > &,\
325  Teuchos::RCP<CrsMatrix< SCALAR, LO, GO, NODE > >& );
326 
327 
328 #endif
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
size_t global_size_t
Global size_t object.
Describes a parallel distribution of objects over processes.