44 #ifndef TPETRA_MATRIX_IO_DEF 45 #define TPETRA_MATRIX_IO_DEF 47 #include "Tpetra_CrsMatrix.hpp" 48 #include "Tpetra_MatrixIO.hpp" 54 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
56 generateMatrix (
const Teuchos::RCP<Teuchos::ParameterList> &plist,
57 const Teuchos::RCP<
const Teuchos::Comm<int> > &comm,
58 const Teuchos::RCP<Node> &node,
62 typedef Teuchos::ScalarTraits<Scalar> ST;
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") {
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));
82 rowMap = Teuchos::rcp (
new map_type (static_cast<global_size_t> (numRows),
83 static_cast<GlobalOrdinal> (0),
84 comm, GloballyDistributed, node));
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) {
93 neighbors.push_back(r);
94 GlobalOrdinal ixy, iz, ix, iy;
98 iy = (ixy - ix)/gridSize;
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()) );
111 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
112 "Tpetra::Utils::generateMatrix(): ParameterList specified unsupported ""mat_type"".");
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 readHBMatrix (
const std::string &filename,
120 const Teuchos::RCP<
const Teuchos::Comm<int> > &comm,
121 const Teuchos::RCP<Node> &node,
124 const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
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;
136 bool isSymmetric=
false;
137 Teuchos::ArrayRCP<double> dvals;
138 Teuchos::ArrayRCP<int> colptrs, rowinds;
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') {
146 if (fail == 0 && numNZ > 0) {
147 if (type[1] ==
'S' || type[1] ==
's') {
154 if (fail == 0 && numNZ > 0) {
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) {
164 for (
int c=0; c < numCols; ++c) {
166 for (
int i=colptrs[c]-1; i != colptrs[c+1]-1; ++i) {
167 if (rowinds[i] != c+1) {
175 svals = Teuchos::arcp<Scalar>(numNZ);
176 colinds = Teuchos::arcp<GlobalOrdinal>(numNZ);
177 rowptrs = Teuchos::arcp<int>(numRows+1);
179 #ifdef HAVE_TPETRA_DEBUG 180 Teuchos::ArrayRCP<size_t> nnzPerRow_debug(nnzPerRow.size());
181 std::copy(nnzPerRow.begin(), nnzPerRow.end(), nnzPerRow_debug.begin());
183 for (
int j=1; j <= numRows; ++j) {
184 rowptrs[j] = rowptrs[j-1] + nnzPerRow[j-1];
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;
192 const size_t entry = rowptrs[row] + nnzPerRow[row];
193 svals[entry] = Teuchos::as<Scalar>(dvals[i]);
194 colinds[entry] = Teuchos::as<GlobalOrdinal>(col);
196 if (isSymmetric && row != col) {
198 const size_t symentry = rowptrs[col] + nnzPerRow[col];
199 svals[symentry] = Teuchos::as<Scalar>(dvals[i]);
200 colinds[symentry] = Teuchos::as<GlobalOrdinal>(row);
205 #ifdef HAVE_TPETRA_DEBUG 208 typename Teuchos::ArrayRCP<size_t>::const_iterator it1, it2;
209 for (it1 = nnzPerRow.begin(), it2 = nnzPerRow_debug.begin(); it1 != nnzPerRow.end(); ++it1, ++it2) {
215 TEUCHOS_TEST_FOR_EXCEPTION(!isequal || nnzPerRow.size() != nnzPerRow_debug.size(), std::logic_error,
216 "Tpetra::Utils::readHBMatrix(): Logic error.");
223 broadcast(*comm,0,&fail);
224 TEUCHOS_TEST_FOR_EXCEPTION(fail == 1, std::runtime_error,
"Tpetra::Utils::readHBMatrix() can only read Real matrices.");
226 broadcast(*comm,0,&numRows);
227 broadcast(*comm,0,&numCols);
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));
236 rowMap = Teuchos::rcp (
new map_type (static_cast<global_size_t> (numRows),
237 static_cast<GlobalOrdinal> (0),
238 comm, GloballyDistributed, node));
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.");
247 Teuchos::ArrayRCP<size_t> myNNZ;
248 if (rowMap->getNodeNumElements()) {
249 myNNZ = Teuchos::arcp<size_t>(rowMap->getNodeNumElements());
252 LocalOrdinal numRowsAlreadyDistributed = rowMap->getNodeNumElements();
253 std::copy(nnzPerRow.begin(), nnzPerRow.begin()+numRowsAlreadyDistributed,myNNZ);
254 for (
int p=1; p < Teuchos::size(*comm); ++p) {
256 Teuchos::receive(*comm,p,&numRowsForP);
258 Teuchos::send<int,size_t>(*comm,numRowsForP,nnzPerRow(numRowsAlreadyDistributed,numRowsForP).getRawPtr(),p);
259 numRowsAlreadyDistributed += numRowsForP;
264 const size_t numMyRows = rowMap->getNodeNumElements();
265 Teuchos::send(*comm,numMyRows,0);
267 Teuchos::receive<int,size_t>(*comm,0,numMyRows,myNNZ(0,numMyRows).getRawPtr());
270 nnzPerRow = Teuchos::null;
272 Teuchos::RCP<const map_type> domMap;
273 if (numRows == numCols) {
277 if (node.is_null ()) {
278 domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm);
281 domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm,node);
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];
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);
298 colinds = Teuchos::null;
299 svals = Teuchos::null;
300 rowptrs = Teuchos::null;
301 A->fillComplete(domMap,rowMap,params);
312 #define TPETRA_MATRIXIO_INSTANT(SCALAR,LO,GO,NODE) \ 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>& ); \ 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 > >& ); 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.