43#include <fei_macros.hpp>
45#include <fei_Include_Trilinos.hpp>
48#include <fei_VectorTraits_Epetra.hpp>
49#include <fei_MatrixTraits_Epetra.hpp>
50#include <fei_LinProbMgr_EpetraBasic.hpp>
53#include <fei_Trilinos_Helpers.hpp>
54#include <fei_ParameterSet.hpp>
55#include <fei_Vector_Impl.hpp>
56#include <fei_VectorReducer.hpp>
57#include <fei_Matrix_Impl.hpp>
58#include <fei_MatrixReducer.hpp>
61namespace Trilinos_Helpers {
66create_Epetra_Map(MPI_Comm comm,
67 const std::vector<int>& local_eqns)
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");
95 MPI_Comm comm = vecspace->getCommunicator();
101 int localSizeBlk = vecspace->getNumBlkIndices_Owned();
102 int globalSizeBlk = vecspace->getGlobalNumBlkIndices();
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]);
112 int errcode = vecspace->getBlkIndices_Owned(localSizeBlk,
113 blkEqnsPtr, blkEqnsPtr+localSizeBlk,
116 FEI_OSTRINGSTREAM osstr;
117 osstr <<
"create_Epetra_BlockMap ERROR, nonzero errcode="<<errcode
118 <<
" returned by vecspace->getBlkIndices_Owned.";
119 throw std::runtime_error(osstr.str());
123 blkEqnsPtr, blkEqnsPtr+localSizeBlk, 0, EComm);
130 bool blockEntries,
bool orderRowsWithLocalColsFirst)
133 matgraph->createGraph(blockEntries);
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;
148 create_Epetra_BlockMap(vecspace) : create_Epetra_Map(comm, local_eqns);
150 if (orderRowsWithLocalColsFirst ==
true &&
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];
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)
235 int localSize = vecSpace->getNumIndices_Owned();
241 Trilinos_Helpers::create_Epetra_CrsGraph(matrixGraph, blockEntryMatrix,
242 orderRowsWithLocalColsFirst);
244 if (blockEntryMatrix) {
248 bool zeroSharedRows =
false;
250 matrixGraph, localSize, zeroSharedRows));
251 zero_Epetra_VbrMatrix(epetraMatrix.get());
258 matrixGraph, localSize));
261 catch(std::runtime_error& exc) {
263 <<
"caught exception: '" << exc.what() <<
"', rethrowing..."
268 if (reducer.
get() != NULL) {
269 feimat.
reset(
new fei::MatrixReducer(reducer, tmpmat));
280 bool blockEntryMatrix,
286 matrixGraph->createGraph(blockEntryMatrix);
287 if (srgraph.
get() == NULL) {
288 throw std::runtime_error(
"create_LPM_EpetraBasic ERROR in fei::MatrixGraph::createGraph");
293 if (reducer.
get() != NULL) {
294 std::vector<int>& reduced_eqns = reducer->getLocalReducedEqns();
295 lpm_epetrabasic->setRowDistribution(reduced_eqns);
296 lpm_epetrabasic->setMatrixGraph(sharedsrg);
297 localSize = reduced_eqns.size();
302 std::vector<int> indices;
303 if (blockEntryMatrix) {
304 localSize = vecSpace->getNumBlkIndices_Owned();
305 indices.resize(localSize*2);
306 err = vecSpace->getBlkIndices_Owned(localSize, &indices[0],
307 &indices[localSize], chkNum);
310 localSize = vecSpace->getNumIndices_Owned();
311 indices.resize(localSize);
312 err = vecSpace->getIndices_Owned(indices);
315 throw std::runtime_error(
"Factory_Trilinos: createMatrix: error in vecSpace->getIndices_Owned");
318 lpm_epetrabasic->setRowDistribution(indices);
319 lpm_epetrabasic->setMatrixGraph(sharedsrg);
326 if (reducer.
get() != NULL) {
327 feimat.
reset(
new fei::MatrixReducer(reducer, tmpmat));
338 Teuchos::ParameterList& paramlist)
341 iter = paramset.
begin(),
342 iter_end = paramset.
end();
344 for(; iter != iter_end; ++iter) {
348 case fei::Param::STRING:
351 case fei::Param::DOUBLE:
354 case fei::Param::INT:
357 case fei::Param::BOOL:
366void copy_parameterlist(
const Teuchos::ParameterList& paramlist,
369 Teuchos::ParameterList::ConstIterator
370 iter = paramlist.begin(),
371 iter_end = paramlist.end();
373 for(; iter != iter_end; ++iter) {
374 const Teuchos::ParameterEntry& param = paramlist.entry(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
393get_Epetra_MultiVector(
fei::Vector* feivec,
bool soln_vec)
396 fei::VectorReducer* feireducer =
dynamic_cast<fei::VectorReducer*
>(feivec);
397 if (feireducer != NULL) vecptr = feireducer->getTargetVector().
get();
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) {
412 LinProbMgr_EpetraBasic* lpm_epetrabasic =
414 if (lpm_epetrabasic == 0) {
415 throw std::runtime_error(
"fei Trilinos_Helpers: ERROR getting LinProbMgr_EpetraBasic");
419 return(lpm_epetrabasic->get_solution_vector().
get());
422 return(lpm_epetrabasic->get_rhs_vector().
get());
429 fei::MatrixReducer* feireducer =
430 dynamic_cast<fei::MatrixReducer*
>(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());
451 fei::MatrixReducer* feireducer =
452 dynamic_cast<fei::MatrixReducer*
>(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) {
467 return(fei_epetra_crs->
getMatrix().get());
470 LinProbMgr_EpetraBasic* lpm_epetrabasic =
471 dynamic_cast<LinProbMgr_EpetraBasic*
>(fei_lpm->
getMatrix().get());
472 if (lpm_epetrabasic == 0) {
473 throw std::runtime_error(
"fei Trilinos_Helpers ERROR getting LinProbMgr_EpetraBasic");
476 return(lpm_epetrabasic->get_A_matrix().
get());
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")) {
497 crsA = get_Epetra_CrsMatrix(feiA.
get());
510 std::vector<double> zeros(maxBlkRowSize*maxBlkColSize, 0);
513 for(
int i=0; i<numMyRows; ++i) {
516 int* colindicesView = NULL;
517 int localrow = rowmap.
LID(row);
527 for(
int j=0; j<rowlength; ++j) {
528 int blkColSize = colmap.
ElementSize(colindicesView[j]);
530 blkRowSize, blkColSize);
int MyGlobalElements(int *MyGlobalElementList) const
const Epetra_Comm & Comm() const
int NumMyElements() const
bool MyGID(int GID_in) const
virtual int NumProc() const=0
virtual int MyPID() const=0
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
const Epetra_BlockMap & RowMap() const
const Epetra_BlockMap & ColMap() const
int SumAll(double *PartialSums, double *GlobalSums, int Count) const
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
int GlobalMaxRowDim() const
const Epetra_CrsGraph & Graph() const
int BeginReplaceMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
int GlobalMaxColDim() const
const Epetra_Map & RowMatrixRowMap() const
fei::SharedPtr< T > getMatrix()
const std::string & getStringValue() const
double getDoubleValue() const
ParamType getType() const
bool getBoolValue() const
const std::string & getName() const
void add(const Param ¶m, bool maintain_unique_keys=true)
const_iterator end() const
const_iterator begin() const
T * getUnderlyingVector()
int localProc(MPI_Comm comm)
std::ostream & console_out()