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,
61 typedef Teuchos::ScalarTraits<Scalar> ST;
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") {
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;
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) {
82 neighbors.push_back(r);
83 GlobalOrdinal ixy, iz, ix, iy;
87 iy = (ixy - ix)/gridSize;
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()) );
100 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
101 "Tpetra::Utils::generateMatrix(): ParameterList specified unsupported ""mat_type"".");
106 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
108 readHBMatrix (
const std::string &filename,
109 const Teuchos::RCP<
const Teuchos::Comm<int> > &comm,
110 const Teuchos::RCP<Node> &node,
113 const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
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;
123 bool isSymmetric=
false;
124 Teuchos::ArrayRCP<double> dvals;
125 Teuchos::ArrayRCP<int> colptrs, rowinds;
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') {
133 if (fail == 0 && numNZ > 0) {
134 if (type[1] ==
'S' || type[1] ==
's') {
141 if (fail == 0 && numNZ > 0) {
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) {
151 for (
int c=0; c < numCols; ++c) {
153 for (
int i=colptrs[c]-1; i != colptrs[c+1]-1; ++i) {
154 if (rowinds[i] != c+1) {
162 svals = Teuchos::arcp<Scalar>(numNZ);
163 colinds = Teuchos::arcp<GlobalOrdinal>(numNZ);
164 rowptrs = Teuchos::arcp<int>(numRows+1);
166 #ifdef HAVE_TPETRA_DEBUG 167 Teuchos::ArrayRCP<size_t> nnzPerRow_debug(nnzPerRow.size());
168 std::copy(nnzPerRow.begin(), nnzPerRow.end(), nnzPerRow_debug.begin());
170 for (
int j=1; j <= numRows; ++j) {
171 rowptrs[j] = rowptrs[j-1] + nnzPerRow[j-1];
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;
179 const size_t entry = rowptrs[row] + nnzPerRow[row];
180 svals[entry] = Teuchos::as<Scalar>(dvals[i]);
181 colinds[entry] = Teuchos::as<GlobalOrdinal>(col);
183 if (isSymmetric && row != col) {
185 const size_t symentry = rowptrs[col] + nnzPerRow[col];
186 svals[symentry] = Teuchos::as<Scalar>(dvals[i]);
187 colinds[symentry] = Teuchos::as<GlobalOrdinal>(row);
192 #ifdef HAVE_TPETRA_DEBUG 195 typename Teuchos::ArrayRCP<size_t>::const_iterator it1, it2;
196 for (it1 = nnzPerRow.begin(), it2 = nnzPerRow_debug.begin(); it1 != nnzPerRow.end(); ++it1, ++it2) {
202 TEUCHOS_TEST_FOR_EXCEPTION(!isequal || nnzPerRow.size() != nnzPerRow_debug.size(), std::logic_error,
203 "Tpetra::Utils::readHBMatrix(): Logic error.");
210 broadcast(*comm,0,&fail);
211 TEUCHOS_TEST_FOR_EXCEPTION(fail == 1, std::runtime_error,
"Tpetra::Utils::readHBMatrix() can only read Real matrices.");
213 broadcast(*comm,0,&numRows);
214 broadcast(*comm,0,&numCols);
216 if (rowMap == Teuchos::null) {
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.");
225 Teuchos::ArrayRCP<size_t> myNNZ;
226 if (rowMap->getNodeNumElements()) {
227 myNNZ = Teuchos::arcp<size_t>(rowMap->getNodeNumElements());
230 LocalOrdinal numRowsAlreadyDistributed = rowMap->getNodeNumElements();
231 std::copy(nnzPerRow.begin(), nnzPerRow.begin()+numRowsAlreadyDistributed,myNNZ);
232 for (
int p=1; p < Teuchos::size(*comm); ++p) {
234 Teuchos::receive(*comm,p,&numRowsForP);
236 Teuchos::send<int,size_t>(*comm,numRowsForP,nnzPerRow(numRowsAlreadyDistributed,numRowsForP).getRawPtr(),p);
237 numRowsAlreadyDistributed += numRowsForP;
242 const size_t numMyRows = rowMap->getNodeNumElements();
243 Teuchos::send(*comm,numMyRows,0);
245 Teuchos::receive<int,size_t>(*comm,0,numMyRows,myNNZ(0,numMyRows).getRawPtr());
248 nnzPerRow = Teuchos::null;
250 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap;
251 if (numRows == numCols) {
255 domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm,node);
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];
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);
271 colinds = Teuchos::null;
272 svals = Teuchos::null;
273 rowptrs = Teuchos::null;
274 A->fillComplete(domMap,rowMap,params);
285 #define TPETRA_MATRIXIO_INSTANT(SCALAR,LO,GO,NODE) \ 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>& ); \ 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 > >& ); 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.