Teko Version of the Day
Teko_MLPreconditionerFactory.cpp
00001 #include "Teko_MLPreconditionerFactory.hpp"
00002 
00003 #include "Teko_MLLinearOp.hpp"
00004 
00005 #include "ml_include.h"
00006 #include "ml_MultiLevelPreconditioner.h"
00007 #include "ml_epetra_utils.h"
00008 #include "ml_RowMatrix.h"
00009 
00010 #include "Thyra_get_Epetra_Operator.hpp"
00011 
00012 namespace Teko {
00013 
00015 // MLPreconditionerState
00017 
00018 MLPreconditionerState::MLPreconditionerState() : isFilled_(false) {}
00019 
00020 void MLPreconditionerState::setMLComm(ML_Comm * comm)
00021 {
00022    mlComm_ = Teuchos::rcpWithDealloc(comm,Teuchos::deallocFunctorDelete<ML_Comm>(cleanup_ML_Comm));
00023 }
00024 
00025 void MLPreconditionerState::setMLOperator(ML_Operator * op)
00026 {
00027    mlOp_ = Teuchos::rcpWithDealloc(op,Teuchos::deallocFunctorDelete<ML_Operator>(cleanup_ML_Operator));
00028 }
00029 
00030 void MLPreconditionerState::setIsFilled(bool value)
00031 {
00032    isFilled_ = value;
00033 }
00034 
00035 bool MLPreconditionerState::isFilled() const
00036 {
00037    return isFilled_;
00038 }
00039 
00040 void MLPreconditionerState::cleanup_ML_Comm(ML_Comm * mlComm)
00041 {
00042    ML_Comm_Destroy(&mlComm);
00043 }
00044 
00045 void MLPreconditionerState::cleanup_ML_Operator(ML_Operator * mlOp)
00046 {
00047    ML_Operator_Destroy(&mlOp);
00048 }
00049 
00050 void MLPreconditionerState::setAggregationMatrices(const std::vector<Epetra_RowMatrix *> & diags)
00051 {
00052    diagonalOps_ = diags;
00053 }
00054 
00055 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> 
00056 MLPreconditionerState::constructMLPreconditioner(const Teuchos::ParameterList & mainList,
00057                                                  const std::vector<Teuchos::RCP<const Teuchos::ParameterList> > & coarseningParams)
00058 {
00059    TEUCHOS_ASSERT(isFilled());
00060    std::vector<Teuchos::ParameterList> cpls(coarseningParams.size());
00061    for(std::size_t i=0;i<coarseningParams.size();i++)
00062       cpls[i] = *coarseningParams[i];
00063  
00064    mlPreconditioner_ = rcp(new ML_Epetra::MultiLevelPreconditioner(&*mlOp_, mainList, &diagonalOps_[0], &cpls[0], diagonalOps_.size()));
00065 
00066    return mlPreconditioner_;
00067 }
00068 
00070 // MLPreconditionerFactory
00072 
00073 MLPreconditionerFactory::MLPreconditionerFactory()
00074 {
00075 }
00076 
00077 LinearOp MLPreconditionerFactory::buildPreconditionerOperator(BlockedLinearOp & blo,BlockPreconditionerState & state) const
00078 {
00079    MLPreconditionerState & mlState = Teuchos::dyn_cast<MLPreconditionerState>(state);
00080 
00081    if(not mlState.isFilled())
00082       fillMLPreconditionerState(blo,mlState);
00083    TEUCHOS_ASSERT(mlState.isFilled());
00084 
00085    Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> mlPrecOp = 
00086          mlState.constructMLPreconditioner(mainParams_,coarseningParams_);
00087 
00088    // return Thyra::epetraLinearOp(mlPrecOp);
00089    return Teuchos::rcp(new MLLinearOp(mlPrecOp));
00090 }
00091 
00092 Teuchos::RCP<PreconditionerState> MLPreconditionerFactory::buildPreconditionerState() const
00093 {
00094    return Teuchos::rcp(new MLPreconditionerState());
00095 }
00096 
00097 void MLPreconditionerFactory::fillMLPreconditionerState(const BlockedLinearOp & blo,MLPreconditionerState & mlState) const
00098 {
00099    TEUCHOS_ASSERT(not mlState.isFilled());
00100 
00101    EpetraExt::CrsMatrix_SolverMap RowMatrixColMapTrans;
00102 
00103    int rowCnt = blockRowCount(blo);
00104    int colCnt = blockRowCount(blo);
00105 
00106    // construct comm object...add to state
00107    ML_Comm * mlComm;
00108    ML_Comm_Create(&mlComm);
00109    mlState.setMLComm(mlComm);
00110 
00111    ML_Operator * mlBlkMat = ML_Operator_Create(mlComm);
00112    ML_Operator_BlkMatInit(mlBlkMat,mlComm,rowCnt,colCnt,ML_DESTROY_EVERYTHING);
00113 
00114    std::vector<Epetra_RowMatrix *> aggMats;
00115    for(int r=0;r<rowCnt;r++) {
00116       for(int c=0;c<colCnt;c++) {
00117          ML_Operator * tmp = 0;
00118          Epetra_RowMatrix * EMat = 0;
00119          std::stringstream ss;
00120 
00121          ss << "fine_"<< r << "_" << c;
00122 
00123          // extract crs matrix
00124          Teuchos::RCP<Epetra_CrsMatrix> crsMat 
00125             = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(Thyra::get_Epetra_Operator(*blo->getNonconstBlock(r,c)));
00126          EMat = ML_Epetra::ModifyEpetraMatrixColMap(*crsMat,RowMatrixColMapTrans, ss.str().c_str());
00127          if(r==c)  // setup diagonal scaling matrices
00128             aggMats.push_back(EMat);
00129 
00130          // extract ml sub operator
00131          tmp = ML_Operator_Create(mlComm);
00132          ML_Operator_WrapEpetraMatrix(EMat, tmp);
00133 
00134          // add it to the block
00135          ML_Operator_BlkMatInsert(mlBlkMat, tmp, r,c);
00136       }
00137    }
00138    ML_Operator_BlkMatFinalize(mlBlkMat);
00139 
00140    // finish setting up state object
00141    mlState.setMLOperator(mlBlkMat); 
00142 
00143    // now set aggregation matrices
00144    mlState.setAggregationMatrices(aggMats);
00145  
00146    mlState.setIsFilled(true); // register state object as filled
00147 }
00148 
00152 void MLPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList & settings)
00153 {
00154    using Teuchos::RCP;
00155    using Teuchos::rcp;
00156 
00157    Teko_DEBUG_SCOPE("MLPreconditionerFactory::initializeFromParameterList",10);
00158 
00159    // clear initial state
00160    coarseningParams_.clear();
00161    blockRowCount_ = 0;
00162 
00163    blockRowCount_ = settings.get<int>("Block Row Count");
00164 
00165    // read in main parameter list: with smoothing information
00167    mainParams_ = settings.sublist("Smoothing Parameters"); 
00168    mainParams_.set<Teuchos::RCP<const Teko::InverseLibrary> >("smoother: teko inverse library",getInverseLibrary());
00169  
00170    // read in aggregation sub lists
00172    const Teuchos::ParameterList & aggSublist = settings.sublist("Block Aggregation"); 
00173 
00174    for(int block=0;block<blockRowCount_;block++) {
00175       // write out sub list name: "Block #"
00176       std::stringstream ss;
00177       ss << "Block " << block;
00178       std::string sublistName = ss.str();
00179 
00180       // grab sublist
00181       RCP<Teuchos::ParameterList> userSpec = rcp(new Teuchos::ParameterList(aggSublist.sublist(sublistName)));
00182       coarseningParams_.push_back(userSpec);
00183    }
00184 }
00185 
00186 } // end namespace Teko
 All Classes Files Functions Variables