dune-fem  2.8-git
spqrsolver.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SPQRSOLVER_HH
2 #define DUNE_FEM_SPQRSOLVER_HH
3 
4 #include <algorithm>
5 #include <iostream>
6 #include <tuple>
7 #include <vector>
8 #include <utility>
9 
10 #include <dune/common/hybridutilities.hh>
14 #include <dune/fem/io/parameter.hh>
20 
21 #if HAVE_SUITESPARSE_SPQR
22 #include <SuiteSparseQR.hpp>
23 
24 namespace Dune
25 {
26 namespace Fem
27 {
28 
46 template<class DiscreteFunction, bool symmetric=false,
47  class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
48 class SPQRInverseOperator;
49 
50 template<class DiscreteFunction, bool symmetric, class Matrix>
51 struct SPQRInverseOperatorTraits
52 {
53  typedef DiscreteFunction DiscreteFunctionType;
54  typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
55 
57  typedef OperatorType PreconditionerType;
58 
59  typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction, Matrix > AssembledOperatorType;
60  typedef ColCompMatrix<typename AssembledOperatorType::MatrixType::MatrixBaseType,int > CCSMatrixType;
61 
62  typedef SPQRInverseOperator< DiscreteFunction, symmetric, Matrix > InverseOperatorType;
63  typedef SolverParameter SolverParameterType;
64 };
65 
66 
67 
68 template< class DF, bool symmetric, class Matrix >
69 class SPQRInverseOperator : public InverseOperatorInterface< SPQRInverseOperatorTraits< DF, symmetric, Matrix > >
70 {
71  typedef SPQRInverseOperatorTraits< DF, symmetric, Matrix > Traits;
72  typedef InverseOperatorInterface< Traits > BaseType;
73  typedef SPQRInverseOperator< DF, symmetric, Matrix > ThisType;
74 public:
75  typedef DF DiscreteFunctionType;
76  typedef typename BaseType :: OperatorType OperatorType;
77  typedef typename BaseType :: SolverDiscreteFunctionType SolverDiscreteFunctionType;
78 
79  // \brief The column-compressed matrix type.
80  typedef typename Traits :: CCSMatrixType CCSMatrixType;
81  typedef typename DiscreteFunctionType::DofType DofType;
82  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
83 
91  SPQRInverseOperator(const double& redEps, const double& absLimit, const int& maxIter, const bool& verbose,
92  const ParameterReader &parameter = Parameter::container() ) :
93  SPQRInverseOperator(parameter)
94  {}
95 
102  SPQRInverseOperator(const double& redEps, const double& absLimit, const int& maxIter,
103  const ParameterReader &parameter = Parameter::container() ) :
104  SPQRInverseOperator(parameter)
105  {}
106 
107  SPQRInverseOperator(const SolverParameter& parameter = SolverParameter( Parameter::container() ) )
108  : BaseType(parameter)
109  , verbose_(BaseType::verbose())
110  , ccsmat_(), cc_(new cholmod_common())
111  {
112  cholmod_l_start(cc_);
113  }
114 
122  SPQRInverseOperator(const OperatorType& op, const double& redEps, const double& absLimit, const int& maxIter, const bool& verbose,
123  const ParameterReader &parameter = Parameter::container() ) :
124  SPQRInverseOperator(parameter)
125  {
126  bind(op);
127  }
128 
135  SPQRInverseOperator(const OperatorType& op, const double& redEps, const double& absLimit, const int& maxIter,
136  const ParameterReader &parameter = Parameter::container() ) :
137  SPQRInverseOperator(parameter)
138  {
139  bind(op);
140  }
141 
142  SPQRInverseOperator(const OperatorType& op, const ParameterReader &parameter = Parameter::container() ) :
143  SPQRInverseOperator(parameter)
144  {
145  bind(op);
146  }
147 
148  SPQRInverseOperator(const SPQRInverseOperator& other) :
149  SPQRInverseOperator(other.parameter())
150  {
151  if( other.operator_ )
152  bind( *(other.operator_) );
153  }
154 
155  // \brief Destructor.
156  ~SPQRInverseOperator()
157  {
158  finalize();
159  cholmod_l_finish(cc_);
160  delete cc_;
161  }
162 
163  using BaseType :: bind;
164 
165  void unbind ()
166  {
167  BaseType::unbind();
168  finalize();
169  }
170 
171  // \brief Decompose matrix.
172  template<typename... A>
173  void prepare(A... ) const
174  {
175  if(assembledOperator_ && !ccsmat_ )
176  {
177  ccsmat_.reset( new CCSMatrixType(assembledOperator_->exportMatrix() ) );
178  decompose();
179  }
180  }
181 
182  // \brief Free allocated memory.
183  virtual void finalize()
184  {
185  if( ccsmat_ )
186  {
187  ccsmat_.reset();
188  cholmod_l_free_sparse(&A_, cc_);
189  cholmod_l_free_dense(&B_, cc_);
190  SuiteSparseQR_free<DofType>(&spqrfactorization_, cc_);
191  }
192  }
193 
200  int apply (const DofType* arg, DofType* dest) const
201  {
202  assert(ccsmat_);
203  const std::size_t dimMat(ccsmat_->N());
204  // fill B
205  for(std::size_t k = 0; k != dimMat; ++k)
206  (static_cast<DofType*>(B_->x))[k] = arg[k];
207  cholmod_dense* BTemp = B_;
208  B_ = SuiteSparseQR_qmult<DofType>(0, spqrfactorization_, B_, cc_);
209  cholmod_dense* X = SuiteSparseQR_solve<DofType>(1, spqrfactorization_, B_, cc_);
210  cholmod_l_free_dense(&BTemp, cc_);
211  // fill x
212  for(std::size_t k = 0; k != dimMat; ++k)
213  dest[k] = (static_cast<DofType*>(X->x))[k];
214  cholmod_l_free_dense(&X, cc_);
215  // output some statistics
216  if(verbose_ > 0)
217  printDecompositionInfo();
218  return 1;
219  }
220 
227  int apply(const SolverDiscreteFunctionType& arg, SolverDiscreteFunctionType& dest) const
228  {
229  prepare();
230  apply(arg.leakPointer(),dest.leakPointer());
231  const_cast<ThisType*>(this)->finalize();
232  return 1;
233  }
234 
241  int apply(const ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpaceType>& arg,
242  ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpaceType>& dest) const
243  {
244  // copy DOF's arg into a consecutive vector
245  std::vector<DofType> vecArg(arg.size());
246  std::copy(arg.dbegin(),arg.dend(),vecArg.begin());
247  std::vector<DofType> vecDest(dest.size());
248  // apply operator
249  apply(vecArg.data(),vecDest.data());
250  // copy back solution into dest
251  std::copy(vecDest.begin(),vecDest.end(),dest.dbegin());
252  return 1;
253  }
254 
261  template<typename... DFs>
262  int apply(const TupleDiscreteFunction<DFs...>& arg,TupleDiscreteFunction<DFs...>& dest) const
263  {
264  // copy DOF's arg into a consecutive vector
265  std::vector<DofType> vecArg(arg.size());
266  auto vecArgIt(vecArg.begin());
267  Hybrid::forEach(std::make_index_sequence<sizeof...(DFs)>{},
268  [&](auto i){vecArgIt=std::copy(std::get<i>(arg).dbegin(),std::get<i>(arg).dend(),vecArgIt);});
269  std::vector<DofType> vecDest(dest.size());
270  // apply operator
271  apply(vecArg.data(),vecDest.data());
272  // copy back solution into dest
273  auto vecDestIt(vecDest.begin());
274  Hybrid::forEach(std::make_index_sequence<sizeof...(DFs)>{},[&](auto i){for(auto& dof:dofs(std::get<i>(dest))) dof=(*(vecDestIt++));});
275  return 1;
276  }
277 
278  void printTexInfo(std::ostream& out) const
279  {
280  out<<"Solver: SPQR direct solver"<<std::endl;
281  }
282 
283  // \brief Print some statistics about the SPQR decomposition.
284  void printDecompositionInfo() const
285  {
286  if( ccsmat_ )
287  {
288  std::cout<<std::endl<<"Solving with SuiteSparseQR"<<std::endl;
289  std::cout<<"Flops Taken: "<<cc_->SPQR_flopcount<<std::endl;
290  std::cout<<"Analysis Time: "<<cc_->SPQR_analyze_time<<" s"<<std::endl;
291  std::cout<<"Factorize Time: "<<cc_->SPQR_factorize_time<<" s"<<std::endl;
292  std::cout<<"Backsolve Time: "<<cc_->SPQR_solve_time<<" s"<<std::endl;
293  std::cout<<"Peak Memory Usage: "<<cc_->memory_usage<<" bytes"<<std::endl;
294  std::cout<<"Rank Estimate: "<<cc_->SPQR_istat[4]<<std::endl<<std::endl;
295  }
296  }
297 
298  void setMaxIterations( int ) {}
299 
303  SuiteSparseQR_factorization<DofType>* getFactorization()
304  {
305  return spqrfactorization_;
306  }
307 
311  CCSMatrixType& getCCSMatrix() const
312  {
313  assert( ccsmat_ );
314  return *ccsmat_;
315  }
316 
317 private:
318  using BaseType :: operator_;
319  using BaseType :: assembledOperator_;
320  const bool verbose_;
321  mutable std::unique_ptr< CCSMatrixType > ccsmat_;
322  mutable cholmod_common* cc_ = nullptr ;
323  mutable cholmod_sparse* A_ = nullptr ;
324  mutable cholmod_dense* B_ = nullptr ;
325  mutable SuiteSparseQR_factorization<DofType>* spqrfactorization_ = nullptr;
326 
327  // \brief Computes the SPQR decomposition.
328  void decompose() const
329  {
330  CCSMatrixType& ccsmat = getCCSMatrix();
331 
332  const std::size_t dimMat(ccsmat.N());
333  const std::size_t nnz(ccsmat.getColStart()[dimMat]);
334  // initialise the matrix A
335  bool sorted(true);
336  bool packed(true);
337  bool real(std::is_same<DofType,double>::value);
338  A_ = cholmod_l_allocate_sparse(dimMat, dimMat, nnz, sorted, packed, symmetric, real, cc_);
339  // copy all the entries of Ap, Ai, Ax
340  for(std::size_t k = 0; k != (dimMat+1); ++k)
341  (static_cast<long int *>(A_->p))[k] = ccsmat.getColStart()[k];
342  for(std::size_t k = 0; k != nnz; ++k)
343  {
344  (static_cast<long int*>(A_->i))[k] = ccsmat.getRowIndex()[k];
345  (static_cast<DofType*>(A_->x))[k] = ccsmat.getValues()[k];
346  }
347  // initialise the vector B
348  B_ = cholmod_l_allocate_dense(dimMat, 1, dimMat, A_->xtype, cc_);
349  // compute factorization of A
350  spqrfactorization_=SuiteSparseQR_factorize<DofType>(SPQR_ORDERING_DEFAULT,SPQR_DEFAULT_TOL,A_,cc_);
351  }
352 };
353 
354 
355 } // end namespace Fem
356 } // end namespace Dune
357 
358 #endif // #if HAVE_SUITESPARSE_SPQR
359 
360 #endif // #ifndef DUNE_FEM_SPQRSOLVER_HH
Definition: bindguard.hh:11
double real(const std::complex< Double > &x)
Definition: double.hh:906
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
static ParameterContainer & container()
Definition: io/parameter.hh:193