47 #ifdef HAVE_FEI_EPETRA 63 #ifdef HAVE_FEI_EPETRA 67 const std::vector<int>& local_eqns)
70 Epetra_MpiComm EComm(comm);
72 Epetra_SerialComm EComm;
75 int localSize = local_eqns.size();
77 EComm.SumAll(&localSize, &globalSize, 1);
79 if (localSize < 0 || globalSize < 0) {
80 throw std::runtime_error(
"Trilinos_Helpers::create_Epetra_Map: negative local or global size.");
83 Epetra_Map EMap(globalSize, localSize, 0, EComm);
90 if (vecspace.
get() == 0) {
91 throw std::runtime_error(
"create_Epetra_Map needs non-null fei::VectorSpace");
96 Epetra_MpiComm EComm(comm);
98 Epetra_SerialComm EComm;
104 if (localSizeBlk < 0 || globalSizeBlk < 0) {
105 throw std::runtime_error(
"Trilinos_Helpers::create_Epetra_BlockMap: fei::VectorSpace has negative local or global size.");
108 std::vector<int> blkEqns(localSizeBlk*2);
109 int* blkEqnsPtr = &(blkEqns[0]);
113 blkEqnsPtr, blkEqnsPtr+localSizeBlk,
117 osstr <<
"create_Epetra_BlockMap ERROR, nonzero errcode="<<errcode
118 <<
" returned by vecspace->getBlkIndices_Owned.";
119 throw std::runtime_error(osstr.str());
122 Epetra_BlockMap EBMap(globalSizeBlk, localSizeBlk,
123 blkEqnsPtr, blkEqnsPtr+localSizeBlk, 0, EComm);
130 bool blockEntries,
bool orderRowsWithLocalColsFirst)
134 if (localSRGraph.
get() == NULL) {
135 throw std::runtime_error(
"create_Epetra_CrsGraph ERROR in fei::MatrixGraph::createGraph");
138 int numLocallyOwnedRows = localSRGraph->
rowNumbers.size();
139 int* rowNumbers = numLocallyOwnedRows>0 ? &(localSRGraph->
rowNumbers[0]) : NULL;
140 int* rowOffsets = &(localSRGraph->
rowOffsets[0]);
141 int* packedColumnIndices = numLocallyOwnedRows>0 ? &(localSRGraph->
packedColumnIndices[0]) : NULL;
144 MPI_Comm comm = vecspace->getCommunicator();
145 std::vector<int>& local_eqns = localSRGraph->
rowNumbers;
147 Epetra_BlockMap emap = blockEntries ?
148 create_Epetra_BlockMap(vecspace) : create_Epetra_Map(comm, local_eqns);
150 if (orderRowsWithLocalColsFirst ==
true &&
151 emap.Comm().NumProc() > 2 && !blockEntries) {
152 bool* used_row =
new bool[local_eqns.size()];
153 for(
unsigned ii=0; ii<local_eqns.size(); ++ii) used_row[ii] =
false;
156 std::vector<int> ordered_local_eqns(local_eqns.size());
157 for(
unsigned ii=0; ii<local_eqns.size(); ++ii) {
158 bool row_has_off_proc_cols =
false;
159 for(
int j=rowOffsets[ii]; j<rowOffsets[ii+1]; ++j) {
160 if (emap.MyGID(packedColumnIndices[j]) ==
false) {
161 row_has_off_proc_cols =
true;
166 if (row_has_off_proc_cols ==
false) {
167 ordered_local_eqns[offset++] = rowNumbers[ii];
172 for(
unsigned ii=0; ii<local_eqns.size(); ++ii) {
173 if (used_row[ii] ==
true)
continue;
174 ordered_local_eqns[offset++] = rowNumbers[ii];
177 emap = create_Epetra_Map(comm, ordered_local_eqns);
183 std::vector<int> rowLengths; rowLengths.reserve(numLocallyOwnedRows);
184 for(
int ii=0; ii<numLocallyOwnedRows; ++ii) {
185 rowLengths.push_back(rowOffsets[ii+1]-rowOffsets[ii]);
188 bool staticProfile =
true;
189 int* rowLengthsPtr = rowLengths.empty() ? NULL : &rowLengths[0];
190 Epetra_CrsGraph egraph(
Copy, emap, rowLengthsPtr, staticProfile);
192 const Epetra_Comm& ecomm = emap.Comm();
195 int firstLocalEqn = numLocallyOwnedRows > 0 ? rowNumbers[0] : -1;
198 for(
int i=0; i<numLocallyOwnedRows; ++i) {
199 int err = egraph.InsertGlobalIndices(firstLocalEqn+i,
201 &(packedColumnIndices[offset]));
204 <<
" inserting row " << firstLocalEqn+i<<
", cols ";
205 for(
int ii=0; ii<rowLengths[i]; ++ii) {
209 throw std::runtime_error(
"... occurred in create_Epetra_CrsGraph");
212 offset += rowLengths[i];
217 egraph.FillComplete();
221 if (!blockEntries) egraph.OptimizeStorage();
228 bool blockEntryMatrix,
230 bool orderRowsWithLocalColsFirst)
240 Epetra_CrsGraph egraph =
241 Trilinos_Helpers::create_Epetra_CrsGraph(matrixGraph, blockEntryMatrix,
242 orderRowsWithLocalColsFirst);
244 if (blockEntryMatrix) {
246 epetraMatrix(
new Epetra_VbrMatrix(
Copy, egraph));
248 bool zeroSharedRows =
false;
250 matrixGraph, localSize, zeroSharedRows));
251 zero_Epetra_VbrMatrix(epetraMatrix.
get());
255 epetraMatrix(
new Epetra_CrsMatrix(
Copy, egraph));
258 matrixGraph, localSize));
261 catch(std::runtime_error& exc) {
263 <<
"caught exception: '" << exc.what() <<
"', rethrowing..." 268 if (reducer.
get() != NULL) {
280 bool blockEntryMatrix,
287 if (srgraph.
get() == NULL) {
288 throw std::runtime_error(
"create_LPM_EpetraBasic ERROR in fei::MatrixGraph::createGraph");
293 if (reducer.
get() != NULL) {
297 localSize = reduced_eqns.size();
302 std::vector<int> indices;
303 if (blockEntryMatrix) {
305 indices.resize(localSize*2);
307 &indices[localSize], chkNum);
311 indices.resize(localSize);
315 throw std::runtime_error(
"Factory_Trilinos: createMatrix: error in vecSpace->getIndices_Owned");
326 if (reducer.
get() != NULL) {
335 #endif //HAVE_FEI_EPETRA 341 iter = paramset.
begin(),
342 iter_end = paramset.
end();
344 for(; iter != iter_end; ++iter) {
370 iter = paramlist.
begin(),
371 iter_end = paramlist.
end();
373 for(; iter != iter_end; ++iter) {
375 if (param.
isType<std::string>()) {
376 paramset.
add(
fei::Param(paramlist.
name(iter).c_str(), Teuchos::getValue<std::string>(param).c_str()));
378 else if (param.
isType<
double>()) {
379 paramset.
add(
fei::Param(paramlist.
name(iter).c_str(), Teuchos::getValue<double>(param)));
381 else if (param.
isType<
int>()) {
382 paramset.
add(
fei::Param(paramlist.
name(iter).c_str(), Teuchos::getValue<int>(param)));
384 else if (param.
isType<
bool>()) {
385 paramset.
add(
fei::Param(paramlist.
name(iter).c_str(), Teuchos::getValue<bool>(param)));
390 #ifdef HAVE_FEI_EPETRA 393 get_Epetra_MultiVector(
fei::Vector* feivec,
bool soln_vec)
404 if (fei_epetra_vec == NULL && fei_lpm == NULL) {
405 throw std::runtime_error(
"failed to obtain Epetra_MultiVector from fei::Vector.");
408 if (fei_epetra_vec != NULL) {
414 if (lpm_epetrabasic == 0) {
415 throw std::runtime_error(
"fei Trilinos_Helpers: ERROR getting LinProbMgr_EpetraBasic");
425 Epetra_VbrMatrix* get_Epetra_VbrMatrix(
fei::Matrix* feimat)
432 if (feireducer != NULL) {
438 if (fei_epetra_vbr == NULL) {
439 throw std::runtime_error(
"failed to obtain Epetra_VbrMatrix from fei::Matrix.");
442 return(fei_epetra_vbr->
getMatrix().get());
445 Epetra_CrsMatrix* get_Epetra_CrsMatrix(
fei::Matrix* feimat)
454 if (feireducer != NULL) {
462 if (fei_epetra_crs == NULL && fei_lpm == NULL) {
463 throw std::runtime_error(
"failed to obtain Epetra_CrsMatrix from fei::Matrix.");
466 if (fei_epetra_crs != NULL) {
472 if (lpm_epetrabasic == 0) {
473 throw std::runtime_error(
"fei Trilinos_Helpers ERROR getting LinProbMgr_EpetraBasic");
483 Epetra_CrsMatrix*& crsA,
484 Epetra_Operator*& opA,
485 Epetra_MultiVector*& x,
486 Epetra_MultiVector*& b)
488 x = get_Epetra_MultiVector(feix.
get(),
true);
489 b = get_Epetra_MultiVector(feib.
get(),
false);
491 const char* matname = feiA->
typeName();
492 if (!strcmp(matname,
"Epetra_VbrMatrix")) {
493 Epetra_VbrMatrix*
A = get_Epetra_VbrMatrix(feiA.
get());
497 crsA = get_Epetra_CrsMatrix(feiA.
get());
502 int zero_Epetra_VbrMatrix(Epetra_VbrMatrix* mat)
504 const Epetra_CrsGraph& crsgraph = mat->Graph();
505 const Epetra_BlockMap& rowmap = crsgraph.RowMap();
506 const Epetra_BlockMap& colmap = crsgraph.ColMap();
507 mat->RowMatrixRowMap();
508 int maxBlkRowSize = mat->GlobalMaxRowDim();
509 int maxBlkColSize = mat->GlobalMaxColDim();
510 std::vector<double> zeros(maxBlkRowSize*maxBlkColSize, 0);
511 int numMyRows = rowmap.NumMyElements();
512 int* myRows = rowmap.MyGlobalElements();
513 for(
int i=0; i<numMyRows; ++i) {
516 int* colindicesView = NULL;
517 int localrow = rowmap.LID(row);
518 int err = crsgraph.ExtractMyRowView(localrow, rowlength, colindicesView);
522 err = mat->BeginReplaceMyValues(localrow, rowlength, colindicesView);
526 int blkRowSize = rowmap.ElementSize(localrow);
527 for(
int j=0; j<rowlength; ++j) {
528 int blkColSize = colmap.ElementSize(colindicesView[j]);
529 err = mat->SubmitBlockEntry(&zeros[0], maxBlkRowSize,
530 blkRowSize, blkColSize);
535 err = mat->EndSubmitEntries();
544 #endif //HAVE_FEI_EPETRA
virtual const char * typeName()=0
int getBlkIndices_Owned(int lenBlkIndices, int *globalBlkIndices, int *blkSizes, int &numBlkIndices)
ParamType getType() const
bool getBoolValue() const
ConstIterator begin() const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
const ParameterEntry & entry(ConstIterator i) const
fei::SharedPtr< Epetra_CrsMatrix > get_A_matrix()
std::vector< int > rowNumbers
const_iterator end() const
int getNumIndices_Owned() const
const std::string & getStringValue() const
virtual fei::SharedPtr< fei::SparseRowGraph > createGraph(bool blockEntryGraph, bool localRowGraph_includeSharedRows=false)=0
const_iterator begin() const
std::vector< int > packedColumnIndices
std::vector< int > rowOffsets
fei::SharedPtr< fei::Vector > getTargetVector()
ConstIterator end() const
virtual fei::SharedPtr< fei::VectorSpace > getRowSpace()=0
void copy_parameterlist(const Teuchos::ParameterList ¶mlist, fei::ParameterSet ¶mset)
void add(const Param ¶m, bool maintain_unique_keys=true)
int getGlobalNumBlkIndices() const
fei::SharedPtr< Epetra_MultiVector > get_rhs_vector()
fei::SharedPtr< T > getMatrix()
const std::string & getName() const
void copy_parameterset(const fei::ParameterSet ¶mset, Teuchos::ParameterList ¶mlist)
MPI_Comm getCommunicator() const
int getNumBlkIndices_Owned() const
const std::string & name() const
std::ostream & console_out()
int getIndices_Owned(std::vector< int > &globalIndices) const
T * getUnderlyingVector()
std::vector< int > & getLocalReducedEqns()
int localProc(MPI_Comm comm)
virtual void setRowDistribution(const std::vector< int > &ownedGlobalRows)=0
double getDoubleValue() const
#define FEI_OSTRINGSTREAM
fei::SharedPtr< Epetra_MultiVector > get_solution_vector()
virtual void setMatrixGraph(fei::SharedPtr< fei::SparseRowGraph > matrixGraph)=0