FEI Version of the Day
Loading...
Searching...
No Matches
fei_Trilinos_Helpers.cpp
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43#include <fei_macros.hpp>
44
45#include <fei_Include_Trilinos.hpp>
46
47#ifdef HAVE_FEI_EPETRA
48#include <fei_VectorTraits_Epetra.hpp>
49#include <fei_MatrixTraits_Epetra.hpp>
50#include <fei_LinProbMgr_EpetraBasic.hpp>
51#endif
52
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>
59//#include <EpetraExt_BlockMapOut.h>
60
61namespace Trilinos_Helpers {
62
63#ifdef HAVE_FEI_EPETRA
64
66create_Epetra_Map(MPI_Comm comm,
67 const std::vector<int>& local_eqns)
68{
69#ifndef FEI_SER
70 Epetra_MpiComm EComm(comm);
71#else
73#endif
74
75 int localSize = local_eqns.size();
76 int globalSize = 0;
77 EComm.SumAll(&localSize, &globalSize, 1);
78
79 if (localSize < 0 || globalSize < 0) {
80 throw std::runtime_error("Trilinos_Helpers::create_Epetra_Map: negative local or global size.");
81 }
82
83 Epetra_Map EMap(globalSize, localSize, 0, EComm);
84 return(EMap);
85}
86
88create_Epetra_BlockMap(const fei::SharedPtr<fei::VectorSpace>& vecspace)
89{
90 if (vecspace.get() == 0) {
91 throw std::runtime_error("create_Epetra_Map needs non-null fei::VectorSpace");
92 }
93
94#ifndef FEI_SER
95 MPI_Comm comm = vecspace->getCommunicator();
96 Epetra_MpiComm EComm(comm);
97#else
99#endif
100
101 int localSizeBlk = vecspace->getNumBlkIndices_Owned();
102 int globalSizeBlk = vecspace->getGlobalNumBlkIndices();
103
104 if (localSizeBlk < 0 || globalSizeBlk < 0) {
105 throw std::runtime_error("Trilinos_Helpers::create_Epetra_BlockMap: fei::VectorSpace has negative local or global size.");
106 }
107
108 std::vector<int> blkEqns(localSizeBlk*2);
109 int* blkEqnsPtr = &(blkEqns[0]);
110
111 int chkNum = 0;
112 int errcode = vecspace->getBlkIndices_Owned(localSizeBlk,
113 blkEqnsPtr, blkEqnsPtr+localSizeBlk,
114 chkNum);
115 if (errcode != 0) {
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());
120 }
121
122 Epetra_BlockMap EBMap(globalSizeBlk, localSizeBlk,
123 blkEqnsPtr, blkEqnsPtr+localSizeBlk, 0, EComm);
124
125 return(EBMap);
126}
127
129create_Epetra_CrsGraph(const fei::SharedPtr<fei::MatrixGraph>& matgraph,
130 bool blockEntries, bool orderRowsWithLocalColsFirst)
131{
133 matgraph->createGraph(blockEntries);
134 if (localSRGraph.get() == NULL) {
135 throw std::runtime_error("create_Epetra_CrsGraph ERROR in fei::MatrixGraph::createGraph");
136 }
137
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;
142
143 fei::SharedPtr<fei::VectorSpace> vecspace = matgraph->getRowSpace();
144 MPI_Comm comm = vecspace->getCommunicator();
145 std::vector<int>& local_eqns = localSRGraph->rowNumbers;
146
147 Epetra_BlockMap emap = blockEntries ?
148 create_Epetra_BlockMap(vecspace) : create_Epetra_Map(comm, local_eqns);
149
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;
154
155 int offset = 0;
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;
162 break;
163 }
164 }
165
166 if (row_has_off_proc_cols == false) {
167 ordered_local_eqns[offset++] = rowNumbers[ii];
168 used_row[ii] = true;
169 }
170 }
171
172 for(unsigned ii=0; ii<local_eqns.size(); ++ii) {
173 if (used_row[ii] == true) continue;
174 ordered_local_eqns[offset++] = rowNumbers[ii];
175 }
176
177 emap = create_Epetra_Map(comm, ordered_local_eqns);
178 delete [] used_row;
179 }
180
181// EpetraExt::BlockMapToMatrixMarketFile("EBMap.np12.mm",emap,"AriaTest");
182
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]);
186 }
187
188 bool staticProfile = true;
189 int* rowLengthsPtr = rowLengths.empty() ? NULL : &rowLengths[0];
190 Epetra_CrsGraph egraph(Copy, emap, rowLengthsPtr, staticProfile);
191
192 const Epetra_Comm& ecomm = emap.Comm();
193 int localProc = ecomm.MyPID();
194
195 int firstLocalEqn = numLocallyOwnedRows > 0 ? rowNumbers[0] : -1;
196
197 int offset = 0;
198 for(int i=0; i<numLocallyOwnedRows; ++i) {
199 int err = egraph.InsertGlobalIndices(firstLocalEqn+i,
200 rowLengths[i],
201 &(packedColumnIndices[offset]));
202 if (err != 0) {
203 fei::console_out() << "proc " << localProc << " err-return " << err
204 << " inserting row " << firstLocalEqn+i<<", cols ";
205 for(int ii=0; ii<rowLengths[i]; ++ii) {
206 fei::console_out() << packedColumnIndices[offset+ii]<<",";
207 }
208 fei::console_out() << FEI_ENDL;
209 throw std::runtime_error("... occurred in create_Epetra_CrsGraph");
210 }
211
212 offset += rowLengths[i];
213 }
214
215 //Epetra_BlockMap* domainmap = const_cast<Epetra_BlockMap*>(&(epetraGraph_->DomainMap()));
216 //Epetra_BlockMap* rangemap = const_cast<Epetra_BlockMap*>(&(epetraGraph_->RangeMap()));
217 egraph.FillComplete();
218
219 //only optimize-storage for graph if we're not in block-matrix mode.
220 //(Epetra_VbrMatrix can't be constructed with an optimize-storage'd graph.)
221 if (!blockEntries) egraph.OptimizeStorage();
222
223 return(egraph);
224}
225
227create_from_Epetra_Matrix(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
228 bool blockEntryMatrix,
230 bool orderRowsWithLocalColsFirst)
231{
232 fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
233 //localSize is in point-equations, because we will only use it for constructing
234 //the fei::Matrix_Impl, and those always deal in point-equations.
235 int localSize = vecSpace->getNumIndices_Owned();
236
239 try {
240 Epetra_CrsGraph egraph =
241 Trilinos_Helpers::create_Epetra_CrsGraph(matrixGraph, blockEntryMatrix,
242 orderRowsWithLocalColsFirst);
243
244 if (blockEntryMatrix) {
246 epetraMatrix(new Epetra_VbrMatrix(Copy, egraph));
247
248 bool zeroSharedRows = false;
249 tmpmat.reset(new fei::Matrix_Impl<Epetra_VbrMatrix>(epetraMatrix,
250 matrixGraph, localSize, zeroSharedRows));
251 zero_Epetra_VbrMatrix(epetraMatrix.get());
252 }
253 else {
255 epetraMatrix(new Epetra_CrsMatrix(Copy, egraph));
256
257 tmpmat.reset(new fei::Matrix_Impl<Epetra_CrsMatrix>(epetraMatrix,
258 matrixGraph, localSize));
259 }
260 }
261 catch(std::runtime_error& exc) {
262 fei::console_out() << "Trilinos_Helpers::create_from_Epetra_Matrix ERROR, "
263 << "caught exception: '" << exc.what() << "', rethrowing..."
264 << FEI_ENDL;
265 throw exc;
266 }
267
268 if (reducer.get() != NULL) {
269 feimat.reset(new fei::MatrixReducer(reducer, tmpmat));
270 }
271 else {
272 feimat = tmpmat;
273 }
274
275 return(feimat);
276}
277
279create_from_LPM_EpetraBasic(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
280 bool blockEntryMatrix,
283 lpm_epetrabasic)
284{
286 matrixGraph->createGraph(blockEntryMatrix);
287 if (srgraph.get() == NULL) {
288 throw std::runtime_error("create_LPM_EpetraBasic ERROR in fei::MatrixGraph::createGraph");
289 }
290
291 fei::SharedPtr<fei::SparseRowGraph> sharedsrg(srgraph);
292 int localSize;
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();
298 }
299 else {
300 fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
301 int err = 0,chkNum;
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);
308 }
309 else {
310 localSize = vecSpace->getNumIndices_Owned();
311 indices.resize(localSize);
312 err = vecSpace->getIndices_Owned(indices);
313 }
314 if (err != 0) {
315 throw std::runtime_error("Factory_Trilinos: createMatrix: error in vecSpace->getIndices_Owned");
316 }
317
318 lpm_epetrabasic->setRowDistribution(indices);
319 lpm_epetrabasic->setMatrixGraph(sharedsrg);
320 }
321
322 fei::SharedPtr<fei::Matrix> tmpmat(new fei::Matrix_Impl<fei::LinearProblemManager>(lpm_epetrabasic, matrixGraph, localSize));
323
325
326 if (reducer.get() != NULL) {
327 feimat.reset(new fei::MatrixReducer(reducer, tmpmat));
328 }
329 else {
330 feimat = tmpmat;
331 }
332
333 return(feimat);
334}
335#endif //HAVE_FEI_EPETRA
336
337void copy_parameterset(const fei::ParameterSet& paramset,
338 Teuchos::ParameterList& paramlist)
339{
341 iter = paramset.begin(),
342 iter_end = paramset.end();
343
344 for(; iter != iter_end; ++iter) {
345 const fei::Param& param = *iter;
346 fei::Param::ParamType ptype = param.getType();
347 switch(ptype) {
348 case fei::Param::STRING:
349 paramlist.set(param.getName(), param.getStringValue());
350 break;
351 case fei::Param::DOUBLE:
352 paramlist.set(param.getName(), param.getDoubleValue());
353 break;
354 case fei::Param::INT:
355 paramlist.set(param.getName(), param.getIntValue());
356 break;
357 case fei::Param::BOOL:
358 paramlist.set(param.getName(), param.getBoolValue());
359 break;
360 default:
361 break;
362 }
363 }
364}
365
366void copy_parameterlist(const Teuchos::ParameterList& paramlist,
367 fei::ParameterSet& paramset)
368{
369 Teuchos::ParameterList::ConstIterator
370 iter = paramlist.begin(),
371 iter_end = paramlist.end();
372
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()));
377 }
378 else if (param.isType<double>()) {
379 paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<double>(param)));
380 }
381 else if (param.isType<int>()) {
382 paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<int>(param)));
383 }
384 else if (param.isType<bool>()) {
385 paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<bool>(param)));
386 }
387 }
388}
389
390#ifdef HAVE_FEI_EPETRA
391
393get_Epetra_MultiVector(fei::Vector* feivec, bool soln_vec)
394{
395 fei::Vector* vecptr = feivec;
396 fei::VectorReducer* feireducer = dynamic_cast<fei::VectorReducer*>(feivec);
397 if (feireducer != NULL) vecptr = feireducer->getTargetVector().get();
398
400 dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>*>(vecptr);
402 dynamic_cast<fei::Vector_Impl<fei::LinearProblemManager>*>(vecptr);
403
404 if (fei_epetra_vec == NULL && fei_lpm == NULL) {
405 throw std::runtime_error("failed to obtain Epetra_MultiVector from fei::Vector.");
406 }
407
408 if (fei_epetra_vec != NULL) {
409 return( fei_epetra_vec->getUnderlyingVector());
410 }
411
412 LinProbMgr_EpetraBasic* lpm_epetrabasic =
413 dynamic_cast<LinProbMgr_EpetraBasic*>(fei_lpm->getUnderlyingVector());
414 if (lpm_epetrabasic == 0) {
415 throw std::runtime_error("fei Trilinos_Helpers: ERROR getting LinProbMgr_EpetraBasic");
416 }
417
418 if (soln_vec) {
419 return(lpm_epetrabasic->get_solution_vector().get());
420 }
421
422 return(lpm_epetrabasic->get_rhs_vector().get());
423}
424
425Epetra_VbrMatrix* get_Epetra_VbrMatrix(fei::Matrix* feimat)
426{
427 fei::Matrix_Impl<Epetra_VbrMatrix>* fei_epetra_vbr =
428 dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(feimat);
429 fei::MatrixReducer* feireducer =
430 dynamic_cast<fei::MatrixReducer*>(feimat);
431
432 if (feireducer != NULL) {
433 fei::SharedPtr<fei::Matrix> feimat2 = feireducer->getTargetMatrix();
434 fei_epetra_vbr =
435 dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(feimat2.get());
436 }
437
438 if (fei_epetra_vbr == NULL) {
439 throw std::runtime_error("failed to obtain Epetra_VbrMatrix from fei::Matrix.");
440 }
441
442 return(fei_epetra_vbr->getMatrix().get());
443}
444
445Epetra_CrsMatrix* get_Epetra_CrsMatrix(fei::Matrix* feimat)
446{
447 fei::Matrix_Impl<Epetra_CrsMatrix>* fei_epetra_crs =
448 dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(feimat);
450 dynamic_cast<fei::Matrix_Impl<fei::LinearProblemManager>*>(feimat);
451 fei::MatrixReducer* feireducer =
452 dynamic_cast<fei::MatrixReducer*>(feimat);
453
454 if (feireducer != NULL) {
455 fei::SharedPtr<fei::Matrix> feimat2 = feireducer->getTargetMatrix();
456 fei_epetra_crs =
457 dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(feimat2.get());
458 fei_lpm =
459 dynamic_cast<fei::Matrix_Impl<fei::LinearProblemManager>*>(feimat2.get());
460 }
461
462 if (fei_epetra_crs == NULL && fei_lpm == NULL) {
463 throw std::runtime_error("failed to obtain Epetra_CrsMatrix from fei::Matrix.");
464 }
465
466 if (fei_epetra_crs != NULL) {
467 return(fei_epetra_crs->getMatrix().get());
468 }
469
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");
474 }
475
476 return(lpm_epetrabasic->get_A_matrix().get());
477}
478
479
480void get_Epetra_pointers(fei::SharedPtr<fei::Matrix> feiA,
483 Epetra_CrsMatrix*& crsA,
484 Epetra_Operator*& opA,
487{
488 x = get_Epetra_MultiVector(feix.get(), true);
489 b = get_Epetra_MultiVector(feib.get(), false);
490
491 const char* matname = feiA->typeName();
492 if (!strcmp(matname, "Epetra_VbrMatrix")) {
493 Epetra_VbrMatrix* A = get_Epetra_VbrMatrix(feiA.get());
494 opA = A;
495 }
496 else {
497 crsA = get_Epetra_CrsMatrix(feiA.get());
498 opA = crsA;
499 }
500}
501
502int zero_Epetra_VbrMatrix(Epetra_VbrMatrix* mat)
503{
504 const Epetra_CrsGraph& crsgraph = mat->Graph();
505 const Epetra_BlockMap& rowmap = crsgraph.RowMap();
506 const Epetra_BlockMap& colmap = crsgraph.ColMap();
507 mat->RowMatrixRowMap();//generate point objects
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) {
514 int row = myRows[i];
515 int rowlength = 0;
516 int* colindicesView = NULL;
517 int localrow = rowmap.LID(row);
518 int err = crsgraph.ExtractMyRowView(localrow, rowlength, colindicesView);
519 if (err != 0) {
520 return err;
521 }
522 err = mat->BeginReplaceMyValues(localrow, rowlength, colindicesView);
523 if (err != 0) {
524 return err;
525 }
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);
531 if (err != 0) {
532 return err;
533 }
534 }
535 err = mat->EndSubmitEntries();
536 if (err != 0) {
537 return err;
538 }
539 }
540
541 return 0;
542}
543
544#endif //HAVE_FEI_EPETRA
545
546}//namespace Trilinos_Helpers
547
int MyGlobalElements(int *MyGlobalElementList) const
int LID(int GID) const
int ElementSize() 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
Definition: fei_Param.hpp:104
double getDoubleValue() const
Definition: fei_Param.hpp:110
ParamType getType() const
Definition: fei_Param.hpp:98
bool getBoolValue() const
Definition: fei_Param.hpp:122
const std::string & getName() const
Definition: fei_Param.hpp:92
int getIntValue() const
Definition: fei_Param.hpp:116
void add(const Param &param, bool maintain_unique_keys=true)
const_iterator end() const
const_iterator begin() const
void reset(T *p=0)
T * get() const
int localProc(MPI_Comm comm)
std::ostream & console_out()