1#ifndef DUNE_FEM_SPQRSOLVER_HH
2#define DUNE_FEM_SPQRSOLVER_HH
10#include <dune/common/hybridutilities.hh>
21#if HAVE_SUITESPARSE_SPQR
22#include <SuiteSparseQR.hpp>
46template<
class DiscreteFunction,
bool symmetric=
false,
47 class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
48class SPQRInverseOperator;
50template<
class DiscreteFunction,
bool symmetric,
class Matrix>
51struct SPQRInverseOperatorTraits
53 typedef DiscreteFunction DiscreteFunctionType;
54 typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
57 typedef OperatorType PreconditionerType;
59 typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction, Matrix > AssembledOperatorType;
60 typedef ColCompMatrix<typename AssembledOperatorType::MatrixType::MatrixBaseType,int > CCSMatrixType;
62 typedef SPQRInverseOperator< DiscreteFunction, symmetric, Matrix > InverseOperatorType;
63 typedef SolverParameter SolverParameterType;
68template<
class DF,
bool symmetric,
class Matrix >
69class SPQRInverseOperator :
public InverseOperatorInterface< SPQRInverseOperatorTraits< DF, symmetric, Matrix > >
71 typedef SPQRInverseOperatorTraits< DF, symmetric, Matrix > Traits;
72 typedef InverseOperatorInterface< Traits > BaseType;
73 typedef SPQRInverseOperator< DF, symmetric, Matrix > ThisType;
75 typedef DF DiscreteFunctionType;
76 typedef typename BaseType :: OperatorType OperatorType;
77 typedef typename BaseType :: SolverDiscreteFunctionType SolverDiscreteFunctionType;
80 typedef typename Traits :: CCSMatrixType CCSMatrixType;
81 typedef typename DiscreteFunctionType::DofType DofType;
82 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
91 SPQRInverseOperator(
const double& redEps,
const double& absLimit,
const int& maxIter,
const bool& verbose,
92 const ParameterReader ¶meter = Parameter::container() ) :
93 SPQRInverseOperator(parameter)
102 SPQRInverseOperator(
const double& redEps,
const double& absLimit,
const int& maxIter,
103 const ParameterReader ¶meter = Parameter::container() ) :
104 SPQRInverseOperator(parameter)
107 SPQRInverseOperator(
const SolverParameter& parameter = SolverParameter( Parameter::container() ) )
108 : BaseType(parameter)
109 , verbose_(BaseType::verbose())
110 , ccsmat_(), cc_(new cholmod_common())
112 cholmod_l_start(cc_);
122 SPQRInverseOperator(
const OperatorType& op,
const double& redEps,
const double& absLimit,
const int& maxIter,
const bool& verbose,
123 const ParameterReader ¶meter = Parameter::container() ) :
124 SPQRInverseOperator(parameter)
135 SPQRInverseOperator(
const OperatorType& op,
const double& redEps,
const double& absLimit,
const int& maxIter,
136 const ParameterReader ¶meter = Parameter::container() ) :
137 SPQRInverseOperator(parameter)
142 SPQRInverseOperator(
const OperatorType& op,
const ParameterReader ¶meter = Parameter::container() ) :
143 SPQRInverseOperator(parameter)
148 SPQRInverseOperator(
const SPQRInverseOperator& other) :
149 SPQRInverseOperator(other.parameter())
151 if( other.operator_ )
152 bind( *(other.operator_) );
156 ~SPQRInverseOperator()
159 cholmod_l_finish(cc_);
163 using BaseType :: bind;
172 template<
typename... A>
173 void prepare(A... )
const
175 if(assembledOperator_ && !ccsmat_ )
177 ccsmat_.reset(
new CCSMatrixType(assembledOperator_->exportMatrix() ) );
183 virtual void finalize()
188 cholmod_l_free_sparse(&A_, cc_);
189 cholmod_l_free_dense(&B_, cc_);
190 SuiteSparseQR_free<DofType>(&spqrfactorization_, cc_);
200 int apply (
const DofType* arg, DofType* dest)
const
203 const std::size_t dimMat(ccsmat_->N());
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_);
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_);
217 printDecompositionInfo();
227 int apply(
const SolverDiscreteFunctionType& arg, SolverDiscreteFunctionType& dest)
const
230 apply(arg.leakPointer(),dest.leakPointer());
231 const_cast<ThisType*
>(
this)->finalize();
241 int apply(
const ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpaceType>& arg,
242 ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpaceType>& dest)
const
245 std::vector<DofType> vecArg(arg.size());
246 std::copy(arg.dbegin(),arg.dend(),vecArg.begin());
247 std::vector<DofType> vecDest(dest.size());
249 apply(vecArg.data(),vecDest.data());
251 std::copy(vecDest.begin(),vecDest.end(),dest.dbegin());
261 template<
typename... DFs>
262 int apply(
const TupleDiscreteFunction<DFs...>& arg,TupleDiscreteFunction<DFs...>& dest)
const
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());
271 apply(vecArg.data(),vecDest.data());
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++));});
278 void printTexInfo(std::ostream& out)
const
280 out<<
"Solver: SPQR direct solver"<<std::endl;
284 void printDecompositionInfo()
const
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;
298 void setMaxIterations(
int ) {}
303 SuiteSparseQR_factorization<DofType>* getFactorization()
305 return spqrfactorization_;
311 CCSMatrixType& getCCSMatrix()
const
318 using BaseType :: operator_;
319 using BaseType :: assembledOperator_;
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;
328 void decompose()
const
330 CCSMatrixType& ccsmat = getCCSMatrix();
332 const std::size_t dimMat(ccsmat.N());
333 const std::size_t nnz(ccsmat.getColStart()[dimMat]);
337 bool real(std::is_same<DofType,double>::value);
338 A_ = cholmod_l_allocate_sparse(dimMat, dimMat, nnz, sorted, packed, symmetric, real, cc_);
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)
344 (
static_cast<long int*
>(A_->i))[k] = ccsmat.getRowIndex()[k];
345 (
static_cast<DofType*
>(A_->x))[k] = ccsmat.getValues()[k];
348 B_ = cholmod_l_allocate_dense(dimMat, 1, dimMat, A_->xtype, cc_);
350 spqrfactorization_=SuiteSparseQR_factorize<DofType>(SPQR_ORDERING_DEFAULT,SPQR_DEFAULT_TOL,A_,cc_);
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