FEI Version of the Day
Loading...
Searching...
No Matches
fei_Matrix_core.cpp
1/*--------------------------------------------------------------------*/
2/* Copyright 2005 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#include <fei_macros.hpp>
9
10#include <fei_utils.hpp>
11
12#include <fei_ParameterSet.hpp>
13#include <fei_LogManager.hpp>
14#include <fei_Matrix_core.hpp>
15#include <fei_CSRMat.hpp>
16#include <fei_TemplateUtils.hpp>
17#include <fei_CommUtils.hpp>
18#include <fei_impl_utils.hpp>
19
20#include <fei_VectorSpace.hpp>
21#include <snl_fei_PointBlockMap.hpp>
22#include <fei_MatrixGraph.hpp>
23#include <snl_fei_Utils.hpp>
24
25#include <fei_MatrixTraits.hpp>
26#include <fei_MatrixTraits_FillableMat.hpp>
27
28#undef fei_file
29#define fei_file "fei_Matrix_core.cpp"
30#include <fei_ErrMacros.hpp>
31
32fei::Matrix_core::Matrix_core(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
33 int numLocalEqns)
34 : name_(),
35 work_indices_(),
36 work_indices2_(),
37 work_ints_(),
38 work_data1D_(),
39 work_data2D_(),
40 eqnComm_(),
41 rhsVector_(),
42 comm_(matrixGraph->getRowSpace()->getCommunicator()),
43 localProc_(0),
44 numProcs_(1),
45 vecSpace_(),
46 matrixGraph_(matrixGraph),
47 remotelyOwned_(),
48 remotelyOwned_last_requested_(NULL),
49 sendProcs_(),
50 recvProcs_(),
51 recv_chars_(),
52 send_chars_(),
53 sendRecvProcsNeedUpdated_(true),
54 proc_last_requested_(-1),
55 haveBlockMatrix_(false),
56 haveFEMatrix_(false),
57 globalOffsets_(),
58 firstLocalOffset_(0),
59 lastLocalOffset_(0)
60{
61 if (matrixGraph.get() == NULL) {
62 throw std::runtime_error("fei::Matrix_core constructed with NULL fei::MatrixGraph");
63 }
64
65 vecSpace_ = matrixGraph->getRowSpace();
66
67 if (vecSpace_->initCompleteAlreadyCalled()) {
68 vecSpace_->getGlobalIndexOffsets(globalOffsets_);
69 eqnComm_.reset(new fei::EqnComm(comm_, numLocalEqns, globalOffsets_));
70 }
71 else {
72 eqnComm_.reset(new fei::EqnComm(comm_, numLocalEqns));
73 }
74
75 localProc_ = fei::localProc(comm_);
76 numProcs_ = fei::numProcs(comm_);
77
78 setName("dbg");
79
80 globalOffsets_ = eqnComm_->getGlobalOffsets();
81
82 firstLocalOffset_ = globalOffsets_[localProc_];
83 lastLocalOffset_ = globalOffsets_[localProc_+1]-1;
84}
85
86std::map<int,fei::FillableMat*>&
87fei::Matrix_core::getRemotelyOwnedMatrices()
88{
89 return remotelyOwned_;
90}
91
92void
93fei::Matrix_core::putScalar_remotelyOwned(double scalar)
94{
95 std::map<int,FillableMat*>::iterator
96 it = remotelyOwned_.begin(),
97 it_end = remotelyOwned_.end();
98 for(; it!=it_end; ++it) {
100 }
101}
102
103void
104fei::Matrix_core::setEqnComm(fei::SharedPtr<fei::EqnComm> eqnComm)
105{
106 eqnComm_ = eqnComm;
107 globalOffsets_ = eqnComm_->getGlobalOffsets();
108 firstLocalOffset_ = globalOffsets_[localProc_];
109 lastLocalOffset_ = globalOffsets_[localProc_+1]-1;
110 sendRecvProcsNeedUpdated_ = true;
111}
112
113fei::Matrix_core::~Matrix_core()
114{
115 std::map<int,FillableMat*>::iterator
116 it = remotelyOwned_.begin(),
117 it_end = remotelyOwned_.end();
118 for(; it!=it_end; ++it) {
119 delete it->second;
120 }
121}
122
123void
124fei::Matrix_core::parameters(const fei::ParameterSet& paramset)
125{
126 const fei::Param* param = paramset.get("name");
127 fei::Param::ParamType ptype = param != NULL ?
128 param->getType() : fei::Param::BAD_TYPE;
129 if (ptype == fei::Param::STRING) {
130 setName(param->getStringValue().c_str());
131 }
132
133 param = paramset.get("FEI_OUTPUT_LEVEL");
134 ptype = param != NULL ? param->getType() : fei::Param::BAD_TYPE;
135 if (ptype == fei::Param::STRING) {
137 log_manager.setOutputLevel(param->getStringValue().c_str());
138 setOutputLevel(fei::utils::string_to_output_level(param->getStringValue()));
139 }
140
141}
142
143void fei::Matrix_core::setName(const char* name)
144{
145 if (name == NULL) return;
146
147 if (name_ == name) return;
148
149 name_ = name;
150}
151
152int fei::Matrix_core::getOwnerProc(int globalEqn) const
153{
154 int len = globalOffsets_.size();
155 if (globalEqn > globalOffsets_[len-1]) return(-1);
156
157 for(int p=len-2; p>=0; --p) {
158 if (globalEqn >= globalOffsets_[p]) {
159 return(p);
160 }
161 }
162
163 return(-1);
164}
165
166void fei::Matrix_core::setRHS(fei::SharedPtr<fei::Vector> rhsvector)
167{
168 rhsVector_ = rhsvector;
169}
170
171void fei::Matrix_core::setCommSizes()
172{
173#ifndef FEI_SER
174 //iterate the remotelyOwned_ map and create a list of processors
175 //that we will be sending data to. (processors which own matrix rows
176 //that we share.)
177 sendProcs_.clear();
178 std::map<int,FillableMat*>::iterator
179 it = remotelyOwned_.begin(),
180 it_end = remotelyOwned_.end();
181 for(; it!=it_end; ++it) {
182 if (it->first == localProc_) continue;
183 if (it->second != NULL) {
184 if (it->second->getNumRows() == 0) {
185 continue;
186 }
187
188 sendProcs_.push_back(it->first);
189 }
190 }
191
192 //now we can find out which procs we'll be receiving from.
193 recvProcs_.clear();
194 fei::mirrorProcs(comm_, sendProcs_, recvProcs_);
195
196 int tag1 = 11111;
197 std::vector<MPI_Request> mpiReqs(recvProcs_.size());
198
199 //next find out the size of each message that will be recvd or sent.
200 std::vector<int> recv_sizes(recvProcs_.size(), 0);
201 for(size_t i=0; i<recvProcs_.size(); ++i) {
202 MPI_Irecv(&recv_sizes[i], 1, MPI_INT, recvProcs_[i],
203 tag1, comm_, &mpiReqs[i]);
204 }
205
206 //pack our to-be-sent data into buffers, and send the
207 //sizes to the receiving procs:
208 send_chars_.resize(sendProcs_.size());
209 recv_chars_.resize(recvProcs_.size());
210
211 for(size_t i=0; i<sendProcs_.size(); ++i) {
212 int proc = sendProcs_[i];
213
214 int num_bytes = fei::impl_utils::num_bytes_FillableMat(*(remotelyOwned_[proc]));
215 send_chars_[i].resize(num_bytes);
216 char* buffer = &(send_chars_[i][0]);
217 fei::impl_utils::pack_FillableMat(*(remotelyOwned_[proc]), buffer);
218
219 int bsize = send_chars_[i].size();
220
221 MPI_Send(&bsize, 1, MPI_INT, proc, tag1, comm_);
222 }
223
224 int numRecvProcs = recvProcs_.size();
225 for(size_t i=0; i<recvProcs_.size(); ++i) {
226 MPI_Status status;
227 int index;
228 MPI_Waitany(numRecvProcs, &mpiReqs[0], &index, &status);
229
230 recv_chars_[index].resize(recv_sizes[index]);
231 }
232
233 sendRecvProcsNeedUpdated_ = false;
234#endif
235}
236
237int fei::Matrix_core::gatherFromOverlap(bool accumulate)
238{
239 if (numProcs() == 1) return(0);
240
241#ifndef FEI_SER
242 //this function gathers shared-but-not-owned data onto the
243 //owning processors.
244
245 if (sendRecvProcsNeedUpdated_) {
246 setCommSizes();
247 }
248
249 std::vector<MPI_Request> mpiReqs(recvProcs_.size());
250
251 int tag1 = 11111;
252 int numRecvProcs = recvProcs_.size();
253
254 //post the recvs.
255 for(size_t i=0; i<recvProcs_.size(); ++i) {
256 int bsize = recv_chars_[i].size();
257
258 MPI_Irecv(&(recv_chars_[i][0]), bsize, MPI_CHAR, recvProcs_[i],
259 tag1, comm_, &mpiReqs[i]);
260 }
261
262 //now pack and send our buffers.
263 for(size_t i=0; i<sendProcs_.size(); ++i) {
264 int proc = sendProcs_[i];
265 fei::FillableMat* remoteMat = remotelyOwned_[proc];
266 char* buffer = &(send_chars_[i][0]);
267 fei::impl_utils::pack_FillableMat(*remoteMat, buffer);
268 remoteMat->setValues(0.0);
269
270 MPI_Send(&(send_chars_[i][0]), send_chars_[i].size(), MPI_CHAR, proc, tag1, comm_);
271 }
272
273 for(size_t ir=0; ir<recvProcs_.size(); ++ir) {
274 int index;
275 MPI_Status status;
276 MPI_Waitany(numRecvProcs, &mpiReqs[0], &index, &status);
277 }
278
279 //and finally, unpack and store the received buffers.
280 CSRMat recvMat;
281 for(size_t ir=0; ir<recvProcs_.size(); ++ir) {
282 size_t len = recv_chars_[ir].size();
283 const char* data = &(recv_chars_[ir][0]);
284 bool all_zeros = fei::impl_utils::unpack_CSRMat(data, data+len, recvMat);
285
286 if (all_zeros) continue;
287
288 int nrows = recvMat.getNumRows();
289
290 for(int i=0; i<nrows; ++i) {
291 int row = recvMat.getGraph().rowNumbers[i];
292 int rowOffset = recvMat.getGraph().rowOffsets[i];
293 int rowLen = recvMat.getGraph().rowOffsets[i+1] - rowOffset;
294
295 int* indices = &(recvMat.getGraph().packedColumnIndices[rowOffset]);
296 double* vals = &(recvMat.getPackedCoefs()[rowOffset]);
297 int err = 0;
298 if (haveBlockMatrix()) {
299 err = giveToBlockMatrix(1, &row, rowLen, indices, &vals, accumulate);
300 }
301 else {
302 err = giveToUnderlyingMatrix(1, &row, rowLen, indices,
303 &vals, accumulate, FEI_DENSE_ROW);
304 }
305 if (err != 0) {
306 FEI_COUT << "fei::Matrix_core::gatherFromOverlap ERROR storing "
307 << "values for row=" << row<<", recv'd from proc " << recvProcs_[i]
308 << FEI_ENDL;
309
310 return(err);
311 }
312 }
313 }
314#endif
315
316 return(0);
317}
318
319void
320fei::Matrix_core::copyTransposeToWorkArrays(int numRows, int numCols,
321 const double*const* values,
322 std::vector<double>& work_1D,
323 std::vector<const double*>& work_2D)
324{
325 //copy the transpose of 'values' into work-arrays.
326
327 int arrayLen = numCols*numRows;
328 if ((int)work_1D.size() != arrayLen) {
329 work_1D.resize(arrayLen);
330 }
331 if ((int)work_2D.size() != numCols) {
332 work_2D.resize(numCols);
333 }
334
335 const double** dataPtr = &work_2D[0];
336 double* data1DPtr = &work_1D[0];
337
338 for(int i=0; i<numCols; ++i) {
339 dataPtr[i] = data1DPtr;
340
341 for(int j=0; j<numRows; ++j) {
342 data1DPtr[j] = values[j][i];
343 }
344
345 data1DPtr += numRows;
346 }
347}
348
349
350int fei::Matrix_core::convertPtToBlk(int numRows,
351 const int* rows,
352 int numCols,
353 const int* cols,
354 int* blkRows,
355 int* blkRowOffsets,
356 int* blkCols,
357 int* blkColOffsets)
358{
359 snl_fei::PointBlockMap* pointBlockMap = vecSpace_->getPointBlockMap();
360
361 int i;
362 for(i=0; i<numRows; ++i) {
363 if (pointBlockMap->getPtEqnInfo(rows[i], blkRows[i], blkRowOffsets[i])!=0){
364 ERReturn(-1);
365 }
366 }
367 for(i=0; i<numCols; ++i) {
368 if(pointBlockMap->getPtEqnInfo(cols[i], blkCols[i], blkColOffsets[i])!=0){
369 ERReturn(-1);
370 }
371 }
372
373 return(0);
374}
375
376int fei::Matrix_core::copyPointRowsToBlockRow(int numPtRows,
377 int numPtCols,
378 const double*const* ptValues,
379 int numBlkCols,
380 const int* blkColDims,
381 double** blkValues)
382{
383 int ptColOffset = 0;
384 for(int blki=0; blki<numBlkCols; ++blki) {
385 int blkvalueOffset = 0;
386 double* blkvalues_i = blkValues[blki];
387 for(int j=0; j<blkColDims[blki]; ++j) {
388 int loc = ptColOffset+j;
389 for(int i=0; i<numPtRows; ++i) {
390 blkvalues_i[blkvalueOffset++] = ptValues[i][loc];
391 }
392 }
393 ptColOffset += blkColDims[blki];
394 }
395
396 return(0);
397}
398
399void fei::Matrix_core::setMatrixGraph(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
400{
401 matrixGraph_ = matrixGraph;
402 vecSpace_ = matrixGraph->getRowSpace();
403}
404
void setOutputLevel(OutputLevel olevel)
static LogManager & getLogManager()
const std::string & getStringValue() const
ParamType getType() const
Definition fei_Param.hpp:98
const Param * get(const char *name) const
void reset(T *p=0)
int getPtEqnInfo(int ptEqn, int &blkEqn, int &blkOffset)
bool unpack_CSRMat(const char *buffer_begin, const char *buffer_end, fei::CSRMat &mat)
fei::OutputLevel string_to_output_level(const std::string &str)
Definition fei_utils.cpp:58
int localProc(MPI_Comm comm)
int mirrorProcs(MPI_Comm comm, std::vector< int > &toProcs, std::vector< int > &fromProcs)
int numProcs(MPI_Comm comm)
static int setValues(T *mat, double scalar)