ShyLU  Version of the Day
shylu_iterativesolver_interface_def.hpp
1 #ifndef SHYLU_ITERATIVESOLVER_INTERFACE_DEF_HPP
2 #define SHYLU_ITERATIVESOLVER_INTERFACE_DEF_HPP
3 
4 #include "ShyLUCore_config.h"
5 #include "shylu_iterativesolver_interface_decl.hpp"
6 
7 #include <Teuchos_XMLParameterListHelpers.hpp>
8 
9 
10 
11 namespace ShyLU
12 {
13  template <class Matrix, class Vector>
14  IterativeSolverInterface<Matrix,Vector>::IterativeSolverInterface()
15  {}
16 
17  template <class Matrix, class Vector>
18  IterativeSolverInterface<Matrix,Vector>::IterativeSolverInterface(Matrix *inA, Teuchos::ParameterList *inpList)
19  {
20  A = inA;
21  pList = inpList;
22  }
23 
24  template <class Matrix, class Vector>
25  IterativeSolverInterface<Matrix,Vector>::~IterativeSolverInterface()
26  {}
27 
28  template <class Matrix, class Vector>
29  int IterativeSolverInterface<Matrix,Vector>::init_matrix(Matrix *inA, Teuchos::ParameterList *inpList)
30  {
31  A=inA;
32  pList = inpList;
33  return 0;
34  }
35 
36  template <class Matrix, class Vector>
37  int IterativeSolverInterface<Matrix,Vector>::solveAztec(Vector *b, Vector *x)
38  {
39  printf("**Error** Aztec only supported by Epetra \n");
40  return 1;
41  }
42 
43 
44  template <>
45  int IterativeSolverInterface<Epetra_CrsMatrix,Epetra_MultiVector>::solveAztec(Epetra_MultiVector *b, Epetra_MultiVector *x)
46  {
47 
48  Teuchos::ParameterList subList = pList->sublist("Aztec Input");
49  //Add additional sublist parameters for aztec
50 
51 
52  e_problem.SetOperator(A);
53  e_problem.SetRHS(x);
54  e_problem.SetLHS(b);
55 
56  //set Aztec paramters
57  solver_aztec = new AztecOO(e_problem);
58  solver_aztec->SetAztecOption(AZ_precond, AZ_none);
59  solver_aztec->SetAztecOption(AZ_solver, AZ_gmres);
60  solver_aztec->SetAztecOption(AZ_max_iter, 1000);
61 
62  //solve
63  solver_aztec->Iterate(1000, 1e-5);
64 
65  return 0;
66  }//end solveAztec
67 
68  template <class Matrix, class Vector>
69  int IterativeSolverInterface<Matrix,Vector>::solveBelos(Vector *b, Vector *x)
70  {
71 
72  typedef typename Matrix::scalar_type scalar_type;
73  typedef typename Matrix::local_ordinal_type local_ordinal_type;
74  typedef typename Matrix::global_ordinal_type global_ordinal_type;
75  typedef typename Matrix::node_type node_type;
76  typedef Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type> OP;
77 
78  //Input list of belos
79  //Teuchos::RCP<Teuchos::ParameterList> belosList = Teuchos::parameterList();
80 
81  //Teuchos::RCP<Teuchos::ParameterList> subList = pList->sublist("Belos Input");
82 
83 
84  Teuchos::RCP<Teuchos::ParameterList> pListRCP = Teuchos::rcp(pList,false);
85  Teuchos::RCP<Teuchos::ParameterList> belosList = Teuchos::sublist(pListRCP, "Belos Input");
86 
87  //pListRCP->sublist("Belos Input");
88 
89  //pList->sublist("Belos Input");
90 
91  string solver_name;
92  if(belosList->isParameter("Solver"))
93  {
94  //solver_name = Teuchos::getParameter<string>(&belosList, "Solver");
95  solver_name = belosList->get<string>("Solver");
96  }
97 
98 
99  //printf("start\n");
100  //Add additional sublipst parameter
101  //belosList.set( "Num blocks", 100);
102  // belosList->set( "Block Size" , 10);
103  // belosList->set( "Maximum Iterations", 1000);
104  //belosList->set( "Maximum Restarts", 10);
105  //belosList->set( "Convergence Tolerance" , 1e-5);
106 
107  printf("load belos parameters\n");
108  Teuchos::RCP<Belos::LinearProblem<scalar_type, Vector, OP> > belos_problem;
109  belos_problem = Teuchos::rcp(new Belos::LinearProblem<scalar_type, Vector, OP>());
110 
111  belos_problem->setOperator(Teuchos::rcp(A,false));
112  belos_problem->setRHS(Teuchos::rcp(x,false));
113  belos_problem->setLHS(Teuchos::rcp(b,false));
114  belos_problem->setProblem();
115 
116  //Based on input
117 
118  Belos::SolverFactory<scalar_type, Vector, OP> factory;
119 
120  Teuchos::RCP<Belos::SolverManager<scalar_type,Vector,OP> > solver_belos;
121 
122  solver_belos = factory.create(solver_name, belosList);
123 
124  solver_belos->setProblem(belos_problem);
125 
126 
127  // FIXME (mfh 23 Aug 2015) Not using this return value results in
128  // a build warning. You should use the return value from Belos to
129  // return something other than 0, in case the solve didn't succeed
130  // (didn't return Belos::Converged). I don't know what ShyLU uses
131  // for a "didn't succeed, but not catastrophically so" error code,
132  // or I would put it in myself.
133  //
134  //Belos::ReturnType ret = solver_belos->solve();
135  (void) solver_belos->solve();
136 
137  return 0;
138  }// end solveBelos
139 
140  template <>
141  int IterativeSolverInterface<Epetra_CrsMatrix,Epetra_MultiVector>::solveBelos(Epetra_MultiVector *b, Epetra_MultiVector *x)
142  {
143  // FIXME (mfh 23 Aug 2015) Maybe you should throw an exception or
144  // something to signal that this routine is unimplemented, instead
145  // of just returning 0.
146 
147  //typedef double scalar_type;
148  //typedef Teuchos::ScalarTraits<scalar_type> scalar_type_traits;
149  //typedef scalar_type_traits::magnitudeType magnitude_type; // unused
150  //typedef Epetra_Operator OP; // unused
151 
152  //Input list of belos
153  //Teuchos::RCP<Teuchos::ParameterList> belosList = Teuchos::parameterList();
154 
155  //Teuchos::RCP<Teuchos::ParameterList> subList = pList->sublist("Belos Input");
156 
157  /*
158  Teuchos::RCP<Teuchos::ParameterList> belosList = Teuchos::rcp(pList)->sublist("Belos Input");
159 
160  string solver_name;
161  if(belosList->isParameter("Solver"))
162  {
163  solver_name = belosList->getEntry("Solver");
164  }
165 
166 
167  //printf("start\n");
168  //Add additional sublipst parameter
169  //belosList.set( "Num blocks", 100);
170  // belosList->set( "Block Size" , 10);
171  // belosList->set( "Maximum Iterations", 1000);
172  //belosList->set( "Maximum Restarts", 10);
173  //belosList->set( "Convergence Tolerance" , 1e-5);
174 
175  printf("load belos parameters\n");
176  Teuchos::RCP<Belos::LinearProblem<scalar_type, Vector, OP> > belos_problem;
177  belos_problem = Teuchos::rcp(new Belos::LinearProblem<scalar_type, Vector, OP>());
178 
179  belos_problem->setOperator(Teuchos::rcp(A,false));
180  belos_problem->setRHS(Teuchos::rcp(x,false));
181  belos_problem->setLHS(Teuchos::rcp(b,false));
182  belos_problem->setProblem();
183 
184  //Based on input
185 
186  Belos::SolverFactory<scalar_type, Vector, OP> factory;
187 
188  Teuchos::RCP<Belos::SolverManager<scalar_type,Vector,OP> > solver_belos;
189 
190  solver_belos = factory.create(solver_name, belosList);
191 
192  solver_belos->setProblem(belos_problem);
193 
194 
195  Belos::ReturnType ret = solver_belos->solve();
196 
197  */
198  return 0;
199  }// end solveBelos
200 
201 
202  template <class Matrix, class Vector>
203  int IterativeSolverInterface<Matrix,Vector>::solve(Vector *b, Vector* x)
204  {
205  string solverpackage = Teuchos::getParameter<string>(*pList, "Iterative Solver Package");
206 
207  cout << "found package " << solverpackage << endl;
208 
209  if(solverpackage.compare("Aztec")==0)
210  {
211  printf("**Warning** Aztec is not preferred \n");
212  solveAztec(b,x);
213  }
214  else if(solverpackage.compare("Belos")==0)
215  {
216  solveBelos(b,x);
217  }
218  else
219  {
220  printf("**Error**Iterative Solver Package Not Provided\n");
221  return 1;
222  }
223  return 0;
224  }//end solve()
225 
226 
227 }//end namespace shylu
228 
229 #endif //end ifdef shylu_interative_fie