dune-fem  2.8-git
viennacl.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SOLVER_VIENNACL_HH
2 #define DUNE_FEM_SOLVER_VIENNACL_HH
3 
4 #if HAVE_VIENNACL
5 
6 #if HAVE_EIGEN
7 #define VIENNACL_WITH_EIGEN
8 #endif
9 
10 // OpenCL overrules OpenMP
11 #if HAVE_OPENCL
12 //#warning "Using OpneCL"
13 //#define VIENNACL_WITH_OPENCL
14 #elif _OPENMP
15 #error
16 #define VIENNACL_WITH_OPENMP
17 #endif
18 
19 #include <viennacl/linalg/bicgstab.hpp>
20 #include <viennacl/linalg/cg.hpp>
21 #include <viennacl/linalg/gmres.hpp>
22 #include <viennacl/linalg/ilu.hpp>
23 #include <viennacl/compressed_matrix.hpp>
24 #include <viennacl/vector.hpp>
25 
29 
31 
32 namespace Dune
33 {
34 
35  namespace Fem
36  {
37 
38  template< class DiscreteFunction, int method = -1 >
39  class ViennaCLInverseOperator;
40 
41  template< class DiscreteFunction, int method >
42  struct ViennaCLInverseOperatorTraits
43  {
44  typedef DiscreteFunction DiscreteFunctionType;
45  typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
46 
48  typedef OperatorType PreconditionerType;
49 
50  typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction > AssembledOperatorType;
51 
52  typedef ViennaCLInverseOperator< DiscreteFunction, method > InverseOperatorType;
53  typedef SolverParameter SolverParameterType;
54  };
55 
56 
57 
58 
59 
60  template< class DiscreteFunction, int method >
61  class ViennaCLInverseOperator
62  : public InverseOperatorInterface< ViennaCLInverseOperatorTraits< DiscreteFunction, method > >
63  {
64  typedef InverseOperatorInterface< ViennaCLInverseOperatorTraits< DiscreteFunction, method > > BaseType;
65 
66  public:
67  typedef typename BaseType::DomainFunctionType DomainFunction;
68  typedef typename BaseType::RangeFunctionType RangeFunction;
69 
70  typedef typename DomainFunction :: DiscreteFunctionSpaceType DomainSpaceType;
71  typedef typename RangeFunction :: DiscreteFunctionSpaceType RangeSpaceType;
72 
73 #if HAVE_EIGEN
74  typedef Fem::EigenLinearOperator< RangeFunction, DomainFunction > EigenOperatorType;
75 #endif
76 
77  typedef Fem::SparseRowLinearOperator< RangeFunction, DomainFunction > OperatorType;
78 
79  typedef typename BaseType :: SolverDiscreteFunctionType SolverDiscreteFunctionType;
80 
81  private:
82  typedef typename OperatorType::MatrixType Matrix;
83 
84  typedef typename Matrix::field_type Field;
85  typedef typename Matrix::size_type size_type;
86 
87  typedef viennacl::compressed_matrix< Field > ViennaCLMatrix;
88  typedef viennacl::vector< Field > ViennaCLVector;
89 
90  //typedef viennacl::linalg::ilu0_precond< ViennaCLMatrix > Preconditioner;
91  typedef viennacl::linalg::no_precond Preconditioner;
92  public:
93  ViennaCLInverseOperator ( const OperatorType &op, double redEps, double absLimit, unsigned int maxIter, bool verbose )
94  : ViennaCLInverseOperator( redEps, absLimit, maxIter, verbose )
95  {
96  bind( op );
97  }
98 
99  ViennaCLInverseOperator ( const OperatorType &op, double redEps, double absLimit, unsigned int maxIter = std::numeric_limits< unsigned int >::max() )
100  : ViennaCLInverseOperator( redEps, absLimit, maxIter, false )
101  {
102  bind( op );
103  }
104 
105  ViennaCLInverseOperator ( double redEps, double absLimit, unsigned int maxIter = std::numeric_limits< unsigned int >::max(),
106  const bool verbose = false,
107  const SolverParameter& parameter = SolverParameter( Parameter::container() ) )
108  : BaseType( parameter ),
109  absLimit_( absLimit ),
110  method_( method < 0 ? parameter.solverMethod({ SolverParameter::gmres, SolverParameter::cg, SolverParameter::bicgstab }) : method )
111  {
112  }
113 
114  ViennaCLInverseOperator ( const SolverParameter& parameter = SolverParameter( Parameter::container() ) )
115  : ViennaCLInverseOperator( -1,
116  parameter.tolerance(), parameter.maxIterations(),
117  parameter.verbose(), parameter )
118  {
119  assert( parameter.errorMeasure() == 0 );
120  }
121 
122  void bind( const OperatorType& op )
123  {
124  BaseType::bind( op );
125 #if HAVE_EIGEN
126  eigenMatrix_ = dynamic_cast<const EigenOperatorType* > (&op);
127  if( eigenMatrix_ )
128  {
129  const auto& matrix = eigenMatrix_->matrix();
130  gupMatrix_.resize( matrix.rows(), matrix.cols() );
131  viennacl::copy( matrix.data(), gpuMatrix_ );
132  }
133  else
134 #endif // HAVE_EIGEN
135  if( assembledOperator_ )
136  {
137  const auto& matrix = assembledOperator_->matrix();
138  gpuMatrix_.resize( matrix.rows(), matrix.cols() );
139  std::vector< std::map< size_type, Field > > cpuMatrix;
140  matrix.fillCSRStorage( cpuMatrix );
141  viennacl::copy( cpuMatrix, gpuMatrix_ );
142  }
143  }
144 
145  void unbind()
146  {
147  BaseType::unbind();
148  gpuMatrix_ = ViennaCLMatrix();
149 #if HAVE_EIGEN
150  eigenOperator_ = nullptr;
151 #endif
152  }
153 
154  int apply( const SolverDiscreteFunctionType& u, SolverDiscreteFunctionType& w ) const
155  {
156  viennacl::vector< Field > vclU( u.size() ), vclW( w.size() );
157  viennacl::copy( u.dbegin(), u.dend(), vclU.begin() );
158 
159  int maxIterations = parameter().maxIterations();
160  int iterations = -1;
161 
162  //std::cout << "Using ViennaCL " << std::endl;
163 
164  if( method_ == SolverParameter::cg )
165  {
166  Preconditioner ilu0; // ( gpuMatrix_, viennacl::linalg::ilu0_tag() );
167  viennacl::linalg::cg_tag tag( absLimit_, maxIterations );
168  vclW = viennacl::linalg::solve( gpuMatrix_, vclU, tag, ilu0 );
169  iterations = tag.iters();
170  }
171  else if ( method_ == SolverParameter::bicgstab )
172  {
173  /*
174  // configuration of preconditioner:
175  viennacl::linalg::chow_patel_tag chow_patel_ilu_config;
176  chow_patel_ilu_config.sweeps(3); // three nonlinear sweeps
177  chow_patel_ilu_config.jacobi_iters(2); // two Jacobi iterations per triangular 'solve' Rx=r
178  // create and compute preconditioner:
179  typedef viennacl::linalg::chow_patel_ilu_precond< ViennaCLMatrix > Preconditioner;
180  //viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix<ScalarType> > chow_patel_ilu(A, chow_patel_ilu_config);
181 
182  typedef viennacl::linalg::block_ilu_precond< ViennaCLMatrix > Preconditioner;
183  Preconditioner ilu0( gpuMatrix_, chow_patel_ilu_config );
184  */
185 
186  //typedef viennacl::linalg::block_ilu_precond< ViennaCLMatrix, viennacl::linalg::ilu0_tag > Preconditioner;
187  Preconditioner ilu0; // ( gpuMatrix_, viennacl::linalg::ilu0_tag() );
188  viennacl::linalg::bicgstab_tag tag( absLimit_, maxIterations );
189  vclW = viennacl::linalg::solve( gpuMatrix_, vclU, tag, ilu0 );
190  iterations = tag.iters();
191  }
192  else if ( method_ == SolverParameter::gmres )
193  {
194  Preconditioner ilu0; // ( gpuMatrix_, viennacl::linalg::ilu0_tag() );
195  viennacl::linalg::gmres_tag tag( absLimit_, maxIterations );
196  vclW = viennacl::linalg::solve( gpuMatrix_, vclU, tag, ilu0 );
197  iterations = tag.iters();
198  }
199  else
200  {
201  DUNE_THROW(NotImplemented,"ViennaCL does not support this solver");
202  }
203 
204  viennacl::copy( vclW.begin(), vclW.end(), w.dbegin() );
205  return iterations ;
206  }
207 
208  protected:
209  ViennaCLMatrix gpuMatrix_;
210 
211  using BaseType :: assembledOperator_;
212  using BaseType :: iterations_;
213  using BaseType :: parameter;
214 #if HAVE_EIGEN
215  const EigenOperatorType* eigenOperator_ = nullptr;
216 #endif
217 
218  double absLimit_;
219 
220 
221  int method_;
222  };
223 
224 
225  // ViennaCLBiCGStabInverseOperator
226  //-----------------------------------
227 
228  template< class DiscreteFunction >
229  using ViennaCLCGInverseOperator = ViennaCLInverseOperator< DiscreteFunction, SolverParameter :: cg >;
230 
231 
232  // ViennaCLBiCGStabInverseOperator
233  //-----------------------------------
234 
235  template< class DiscreteFunction >
236  using ViennaCLBiCGStabInverseOperator = ViennaCLInverseOperator< DiscreteFunction, SolverParameter :: bicgstab >;
237 
238 
239  // ViennaCLGMResInverseOperator
240  //-----------------------------------
241 
242  template< class DiscreteFunction >
243  using ViennaCLGMResInverseOperator = ViennaCLInverseOperator< DiscreteFunction, SolverParameter :: gmres >;
244 
245  } // namespace Fem
246 
247 } // namespace Dune
248 
249 #endif
250 
251 #endif // #ifndef DUNE_FEM_SOLVER_VIENNACL_HH
Definition: bindguard.hh:11
static constexpr T max(T a)
Definition: utility.hh:77
static ParameterContainer & container()
Definition: io/parameter.hh:193
static const int cg
Definition: solver/parameter.hh:24
static const int gmres
Definition: solver/parameter.hh:26
static const int bicgstab
Definition: solver/parameter.hh:25