FEI Version of the Day
Loading...
Searching...
No Matches
fei_Matrix_Local.cpp
1/*--------------------------------------------------------------------*/
2/* Copyright 2007 Sandia Corporation. */
3/* Under the terms of Contract DE-AC04-94AL85000, there is a */
4/* non-exclusive license for use of this work by or on behalf */
5/* of the U.S. Government. Export of this program may require */
6/* a license from the United States Government. */
7/*--------------------------------------------------------------------*/
8
9#include <fei_ParameterSet.hpp>
10#include "fei_Matrix_Local.hpp"
11#include "fei_Matrix_core.hpp"
12#include "fei_sstream.hpp"
13#include "fei_fstream.hpp"
14
15namespace fei {
16
17Matrix_Local::Matrix_Local(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
19 : matrixGraph_(matrixGraph),
20 sparseRowGraph_(sparseRowGraph),
21 coefs_(sparseRowGraph->packedColumnIndices.size(), 0.0),
22 stateChanged_(false),
23 work_data1D_(),
24 work_data2D_()
25{
26}
27
28Matrix_Local::~Matrix_Local()
29{
30}
31
32//----------------------------------------------------------------------------
34Matrix_Local::create_Matrix_Local(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
35 bool blockEntry)
36{
38 matrixGraph->createGraph(blockEntry, true);
39 fei::SharedPtr<fei::Matrix> mat(new fei::Matrix_Local(matrixGraph, srg));
40 return(mat);
41}
42
43//----------------------------------------------------------------------------
44const char*
45Matrix_Local::typeName()
46{ return( "fei::Matrix_Local" ); }
47
48//----------------------------------------------------------------------------
49int
50Matrix_Local::parameters(const fei::ParameterSet& /*paramset*/)
51{
52 return(0);
53}
54
55//----------------------------------------------------------------------------
56int
57Matrix_Local::parameters(int /*numParams*/, const char* const* /*paramStrings*/)
58{
59 return(0);
60}
61
63Matrix_Local::getMatrixGraph() const
64{ return( matrixGraph_ ); }
65
66void
67Matrix_Local::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
68{ matrixGraph_ = matrixGraph; }
69
70int
71Matrix_Local::getGlobalNumRows() const
72{ return( sparseRowGraph_->rowNumbers.size() ); }
73
74int
75Matrix_Local::getLocalNumRows() const
76{ return( getGlobalNumRows() ); }
77
78int
79Matrix_Local::getRowIndex(int rowNumber) const
80{
81 int* rows = &(sparseRowGraph_->rowNumbers[0]);
82 int numRows = getLocalNumRows();
83 return( fei::binarySearch(rowNumber, rows, numRows) );
84}
85
86int
87Matrix_Local::getRowLength(int row, int& length) const
88{
89 int idx = getRowIndex(row);
90 if (idx < 0) return(idx);
91
92 length = sparseRowGraph_->rowOffsets[idx+1] -
93 sparseRowGraph_->rowOffsets[idx];
94 return(0);
95}
96
97int
98Matrix_Local::putScalar(double scalar)
99{
100 for(unsigned i=0; i<coefs_.size(); ++i) coefs_[i] = scalar;
101 stateChanged_ = true;
102 return(0);
103}
104
105int
106Matrix_Local::copyOutRow(int row, int len, double* coefs, int* indices) const
107{
108 int idx = getRowIndex(row);
109 if (idx < 0) return(idx);
110
111 int offset = sparseRowGraph_->rowOffsets[idx];
112 int length = sparseRowGraph_->rowOffsets[idx+1]-offset;
113 if (length > len) length = len;
114
115 for(int i=0; i<length; ++i) {
116 indices[i] = sparseRowGraph_->packedColumnIndices[offset+i];
117 coefs[i] = coefs_[offset+i];
118 }
119
120 return(0);
121}
122
123int
124Matrix_Local::giveToMatrix(int numRows, const int* rows,
125 int numCols, const int* cols,
126 const double* const* values,
127 bool sumInto,
128 int format)
129{
130 if (numRows == 0 || numCols == 0) {
131 return(0);
132 }
133
134 if (format != FEI_DENSE_ROW && format != FEI_DENSE_COL) {
135 return(-1);
136 }
137
138 const double** myvalues = const_cast<const double**>(values);
139 if (format != FEI_DENSE_ROW) {
140 fei::Matrix_core::copyTransposeToWorkArrays(numRows, numCols, values,
141 work_data1D_, work_data2D_);
142 myvalues = &work_data2D_[0];
143 }
144
145 for(int i=0; i<numRows; ++i) {
146 int idx = getRowIndex(rows[i]);
147 if (idx < 0) {
148 throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, row not found.");
149 }
150
151 int offset = sparseRowGraph_->rowOffsets[idx];
152 int len = sparseRowGraph_->rowOffsets[idx+1] - offset;
153
154 int* colInds = &(sparseRowGraph_->packedColumnIndices[offset]);
155 double* coefs = &(coefs_[offset]);
156
157 for(int j=0; j<numCols; ++j) {
158 int idx2 = fei::binarySearch(cols[j], colInds, len);
159 if (idx2 < 0) {
160 throw std::runtime_error("fei::Matrix_Local::sumIn ERROR, col not found.");
161 }
162
163 if (sumInto) {
164 coefs[idx2] += myvalues[i][j];
165 }
166 else {
167 coefs[idx2] = myvalues[i][j];
168 }
169 }
170 }
171
172 stateChanged_ = true;
173 return(0);
174}
175
176int
177Matrix_Local::sumIn(int numRows, const int* rows,
178 int numCols, const int* cols,
179 const double* const* values,
180 int format)
181{
182 return( giveToMatrix(numRows, rows, numCols, cols, values,
183 true, format) );
184}
185
186int
187Matrix_Local::copyIn(int numRows, const int* rows,
188 int numCols, const int* cols,
189 const double* const* values,
190 int format)
191{
192 return( giveToMatrix(numRows, rows, numCols, cols, values,
193 false, format) );
194}
195
196int
197Matrix_Local::sumInFieldData(int fieldID,
198 int idType,
199 int rowID,
200 int colID,
201 const double* const* data,
202 int format)
203{
204 fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
205 fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
206
207 int fieldSize = (int)rspace->getFieldSize(fieldID);
208 std::vector<int> indices(2*fieldSize);
209
210 rspace->getGlobalIndex(idType, rowID, fieldID, indices[0]);
211 for(int i=1; i<fieldSize; ++i) {
212 indices[i] = indices[0]+i;
213 }
214
215 cspace->getGlobalIndex(idType, colID, fieldID, indices[fieldSize]);
216 for(int i=1; i<fieldSize; ++i) {
217 indices[fieldSize+i] = indices[fieldSize]+i;
218 }
219
220 return( giveToMatrix(fieldSize, &indices[0], fieldSize, &indices[fieldSize],
221 data, true, format) );
222}
223
224int
225Matrix_Local::sumInFieldData(int fieldID,
226 int idType,
227 int rowID,
228 int colID,
229 const double* data,
230 int format)
231{
232 fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
233
234 int fieldSize = (int)rspace->getFieldSize(fieldID);
235 std::vector<const double*> data2D(fieldSize);
236
237 int offset = 0;
238 for(int i=0; i<fieldSize; ++i) {
239 data2D[i] = &data[offset];
240 offset += fieldSize;
241 }
242
243 return( sumInFieldData(fieldID, idType, rowID, colID,
244 &data2D[0], format) );
245}
246
247int
248Matrix_Local::sumIn(int blockID, int connectivityID,
249 const double* const* values,
250 int format)
251{
252 int numIndices = matrixGraph_->getConnectivityNumIndices(blockID);
253 std::vector<int> indices(numIndices);
254
255 matrixGraph_->getConnectivityIndices(blockID, connectivityID,
256 numIndices, &indices[0], numIndices);
257
258 return( giveToMatrix(numIndices, &indices[0], numIndices, &indices[0],
259 values, true, format) );
260}
261
262int
263Matrix_Local::globalAssemble()
264{ return(0); }
265
266int
267Matrix_Local::multiply(fei::Vector* x,
268 fei::Vector* y)
269{
270 FEI_COUT << "fei::Matrix_Local::multiply NOT IMPLEMENTED."<<FEI_ENDL;
271 return(-1);
272}
273
274void
275Matrix_Local::setCommSizes()
276{
277}
278
279int
280Matrix_Local::gatherFromOverlap(bool accumulate)
281{
282 (void)accumulate;
283 return(0);
284}
285
286int
287Matrix_Local::writeToFile(const char* filename,
288 bool matrixMarketFormat)
289{
290 fei::SharedPtr<fei::VectorSpace> vspace = matrixGraph_->getRowSpace();
291
292 MPI_Comm comm = vspace->getCommunicator();
293
294 int localProc = fei::localProc(comm);
295
296 FEI_OSTRINGSTREAM osstr;
297 osstr << filename << "." << localProc << ".mtx";
298 std::string fullname = osstr.str();
299
300 FEI_OFSTREAM ofstr(fullname.c_str(), IOS_OUT);
301
302 return( writeToStream(ofstr, matrixMarketFormat) );
303}
304
305int
306Matrix_Local::writeToStream(FEI_OSTREAM& ostrm,
307 bool matrixMarketFormat)
308{
309 static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
310
311 fei::SharedPtr<fei::VectorSpace> rspace = matrixGraph_->getRowSpace();
312 fei::SharedPtr<fei::VectorSpace> cspace = matrixGraph_->getColSpace();
313
314 int numRows = getLocalNumRows();
315 int numCols = cspace->getEqnNumbers().size();
316 int nnz = coefs_.size();
317
318 if (matrixMarketFormat) {
319 ostrm << mmbanner << FEI_ENDL;
320 ostrm << numRows << " " << numCols << " " << nnz << FEI_ENDL;
321 }
322 else {
323 ostrm << numRows << " " << numCols << " "<< FEI_ENDL;
324 }
325
326 std::vector<int>& rowNumbers = sparseRowGraph_->rowNumbers;
327 std::vector<int>& rowOffsets = sparseRowGraph_->rowOffsets;
328 std::vector<int>& colIndices = sparseRowGraph_->packedColumnIndices;
329
330 ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
331 ostrm.precision(13);
332
333 int offset = 0;
334 for(unsigned i=0; i<rowNumbers.size(); ++i) {
335 int rowlen = rowOffsets[i+1]-rowOffsets[i];
336
337 for(int j=0; j<rowlen; ++j) {
338 if (matrixMarketFormat) {
339 ostrm << rowNumbers[i]+1 << " " << colIndices[offset]+1
340 << " " << coefs_[offset] << FEI_ENDL;
341 }
342 else {
343 ostrm << rowNumbers[i] << " " << colIndices[offset]
344 << " " << coefs_[offset] << FEI_ENDL;
345 }
346 ++offset;
347 }
348 }
349
350 return(0);
351}
352
353bool
354Matrix_Local::usingBlockEntryStorage()
355{ return( false ); }
356
357void
358Matrix_Local::markState()
359{
360 stateChanged_ = false;
361}
362
363bool
364Matrix_Local::changedSinceMark()
365{ return(stateChanged_); }
366
367const std::vector<int>&
368Matrix_Local::getRowNumbers() const
369{ return( sparseRowGraph_->rowNumbers ); }
370
371const std::vector<int>&
372Matrix_Local::getRowOffsets() const
373{ return( sparseRowGraph_->rowOffsets ); }
374
375const std::vector<int>&
376Matrix_Local::getColumnIndices() const
377{ return( sparseRowGraph_->packedColumnIndices ); }
378
379const std::vector<double>&
380Matrix_Local::getCoefs() const
381{ return( coefs_ ); }
382
383}//namespace fei
384
int localProc(MPI_Comm comm)
int binarySearch(const T &item, const T *list, int len)
void writeToStream(snl_fei::RaggedTable< MAP_TYPE, SET_TYPE > &table, FEI_OSTREAM &os, const char *lineprefix=NULL)