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