FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
fei_MatrixReducer.cpp
Go to the documentation of this file.
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
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
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
43
44int
46{
47 return(target_->parameters(paramset));
48}
49
50void
55
56int
61
62int
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)) {
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)) {
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
126 int idType,
127 int rowID,
128 int colID,
129 const double* const* data,
130 int format)
131{
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
157 int idType,
158 int rowID,
159 int colID,
160 const double* data,
161 int format)
162{
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{
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
203
206
208{
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
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
293 bool matrixMarketFormat)
294{
295 return(target_->writeToStream(ostrm, matrixMarketFormat));
296}
297
300
303
304}//namespace fei
305
virtual fei::SharedPtr< fei::VectorSpace > getColSpace()=0
virtual fei::SharedPtr< fei::VectorSpace > getRowSpace()=0
int writeToStream(FEI_OSTREAM &ostrm, bool matrixMarketFormat=true)
int putScalar(double scalar)
int getRowLength(int row, int &length) const
MatrixReducer(fei::SharedPtr< fei::Reducer > reducer, fei::SharedPtr< fei::Matrix > target)
int copyIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)
void setMatrixGraph(fei::SharedPtr< fei::MatrixGraph > matrixGraph)
int multiply(fei::Vector *x, fei::Vector *y)
int gatherFromOverlap(bool accumulate=true)
int writeToFile(const char *filename, bool matrixMarketFormat=true)
int sumIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)
int copyOutRow(int row, int len, double *coefs, int *indices) const
int sumInFieldData(int fieldID, int idType, int rowID, int colID, const double *const *data, int format=0)
fei::SharedPtr< fei::Reducer > reducer_
fei::SharedPtr< fei::MatrixGraph > getMatrixGraph() const
fei::SharedPtr< fei::Matrix > target_
int parameters(const fei::ParameterSet &paramset)
void setEqnComm(fei::SharedPtr< fei::EqnComm > eqnComm)
virtual int copyOutRow(int row, int len, double *coefs, int *indices) const =0
virtual int writeToStream(FEI_OSTREAM &ostrm, bool matrixMarketFormat=true)=0
virtual fei::SharedPtr< fei::MatrixGraph > getMatrixGraph() const =0
virtual void setCommSizes()=0
virtual int gatherFromOverlap(bool accumulate=true)=0
virtual bool changedSinceMark()=0
virtual int parameters(const fei::ParameterSet &paramset)=0
virtual void setMatrixGraph(fei::SharedPtr< fei::MatrixGraph > matrixGraph)=0
virtual int putScalar(double scalar)=0
virtual int getRowLength(int row, int &length) const =0
virtual int globalAssemble()=0
virtual int multiply(fei::Vector *x, fei::Vector *y)=0
virtual void markState()=0
int translateToReducedEqn(int unreducedEqn) const
bool isSlaveEqn(int unreducedEqn) const
int addMatrixValues(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, bool sum_into, fei::Matrix &feimat, int format)
std::vector< int > & getLocalReducedEqns()
void assembleReducedMatrix(fei::Matrix &matrix)
int translateFromReducedEqn(int reduced_eqn) const
int getNumIndices_Owned() const
int getGlobalNumIndices() const
MPI_Comm getCommunicator() const
#define CHK_ERR(a)
#define CHK_MPI(a)
#define FEI_DENSE_ROW
Definition fei_defs.h:77
#define FEI_OFSTREAM
#define FEI_OSTREAM
#define IOS_APP
#define FEI_ENDL
#define IOS_FLOATFIELD
#define IOS_SCIENTIFIC
#define IOS_OUT
#define MPI_Comm
Definition fei_mpi.h:56
#define FEI_OSTRINGSTREAM
int localProc(MPI_Comm comm)
int numProcs(MPI_Comm comm)
int GlobalSum(MPI_Comm comm, std::vector< T > &local, std::vector< T > &global)
void Barrier(MPI_Comm comm)