dune-fem  2.8-git
umfpacksolver.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_UMFPACKSOLVER_HH
2 #define DUNE_FEM_UMFPACKSOLVER_HH
3 
4 #include <algorithm>
5 #include <iostream>
6 #include <tuple>
7 #include <utility>
8 #include <vector>
9 
10 #include <dune/common/hybridutilities.hh>
14 #include <dune/fem/io/parameter.hh>
20 
21 #if HAVE_SUITESPARSE_UMFPACK
22 #include <dune/fem/misc/umfpack.hh>
23 
24 namespace Dune
25 {
26 namespace Fem
27 {
28 
38 template<class DiscreteFunction, class Matrix>
39 class UMFPACKInverseOperator;
40 
41 template<class DiscreteFunction, class Matrix>
42 struct UMFPACKInverseOperatorTraits
43 {
44  typedef DiscreteFunction DiscreteFunctionType;
45  typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
46 
48  typedef OperatorType PreconditionerType;
49 
50  typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction, Matrix > AssembledOperatorType;
51 
52  typedef UMFPACKInverseOperator< DiscreteFunction,Matrix> InverseOperatorType;
53  typedef SolverParameter SolverParameterType;
54 };
55 
63 template<class DiscreteFunction,
64  class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
65 class UMFPACKInverseOperator :
66  public InverseOperatorInterface< UMFPACKInverseOperatorTraits< DiscreteFunction, Matrix > >
67 {
68 public:
69  typedef UMFPACKInverseOperatorTraits< DiscreteFunction, Matrix > Traits;
70  typedef InverseOperatorInterface< Traits > BaseType;
71 
72  typedef typename BaseType :: SolverDiscreteFunctionType
73  SolverDiscreteFunctionType;
74 
75  typedef typename BaseType :: OperatorType OperatorType;
76  typedef typename BaseType :: AssembledOperatorType AssembledOperatorType;
77 
78  typedef UMFPACKInverseOperator< DiscreteFunction,Matrix> ThisType;
79 
80  typedef DiscreteFunction DiscreteFunctionType;
81 
82  // \brief The column-compressed matrix type.
83  typedef ColCompMatrix<typename AssembledOperatorType::MatrixType::MatrixBaseType,long int> CCSMatrixType;
84  typedef typename DiscreteFunctionType::DofType DofType;
85  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
86 
87  using BaseType :: parameter_;
88 
92  UMFPACKInverseOperator(const SolverParameter &parameter = SolverParameter(Parameter::container()) )
93  : BaseType(parameter),
94  ccsmat_(),
95  UMF_Symbolic( nullptr ),
96  UMF_Numeric( nullptr )
97  {
98  Caller::defaults(UMF_Control);
99  UMF_Control[UMFPACK_PRL] = 4;
100  }
101 
102  // \brief Destructor.
103  ~UMFPACKInverseOperator()
104  {
105  unbind();
106  }
107 
108  void bind( const OperatorType& op )
109  {
110  // clear old storage
111  finalize();
112  BaseType::bind( op );
113  assert( assembledOperator_ );
114  }
115 
116  void unbind ()
117  {
118  finalize();
119  BaseType::unbind();
120  }
121 
122  // \brief Decompose matrix.
123  template<typename... A>
124  void prepare(A... ) const
125  {
126  if(assembledOperator_ && !ccsmat_)
127  {
128  ccsmat_ = std::make_unique<CCSMatrixType>(assembledOperator_->exportMatrix());
129  decompose();
130  }
131  }
132 
133  // \brief Free allocated memory.
134  virtual void finalize()
135  {
136  if( ccsmat_ )
137  {
138  getCCSMatrix().free();
139  ccsmat_.reset();
140  Caller::free_symbolic(&UMF_Symbolic); UMF_Symbolic = nullptr;
141  Caller::free_numeric(&UMF_Numeric); UMF_Numeric = nullptr;
142  }
143  }
144 
151  int apply(const DofType* arg, DofType* dest) const
152  {
153  prepare();
154 
155  double UMF_Apply_Info[UMFPACK_INFO];
156  Caller::solve(UMFPACK_A, getCCSMatrix().getColStart(), getCCSMatrix().getRowIndex(), getCCSMatrix().getValues(),
157  dest, const_cast<DofType*>(arg), UMF_Numeric, UMF_Control, UMF_Apply_Info);
158  if( Parameter::verbose() && parameter_->verbose() )
159  {
160  Caller::report_status(UMF_Control, UMF_Apply_Info[UMFPACK_STATUS]);
161  std::cout <<"[UMFPack Solve]" << std::endl;
162  std::cout << "Wallclock Time: " << UMF_Apply_Info[UMFPACK_SOLVE_WALLTIME]
163  << " (CPU Time: " << UMF_Apply_Info[UMFPACK_SOLVE_TIME] << ")" << std::endl;
164  std::cout << "Flops Taken: " << UMF_Apply_Info[UMFPACK_SOLVE_FLOPS] << std::endl;
165  std::cout << "Iterative Refinement steps taken: " << UMF_Apply_Info[UMFPACK_IR_TAKEN] << std::endl;
166  std::cout << "Error Estimate: " << UMF_Apply_Info[UMFPACK_OMEGA1] << " resp. " << UMF_Apply_Info[UMFPACK_OMEGA2] << std::endl;
167  }
168 
169  const_cast<ThisType*>(this)->finalize();
170 
171  return 1;
172  }
173 
180  int apply(const SolverDiscreteFunctionType& arg, SolverDiscreteFunctionType& dest) const
181  {
182  return apply(arg.leakPointer(), dest.leakPointer());
183  }
184 
191  template<typename... DFs>
192  int apply(const TupleDiscreteFunction<DFs...>& arg,TupleDiscreteFunction<DFs...>& dest) const
193  {
194  // copy DOF's arg into a consecutive vector
195  std::vector<DofType> vecArg(arg.size());
196  auto vecArgIt(vecArg.begin());
197  Hybrid::forEach(std::make_index_sequence<sizeof...(DFs)>{},
198  [&](auto i){vecArgIt=std::copy(std::get<i>(arg).dbegin(),std::get<i>(arg).dend(),vecArgIt);});
199  std::vector<DofType> vecDest(dest.size());
200  // apply operator
201  apply(vecArg.data(),vecDest.data());
202  // copy back solution into dest
203  auto vecDestIt(vecDest.begin());
204  Hybrid::forEach(std::make_index_sequence<sizeof...(DFs)>{},[&](auto i){for(auto& dof:dofs(std::get<i>(dest))) dof=(*(vecDestIt++));});
205  return 1;
206  }
207 
208  void printTexInfo(std::ostream& out) const
209  {
210  out<<"Solver: UMFPACK direct solver"<<std::endl;
211  }
212 
213  void setMaxIterations ( int ) {}
214 
218  CCSMatrixType& getCCSMatrix() const
219  {
220  assert( ccsmat_ );
221  return *ccsmat_;
222  }
223 
224  // \brief Print some statistics about the UMFPACK decomposition.
225  void printDecompositionInfo() const
226  {
227  Caller::report_info(UMF_Control,UMF_Decomposition_Info);
228  }
229 
230  UMFPACKInverseOperator(const UMFPACKInverseOperator &other)
231  : BaseType(other),
232  ccsmat_(),
233  UMF_Symbolic(other.UMF_Symbolic),
234  UMF_Numeric(other.UMF_Numeric)
235  {
236  for (int i=0;i<UMFPACK_CONTROL;++i) UMF_Control[i] = other.UMF_Control[i];
237  for (int i=0;i<UMFPACK_INFO;++i) UMF_Decomposition_Info[i] = other.UMF_Decomposition_Info[i];
238  }
239 
240 protected:
241  using BaseType::assembledOperator_;
242  mutable std::unique_ptr<CCSMatrixType> ccsmat_;
243  mutable void *UMF_Symbolic;
244  mutable void *UMF_Numeric;
245  mutable double UMF_Control[UMFPACK_CONTROL];
246  mutable double UMF_Decomposition_Info[UMFPACK_INFO];
247 
248  typedef typename Dune::UMFPackMethodChooser<DofType> Caller;
249 
250  // \brief Computes the UMFPACK decomposition.
251  void decompose() const
252  {
253  const std::size_t dimMat(getCCSMatrix().N());
254  Caller::symbolic(static_cast<int>(dimMat), static_cast<int>(dimMat), getCCSMatrix().getColStart(), getCCSMatrix().getRowIndex(),
255  reinterpret_cast<double*>(getCCSMatrix().getValues()), &UMF_Symbolic, UMF_Control, UMF_Decomposition_Info);
256  Caller::numeric(getCCSMatrix().getColStart(), getCCSMatrix().getRowIndex(), reinterpret_cast<double*>(getCCSMatrix().getValues()),
257  UMF_Symbolic, &UMF_Numeric, UMF_Control, UMF_Decomposition_Info);
258  if( Parameter::verbose() && parameter_->verbose() )
259  {
260  Caller::report_status(UMF_Control,UMF_Decomposition_Info[UMFPACK_STATUS]);
261  std::cout << "[UMFPack Decomposition]" << std::endl;
262  std::cout << "Wallclock Time taken: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_WALLTIME]
263  << " (CPU Time: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_TIME] << ")" << std::endl;
264  std::cout << "Flops taken: " << UMF_Decomposition_Info[UMFPACK_FLOPS] << std::endl;
265  std::cout << "Peak Memory Usage: " << UMF_Decomposition_Info[UMFPACK_PEAK_MEMORY]*UMF_Decomposition_Info[UMFPACK_SIZE_OF_UNIT]
266  << " bytes" << std::endl;
267  std::cout << "Condition number estimate: " << 1./UMF_Decomposition_Info[UMFPACK_RCOND] << std::endl;
268  std::cout << "Numbers of non-zeroes in decomposition: L: " << UMF_Decomposition_Info[UMFPACK_LNZ]
269  << " U: " << UMF_Decomposition_Info[UMFPACK_UNZ] << std::endl;
270  }
271  }
272 };
273 
274 }
275 }
276 
277 #endif // #if HAVE_SUITESPARSE_UMFPACK
278 
279 #endif // #ifndef DUNE_FEM_UMFPACKSOLVER_HH
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
static ParameterContainer & container()
Definition: io/parameter.hh:193