Teko Version of the Day
Teko_mlutils.cpp
00001 #include "Teko_mlutils.hpp"
00002 
00003 #include <vector>
00004 
00005 #ifdef HAVE_MPI
00006    #include "mpi.h"
00007    #include "Epetra_MpiComm.h"
00008 #else
00009    #include "Epetra_SerialComm.h"
00010 #endif
00011 #include "Epetra_Vector.h"
00012 
00013 #include "ml_epetra_utils.h"
00014 #include "ml_op_utils.h"
00015 
00016 #include "Thyra_EpetraLinearOp.hpp"
00017 
00018 #include "EpetraExt_RowMatrixOut.h"
00019 
00020 #include "Teuchos_XMLParameterListHelpers.hpp"
00021 
00022 #include "Teko_InverseFactory.hpp"
00023 #include "Teko_InverseLibrary.hpp"
00024 #include "Teko_Utilities.hpp"
00025 #include "Teko_PreconditionerFactory.hpp"
00026 #include "Teko_EpetraOperatorWrapper.hpp"
00027 #include "Teko_SmootherPreconditionerFactory.hpp"
00028 
00029 namespace Teko {
00030 namespace mlutils {
00031 
00032 // build a very simple row map from the ML_Operator
00033 Teuchos::RCP<Epetra_Map> buildRowMap(ML_Operator * mlOp)
00034 {
00035    ML_Comm *comm = mlOp->comm;
00036 #ifdef ML_MPI
00037    MPI_Comm mpi_comm;
00038    mpi_comm = comm->USR_comm;
00039    Epetra_MpiComm eComm( mpi_comm ) ;
00040 #else
00041    Epetra_SerialComm eComm ;
00042 #endif
00043 
00044    // get operators row count, build map...who cares what GIDs are
00045    int rowCnt = mlOp->outvec_leng;
00046    return Teuchos::rcp(new Epetra_Map(-1,rowCnt,0,eComm));
00047 }
00048 
00052 Teuchos::RCP<Epetra_CrsMatrix> convertToCrsMatrix(ML_Operator * mlOp,
00053                                                   const Teuchos::RCP<Epetra_Map> & rowMapArg)
00054 {
00055    ML_Comm *comm = mlOp->comm;
00056 #ifdef ML_MPI
00057    MPI_Comm mpi_comm;
00058    mpi_comm = comm->USR_comm;
00059    Epetra_MpiComm eComm( mpi_comm ) ;
00060 #else
00061    Epetra_SerialComm eComm ;
00062 #endif
00063 
00064    // build a default row map
00065    Teuchos::RCP<Epetra_Map> rowMap = rowMapArg;
00066    if(rowMapArg==Teuchos::null) 
00067       rowMap = buildRowMap(mlOp);
00068 
00069    // build lightweight CrsMatrix wrapper 
00070    Epetra_CrsMatrix * crsMatrixPtr = 0;
00071    int maxNumNonzeros = 0;
00072    double cpuTime = 0.0;
00073    bool verbose = false;
00074    ML_Operator2EpetraCrsMatrix(mlOp,crsMatrixPtr,maxNumNonzeros,false,cpuTime,verbose);
00075 
00076    return Teuchos::rcp(crsMatrixPtr);
00077 }
00078 
00079 Teko::LinearOp buildTekoBlockOp(ML_Operator * mlOp,int level)
00080 {
00081    Teko_DEBUG_SCOPE("Teko::mlutils::buildTekoBlockOp",0);
00082   
00083    using Teuchos::RCP;
00084 
00085    int numRows = ML_Operator_BlkMatNumBlockRows(mlOp);  
00086    int numCols = ML_Operator_BlkMatNumBlockCols(mlOp);  
00087 
00088    Teko::BlockedLinearOp tekoOp = Teko::createBlockedOp();
00089    Teko::beginBlockFill(tekoOp,numRows,numCols);
00090       for(int i=0;i<numRows;i++) {
00091          for(int j=0;j<numCols;j++) {
00092             ML_Operator * subBlock = ML_Operator_BlkMatExtract(mlOp,i,j);
00093  
00094             // if required construct and add a block to tekoOp
00095             if(subBlock!=0) {
00096                // create a CRS matrix from ML operator
00097                RCP<const Epetra_CrsMatrix> eCrsMat = convertToCrsMatrix(subBlock);
00098 
00099                #if 0
00100                std::stringstream ss;
00101                ss << "A(" << level << ")_" << i << "_" << j << ".mm";
00102                EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*eCrsMat);
00103                #endif
00104  
00105                // build Teko operator, add to Teko blocked operator
00106                Teko::LinearOp tekoSubBlock = Thyra::epetraLinearOp(eCrsMat); 
00107                Teko::setBlock(i,j,tekoOp,tekoSubBlock);
00108             }
00109          }
00110       }
00111    Teko::endBlockFill(tekoOp);
00112 
00113    return tekoOp.getConst();
00114 }
00115 
00116 int smoother(ML_Smoother *mydata, int leng1, double x[], int leng2, 
00117              double rhs[]) 
00118 {
00119    Teko_DEBUG_SCOPE("Teko::mlutils::smoother",10);
00120 
00121    // std::cout << "init guess = " << mydata->init_guess << std::endl;
00122 
00123    // grab data object
00124    SmootherData * smootherData = (struct SmootherData *) ML_Get_MySmootherData(mydata);
00125 
00126    Epetra_Vector X(View, smootherData->Amat->OperatorDomainMap(), x);
00127    Epetra_Vector Y(View, smootherData->Amat->OperatorRangeMap(),rhs);
00128 
00129    smootherData->smootherOperator->Apply(Y,X);
00130 
00131    return 0;
00132 }
00133 
00134 Teuchos::RCP<Teko::InverseFactory> ML_Gen_Init_Teko_Prec(const std::string smoothName,const Teko::InverseLibrary & invLib)
00135 {
00136    Teko_DEBUG_SCOPE("Teko::mlutils::ML_Gen_Init_Teko_Prec",10);
00137 
00138    // Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler());
00139    // Teuchos::RCP<Teko::InverseLibrary> invLib
00140    //       = Teko::InverseLibrary::buildFromParameterList(pl);
00141    // invLib->setRequestHandler(rh);
00142 
00143    Teuchos::RCP<Teko::PreconditionerFactory> precFact;
00144    Teuchos::RCP<Teko::InverseFactory> invFact = invLib.getInverseFactory(smoothName);
00145 
00146    return invFact;
00147 }
00148 
00149 extern "C" 
00150 void ML_Destroy_Smoother_Teko(void * data) 
00151 {
00152    Teko::mlutils::SmootherData * sData = (Teko::mlutils::SmootherData*) data;
00153    delete sData;
00154 }
00155 
00156 extern "C" 
00157 int ML_Gen_Smoother_Teko(ML *ml, int level, int pre_or_post, int ntimes,const Teuchos::RCP<const Teuchos::ParameterList> & tekoPL,
00158                          const Teuchos::RCP<const Teko::InverseLibrary> & invLib_in,const std::string & inverse, bool isBlocked)
00159 {
00160    Teko_DEBUG_SCOPE("Teko::mlutils::ML_Gen_Smoother_Teko",10);
00161    ML_Operator * BlkMat = &(ml->Amat[level]);
00162 
00163    Teuchos::RCP<const Teko::InverseLibrary> invLib;
00164    if(invLib_in ==Teuchos::null) {
00165       Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler());
00166       Teuchos::RCP<Teko::InverseLibrary> ncInvLib = Teko::InverseLibrary::buildFromParameterList(*tekoPL);
00167       ncInvLib->setRequestHandler(rh);
00168 
00169       invLib = ncInvLib.getConst();
00170    }
00171    else {
00172       // this assumes a request handler has already been set
00173       invLib = invLib_in;
00174    }
00175 
00176    Teko::LinearOp tekoAmat;
00177    if(isBlocked) {
00178       // build teko blocked operator
00179       tekoAmat = Teko::mlutils::buildTekoBlockOp(BlkMat,level);
00180    }
00181    else {
00182       // build unblocked operator
00183       Teuchos::RCP<const Epetra_CrsMatrix> eCrsMat = convertToCrsMatrix(BlkMat);
00184       tekoAmat = Thyra::epetraLinearOp(eCrsMat); 
00185    }
00186 
00187    // Teuchos::RCP<Teko::mlutils::SmootherData> data = Teuchos::rcp(new Teko::mlutils::SmootherData);
00188    Teko::mlutils::SmootherData * data = new Teko::mlutils::SmootherData;
00189    data->Amat = Teuchos::rcp(new Teko::Epetra::EpetraOperatorWrapper(tekoAmat));
00190 
00191    // Teuchos::ParameterList pl = *Teuchos::getParametersFromXmlFile(filename);
00192    Teuchos::RCP<Teko::InverseFactory> invFact = ML_Gen_Init_Teko_Prec(inverse,*invLib);
00193    Teuchos::RCP<Teko::RequestHandler> rh = invFact->getRequestHandler();
00194 
00195    // build smoother operator
00196    Teko::LinearOp precOp = Teko::buildInverse(*invFact,tekoAmat);
00197    Teko::LinearOp smootherOp = Teko::buildSmootherLinearOp(tekoAmat,precOp,ntimes,true);
00198 
00199    // get an epetra operator wrapper
00200    data->smootherOperator = Teuchos::rcp(new Teko::Epetra::EpetraOperatorWrapper(smootherOp));
00201 
00202    int ret_val = ML_Set_Smoother(ml,level, pre_or_post, (void*) data, Teko::mlutils::smoother, NULL);
00203    ml->post_smoother[level].data_destroy = ML_Destroy_Smoother_Teko;
00204  
00205    return ret_val;
00206 } 
00207 
00208 }
00209 }
 All Classes Files Functions Variables