FEI Version of the Day
fei_Trilinos_Helpers.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2006 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 #include <fei_macros.hpp>
00009 
00010 #include <fei_Include_Trilinos.hpp>
00011 
00012 #ifdef HAVE_FEI_EPETRA
00013 #include <fei_VectorTraits_Epetra.hpp>
00014 #include <fei_MatrixTraits_Epetra.hpp>
00015 #include <fei_LinProbMgr_EpetraBasic.hpp>
00016 #endif
00017 
00018 #include <fei_Trilinos_Helpers.hpp>
00019 #include <fei_ParameterSet.hpp>
00020 #include <fei_Vector_Impl.hpp>
00021 #include <fei_VectorReducer.hpp>
00022 #include <fei_Matrix_Impl.hpp>
00023 #include <fei_MatrixReducer.hpp>
00024 //#include <EpetraExt_BlockMapOut.h>
00025 
00026 namespace Trilinos_Helpers {
00027 
00028 #ifdef HAVE_FEI_EPETRA
00029 
00030 Epetra_Map
00031 create_Epetra_Map(MPI_Comm comm,
00032                   const std::vector<int>& local_eqns)
00033 {
00034 #ifndef FEI_SER
00035   Epetra_MpiComm EComm(comm);
00036 #else
00037   Epetra_SerialComm EComm;
00038 #endif
00039 
00040   int localSize = local_eqns.size();
00041   int globalSize = 0;
00042   EComm.SumAll(&localSize, &globalSize, 1);
00043 
00044   if (localSize < 0 || globalSize < 0) {
00045     throw std::runtime_error("Trilinos_Helpers::create_Epetra_Map: negative local or global size.");
00046   }
00047 
00048   Epetra_Map EMap(globalSize, localSize, 0, EComm);
00049   return(EMap);
00050 }
00051 
00052 Epetra_BlockMap
00053 create_Epetra_BlockMap(const fei::SharedPtr<fei::VectorSpace>& vecspace)
00054 {
00055   if (vecspace.get() == 0) {
00056     throw std::runtime_error("create_Epetra_Map needs non-null fei::VectorSpace");
00057   }
00058 
00059 #ifndef FEI_SER
00060   MPI_Comm comm = vecspace->getCommunicator();
00061   Epetra_MpiComm EComm(comm);
00062 #else
00063   Epetra_SerialComm EComm;
00064 #endif
00065 
00066   int localSizeBlk = vecspace->getNumBlkIndices_Owned();
00067   int globalSizeBlk = vecspace->getGlobalNumBlkIndices();
00068 
00069   if (localSizeBlk < 0 || globalSizeBlk < 0) {
00070     throw std::runtime_error("Trilinos_Helpers::create_Epetra_BlockMap: fei::VectorSpace has negative local or global size.");
00071   }
00072 
00073   std::vector<int> blkEqns(localSizeBlk*2);
00074   int* blkEqnsPtr = &(blkEqns[0]);
00075 
00076   int chkNum = 0;
00077   int errcode = vecspace->getBlkIndices_Owned(localSizeBlk,
00078                 blkEqnsPtr, blkEqnsPtr+localSizeBlk,
00079                 chkNum);
00080   if (errcode != 0) {
00081     FEI_OSTRINGSTREAM osstr;
00082     osstr << "create_Epetra_BlockMap ERROR, nonzero errcode="<<errcode
00083     << " returned by vecspace->getBlkIndices_Owned.";
00084     throw std::runtime_error(osstr.str());
00085   }
00086 
00087   Epetra_BlockMap EBMap(globalSizeBlk, localSizeBlk,
00088       blkEqnsPtr, blkEqnsPtr+localSizeBlk, 0, EComm);
00089 
00090   return(EBMap);
00091 }
00092 
00093 Epetra_CrsGraph
00094 create_Epetra_CrsGraph(const fei::SharedPtr<fei::MatrixGraph>& matgraph,
00095            bool blockEntries, bool orderRowsWithLocalColsFirst)
00096 {
00097   fei::SharedPtr<fei::SparseRowGraph> localSRGraph =
00098     matgraph->createGraph(blockEntries);
00099   if (localSRGraph.get() == NULL) {
00100     throw std::runtime_error("create_Epetra_CrsGraph ERROR in fei::MatrixGraph::createGraph");
00101   }
00102 
00103   int numLocallyOwnedRows = localSRGraph->rowNumbers.size();
00104   int* rowNumbers = &(localSRGraph->rowNumbers[0]);
00105   int* rowOffsets = &(localSRGraph->rowOffsets[0]);
00106   int* packedColumnIndices = &(localSRGraph->packedColumnIndices[0]);
00107 
00108   fei::SharedPtr<fei::VectorSpace> vecspace = matgraph->getRowSpace();
00109   MPI_Comm comm = vecspace->getCommunicator();
00110   std::vector<int>& local_eqns = localSRGraph->rowNumbers;
00111 
00112   Epetra_BlockMap emap = blockEntries ?
00113     create_Epetra_BlockMap(vecspace) : create_Epetra_Map(comm, local_eqns);
00114 
00115   if (orderRowsWithLocalColsFirst == true &&
00116       emap.Comm().NumProc() > 2 && !blockEntries) {
00117     bool* used_row = new bool[local_eqns.size()];
00118     for(unsigned ii=0; ii<local_eqns.size(); ++ii) used_row[ii] = false;
00119 
00120     int offset = 0;
00121     std::vector<int> ordered_local_eqns(local_eqns.size());
00122     for(unsigned ii=0; ii<local_eqns.size(); ++ii) {
00123       bool row_has_off_proc_cols = false;
00124       for(int j=rowOffsets[ii]; j<rowOffsets[ii+1]; ++j) {
00125         if (emap.MyGID(packedColumnIndices[j]) == false) {
00126           row_has_off_proc_cols = true;
00127           break;
00128         }
00129       }
00130 
00131       if (row_has_off_proc_cols == false) {
00132         ordered_local_eqns[offset++] = rowNumbers[ii];
00133         used_row[ii] = true;
00134       }
00135     }
00136 
00137     for(unsigned ii=0; ii<local_eqns.size(); ++ii) {
00138       if (used_row[ii] == true) continue;
00139       ordered_local_eqns[offset++] = rowNumbers[ii];
00140     }
00141 
00142     emap = create_Epetra_Map(comm, ordered_local_eqns);
00143   }
00144 
00145 //  EpetraExt::BlockMapToMatrixMarketFile("EBMap.np12.mm",emap,"AriaTest");
00146 
00147   std::vector<int> rowLengths(numLocallyOwnedRows);
00148   for(int ii=0; ii<numLocallyOwnedRows; ++ii) {
00149     rowLengths[ii] = rowOffsets[ii+1]-rowOffsets[ii];
00150   }
00151 
00152   bool staticProfile = true;
00153   Epetra_CrsGraph egraph(Copy, emap, &(rowLengths[0]), staticProfile);
00154 
00155   const Epetra_Comm& ecomm = emap.Comm();
00156   int localProc = ecomm.MyPID();
00157 
00158   int firstLocalEqn = numLocallyOwnedRows > 0 ? rowNumbers[0] : -1;
00159 
00160   int offset = 0;
00161   for(int i=0; i<numLocallyOwnedRows; ++i) {
00162     int err = egraph.InsertGlobalIndices(firstLocalEqn+i,
00163            rowLengths[i],
00164            &(packedColumnIndices[offset]));
00165     if (err != 0) {
00166       FEI_CERR << "proc " << localProc << " err-return " << err
00167                << " inserting row " << firstLocalEqn+i<<", cols ";
00168       for(int ii=0; ii<rowLengths[i]; ++ii) {
00169   FEI_CERR << packedColumnIndices[offset+ii]<<",";
00170       }
00171       FEI_CERR << FEI_ENDL;
00172       throw std::runtime_error("... occurred in create_Epetra_CrsGraph");
00173     }
00174 
00175     offset += rowLengths[i];
00176   }
00177 
00178   //Epetra_BlockMap* domainmap = const_cast<Epetra_BlockMap*>(&(epetraGraph_->DomainMap()));
00179   //Epetra_BlockMap* rangemap = const_cast<Epetra_BlockMap*>(&(epetraGraph_->RangeMap()));
00180   egraph.FillComplete();
00181 
00182   return(egraph);
00183 }
00184 
00185 fei::SharedPtr<fei::Matrix>
00186 create_from_Epetra_Matrix(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00187                           bool blockEntryMatrix,
00188                           fei::SharedPtr<fei::Reducer> reducer,
00189                           bool orderRowsWithLocalColsFirst)
00190 {
00191   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
00192   //localSize is in point-equations, because we will only use it for constructing
00193   //the fei::Matrix_Impl, and those always deal in point-equations.
00194   int localSize = vecSpace->getNumIndices_Owned();
00195 
00196   fei::SharedPtr<fei::Matrix> tmpmat;
00197   fei::SharedPtr<fei::Matrix> feimat;
00198   try {
00199     Epetra_CrsGraph egraph =
00200       Trilinos_Helpers::create_Epetra_CrsGraph(matrixGraph, blockEntryMatrix,
00201                                                orderRowsWithLocalColsFirst);
00202 
00203     if (blockEntryMatrix) {
00204       fei::SharedPtr<Epetra_VbrMatrix>
00205         epetraMatrix(new Epetra_VbrMatrix(Copy, egraph));
00206 
00207       tmpmat.reset(new fei::Matrix_Impl<Epetra_VbrMatrix>(epetraMatrix,
00208                                                    matrixGraph, localSize));
00209       zero_Epetra_VbrMatrix(epetraMatrix.get());
00210     }
00211     else {
00212       fei::SharedPtr<Epetra_CrsMatrix>
00213         epetraMatrix(new Epetra_CrsMatrix(Copy, egraph));
00214 
00215       tmpmat.reset(new fei::Matrix_Impl<Epetra_CrsMatrix>(epetraMatrix,
00216                                           matrixGraph, localSize));
00217     }
00218   }
00219   catch(std::runtime_error& exc) {
00220     FEI_CERR << "Trilinos_Helpers::create_from_Epetra_Matrix ERROR, "
00221            << "caught exception: '" << exc.what() << "', rethrowing..."
00222            << FEI_ENDL;
00223     throw exc;
00224   }
00225 
00226   if (reducer.get() != NULL) {
00227     feimat.reset(new fei::MatrixReducer(reducer, tmpmat));
00228   }
00229   else {
00230     feimat = tmpmat;
00231   }
00232 
00233   return(feimat);
00234 }
00235 
00236 fei::SharedPtr<fei::Matrix>
00237 create_from_LPM_EpetraBasic(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00238                             bool blockEntryMatrix,
00239                             fei::SharedPtr<fei::Reducer> reducer,
00240                             fei::SharedPtr<fei::LinearProblemManager>
00241                               lpm_epetrabasic)
00242 {
00243   fei::SharedPtr<fei::SparseRowGraph> srgraph =
00244     matrixGraph->createGraph(blockEntryMatrix);
00245   if (srgraph.get() == NULL) {
00246     throw std::runtime_error("create_LPM_EpetraBasic ERROR in fei::MatrixGraph::createGraph");
00247   }
00248 
00249   fei::SharedPtr<fei::SparseRowGraph> sharedsrg(srgraph);
00250   int localSize;
00251   if (reducer.get() != NULL) {
00252     std::vector<int>& reduced_eqns = reducer->getLocalReducedEqns();
00253     lpm_epetrabasic->setRowDistribution(reduced_eqns);
00254     lpm_epetrabasic->setMatrixGraph(sharedsrg);
00255     localSize = reduced_eqns.size();
00256   }
00257   else {
00258     fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
00259     int err = 0,chkNum;
00260     std::vector<int> indices;
00261     if (blockEntryMatrix) {
00262       localSize = vecSpace->getNumBlkIndices_Owned();
00263       indices.resize(localSize*2);
00264       err = vecSpace->getBlkIndices_Owned(localSize, &indices[0],
00265                                   &indices[localSize], chkNum);
00266     }
00267     else {
00268       localSize = vecSpace->getNumIndices_Owned();
00269       indices.resize(localSize);
00270       err = vecSpace->getIndices_Owned(indices);
00271     }
00272     if (err != 0) {
00273       throw std::runtime_error("Factory_Trilinos: createMatrix: error in vecSpace->getIndices_Owned");
00274     }
00275 
00276     lpm_epetrabasic->setRowDistribution(indices);
00277     lpm_epetrabasic->setMatrixGraph(sharedsrg);
00278   }
00279 
00280   fei::SharedPtr<fei::Matrix> tmpmat(new fei::Matrix_Impl<fei::LinearProblemManager>(lpm_epetrabasic, matrixGraph, localSize));
00281 
00282   fei::SharedPtr<fei::Matrix> feimat;
00283 
00284   if (reducer.get() != NULL) {
00285     feimat.reset(new fei::MatrixReducer(reducer, tmpmat));
00286   }
00287   else {
00288     feimat = tmpmat;
00289   }
00290 
00291   return(feimat);
00292 }
00293 #endif //HAVE_FEI_EPETRA
00294 
00295 void copy_parameterset(const fei::ParameterSet& paramset,
00296                        Teuchos::ParameterList& paramlist)
00297 {
00298   fei::ParameterSet::const_iterator
00299     iter = paramset.begin(),
00300     iter_end = paramset.end();
00301 
00302   for(; iter != iter_end; ++iter) {
00303     const fei::Param& param = *iter;
00304     fei::Param::ParamType ptype = param.getType();
00305     switch(ptype) {
00306     case fei::Param::STRING:
00307       paramlist.set(param.getName(), param.getStringValue());
00308       break;
00309     case fei::Param::DOUBLE:
00310       paramlist.set(param.getName(), param.getDoubleValue());
00311       break;
00312     case fei::Param::INT:
00313       paramlist.set(param.getName(), param.getIntValue());
00314       break;
00315     case fei::Param::BOOL:
00316       paramlist.set(param.getName(), param.getBoolValue());
00317       break;
00318     default:
00319       break;
00320     }
00321   }
00322 }
00323 
00324 void copy_parameterlist(const Teuchos::ParameterList& paramlist,
00325                         fei::ParameterSet& paramset)
00326 {
00327   Teuchos::ParameterList::ConstIterator
00328     iter = paramlist.begin(),
00329     iter_end = paramlist.end();
00330 
00331   for(; iter != iter_end; ++iter) {
00332     const Teuchos::ParameterEntry& param = paramlist.entry(iter);
00333     if (param.isType<std::string>()) {
00334       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<std::string>(param).c_str()));
00335     }
00336     else if (param.isType<double>()) {
00337       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<double>(param)));
00338     }
00339     else if (param.isType<int>()) {
00340       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<int>(param)));
00341     }
00342     else if (param.isType<bool>()) {
00343       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<bool>(param)));
00344     }
00345   }
00346 }
00347 
00348 #ifdef HAVE_FEI_EPETRA
00349 
00350 Epetra_MultiVector*
00351 get_Epetra_MultiVector(fei::Vector* feivec, bool soln_vec)
00352 {
00353   fei::Vector* vecptr = feivec;
00354   fei::VectorReducer* feireducer = dynamic_cast<fei::VectorReducer*>(feivec);
00355   if (feireducer != NULL) vecptr = feireducer->getTargetVector().get();
00356 
00357   fei::Vector_Impl<Epetra_MultiVector>* fei_epetra_vec =
00358     dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>*>(vecptr);
00359   fei::Vector_Impl<fei::LinearProblemManager>* fei_lpm =
00360     dynamic_cast<fei::Vector_Impl<fei::LinearProblemManager>*>(vecptr);
00361 
00362   if (fei_epetra_vec == NULL && fei_lpm == NULL) {
00363     throw std::runtime_error("failed to obtain Epetra_MultiVector from fei::Vector.");
00364   }
00365 
00366   if (fei_epetra_vec != NULL) {
00367     return( fei_epetra_vec->getUnderlyingVector());
00368   }
00369 
00370   LinProbMgr_EpetraBasic* lpm_epetrabasic =
00371     dynamic_cast<LinProbMgr_EpetraBasic*>(fei_lpm->getUnderlyingVector());
00372   if (lpm_epetrabasic == 0) {
00373     throw std::runtime_error("fei Trilinos_Helpers: ERROR getting LinProbMgr_EpetraBasic");
00374   }
00375 
00376   if (soln_vec) {
00377     return(lpm_epetrabasic->get_solution_vector().get());
00378   }
00379 
00380   return(lpm_epetrabasic->get_rhs_vector().get());
00381 }
00382 
00383 Epetra_VbrMatrix* get_Epetra_VbrMatrix(fei::Matrix* feimat)
00384 {
00385   fei::Matrix_Impl<Epetra_VbrMatrix>* fei_epetra_vbr =
00386     dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(feimat);
00387   fei::MatrixReducer* feireducer =
00388     dynamic_cast<fei::MatrixReducer*>(feimat);
00389 
00390   if (feireducer != NULL) {
00391     fei::SharedPtr<fei::Matrix> feimat2 = feireducer->getTargetMatrix();
00392     fei_epetra_vbr =
00393       dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(feimat2.get());
00394   }
00395 
00396   if (fei_epetra_vbr == NULL) {
00397     throw std::runtime_error("failed to obtain Epetra_VbrMatrix from fei::Matrix.");
00398   }
00399 
00400   return(fei_epetra_vbr->getMatrix().get());
00401 }
00402 
00403 Epetra_CrsMatrix* get_Epetra_CrsMatrix(fei::Matrix* feimat)
00404 {
00405   fei::Matrix_Impl<Epetra_CrsMatrix>* fei_epetra_crs =
00406     dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(feimat);
00407   fei::Matrix_Impl<fei::LinearProblemManager>* fei_lpm =
00408     dynamic_cast<fei::Matrix_Impl<fei::LinearProblemManager>*>(feimat);
00409   fei::MatrixReducer* feireducer =
00410     dynamic_cast<fei::MatrixReducer*>(feimat);
00411 
00412   if (feireducer != NULL) {
00413     fei::SharedPtr<fei::Matrix> feimat2 = feireducer->getTargetMatrix();
00414     fei_epetra_crs =
00415       dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(feimat2.get());
00416     fei_lpm =
00417       dynamic_cast<fei::Matrix_Impl<fei::LinearProblemManager>*>(feimat2.get());
00418   }
00419 
00420   if (fei_epetra_crs == NULL && fei_lpm == NULL) {
00421     throw std::runtime_error("failed to obtain Epetra_CrsMatrix from fei::Matrix.");
00422   }
00423 
00424   if (fei_epetra_crs != NULL) {
00425     return(fei_epetra_crs->getMatrix().get());
00426   }
00427 
00428   LinProbMgr_EpetraBasic* lpm_epetrabasic =
00429     dynamic_cast<LinProbMgr_EpetraBasic*>(fei_lpm->getMatrix().get());
00430   if (lpm_epetrabasic == 0) {
00431     throw std::runtime_error("fei Trilinos_Helpers ERROR getting LinProbMgr_EpetraBasic");
00432   }
00433 
00434   return(lpm_epetrabasic->get_A_matrix().get());
00435 }
00436 
00437 
00438 void get_Epetra_pointers(fei::SharedPtr<fei::Matrix> feiA,
00439                          fei::SharedPtr<fei::Vector> feix,
00440                          fei::SharedPtr<fei::Vector> feib,
00441                          Epetra_CrsMatrix*& crsA,
00442                          Epetra_Operator*& opA,
00443                          Epetra_MultiVector*& x,
00444                          Epetra_MultiVector*& b)
00445 {
00446   x = get_Epetra_MultiVector(feix.get(), true);
00447   b = get_Epetra_MultiVector(feib.get(), false);
00448 
00449   const char* matname = feiA->typeName();
00450   if (!strcmp(matname, "Epetra_VbrMatrix")) {
00451     Epetra_VbrMatrix* A = get_Epetra_VbrMatrix(feiA.get());
00452     opA = A;
00453   }
00454   else {
00455     crsA = get_Epetra_CrsMatrix(feiA.get());
00456     opA = crsA;
00457   }
00458 }
00459 
00460 int zero_Epetra_VbrMatrix(Epetra_VbrMatrix* mat)
00461 {
00462   const Epetra_CrsGraph& crsgraph = mat->Graph();
00463   const Epetra_BlockMap& rowmap = crsgraph.RowMap();
00464   const Epetra_BlockMap& colmap = crsgraph.ColMap();
00465   int maxBlkRowSize = mat->GlobalMaxRowDim();
00466   int maxBlkColSize = mat->GlobalMaxColDim();
00467   std::vector<double> zeros(maxBlkRowSize*maxBlkColSize, 0);
00468   int numMyRows = rowmap.NumMyElements();
00469   int* myRows = rowmap.MyGlobalElements();
00470   for(int i=0; i<numMyRows; ++i) {
00471     int row = myRows[i];
00472     int rowlength = 0;
00473     int* colindicesView = NULL;
00474     int localrow = rowmap.LID(row);
00475     int err = crsgraph.ExtractMyRowView(localrow, rowlength, colindicesView);
00476     if (err != 0) {
00477       return err;
00478     }
00479     err = mat->BeginReplaceMyValues(localrow, rowlength, colindicesView);
00480     if (err != 0) {
00481       return err;
00482     }
00483     int blkRowSize = rowmap.ElementSize(localrow);
00484     for(int j=0; j<rowlength; ++j) {
00485       int blkColSize = colmap.ElementSize(colindicesView[j]);
00486       err = mat->SubmitBlockEntry(&zeros[0], maxBlkRowSize,
00487                             blkRowSize, blkColSize);
00488       if (err != 0) {
00489         return err;
00490       }
00491     }
00492     err = mat->EndSubmitEntries();
00493     if (err != 0) {
00494       return err;
00495     }
00496   }
00497 
00498   return 0;
00499 }
00500 
00501 #endif //HAVE_FEI_EPETRA
00502 
00503 }//namespace Trilinos_Helpers
00504 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends