dune-fem 2.8.0
Loading...
Searching...
No Matches
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>
20
21#if HAVE_SUITESPARSE_SPQR
22#include <SuiteSparseQR.hpp>
23
24namespace Dune
25{
26namespace Fem
27{
28
46template<class DiscreteFunction, bool symmetric=false,
47 class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
48class SPQRInverseOperator;
49
50template<class DiscreteFunction, bool symmetric, class Matrix>
51struct 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
68template< class DF, bool symmetric, class Matrix >
69class 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;
74public:
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
317private:
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