Teko_mlutils.cpp

#include "Teko_mlutils.hpp"

#include <vector>

#ifdef HAVE_MPI
   #include "mpi.h"
   #include "Epetra_MpiComm.h"
#else
   #include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"

#include "ml_epetra_utils.h"
#include "ml_op_utils.h"

#include "Thyra_EpetraLinearOp.hpp"

#include "EpetraExt_RowMatrixOut.h"

#include "Teuchos_XMLParameterListHelpers.hpp"

#include "Teko_InverseFactory.hpp"
#include "Teko_InverseLibrary.hpp"
#include "Teko_Utilities.hpp"
#include "Teko_PreconditionerFactory.hpp"
#include "Teko_EpetraOperatorWrapper.hpp"
#include "Teko_SmootherPreconditionerFactory.hpp"

namespace Teko {
namespace mlutils {

// build a very simple row map from the ML_Operator
Teuchos::RCP<Epetra_Map> buildRowMap(ML_Operator * mlOp)
{
   ML_Comm *comm = mlOp->comm;
#ifdef ML_MPI
   MPI_Comm mpi_comm;
   mpi_comm = comm->USR_comm;
   Epetra_MpiComm eComm( mpi_comm ) ;
#else
   Epetra_SerialComm eComm ;
#endif

   // get operators row count, build map...who cares what GIDs are
   int rowCnt = mlOp->outvec_leng;
   return Teuchos::rcp(new Epetra_Map(-1,rowCnt,0,eComm));
}

Teuchos::RCP<Epetra_CrsMatrix> convertToCrsMatrix(ML_Operator * mlOp,
                                                  const Teuchos::RCP<Epetra_Map> & rowMapArg)
{
   ML_Comm *comm = mlOp->comm;
#ifdef ML_MPI
   MPI_Comm mpi_comm;
   mpi_comm = comm->USR_comm;
   Epetra_MpiComm eComm( mpi_comm ) ;
#else
   Epetra_SerialComm eComm ;
#endif

   // build a default row map
   Teuchos::RCP<Epetra_Map> rowMap = rowMapArg;
   if(rowMapArg==Teuchos::null) 
      rowMap = buildRowMap(mlOp);

   // build lightweight CrsMatrix wrapper 
   Epetra_CrsMatrix * crsMatrixPtr = 0;
   int maxNumNonzeros = 0;
   double cpuTime = 0.0;
   bool verbose = false;
   ML_Operator2EpetraCrsMatrix(mlOp,crsMatrixPtr,maxNumNonzeros,false,cpuTime,verbose);

   return Teuchos::rcp(crsMatrixPtr);
}

Teko::LinearOp buildTekoBlockOp(ML_Operator * mlOp,int level)
{
   Teko_DEBUG_SCOPE("Teko::mlutils::buildTekoBlockOp",0);
  
   using Teuchos::RCP;

   int numRows = ML_Operator_BlkMatNumBlockRows(mlOp);  
   int numCols = ML_Operator_BlkMatNumBlockCols(mlOp);  

   Teko::BlockedLinearOp tekoOp = Teko::createBlockedOp();
   Teko::beginBlockFill(tekoOp,numRows,numCols);
      for(int i=0;i<numRows;i++) {
         for(int j=0;j<numCols;j++) {
            ML_Operator * subBlock = ML_Operator_BlkMatExtract(mlOp,i,j);
 
            // if required construct and add a block to tekoOp
            if(subBlock!=0) {
               // create a CRS matrix from ML operator
               RCP<const Epetra_CrsMatrix> eCrsMat = convertToCrsMatrix(subBlock);

               #if 0
               std::stringstream ss;
               ss << "A(" << level << ")_" << i << "_" << j << ".mm";
               EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*eCrsMat);
               #endif
 
               // build Teko operator, add to Teko blocked operator
               Teko::LinearOp tekoSubBlock = Thyra::epetraLinearOp(eCrsMat); 
               Teko::setBlock(i,j,tekoOp,tekoSubBlock);
            }
         }
      }
   Teko::endBlockFill(tekoOp);

   return tekoOp.getConst();
}

int smoother(ML_Smoother *mydata, int leng1, double x[], int leng2, 
             double rhs[]) 
{
   Teko_DEBUG_SCOPE("Teko::mlutils::smoother",10);

   // std::cout << "init guess = " << mydata->init_guess << std::endl;

   // grab data object
   SmootherData * smootherData = (struct SmootherData *) ML_Get_MySmootherData(mydata);

   Epetra_Vector X(View, smootherData->Amat->OperatorDomainMap(), x);
   Epetra_Vector Y(View, smootherData->Amat->OperatorRangeMap(),rhs);

   smootherData->smootherOperator->Apply(Y,X);

   return 0;
}

Teuchos::RCP<Teko::InverseFactory> ML_Gen_Init_Teko_Prec(const std::string smoothName,const Teko::InverseLibrary & invLib)
{
   Teko_DEBUG_SCOPE("Teko::mlutils::ML_Gen_Init_Teko_Prec",10);

   // Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler());
   // Teuchos::RCP<Teko::InverseLibrary> invLib
   //       = Teko::InverseLibrary::buildFromParameterList(pl);
   // invLib->setRequestHandler(rh);

   Teuchos::RCP<Teko::PreconditionerFactory> precFact;
   Teuchos::RCP<Teko::InverseFactory> invFact = invLib.getInverseFactory(smoothName);

   return invFact;
}

extern "C" 
void ML_Destroy_Smoother_Teko(void * data) 
{
   Teko::mlutils::SmootherData * sData = (Teko::mlutils::SmootherData*) data;
   delete sData;
}

extern "C" 
int ML_Gen_Smoother_Teko(ML *ml, int level, int pre_or_post, int ntimes,const Teuchos::RCP<const Teuchos::ParameterList> & tekoPL,
                         const Teuchos::RCP<const Teko::InverseLibrary> & invLib_in,const std::string & inverse, bool isBlocked)
{
   Teko_DEBUG_SCOPE("Teko::mlutils::ML_Gen_Smoother_Teko",10);
   ML_Operator * BlkMat = &(ml->Amat[level]);

   Teuchos::RCP<const Teko::InverseLibrary> invLib;
   if(invLib_in ==Teuchos::null) {
      Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler());
      Teuchos::RCP<Teko::InverseLibrary> ncInvLib = Teko::InverseLibrary::buildFromParameterList(*tekoPL);
      ncInvLib->setRequestHandler(rh);

      invLib = ncInvLib.getConst();
   }
   else {
      // this assumes a request handler has already been set
      invLib = invLib_in;
   }

   Teko::LinearOp tekoAmat;
   if(isBlocked) {
      // build teko blocked operator
      tekoAmat = Teko::mlutils::buildTekoBlockOp(BlkMat,level);
   }
   else {
      // build unblocked operator
      Teuchos::RCP<const Epetra_CrsMatrix> eCrsMat = convertToCrsMatrix(BlkMat);
      tekoAmat = Thyra::epetraLinearOp(eCrsMat); 
   }

   // Teuchos::RCP<Teko::mlutils::SmootherData> data = Teuchos::rcp(new Teko::mlutils::SmootherData);
   Teko::mlutils::SmootherData * data = new Teko::mlutils::SmootherData;
   data->Amat = Teuchos::rcp(new Teko::Epetra::EpetraOperatorWrapper(tekoAmat));

   // Teuchos::ParameterList pl = *Teuchos::getParametersFromXmlFile(filename);
   Teuchos::RCP<Teko::InverseFactory> invFact = ML_Gen_Init_Teko_Prec(inverse,*invLib);
   Teuchos::RCP<Teko::RequestHandler> rh = invFact->getRequestHandler();

   // build smoother operator
   Teko::LinearOp precOp = Teko::buildInverse(*invFact,tekoAmat);
   Teko::LinearOp smootherOp = Teko::buildSmootherLinearOp(tekoAmat,precOp,ntimes,true);

   // get an epetra operator wrapper
   data->smootherOperator = Teuchos::rcp(new Teko::Epetra::EpetraOperatorWrapper(smootherOp));

   int ret_val = ML_Set_Smoother(ml,level, pre_or_post, (void*) data, Teko::mlutils::smoother, NULL);
   ml->post_smoother[level].data_destroy = ML_Destroy_Smoother_Teko;
 
   return ret_val;
} 

}
}