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  typedef Teuchos::ScalarTraits<Scalar> ST;
62  using Teuchos::as;
63  TEUCHOS_TEST_FOR_EXCEPTION( plist == Teuchos::null, std::runtime_error,
64  "Tpetra::Utils::generateMatrix(): ParameterList is null.");
65  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::isParameterType<std::string>(*plist,"mat_type") == false, std::runtime_error,
66  "Tpetra::Utils::generateMatrix(): ParameterList did not contain string parameter ""mat_type"".");
67  std::string mat_type = plist->get<std::string>("mat_type");
68  if (mat_type == "Lap3D") {
69  // 3D Laplacian, grid is a cube with dimension gridSize x gridSize x gridSize
70  const GlobalOrdinal gridSize = as<GlobalOrdinal>(plist->get<int>("gridSize",100));
71  const GlobalOrdinal gS2 = gridSize*gridSize;
72  const GlobalOrdinal numRows = gS2*gridSize;
73  Teuchos::RCP<Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap;
74  rowMap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(as<global_size_t>(numRows),as<GlobalOrdinal>(0),comm,GloballyDistributed,node));
76  // fill matrix, one row at a time
77  Teuchos::Array<GlobalOrdinal> neighbors;
78  Teuchos::Array<Scalar> values(7, -ST::one());
79  values[0] = (Scalar)6;
80  for (GlobalOrdinal r = rowMap->getMinGlobalIndex(); r <= rowMap->getMaxGlobalIndex(); ++r) {
81  neighbors.clear();
82  neighbors.push_back(r); // add diagonal
83  GlobalOrdinal ixy, iz, ix, iy; // (x,y,z) coords and index in xy plane
84  ixy = r%gS2;
85  iz = (r - ixy)/gS2;
86  ix = ixy%gridSize;
87  iy = (ixy - ix)/gridSize;
88  //
89  if ( ix != 0 ) neighbors.push_back( r-1 );
90  if ( ix != gridSize-1 ) neighbors.push_back( r+1 );
91  if ( iy != 0 ) neighbors.push_back( r-gridSize );
92  if ( iy != gridSize-1 ) neighbors.push_back( r+gridSize );
93  if ( iz != 0 ) neighbors.push_back( r-gS2 );
94  if ( iz != gridSize-1 ) neighbors.push_back( r+gS2 );
95  A->insertGlobalValues( r, neighbors(), values(0,neighbors.size()) );
96  }
97  A->fillComplete();
98  }
99  else {
100  TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,
101  "Tpetra::Utils::generateMatrix(): ParameterList specified unsupported ""mat_type"".");
102  }
103 }
104 
105 
106 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
107 void
108 readHBMatrix (const std::string &filename,
109  const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
110  const Teuchos::RCP<Node> &node,
112  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap,
113  const Teuchos::RCP<Teuchos::ParameterList> &params)
114 {
115  const int myRank = comm->getRank();
116  int numRows,numCols,numNZ;
117  Teuchos::ArrayRCP<Scalar> svals;
118  Teuchos::ArrayRCP<GlobalOrdinal> colinds;
119  Teuchos::ArrayRCP<int> rowptrs;
120  Teuchos::ArrayRCP<size_t> nnzPerRow;
121  int fail = 0;
122  if (myRank == 0) {
123  bool isSymmetric=false;
124  Teuchos::ArrayRCP<double> dvals;
125  Teuchos::ArrayRCP<int> colptrs, rowinds;
126  std::string type;
127  Tpetra::Utils::readHBMatDouble(filename,numRows,numCols,numNZ,type,colptrs,rowinds,dvals);
128  TEUCHOS_TEST_FOR_EXCEPT(type.size() != 3);
129  if (type[0] != 'R' && type[0] != 'r') {
130  // only real matrices right now
131  fail = 1;
132  }
133  if (fail == 0 && numNZ > 0) {
134  if (type[1] == 'S' || type[1] == 's') {
135  isSymmetric = true;
136  }
137  else {
138  isSymmetric = false;
139  }
140  }
141  if (fail == 0 && numNZ > 0) {
142  // find num non-zero per row
143  nnzPerRow = Teuchos::arcp<size_t>(numRows);
144  std::fill(nnzPerRow.begin(), nnzPerRow.end(), 0);
145  for (Teuchos::ArrayRCP<int>::const_iterator ri=rowinds.begin(); ri != rowinds.end(); ++ri) {
146  // count each row index towards its row
147  ++nnzPerRow[*ri-1];
148  }
149  if (isSymmetric) {
150  // count each column toward the corresponding row as well
151  for (int c=0; c < numCols; ++c) {
152  // the diagonal was already counted; neglect it, if it exists
153  for (int i=colptrs[c]-1; i != colptrs[c+1]-1; ++i) {
154  if (rowinds[i] != c+1) {
155  ++nnzPerRow[c];
156  ++numNZ;
157  }
158  }
159  }
160  }
161  // allocate/set new matrix data
162  svals = Teuchos::arcp<Scalar>(numNZ);
163  colinds = Teuchos::arcp<GlobalOrdinal>(numNZ);
164  rowptrs = Teuchos::arcp<int>(numRows+1);
165  rowptrs[0] = 0;
166 #ifdef HAVE_TPETRA_DEBUG
167  Teuchos::ArrayRCP<size_t> nnzPerRow_debug(nnzPerRow.size());
168  std::copy(nnzPerRow.begin(), nnzPerRow.end(), nnzPerRow_debug.begin());
169 #endif
170  for (int j=1; j <= numRows; ++j) {
171  rowptrs[j] = rowptrs[j-1] + nnzPerRow[j-1];
172  nnzPerRow[j-1] = 0;
173  }
174  // translate from column-oriented to row-oriented
175  for (int col=0; col<numCols; ++col) {
176  for (int i=colptrs[col]-1; i != colptrs[col+1]-1; ++i) {
177  const int row = rowinds[i]-1;
178  // add entry to (row,col), with value dvals[i]
179  const size_t entry = rowptrs[row] + nnzPerRow[row];
180  svals[entry] = Teuchos::as<Scalar>(dvals[i]);
181  colinds[entry] = Teuchos::as<GlobalOrdinal>(col);
182  ++nnzPerRow[row];
183  if (isSymmetric && row != col) {
184  // add entry to (col,row), with value dvals[i]
185  const size_t symentry = rowptrs[col] + nnzPerRow[col];
186  svals[symentry] = Teuchos::as<Scalar>(dvals[i]);
187  colinds[symentry] = Teuchos::as<GlobalOrdinal>(row);
188  ++nnzPerRow[col];
189  }
190  }
191  }
192 #ifdef HAVE_TPETRA_DEBUG
193  {
194  bool isequal = true;
195  typename Teuchos::ArrayRCP<size_t>::const_iterator it1, it2;
196  for (it1 = nnzPerRow.begin(), it2 = nnzPerRow_debug.begin(); it1 != nnzPerRow.end(); ++it1, ++it2) {
197  if (*it1 != *it2) {
198  isequal = false;
199  break;
200  }
201  }
202  TEUCHOS_TEST_FOR_EXCEPTION(!isequal || nnzPerRow.size() != nnzPerRow_debug.size(), std::logic_error,
203  "Tpetra::Utils::readHBMatrix(): Logic error.");
204  }
205 #endif
206  }
207  // std::cout << "Matrix " << filename << " of type " << type << ": " << numRows << " by " << numCols << ", " << numNZ << " nonzeros" << std::endl;
208  }
209  // check for read errors
210  broadcast(*comm,0,&fail);
211  TEUCHOS_TEST_FOR_EXCEPTION(fail == 1, std::runtime_error, "Tpetra::Utils::readHBMatrix() can only read Real matrices.");
212  // distribute global matrix info
213  broadcast(*comm,0,&numRows);
214  broadcast(*comm,0,&numCols);
215  // create map with uniform partitioning
216  if (rowMap == Teuchos::null) {
217  rowMap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>((global_size_t)numRows,(GlobalOrdinal)0,comm,GloballyDistributed,node));
218  }
219  else {
220  TEUCHOS_TEST_FOR_EXCEPTION( rowMap->getGlobalNumElements() != (global_size_t)numRows, std::runtime_error,
221  "Tpetra::Utils::readHBMatrix(): specified map has incorrect number of elements.");
222  TEUCHOS_TEST_FOR_EXCEPTION( rowMap->isDistributed() == false && comm->getSize() > 1, std::runtime_error,
223  "Tpetra::Utils::readHBMatrix(): specified map is not distributed.");
224  }
225  Teuchos::ArrayRCP<size_t> myNNZ;
226  if (rowMap->getNodeNumElements()) {
227  myNNZ = Teuchos::arcp<size_t>(rowMap->getNodeNumElements());
228  }
229  if (myRank == 0) {
230  LocalOrdinal numRowsAlreadyDistributed = rowMap->getNodeNumElements();
231  std::copy(nnzPerRow.begin(), nnzPerRow.begin()+numRowsAlreadyDistributed,myNNZ);
232  for (int p=1; p < Teuchos::size(*comm); ++p) {
233  size_t numRowsForP;
234  Teuchos::receive(*comm,p,&numRowsForP);
235  if (numRowsForP) {
236  Teuchos::send<int,size_t>(*comm,numRowsForP,nnzPerRow(numRowsAlreadyDistributed,numRowsForP).getRawPtr(),p);
237  numRowsAlreadyDistributed += numRowsForP;
238  }
239  }
240  }
241  else {
242  const size_t numMyRows = rowMap->getNodeNumElements();
243  Teuchos::send(*comm,numMyRows,0);
244  if (numMyRows) {
245  Teuchos::receive<int,size_t>(*comm,0,numMyRows,myNNZ(0,numMyRows).getRawPtr());
246  }
247  }
248  nnzPerRow = Teuchos::null;
249  // create column map
250  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap;
251  if (numRows == numCols) {
252  domMap = rowMap;
253  }
254  else {
255  domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm,node);
256  }
258  // free this locally, A will keep it allocated as long as it is needed by A (up until allocation of nonzeros)
259  myNNZ = Teuchos::null;
260  if (myRank == 0 && numNZ > 0) {
261  for (int r=0; r < numRows; ++r) {
262  const LocalOrdinal nnz = rowptrs[r+1] - rowptrs[r];
263  if (nnz > 0) {
264  Teuchos::ArrayView<const GlobalOrdinal> inds = colinds(rowptrs[r],nnz);
265  Teuchos::ArrayView<const Scalar> vals = svals( rowptrs[r],nnz);
266  A->insertGlobalValues(r, inds, vals);
267  }
268  }
269  }
270  // don't need these anymore
271  colinds = Teuchos::null;
272  svals = Teuchos::null;
273  rowptrs = Teuchos::null;
274  A->fillComplete(domMap,rowMap,params);
275 }
276 } // namespace Utils
277 } // namespace Tpetra
278 
279 //
280 // Explicit instantiation macro
281 //
282 // Must be expanded from within the Tpetra::Utils namespace!
283 //
284 
285 #define TPETRA_MATRIXIO_INSTANT(SCALAR,LO,GO,NODE) \
286  template void \
287  readHBMatrix< SCALAR, LO, GO, NODE > (const std::string&, \
288  const Teuchos::RCP<const Teuchos::Comm<int> > &, \
289  const Teuchos::RCP< NODE > &, \
290  Teuchos::RCP<CrsMatrix< SCALAR, LO, GO, NODE > >&, \
291  Teuchos::RCP<const Tpetra::Map< LO, GO, NODE> >, \
292  const Teuchos::RCP<Teuchos::ParameterList>& ); \
293  \
294  template void \
295  generateMatrix< SCALAR, LO, GO, NODE> (const Teuchos::RCP<Teuchos::ParameterList>&, \
296  const Teuchos::RCP<const Teuchos::Comm<int> > &, \
297  const Teuchos::RCP< NODE > &,\
298  Teuchos::RCP<CrsMatrix< SCALAR, LO, GO, NODE > >& );
299 
300 
301 #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.