9 #include "fei_fstream.hpp"
11 #include "fei_Vector_core.hpp"
12 #include "fei_VectorSpace.hpp"
13 #include "fei_CSVec.hpp"
14 #include "snl_fei_RecordCollection.hpp"
15 #include "fei_TemplateUtils.hpp"
16 #include "fei_impl_utils.hpp"
21 #define fei_file "fei_Vector_core.cpp"
23 #include "fei_ErrMacros.hpp"
29 comm_(vecSpace->getCommunicator()),
36 remotelyOwnedProcs_(),
43 sendRecvProcsNeedUpdated_(true),
44 overlapAlreadySet_(false),
47 eqnComm_.
reset(
new fei::EqnComm(comm_,numLocalEqns));
49 const std::vector<int>& offsets = eqnComm_->getGlobalOffsets();
52 numLocal_ = lastLocalOffset_ - firstLocalOffset_ + 1;
56 os<<dbgprefix_<<
" ctor firstLocal="<<firstLocalOffset_<<
", lastLocal="
57 <<lastLocalOffset_<<FEI_ENDL;
63 for(
unsigned i=0; i<remotelyOwned_.size(); ++i) {
64 delete remotelyOwned_[i];
69 const int* remoteEqns)
71 if (numRemoteEqns == 0 && remoteEqns == NULL) {
72 if (overlapAlreadySet_)
return;
75 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
76 FEI_OSTREAM& os = *output_stream_;
77 os << dbgprefix_<<
"setOverlap"<<FEI_ENDL;
82 if (numRemoteEqns != 0 && remoteEqns != NULL) {
83 for(
int i=0; i<numRemoteEqns; ++i) {
84 int proc = eqnComm_->getOwnerProc(remoteEqns[i]);
85 if (proc == local_proc)
continue;
86 fei::CSVec* remoteVec = getRemotelyOwned(proc);
87 fei::add_entry(*remoteVec, remoteEqns[i], 0.0);
91 std::vector<int> eqns;
92 vecSpace_->getIndices_SharedAndOwned(eqns);
94 for(
size_t i=0; i<eqns.size(); ++i) {
95 int proc = eqnComm_->getOwnerProc(eqns[i]);
96 if (proc == local_proc)
continue;
97 fei::CSVec* remoteVec = getRemotelyOwned(proc);
99 fei::add_entry(*remoteVec, eqns[i], 0.0);
103 overlapAlreadySet_ =
true;
104 sendRecvProcsNeedUpdated_ =
true;
114 if (!overlapAlreadySet_) {
121 std::vector<int> recvProcs;
122 for(
unsigned i=0; i<remotelyOwned_.size(); ++i) {
124 if (remotelyOwned_[i]->size() == 0)
continue;
126 recvProcs.push_back(remotelyOwnedProcs_[i]);
130 std::vector<int> sendProcs;
134 std::vector<std::vector<int> > send_ints(sendProcs.size());
135 std::vector<std::vector<double> > send_doubles(sendProcs.size());
136 std::vector<int> send_sizes(sendProcs.size());
138 std::vector<MPI_Request> mpiReqs(sendProcs.size()+recvProcs.size());
139 std::vector<MPI_Status> mpiStatuses(sendProcs.size()+recvProcs.size());
146 for(
unsigned i=0; i<sendProcs.size(); ++i) {
147 MPI_Irecv(&send_sizes[i], 1, MPI_INT, sendProcs[i],
148 tag1, comm_, &mpiReqs[i]);
154 for(
unsigned i=0; i<recvProcs.size(); ++i) {
155 int proc = recvProcs[i];
157 fei::CSVec* remoteVec = getRemotelyOwned(proc);
158 int size = remoteVec->size();
159 MPI_Send(&size, 1, MPI_INT, proc, tag1, comm_);
162 MPI_Waitall(sendProcs.size(), &mpiReqs[0], &mpiStatuses[0]);
166 for(
unsigned i=0; i<sendProcs.size(); ++i) {
167 int proc = sendProcs[i];
168 int size = send_sizes[i];
169 send_ints[i].resize(size);
170 MPI_Irecv(&(send_ints[i][0]), size, MPI_INT, proc, tag1,
172 send_doubles[i].resize(size);
176 for(
unsigned i=0; i<recvProcs.size(); ++i) {
177 int proc = recvProcs[i];
178 fei::CSVec* remoteVec = getRemotelyOwned(proc);
179 int size = remoteVec->size();
180 int* indices = &(remoteVec->indices())[0];
181 MPI_Send(indices, size, MPI_INT, proc, tag1, comm_);
184 MPI_Waitall(sendProcs.size(), &mpiReqs[0], &mpiStatuses[0]);
187 for(
unsigned i=0; i<recvProcs.size(); ++i) {
188 int proc = recvProcs[i];
189 fei::CSVec* remoteVec = getRemotelyOwned(proc);
190 int size = remoteVec->size();
191 double* coefs = &(remoteVec->coefs())[0];
192 MPI_Irecv(coefs, size, MPI_DOUBLE, proc, tag2, comm_, &mpiReqs[i]);
196 for(
unsigned i=0; i<sendProcs.size(); ++i) {
197 int proc = sendProcs[i];
199 int num = send_sizes[i];
200 int err = copyOutOfUnderlyingVector(num, &(send_ints[i][0]),
201 &(send_doubles[i][0]), 0);
203 FEI_COUT <<
"fei::Vector_core::scatterToOverlap ERROR getting data to send."<<FEI_ENDL;
207 MPI_Send(&(send_doubles[i][0]), num, MPI_DOUBLE, proc, tag2, comm_);
210 MPI_Waitall(recvProcs.size(), &mpiReqs[0], &mpiStatuses[0]);
212 #endif //#ifndef FEI_SER
220 int vectorIndex)
const
222 for(
int i=0; i<numValues; ++i) {
223 int ind = indices[i];
225 int local = ind - firstLocalOffset_;
226 if (local < 0 || local >= numLocal_) {
231 int proc = eqnComm_->getOwnerProc(ind);
232 const fei::CSVec* remoteVec = getRemotelyOwned(proc);
234 int insertPoint = -1;
238 <<
", index " << ind <<
" not in remotelyOwned_ vec object for proc "
243 values[i] = remoteVec->coefs()[idx];
247 CHK_ERR( copyOutOfUnderlyingVector(1, &ind, &(values[i]), vectorIndex) );
256 const double* values,
262 for(
int i=0; i<numValues; ++i) {
263 int ind = indices[i];
264 double val = values[i];
273 int local = ind - firstLocalOffset_;
274 if (local < 0 || local >= numLocal_) {
275 int proc = eqnComm_->getOwnerProc(ind);
276 if (output_level_ >= fei::BRIEF_LOGS && output_stream_ != NULL) {
277 FEI_OSTREAM& os = *output_stream_;
278 os << dbgprefix_<<
"giveToVector remote["<<proc<<
"]("
279 <<ind<<
","<<val<<
")"<<FEI_ENDL;
282 if (proc != prev_proc) {
283 remoteVec = getRemotelyOwned(proc);
284 prev_vec = remoteVec;
289 fei::add_entry( *remoteVec, ind, val);
292 fei::put_entry( *remoteVec, ind, val);
296 int err = giveToUnderlyingVector(1, &ind, &val, sumInto, vectorIndex);
299 <<
", val: " << val << FEI_ENDL;
316 if (vecSpace_.get() == NULL) ERReturn(-1);
318 int fieldSize = vecSpace_->getFieldSize(fieldID);
320 work_indices_.resize(numIDs*fieldSize);
321 int* indicesPtr = &work_indices_[0];
323 CHK_ERR( vecSpace_->getGlobalIndices(numIDs, IDs, idType, fieldID,
326 CHK_ERR( giveToVector(numIDs*fieldSize, indicesPtr, data, sumInto, vectorIndex) );
331 int fei::Vector_core::assembleFieldDataLocalIDs(
int fieldID,
339 if (vecSpace_.get() == NULL) ERReturn(-1);
341 int fieldSize = vecSpace_->getFieldSize(fieldID);
343 work_indices_.resize(numIDs*fieldSize);
344 int* indicesPtr = &work_indices_[0];
346 CHK_ERR( vecSpace_->getGlobalIndicesLocalIDs(numIDs, localIDs, idType, fieldID,
349 CHK_ERR( giveToVector(numIDs*fieldSize, indicesPtr, data, sumInto, vectorIndex) );
354 void fei::Vector_core::pack_send_buffers(
const std::vector<int>& sendProcs,
355 const std::vector<fei::CSVec*>& remotelyOwned,
356 std::vector<std::vector<char> >& send_chars,
358 bool zeroRemotelyOwnedAfterPacking)
360 for(
size_t i=0; i<sendProcs.size(); ++i) {
361 int proc = sendProcs[i];
362 fei::CSVec* remoteVec = getRemotelyOwned(proc);
363 fei::impl_utils::pack_indices_coefs(remoteVec->indices(),
364 remoteVec->coefs(), send_chars[i], resize_buffer);
366 if (zeroRemotelyOwnedAfterPacking) {
367 fei::set_values(*remoteVec, 0.0);
372 void fei::Vector_core::setCommSizes()
377 for(
unsigned i=0; i<remotelyOwned_.size(); ++i) {
379 if (remotelyOwned_[i]->size() == 0)
continue;
381 sendProcs_.push_back(remotelyOwnedProcs_[i]);
384 std::vector<int> tmpSendProcs;
385 vecSpace_->getSendProcs(tmpSendProcs);
386 for(
size_t i=0; i<tmpSendProcs.size(); ++i) {
388 for(
size_t j=0; j<sendProcs_.size(); ++j) {
389 if (sendProcs_[j] == tmpSendProcs[i]) {
393 if (sendProcs_[j] > tmpSendProcs[i]) {
394 sendProcs_.insert(sendProcs_.begin()+j, tmpSendProcs[i]);
399 if (!found) sendProcs_.push_back(tmpSendProcs[i]);
405 std::vector<MPI_Request> mpiReqs;
408 send_chars_.resize(sendProcs_.size());
409 recv_chars_.resize(recvProcs_.size());
411 bool resize_buffer =
true;
412 bool zero_remotely_owned_after_packing =
false;
413 pack_send_buffers(sendProcs_, remotelyOwned_, send_chars_,
414 resize_buffer, zero_remotely_owned_after_packing);
416 recv_sizes_.resize(recvProcs_.size());
417 mpiReqs.resize(recvProcs_.size());
420 for(
size_t i=0; i<recvProcs_.size(); ++i) {
421 int proc = recvProcs_[i];
422 MPI_Irecv(&recv_sizes_[i], 1, MPI_INT, proc,
423 tag1, comm_, &mpiReqs[i]);
427 for(
unsigned i=0; i<sendProcs_.size(); ++i) {
428 int proc = sendProcs_[i];
429 int size = send_chars_[i].size();
430 MPI_Send(&size, 1, MPI_INT, proc, tag1, comm_);
433 for(
size_t i=0; i<recvProcs_.size(); ++i) {
436 MPI_Waitany(mpiReqs.size(), &mpiReqs[0], &index, &status);
438 recv_chars_[index].resize(recv_sizes_[index]);
441 sendRecvProcsNeedUpdated_ =
false;
452 std::vector<MPI_Request> mpiReqs;
455 if (sendRecvProcsNeedUpdated_) {
459 mpiReqs.resize(recvProcs_.size());
462 for(
size_t i=0; i<recvProcs_.size(); ++i) {
463 MPI_Irecv(&(recv_chars_[i][0]), recv_sizes_[i], MPI_CHAR, recvProcs_[i],
464 tag1, comm_, &mpiReqs[i]);
467 bool resize_buffer =
false;
468 bool zero_remotely_owned_after_packing =
true;
469 pack_send_buffers(sendProcs_, remotelyOwned_, send_chars_,
470 resize_buffer, zero_remotely_owned_after_packing);
473 for(
size_t i=0; i<sendProcs_.size(); ++i) {
474 int proc = sendProcs_[i];
476 int size = send_chars_[i].size();
477 MPI_Send(&(send_chars_[i][0]), size, MPI_CHAR, proc, tag1, comm_);
480 int numRecvProcs = recvProcs_.size();
481 for(
size_t i=0; i<recvProcs_.size(); ++i) {
484 MPI_Waitany(numRecvProcs, &mpiReqs[0], &index, &status);
487 std::vector<int> indices;
488 std::vector<double> coefs;
490 for(
size_t i=0; i<recvProcs_.size(); ++i) {
491 fei::impl_utils::unpack_indices_coefs(recv_chars_[i], indices, coefs);
492 int num = indices.size();
493 if (num == 0)
continue;
494 int err = giveToUnderlyingVector(num, &(indices[0]),
495 &(coefs[0]), accumulate, 0);
502 #endif //#ifndef FEI_SER
514 if (vecSpace_.get() == NULL) ERReturn(-1);
516 int fieldSize = vecSpace_->getFieldSize(fieldID);
520 CHK_ERR( vecSpace_->getRecordCollection(idType, collection) );
524 std::vector<int>& eqnNums = vecSpace_->getEqnNumbers();
525 int* vspcEqnPtr = eqnNums.size() > 0 ? &eqnNums[0] : NULL;
528 for(
int i=0; i<numIDs; ++i) {
541 dofOffset = eqnNumbers[foffset] - eqnNumbers[0];
542 for(
int j=0; j<fieldSize; ++j) {
543 CHK_ERR( copyOut_FE(nodeNumber, dofOffset+j, data[offset++]));
548 work_indices_.resize(numIDs*fieldSize*2);
549 int* indicesPtr = &work_indices_[0];
551 CHK_ERR( vecSpace_->getGlobalIndices(numIDs, IDs, idType,
552 fieldID, indicesPtr) );
554 CHK_ERR( copyOut(numIDs*fieldSize, indicesPtr, data) );
561 bool matrixMarketFormat)
568 static char mmbanner[] =
"%%MatrixMarket matrix array real general";
572 if (p != localProc)
continue;
574 FEI_OFSTREAM* outFile = NULL;
576 outFile =
new FEI_OFSTREAM(filename, IOS_OUT);
577 FEI_OFSTREAM& ofref = *outFile;
578 if (matrixMarketFormat) {
579 ofref << mmbanner << FEI_ENDL;
580 ofref << eqnComm_->getGlobalOffsets()[
numProcs] <<
" 1" << FEI_ENDL;
583 ofref << eqnComm_->getGlobalOffsets()[
numProcs] << FEI_ENDL;
586 else outFile =
new FEI_OFSTREAM(filename, IOS_APP);
587 FEI_OFSTREAM& ofref = *outFile;
588 ofref.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
591 for(
int i=firstLocalOffset_; i<=lastLocalOffset_; ++i) {
592 CHK_ERR( copyOut(1, &i, &coef) );
593 if (matrixMarketFormat) {
594 ofref <<
" " << coef << FEI_ENDL;
597 ofref << i <<
" " << coef << FEI_ENDL;
608 bool matrixMarketFormat)
615 static char mmbanner[] =
"%%MatrixMarket matrix array real general";
617 IOS_FMTFLAGS oldf = ostrm.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
620 for(
int proc=0; proc<
numProcs; ++proc) {
622 if (proc != local_proc)
continue;
625 if (matrixMarketFormat) {
626 ostrm << mmbanner << FEI_ENDL;
627 ostrm << eqnComm_->getGlobalOffsets()[
numProcs] <<
" 1" << FEI_ENDL;
630 ostrm << eqnComm_->getGlobalOffsets()[
numProcs] << FEI_ENDL;
634 for(
size_t p=0; p<remotelyOwned_.size(); ++p) {
635 if (remotelyOwnedProcs_[p] > local_proc)
continue;
636 for(
size_t ii=0; ii<remotelyOwned_[p]->size(); ++ii) {
637 if (matrixMarketFormat) {
638 ostrm <<
" " << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
641 ostrm <<
" " << remotelyOwned_[p]->indices()[ii] <<
" "
642 << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
647 for(
int i=firstLocalOffset_; i<=lastLocalOffset_; ++i) {
648 CHK_ERR( copyOut(1, &i, &coef) );
649 if (matrixMarketFormat) {
650 ostrm <<
" " << coef << FEI_ENDL;
653 ostrm <<
" " << i <<
" " << coef << FEI_ENDL;
657 for(
size_t p=0; p<remotelyOwned_.size(); ++p) {
658 if (remotelyOwnedProcs_[p] < local_proc)
continue;
659 for(
size_t ii=0; ii<remotelyOwned_[p]->size(); ++ii) {
660 if (matrixMarketFormat) {
661 ostrm <<
" " << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
664 ostrm <<
" " << remotelyOwned_[p]->indices()[ii] <<
" "
665 << remotelyOwned_[p]->coefs()[ii] << FEI_ENDL;
671 ostrm.setf(oldf, IOS_FLOATFIELD);
GlobalIDType getNumber() const
void setOverlap(int numRemoteEqns=0, const int *remoteEqns=NULL)
int giveToVector(int numValues, const int *indices, const double *values, bool sumInto=true, int vectorIndex=0)
virtual int writeToStream(FEI_OSTREAM &ostrm, bool matrixMarketFormat=true)
int copyOut(int numValues, const int *indices, double *values, int vectorIndex=0) const
int binarySearch(const T &item, const T *list, int len)
virtual int copyOutFieldData(int fieldID, int idType, int numIDs, const int *IDs, double *data, int vectorIndex=0)
int mirrorProcs(MPI_Comm comm, std::vector< int > &toProcs, std::vector< int > &fromProcs)
OutputLevel output_level_
std::ostream & console_out()
FEI_OSTREAM * output_stream_
int localProc(MPI_Comm comm)
Vector_core(fei::SharedPtr< fei::VectorSpace > vecSpace, int numLocalEqns)
int getFieldEqnOffset(int fieldID, int &offset) const
virtual int writeToFile(const char *filename, bool matrixMarketFormat=true)
virtual int scatterToOverlap()
int getOffsetIntoEqnNumbers() const
fei::FieldMask * getFieldMask()
fei::Record< int > * getRecordWithID(int ID)
virtual int gatherFromOverlap(bool accumulate=true)
int numProcs(MPI_Comm comm)
int assembleFieldData(int fieldID, int idType, int numIDs, const int *IDs, const double *data, bool sumInto=true, int vectorIndex=0)