/[escript]/branches/trilinos_from_5897/trilinoswrap/src/TrilinosMatrixAdapter.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/trilinoswrap/src/TrilinosMatrixAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5913 by caltinay, Tue Feb 9 03:24:30 2016 UTC revision 5914 by caltinay, Wed Feb 10 06:15:12 2016 UTC
# Line 24  Line 24 
24    
25  #include <BelosSolverFactory.hpp>  #include <BelosSolverFactory.hpp>
26  #include <BelosTpetraAdapter.hpp>  #include <BelosTpetraAdapter.hpp>
27    #include <Ifpack2_Factory.hpp>
28  #include <Kokkos_DefaultNode.hpp>  #include <Kokkos_DefaultNode.hpp>
29  #include <MatrixMarket_Tpetra.hpp>  #include <MatrixMarket_Tpetra.hpp>
30  #include <Tpetra_DefaultPlatform.hpp>  #include <Tpetra_DefaultPlatform.hpp>
31  #include <Tpetra_Vector.hpp>  #include <Tpetra_Vector.hpp>
32    #include <MueLu_CreateTpetraPreconditioner.hpp>
33    
34  namespace bp = boost::python;  namespace bp = boost::python;
35  using Teuchos::RCP;  using Teuchos::RCP;
# Line 53  void TrilinosMatrixAdapter::fillComplete Line 55  void TrilinosMatrixAdapter::fillComplete
55      mat->fillComplete(mat->getDomainMap(), mat->getRangeMap(), params);      mat->fillComplete(mat->getDomainMap(), mat->getRangeMap(), params);
56  }  }
57    
58  void TrilinosMatrixAdapter::add(const IndexVector& rowIdx,  void TrilinosMatrixAdapter::add(const std::vector<LO>& rowIdx,
59                                  const std::vector<double>& array)                                  const std::vector<double>& array)
60  {  {
61      const int blockSize = getBlockSize();      const int blockSize = getBlockSize();
# Line 120  void TrilinosMatrixAdapter::ypAx(escript Line 122  void TrilinosMatrixAdapter::ypAx(escript
122      Y->get1dCopy(yView, yView.size());      Y->get1dCopy(yView, yView.size());
123  }  }
124    
125    RCP<SolverType> TrilinosMatrixAdapter::createSolver(
126                                             const escript::SolverBuddy& sb) const
127    {
128        Belos::SolverFactory<ST, VectorType, OpType> factory;
129        RCP<SolverType> solver;
130        RCP<Teuchos::ParameterList> solverParams = Teuchos::parameterList();
131    
132        solverParams->set("Convergence Tolerance", sb.getTolerance());
133        solverParams->set("Maximum Iterations", sb.getIterMax());
134        if (sb.isVerbose()) {
135            solverParams->set("Verbosity", Belos::Errors + Belos::Warnings +
136                    Belos::TimingDetails + Belos::StatusTestDetails);
137        }
138        switch (sb.getSolverMethod()) {
139            case escript::SO_DEFAULT:
140            case escript::SO_METHOD_PCG:
141                solver = factory.create("CG", solverParams);
142                break;
143            case escript::SO_METHOD_GMRES:
144                solver = factory.create("GMRES", solverParams);
145                break;
146            case escript::SO_METHOD_LSQR:
147                solver = factory.create("LSQR", solverParams);
148                break;
149            case escript::SO_METHOD_MINRES:
150                solver = factory.create("MINRES", solverParams);
151                break;
152            case escript::SO_METHOD_TFQMR:
153                solver = factory.create("TFQMR", solverParams);
154                break;
155            default:
156                throw TrilinosAdapterException("Unsupported solver type requested.");
157        }
158        return solver;
159    }
160    
161    RCP<OpType> TrilinosMatrixAdapter::createPreconditioner(
162                                             const escript::SolverBuddy& sb) const
163    {
164        RCP<Teuchos::ParameterList> params = Teuchos::parameterList();
165        Ifpack2::Factory factory;
166        RCP<OpType> prec;
167        RCP<Ifpack2::Preconditioner<ST,LO,GO,NT> > ifprec;
168    
169        switch (sb.getPreconditioner()) {
170            case escript::SO_PRECONDITIONER_NONE:
171                break;
172            case escript::SO_PRECONDITIONER_AMG:
173                {
174                    params->set("max levels", sb.getLevelMax());
175                    params->set("number of equations", getBlockSize());
176                    params->set("cycle type", sb.getCycleType()==1 ? "V" : "W");
177                    params->set("problem: symmetric", sb.isSymmetric());
178                    params->set("verbosity", sb.isVerbose()? "high":"low");
179                    RCP<OpType> A(mat);
180                    prec = MueLu::CreateTpetraPreconditioner(A, *params);
181                }
182                break;
183            case escript::SO_PRECONDITIONER_ILUT:
184                ifprec = factory.create<MatrixType>("ILUT", mat);
185                params->set("fact: drop tolerance", sb.getDropTolerance());
186                break;
187            case escript::SO_PRECONDITIONER_GAUSS_SEIDEL:
188                ifprec = factory.create<MatrixType>("RELAXATION", mat);
189                params->set("relaxation: type", "Gauss-Seidel");
190                params->set("relaxation: sweeps", sb.getNumSweeps());
191                params->set("relaxation: damping factor", sb.getRelaxationFactor());
192                break;
193            case escript::SO_PRECONDITIONER_JACOBI:
194                ifprec = factory.create<MatrixType>("RELAXATION", mat);
195                params->set("relaxation: type", "Jacobi");
196                params->set("relaxation: sweeps", sb.getNumSweeps());
197                params->set("relaxation: damping factor", sb.getRelaxationFactor());
198                break;
199            case escript::SO_PRECONDITIONER_RILU:
200                ifprec = factory.create<MatrixType>("RILUK", mat);
201                params->set("fact: relax value", sb.getRelaxationFactor());
202                break;
203            default:
204                throw TrilinosAdapterException("Unsupported preconditioner requested.");
205        }
206        if (!ifprec.is_null()) {
207            ifprec->setParameters(*params);
208            ifprec->initialize();
209            ifprec->compute();
210            prec = ifprec;
211        }
212        return prec;
213    }
214    
215  void TrilinosMatrixAdapter::setToSolution(escript::Data& out, escript::Data& in,  void TrilinosMatrixAdapter::setToSolution(escript::Data& out, escript::Data& in,
216                                   bp::object& options) const                                   bp::object& options) const
217  {  {
# Line 147  void TrilinosMatrixAdapter::setToSolutio Line 239  void TrilinosMatrixAdapter::setToSolutio
239      RCP<VectorType> X = rcp(new VectorType(mat->getDomainMap(), 1));      RCP<VectorType> X = rcp(new VectorType(mat->getDomainMap(), 1));
240    
241      //solve...      //solve...
242      typedef Tpetra::Operator<ST, LO, GO, MatrixType::node_type> OpType;      RCP<SolverType> solver = createSolver(sb);
243      typedef Belos::LinearProblem<ST, VectorType, OpType> problem_type;      RCP<OpType> prec = createPreconditioner(sb);
244      RCP<Teuchos::ParameterList> solverParams = Teuchos::parameterList();      RCP<ProblemType> problem = rcp(new ProblemType(mat, X, B));
245      solverParams->set("Convergence Tolerance", sb.getTolerance());  
246      solverParams->set("Maximum Iterations", sb.getIterMax());      if (!prec.is_null()) {
247      if (sb.isVerbose()) {          problem->setLeftPrec(prec);
         solverParams->set("Verbosity", Belos::Errors + Belos::Warnings +  
                 Belos::TimingDetails + Belos::StatusTestDetails);  
     }  
     Belos::SolverFactory<ST, VectorType, OpType> factory;  
     RCP<Belos::SolverManager<ST, VectorType, OpType> > solver =  
         factory.create("CG", solverParams);  
     RCP<problem_type> problem = rcp(new problem_type(mat, X, B));  
   
 /*  
     if (sb.getPreconditioner() == escript::SO_PRECONDITIONER_NONE) {  
         if (sb.isVerbose())  
             std::cout << "No preconditioner applied" << std::endl;  
     } else if (sb.getPreconditioner() == escript::SO_PRECONDITIONER_JACOBI) {  
         if (sb.isVerbose())  
             std::cout << "Using Jacobi preconditioner" << std::endl;  
     } else {  
         throw RipleyException("Unsupported preconditioner requested.");  
248      }      }
 */  
     //problem->setRightPrec(M);  
249      problem->setProblem();      problem->setProblem();
250      solver->setProblem(problem);      solver->setProblem(problem);
251      Belos::ReturnType result = solver->solve();      Belos::ReturnType result = solver->solve();

Legend:
Removed from v.5913  
changed lines
  Added in v.5914

  ViewVC Help
Powered by ViewVC 1.1.26