FEI Version of the Day
Loading...
Searching...
No Matches
fei_MatrixReducer.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_MatrixReducer.hpp"
11#include "fei_EqnComm.hpp"
12#include "fei_Matrix_core.hpp"
13#include "fei_sstream.hpp"
14#include "fei_fstream.hpp"
15#include <fei_CommUtils.hpp>
16
17namespace fei {
18
19MatrixReducer::MatrixReducer(fei::SharedPtr<fei::Reducer> reducer,
21 : reducer_(reducer),
22 target_(target),
23 globalAssembleCalled_(false),
24 changedSinceMark_(false)
25{
27 target->getMatrixGraph()->getRowSpace();
28 MPI_Comm comm = vspace->getCommunicator();
29 int numLocal = reducer_->getLocalReducedEqns().size();
30 fei::SharedPtr<fei::EqnComm> eqnComm(new fei::EqnComm(comm, numLocal));
31 fei::Matrix_core* target_core =
32 dynamic_cast<fei::Matrix_core*>(target_.get());
33 if (target_core == NULL) {
34 throw std::runtime_error("fei::MatrixReducer ERROR, target matrix not dynamic_cast-able to fei::Matrix_core.");
35 }
36
37 target_core->setEqnComm(eqnComm);
38}
39
40MatrixReducer::~MatrixReducer()
41{
42}
43
44int
45MatrixReducer::parameters(const fei::ParameterSet& paramset)
46{
47 return(target_->parameters(paramset));
48}
49
50void
51MatrixReducer::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
52{
53 target_->setMatrixGraph(matrixGraph);
54}
55
56int
57MatrixReducer::getGlobalNumRows() const
58{
59 return(target_->getMatrixGraph()->getRowSpace()->getGlobalNumIndices());
60}
61
62int
63MatrixReducer::getLocalNumRows() const
64{
65 return(target_->getMatrixGraph()->getRowSpace()->getNumIndices_Owned());
66}
67
68int MatrixReducer::putScalar(double scalar)
69{ return(target_->putScalar(scalar)); }
70
71int
72MatrixReducer::getRowLength(int row, int& length) const
73{
74 if (reducer_->isSlaveEqn(row)) {
75 FEI_OSTRINGSTREAM osstr;
76 osstr << "fei::MatrixReducer::getRowLength ERROR, row="<<row<<" is a slave eqn. You can't get a slave row from the reduced matrix.";
77 throw std::runtime_error(osstr.str());
78 }
79
80 int reducedrow = reducer_->translateToReducedEqn(row);
81 return(target_->getRowLength(reducedrow, length));
82}
83
84int
85MatrixReducer::copyOutRow(int row, int len, double* coefs, int* indices) const
86{
87 if (reducer_->isSlaveEqn(row)) {
88 FEI_OSTRINGSTREAM osstr;
89 osstr << "fei::MatrixReducer::copyOutRow ERROR, requested row ("<<row
90 <<") is a slave eqn. You can't get a slave row from the reduced matrix.";
91 throw std::runtime_error(osstr.str());
92 }
93
94 int reducedrow = reducer_->translateToReducedEqn(row);
95 int err = target_->copyOutRow(reducedrow, len, coefs, indices);
96 for(int i=0; i<len; ++i) {
97 indices[i] = reducer_->translateFromReducedEqn(indices[i]);
98 }
99 return(err);
100}
101
102int
103MatrixReducer::sumIn(int numRows, const int* rows,
104 int numCols, const int* cols,
105 const double* const* values,
106 int format)
107{
108 int err = reducer_->addMatrixValues(numRows, rows, numCols, cols,
109 values, true, *target_, format);
110 return(err);
111}
112
113int
114MatrixReducer::copyIn(int numRows, const int* rows,
115 int numCols, const int* cols,
116 const double* const* values,
117 int format)
118{
119 int err = reducer_->addMatrixValues(numRows, rows, numCols, cols,
120 values, false, *target_, format);
121 return(err);
122}
123
124int
125MatrixReducer::sumInFieldData(int fieldID,
126 int idType,
127 int rowID,
128 int colID,
129 const double* const* data,
130 int format)
131{
133 target_->getMatrixGraph()->getRowSpace();
135 target_->getMatrixGraph()->getColSpace();
136
137 unsigned fieldSize = rowSpace->getFieldSize(fieldID);
138 std::vector<int> indices(fieldSize*2);
139 int* rowIndices = &indices[0];
140 int* colIndices = rowIndices+fieldSize;
141
142 rowSpace->getGlobalIndices(1, &rowID, idType, fieldID, rowIndices);
143 colSpace->getGlobalIndices(1, &colID, idType, fieldID, colIndices);
144
145 if (format != FEI_DENSE_ROW) {
146 throw std::runtime_error("MatrixReducer: bad format");
147 }
148
149 int err = reducer_->addMatrixValues(fieldSize, rowIndices,
150 fieldSize, colIndices,
151 data, true, *target_, format);
152 return(err);
153}
154
155int
156MatrixReducer::sumInFieldData(int fieldID,
157 int idType,
158 int rowID,
159 int colID,
160 const double* data,
161 int format)
162{
164 target_->getMatrixGraph()->getRowSpace();
165
166 unsigned fieldSize = rowSpace->getFieldSize(fieldID);
167
168 std::vector<const double*> data_2d(fieldSize);
169 for(unsigned i=0; i<fieldSize; ++i) {
170 data_2d[i] = &data[i*fieldSize];
171 }
172
173 return(sumInFieldData(fieldID, idType, rowID, colID, &data_2d[0], format));
174}
175
176int
177MatrixReducer::sumIn(int blockID, int connectivityID,
178 const double* const* values,
179 int format)
180{
181 fei::SharedPtr<fei::MatrixGraph> matGraph = getMatrixGraph();
182 int numRowIndices, numColIndices, dummy;
183 matGraph->getConnectivityNumIndices(blockID, numRowIndices, numColIndices);
184
185 std::vector<int> indices(numRowIndices+numColIndices);
186 int* rowIndices = &indices[0];
187 int* colIndices = rowIndices+numRowIndices;
188
189 matGraph->getConnectivityIndices(blockID, connectivityID,
190 numRowIndices, rowIndices, dummy,
191 numColIndices, colIndices, dummy);
192
193 return(sumIn(numRowIndices, rowIndices, numColIndices, colIndices,
194 values, format));
195}
196
197int
198MatrixReducer::globalAssemble()
199{
200 reducer_->assembleReducedMatrix(*target_);
201 return(target_->globalAssemble());
202}
203
204int MatrixReducer::multiply(fei::Vector* x, fei::Vector* y)
205{ return(target_->multiply(x, y)); }
206
207int MatrixReducer::gatherFromOverlap(bool accumulate)
208{
209 reducer_->assembleReducedMatrix(*target_);
210 target_->setCommSizes();
211 return(target_->gatherFromOverlap(accumulate));
212}
213
214int MatrixReducer::writeToFile(const char* filename,
215 bool matrixMarketFormat)
216{
217 static char mmbanner[] = "%%MatrixMarket matrix coordinate real general";
218 std::vector<int>& localrows = reducer_->getLocalReducedEqns();
219 int localNumRows = localrows.size();
220
221 int globalNNZ = 0;
222 int localNNZ = 0;
223
224 for(int i=0; i<localNumRows; ++i) {
225 int len;
226 CHK_ERR( target_->getRowLength(localrows[i], len) );
227 localNNZ += len;
228 }
229
230 MPI_Comm comm = getMatrixGraph()->getRowSpace()->getCommunicator();
231
232 CHK_MPI( fei::GlobalSum(comm, localNNZ, globalNNZ) );
233 int globalNumRows = 0;
234 CHK_MPI( fei::GlobalSum(comm, localNumRows, globalNumRows) );
235
236 int globalNumCols = globalNumRows;
237
238 for(int p=0; p<fei::numProcs(comm); ++p) {
239 fei::Barrier(comm);
240 if (p != fei::localProc(comm)) continue;
241
242 FEI_OFSTREAM* outFile = NULL;
243 if (p==0) {
244 outFile = new FEI_OFSTREAM(filename, IOS_OUT);
245 FEI_OFSTREAM& ofs = *outFile;
246 if (matrixMarketFormat) {
247 ofs << mmbanner << FEI_ENDL;
248 ofs <<globalNumRows<< " " <<globalNumCols<< " " <<globalNNZ<<FEI_ENDL;
249 }
250 else {
251 ofs <<globalNumRows<< " " <<globalNumCols<<FEI_ENDL;
252 }
253 }
254 else outFile = new FEI_OFSTREAM(filename, IOS_APP);
255
256 outFile->setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
257 outFile->precision(13);
258 FEI_OFSTREAM& ofs = *outFile;
259
260 int rowLength;
261 std::vector<int> work_indices;
262 std::vector<double> work_data1D;
263
264 for(int i=0; i<localNumRows; ++i) {
265 int row = localrows[i];
266 CHK_ERR( target_->getRowLength(row, rowLength) );
267
268 work_indices.resize(rowLength);
269 work_data1D.resize(rowLength);
270
271 int* indPtr = &work_indices[0];
272 double* coefPtr = &work_data1D[0];
273
274 CHK_ERR( target_->copyOutRow(row, rowLength, coefPtr, indPtr) );
275
276 for(int j=0; j<rowLength; ++j) {
277 if (matrixMarketFormat) {
278 ofs << row+1 <<" "<<indPtr[j]+1<<" "<<coefPtr[j]<<FEI_ENDL;
279 }
280 else {
281 ofs << row <<" "<<indPtr[j]<<" "<<coefPtr[j]<<FEI_ENDL;
282 }
283 }
284 }
285
286 delete outFile;
287 }
288
289 return(0);
290}
291
292int MatrixReducer::writeToStream(FEI_OSTREAM& ostrm,
293 bool matrixMarketFormat)
294{
295 return(target_->writeToStream(ostrm, matrixMarketFormat));
296}
297
298void MatrixReducer::markState()
299{ target_->markState(); }
300
301bool MatrixReducer::changedSinceMark()
302{ return(target_->changedSinceMark()); }
303
304}//namespace fei
305
int localProc(MPI_Comm comm)
int numProcs(MPI_Comm comm)
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)