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     delete [] used_row;
00144   }
00145 
00146 //  EpetraExt::BlockMapToMatrixMarketFile("EBMap.np12.mm",emap,"AriaTest");
00147 
00148   std::vector<int> rowLengths; rowLengths.reserve(numLocallyOwnedRows);
00149   for(int ii=0; ii<numLocallyOwnedRows; ++ii) {
00150     rowLengths.push_back(rowOffsets[ii+1]-rowOffsets[ii]);
00151   }
00152 
00153   bool staticProfile = true;
00154   Epetra_CrsGraph egraph(Copy, emap, &(rowLengths[0]), staticProfile);
00155 
00156   const Epetra_Comm& ecomm = emap.Comm();
00157   int localProc = ecomm.MyPID();
00158 
00159   int firstLocalEqn = numLocallyOwnedRows > 0 ? rowNumbers[0] : -1;
00160 
00161   int offset = 0;
00162   for(int i=0; i<numLocallyOwnedRows; ++i) {
00163     int err = egraph.InsertGlobalIndices(firstLocalEqn+i,
00164            rowLengths[i],
00165            &(packedColumnIndices[offset]));
00166     if (err != 0) {
00167       fei::console_out() << "proc " << localProc << " err-return " << err
00168                << " inserting row " << firstLocalEqn+i<<", cols ";
00169       for(int ii=0; ii<rowLengths[i]; ++ii) {
00170   fei::console_out() << packedColumnIndices[offset+ii]<<",";
00171       }
00172       fei::console_out() << FEI_ENDL;
00173       throw std::runtime_error("... occurred in create_Epetra_CrsGraph");
00174     }
00175 
00176     offset += rowLengths[i];
00177   }
00178 
00179   //Epetra_BlockMap* domainmap = const_cast<Epetra_BlockMap*>(&(epetraGraph_->DomainMap()));
00180   //Epetra_BlockMap* rangemap = const_cast<Epetra_BlockMap*>(&(epetraGraph_->RangeMap()));
00181   egraph.FillComplete();
00182 
00183   return(egraph);
00184 }
00185 
00186 fei::SharedPtr<fei::Matrix>
00187 create_from_Epetra_Matrix(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00188                           bool blockEntryMatrix,
00189                           fei::SharedPtr<fei::Reducer> reducer,
00190                           bool orderRowsWithLocalColsFirst)
00191 {
00192   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
00193   //localSize is in point-equations, because we will only use it for constructing
00194   //the fei::Matrix_Impl, and those always deal in point-equations.
00195   int localSize = vecSpace->getNumIndices_Owned();
00196 
00197   fei::SharedPtr<fei::Matrix> tmpmat;
00198   fei::SharedPtr<fei::Matrix> feimat;
00199   try {
00200     Epetra_CrsGraph egraph =
00201       Trilinos_Helpers::create_Epetra_CrsGraph(matrixGraph, blockEntryMatrix,
00202                                                orderRowsWithLocalColsFirst);
00203 
00204     if (blockEntryMatrix) {
00205       fei::SharedPtr<Epetra_VbrMatrix>
00206         epetraMatrix(new Epetra_VbrMatrix(Copy, egraph));
00207 
00208       tmpmat.reset(new fei::Matrix_Impl<Epetra_VbrMatrix>(epetraMatrix,
00209                                                    matrixGraph, localSize));
00210       zero_Epetra_VbrMatrix(epetraMatrix.get());
00211     }
00212     else {
00213       fei::SharedPtr<Epetra_CrsMatrix>
00214         epetraMatrix(new Epetra_CrsMatrix(Copy, egraph));
00215 
00216       tmpmat.reset(new fei::Matrix_Impl<Epetra_CrsMatrix>(epetraMatrix,
00217                                           matrixGraph, localSize));
00218     }
00219   }
00220   catch(std::runtime_error& exc) {
00221     fei::console_out() << "Trilinos_Helpers::create_from_Epetra_Matrix ERROR, "
00222            << "caught exception: '" << exc.what() << "', rethrowing..."
00223            << FEI_ENDL;
00224     throw exc;
00225   }
00226 
00227   if (reducer.get() != NULL) {
00228     feimat.reset(new fei::MatrixReducer(reducer, tmpmat));
00229   }
00230   else {
00231     feimat = tmpmat;
00232   }
00233 
00234   return(feimat);
00235 }
00236 
00237 fei::SharedPtr<fei::Matrix>
00238 create_from_LPM_EpetraBasic(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00239                             bool blockEntryMatrix,
00240                             fei::SharedPtr<fei::Reducer> reducer,
00241                             fei::SharedPtr<fei::LinearProblemManager>
00242                               lpm_epetrabasic)
00243 {
00244   fei::SharedPtr<fei::SparseRowGraph> srgraph =
00245     matrixGraph->createGraph(blockEntryMatrix);
00246   if (srgraph.get() == NULL) {
00247     throw std::runtime_error("create_LPM_EpetraBasic ERROR in fei::MatrixGraph::createGraph");
00248   }
00249 
00250   fei::SharedPtr<fei::SparseRowGraph> sharedsrg(srgraph);
00251   int localSize;
00252   if (reducer.get() != NULL) {
00253     std::vector<int>& reduced_eqns = reducer->getLocalReducedEqns();
00254     lpm_epetrabasic->setRowDistribution(reduced_eqns);
00255     lpm_epetrabasic->setMatrixGraph(sharedsrg);
00256     localSize = reduced_eqns.size();
00257   }
00258   else {
00259     fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
00260     int err = 0,chkNum;
00261     std::vector<int> indices;
00262     if (blockEntryMatrix) {
00263       localSize = vecSpace->getNumBlkIndices_Owned();
00264       indices.resize(localSize*2);
00265       err = vecSpace->getBlkIndices_Owned(localSize, &indices[0],
00266                                   &indices[localSize], chkNum);
00267     }
00268     else {
00269       localSize = vecSpace->getNumIndices_Owned();
00270       indices.resize(localSize);
00271       err = vecSpace->getIndices_Owned(indices);
00272     }
00273     if (err != 0) {
00274       throw std::runtime_error("Factory_Trilinos: createMatrix: error in vecSpace->getIndices_Owned");
00275     }
00276 
00277     lpm_epetrabasic->setRowDistribution(indices);
00278     lpm_epetrabasic->setMatrixGraph(sharedsrg);
00279   }
00280 
00281   fei::SharedPtr<fei::Matrix> tmpmat(new fei::Matrix_Impl<fei::LinearProblemManager>(lpm_epetrabasic, matrixGraph, localSize));
00282 
00283   fei::SharedPtr<fei::Matrix> feimat;
00284 
00285   if (reducer.get() != NULL) {
00286     feimat.reset(new fei::MatrixReducer(reducer, tmpmat));
00287   }
00288   else {
00289     feimat = tmpmat;
00290   }
00291 
00292   return(feimat);
00293 }
00294 #endif //HAVE_FEI_EPETRA
00295 
00296 void copy_parameterset(const fei::ParameterSet& paramset,
00297                        Teuchos::ParameterList& paramlist)
00298 {
00299   fei::ParameterSet::const_iterator
00300     iter = paramset.begin(),
00301     iter_end = paramset.end();
00302 
00303   for(; iter != iter_end; ++iter) {
00304     const fei::Param& param = *iter;
00305     fei::Param::ParamType ptype = param.getType();
00306     switch(ptype) {
00307     case fei::Param::STRING:
00308       paramlist.set(param.getName(), param.getStringValue());
00309       break;
00310     case fei::Param::DOUBLE:
00311       paramlist.set(param.getName(), param.getDoubleValue());
00312       break;
00313     case fei::Param::INT:
00314       paramlist.set(param.getName(), param.getIntValue());
00315       break;
00316     case fei::Param::BOOL:
00317       paramlist.set(param.getName(), param.getBoolValue());
00318       break;
00319     default:
00320       break;
00321     }
00322   }
00323 }
00324 
00325 void copy_parameterlist(const Teuchos::ParameterList& paramlist,
00326                         fei::ParameterSet& paramset)
00327 {
00328   Teuchos::ParameterList::ConstIterator
00329     iter = paramlist.begin(),
00330     iter_end = paramlist.end();
00331 
00332   for(; iter != iter_end; ++iter) {
00333     const Teuchos::ParameterEntry& param = paramlist.entry(iter);
00334     if (param.isType<std::string>()) {
00335       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<std::string>(param).c_str()));
00336     }
00337     else if (param.isType<double>()) {
00338       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<double>(param)));
00339     }
00340     else if (param.isType<int>()) {
00341       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<int>(param)));
00342     }
00343     else if (param.isType<bool>()) {
00344       paramset.add(fei::Param(paramlist.name(iter).c_str(), Teuchos::getValue<bool>(param)));
00345     }
00346   }
00347 }
00348 
00349 #ifdef HAVE_FEI_EPETRA
00350 
00351 Epetra_MultiVector*
00352 get_Epetra_MultiVector(fei::Vector* feivec, bool soln_vec)
00353 {
00354   fei::Vector* vecptr = feivec;
00355   fei::VectorReducer* feireducer = dynamic_cast<fei::VectorReducer*>(feivec);
00356   if (feireducer != NULL) vecptr = feireducer->getTargetVector().get();
00357 
00358   fei::Vector_Impl<Epetra_MultiVector>* fei_epetra_vec =
00359     dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>*>(vecptr);
00360   fei::Vector_Impl<fei::LinearProblemManager>* fei_lpm =
00361     dynamic_cast<fei::Vector_Impl<fei::LinearProblemManager>*>(vecptr);
00362 
00363   if (fei_epetra_vec == NULL && fei_lpm == NULL) {
00364     throw std::runtime_error("failed to obtain Epetra_MultiVector from fei::Vector.");
00365   }
00366 
00367   if (fei_epetra_vec != NULL) {
00368     return( fei_epetra_vec->getUnderlyingVector());
00369   }
00370 
00371   LinProbMgr_EpetraBasic* lpm_epetrabasic =
00372     dynamic_cast<LinProbMgr_EpetraBasic*>(fei_lpm->getUnderlyingVector());
00373   if (lpm_epetrabasic == 0) {
00374     throw std::runtime_error("fei Trilinos_Helpers: ERROR getting LinProbMgr_EpetraBasic");
00375   }
00376 
00377   if (soln_vec) {
00378     return(lpm_epetrabasic->get_solution_vector().get());
00379   }
00380 
00381   return(lpm_epetrabasic->get_rhs_vector().get());
00382 }
00383 
00384 Epetra_VbrMatrix* get_Epetra_VbrMatrix(fei::Matrix* feimat)
00385 {
00386   fei::Matrix_Impl<Epetra_VbrMatrix>* fei_epetra_vbr =
00387     dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(feimat);
00388   fei::MatrixReducer* feireducer =
00389     dynamic_cast<fei::MatrixReducer*>(feimat);
00390 
00391   if (feireducer != NULL) {
00392     fei::SharedPtr<fei::Matrix> feimat2 = feireducer->getTargetMatrix();
00393     fei_epetra_vbr =
00394       dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(feimat2.get());
00395   }
00396 
00397   if (fei_epetra_vbr == NULL) {
00398     throw std::runtime_error("failed to obtain Epetra_VbrMatrix from fei::Matrix.");
00399   }
00400 
00401   return(fei_epetra_vbr->getMatrix().get());
00402 }
00403 
00404 Epetra_CrsMatrix* get_Epetra_CrsMatrix(fei::Matrix* feimat)
00405 {
00406   fei::Matrix_Impl<Epetra_CrsMatrix>* fei_epetra_crs =
00407     dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(feimat);
00408   fei::Matrix_Impl<fei::LinearProblemManager>* fei_lpm =
00409     dynamic_cast<fei::Matrix_Impl<fei::LinearProblemManager>*>(feimat);
00410   fei::MatrixReducer* feireducer =
00411     dynamic_cast<fei::MatrixReducer*>(feimat);
00412 
00413   if (feireducer != NULL) {
00414     fei::SharedPtr<fei::Matrix> feimat2 = feireducer->getTargetMatrix();
00415     fei_epetra_crs =
00416       dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(feimat2.get());
00417     fei_lpm =
00418       dynamic_cast<fei::Matrix_Impl<fei::LinearProblemManager>*>(feimat2.get());
00419   }
00420 
00421   if (fei_epetra_crs == NULL && fei_lpm == NULL) {
00422     throw std::runtime_error("failed to obtain Epetra_CrsMatrix from fei::Matrix.");
00423   }
00424 
00425   if (fei_epetra_crs != NULL) {
00426     return(fei_epetra_crs->getMatrix().get());
00427   }
00428 
00429   LinProbMgr_EpetraBasic* lpm_epetrabasic =
00430     dynamic_cast<LinProbMgr_EpetraBasic*>(fei_lpm->getMatrix().get());
00431   if (lpm_epetrabasic == 0) {
00432     throw std::runtime_error("fei Trilinos_Helpers ERROR getting LinProbMgr_EpetraBasic");
00433   }
00434 
00435   return(lpm_epetrabasic->get_A_matrix().get());
00436 }
00437 
00438 
00439 void get_Epetra_pointers(fei::SharedPtr<fei::Matrix> feiA,
00440                          fei::SharedPtr<fei::Vector> feix,
00441                          fei::SharedPtr<fei::Vector> feib,
00442                          Epetra_CrsMatrix*& crsA,
00443                          Epetra_Operator*& opA,
00444                          Epetra_MultiVector*& x,
00445                          Epetra_MultiVector*& b)
00446 {
00447   x = get_Epetra_MultiVector(feix.get(), true);
00448   b = get_Epetra_MultiVector(feib.get(), false);
00449 
00450   const char* matname = feiA->typeName();
00451   if (!strcmp(matname, "Epetra_VbrMatrix")) {
00452     Epetra_VbrMatrix* A = get_Epetra_VbrMatrix(feiA.get());
00453     opA = A;
00454   }
00455   else {
00456     crsA = get_Epetra_CrsMatrix(feiA.get());
00457     opA = crsA;
00458   }
00459 }
00460 
00461 int zero_Epetra_VbrMatrix(Epetra_VbrMatrix* mat)
00462 {
00463   const Epetra_CrsGraph& crsgraph = mat->Graph();
00464   const Epetra_BlockMap& rowmap = crsgraph.RowMap();
00465   const Epetra_BlockMap& colmap = crsgraph.ColMap();
00466   mat->RowMatrixRowMap();//generate point objects
00467   int maxBlkRowSize = mat->GlobalMaxRowDim();
00468   int maxBlkColSize = mat->GlobalMaxColDim();
00469   std::vector<double> zeros(maxBlkRowSize*maxBlkColSize, 0);
00470   int numMyRows = rowmap.NumMyElements();
00471   int* myRows = rowmap.MyGlobalElements();
00472   for(int i=0; i<numMyRows; ++i) {
00473     int row = myRows[i];
00474     int rowlength = 0;
00475     int* colindicesView = NULL;
00476     int localrow = rowmap.LID(row);
00477     int err = crsgraph.ExtractMyRowView(localrow, rowlength, colindicesView);
00478     if (err != 0) {
00479       return err;
00480     }
00481     err = mat->BeginReplaceMyValues(localrow, rowlength, colindicesView);
00482     if (err != 0) {
00483       return err;
00484     }
00485     int blkRowSize = rowmap.ElementSize(localrow);
00486     for(int j=0; j<rowlength; ++j) {
00487       int blkColSize = colmap.ElementSize(colindicesView[j]);
00488       err = mat->SubmitBlockEntry(&zeros[0], maxBlkRowSize,
00489                             blkRowSize, blkColSize);
00490       if (err != 0) {
00491         return err;
00492       }
00493     }
00494     err = mat->EndSubmitEntries();
00495     if (err != 0) {
00496       return err;
00497     }
00498   }
00499 
00500   return 0;
00501 }
00502 
00503 #endif //HAVE_FEI_EPETRA
00504 
00505 }//namespace Trilinos_Helpers
00506 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends