Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
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
51namespace Tpetra {
52namespace Utils {
53
54
55template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
56void
57readHBMatrix (const std::string &filename,
58 const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
60 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap,
61 const Teuchos::RCP<Teuchos::ParameterList> &params)
62{
64
65 const int myRank = comm->getRank();
66 int numRows,numCols,numNZ;
67 Teuchos::ArrayRCP<Scalar> svals;
68 Teuchos::ArrayRCP<GlobalOrdinal> colinds;
69 Teuchos::ArrayRCP<int> rowptrs;
70 Teuchos::ArrayRCP<size_t> nnzPerRow;
71 int fail = 0;
72 if (myRank == 0) {
73 bool isSymmetric=false;
74 Teuchos::ArrayRCP<double> dvals;
75 Teuchos::ArrayRCP<int> colptrs, rowinds;
76 std::string type;
77 Tpetra::Utils::readHBMatDouble(filename,numRows,numCols,numNZ,type,colptrs,rowinds,dvals);
78 TEUCHOS_TEST_FOR_EXCEPT(type.size() != 3);
79 if (type[0] != 'R' && type[0] != 'r') {
80 // only real matrices right now
81 fail = 1;
82 }
83 if (fail == 0 && numNZ > 0) {
84 if (type[1] == 'S' || type[1] == 's') {
85 isSymmetric = true;
86 }
87 else {
88 isSymmetric = false;
89 }
90 }
91 if (fail == 0 && numNZ > 0) {
92 // find num non-zero per row
93 nnzPerRow = Teuchos::arcp<size_t>(numRows);
94 std::fill(nnzPerRow.begin(), nnzPerRow.end(), 0);
95 for (Teuchos::ArrayRCP<int>::const_iterator ri=rowinds.begin(); ri != rowinds.end(); ++ri) {
96 // count each row index towards its row
97 ++nnzPerRow[*ri-1];
98 }
99 if (isSymmetric) {
100 // count each column toward the corresponding row as well
101 for (int c=0; c < numCols; ++c) {
102 // the diagonal was already counted; neglect it, if it exists
103 for (int i=colptrs[c]-1; i != colptrs[c+1]-1; ++i) {
104 if (rowinds[i] != c+1) {
105 ++nnzPerRow[c];
106 ++numNZ;
107 }
108 }
109 }
110 }
111 // allocate/set new matrix data
112 svals = Teuchos::arcp<Scalar>(numNZ);
113 colinds = Teuchos::arcp<GlobalOrdinal>(numNZ);
114 rowptrs = Teuchos::arcp<int>(numRows+1);
115 rowptrs[0] = 0;
116#ifdef HAVE_TPETRA_DEBUG
117 Teuchos::ArrayRCP<size_t> nnzPerRow_debug(nnzPerRow.size());
118 std::copy(nnzPerRow.begin(), nnzPerRow.end(), nnzPerRow_debug.begin());
119#endif
120 for (int j=1; j <= numRows; ++j) {
121 rowptrs[j] = rowptrs[j-1] + nnzPerRow[j-1];
122 nnzPerRow[j-1] = 0;
123 }
124 // translate from column-oriented to row-oriented
125 for (int col=0; col<numCols; ++col) {
126 for (int i=colptrs[col]-1; i != colptrs[col+1]-1; ++i) {
127 const int row = rowinds[i]-1;
128 // add entry to (row,col), with value dvals[i]
129 const size_t entry = rowptrs[row] + nnzPerRow[row];
130 svals[entry] = Teuchos::as<Scalar>(dvals[i]);
131 colinds[entry] = Teuchos::as<GlobalOrdinal>(col);
132 ++nnzPerRow[row];
133 if (isSymmetric && row != col) {
134 // add entry to (col,row), with value dvals[i]
135 const size_t symentry = rowptrs[col] + nnzPerRow[col];
136 svals[symentry] = Teuchos::as<Scalar>(dvals[i]);
137 colinds[symentry] = Teuchos::as<GlobalOrdinal>(row);
138 ++nnzPerRow[col];
139 }
140 }
141 }
142#ifdef HAVE_TPETRA_DEBUG
143 {
144 bool isequal = true;
145 typename Teuchos::ArrayRCP<size_t>::const_iterator it1, it2;
146 for (it1 = nnzPerRow.begin(), it2 = nnzPerRow_debug.begin(); it1 != nnzPerRow.end(); ++it1, ++it2) {
147 if (*it1 != *it2) {
148 isequal = false;
149 break;
150 }
151 }
152 TEUCHOS_TEST_FOR_EXCEPTION(!isequal || nnzPerRow.size() != nnzPerRow_debug.size(), std::logic_error,
153 "Tpetra::Utils::readHBMatrix(): Logic error.");
154 }
155#endif
156 }
157 // std::cout << "Matrix " << filename << " of type " << type << ": " << numRows << " by " << numCols << ", " << numNZ << " nonzeros" << std::endl;
158 }
159 // check for read errors
160 broadcast(*comm,0,&fail);
161 TEUCHOS_TEST_FOR_EXCEPTION(fail == 1, std::runtime_error, "Tpetra::Utils::readHBMatrix() can only read Real matrices.");
162 // distribute global matrix info
163 broadcast(*comm,0,&numRows);
164 broadcast(*comm,0,&numCols);
165 // create map with uniform partitioning
166 if (rowMap == Teuchos::null) {
167 rowMap = Teuchos::rcp (new map_type (static_cast<global_size_t> (numRows),
168 static_cast<GlobalOrdinal> (0),
169 comm, GloballyDistributed));
170 }
171 else {
172 TEUCHOS_TEST_FOR_EXCEPTION( rowMap->getGlobalNumElements() != (global_size_t)numRows, std::runtime_error,
173 "Tpetra::Utils::readHBMatrix(): specified map has incorrect number of elements.");
174 TEUCHOS_TEST_FOR_EXCEPTION( rowMap->isDistributed() == false && comm->getSize() > 1, std::runtime_error,
175 "Tpetra::Utils::readHBMatrix(): specified map is not distributed.");
176 }
177 Teuchos::Array<size_t> myNNZ (rowMap->getLocalNumElements ());
178 if (myRank == 0) {
179 LocalOrdinal numRowsAlreadyDistributed = rowMap->getLocalNumElements();
180 std::copy(nnzPerRow.begin(), nnzPerRow.begin()+numRowsAlreadyDistributed, myNNZ.begin());
181 for (int p=1; p < Teuchos::size(*comm); ++p) {
182 size_t numRowsForP;
183 Teuchos::receive(*comm,p,&numRowsForP);
184 if (numRowsForP) {
185 Teuchos::send<int,size_t>(*comm,numRowsForP,nnzPerRow(numRowsAlreadyDistributed,numRowsForP).getRawPtr(),p);
186 numRowsAlreadyDistributed += numRowsForP;
187 }
188 }
189 }
190 else {
191 const size_t numMyRows = rowMap->getLocalNumElements();
192 Teuchos::send(*comm,numMyRows,0);
193 if (numMyRows) {
194 Teuchos::receive<int,size_t>(*comm,0,numMyRows,myNNZ(0,numMyRows).getRawPtr());
195 }
196 }
197 nnzPerRow = Teuchos::null;
198 // create column map
199 Teuchos::RCP<const map_type> domMap;
200 if (numRows == numCols) {
201 domMap = rowMap;
202 }
203 else {
204 domMap = createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numCols,comm);
205 }
207 // free this locally, A will keep it allocated as long as it is needed by A (up until allocation of nonzeros)
208 {
209 // Classic idiom for freeing an std::vector; resize doesn't
210 // promise deallocation.
211 Teuchos::Array<size_t> empty;
212 std::swap (myNNZ, empty);
213 }
214 if (myRank == 0 && numNZ > 0) {
215 for (int r=0; r < numRows; ++r) {
216 const LocalOrdinal nnz = rowptrs[r+1] - rowptrs[r];
217 if (nnz > 0) {
218 Teuchos::ArrayView<const GlobalOrdinal> inds = colinds(rowptrs[r],nnz);
219 Teuchos::ArrayView<const Scalar> vals = svals( rowptrs[r],nnz);
220 A->insertGlobalValues(r, inds, vals);
221 }
222 }
223 }
224 // don't need these anymore
225 colinds = Teuchos::null;
226 svals = Teuchos::null;
227 rowptrs = Teuchos::null;
228 A->fillComplete(domMap,rowMap,params);
229}
230
231
232} // namespace Utils
233} // namespace Tpetra
234
235//
236// Explicit instantiation macro
237//
238// Must be expanded from within the Tpetra::Utils namespace!
239//
240
241
242#define TPETRA_MATRIXIO_INSTANT(SCALAR,LO,GO,NODE) \
243 template void \
244 readHBMatrix< SCALAR, LO, GO, NODE > (const std::string&, \
245 const Teuchos::RCP<const Teuchos::Comm<int> > &, \
246 Teuchos::RCP<CrsMatrix< SCALAR, LO, GO, NODE > >&, \
247 Teuchos::RCP<const Tpetra::Map< LO, GO, NODE> >, \
248 const Teuchos::RCP<Teuchos::ParameterList>& );
249
250
251
252#endif
Sparse matrix that presents a row-oriented interface that lets users read or modify entries.
A parallel distribution of indices over processes.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
size_t global_size_t
Global size_t object.