1#ifndef DUNE_FEM_UMFPACKSOLVER_HH
2#define DUNE_FEM_UMFPACKSOLVER_HH
10#include <dune/common/hybridutilities.hh>
21#if HAVE_SUITESPARSE_UMFPACK
38template<
class DiscreteFunction,
class Matrix>
39class UMFPACKInverseOperator;
41template<
class DiscreteFunction,
class Matrix>
42struct UMFPACKInverseOperatorTraits
44 typedef DiscreteFunction DiscreteFunctionType;
45 typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
48 typedef OperatorType PreconditionerType;
50 typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction, Matrix > AssembledOperatorType;
52 typedef UMFPACKInverseOperator< DiscreteFunction,Matrix> InverseOperatorType;
53 typedef SolverParameter SolverParameterType;
63template<
class DiscreteFunction,
64 class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
65class UMFPACKInverseOperator :
66 public InverseOperatorInterface< UMFPACKInverseOperatorTraits< DiscreteFunction, Matrix > >
69 typedef UMFPACKInverseOperatorTraits< DiscreteFunction, Matrix > Traits;
70 typedef InverseOperatorInterface< Traits > BaseType;
72 typedef typename BaseType :: SolverDiscreteFunctionType
73 SolverDiscreteFunctionType;
75 typedef typename BaseType :: OperatorType OperatorType;
76 typedef typename BaseType :: AssembledOperatorType AssembledOperatorType;
78 typedef UMFPACKInverseOperator< DiscreteFunction,Matrix> ThisType;
80 typedef DiscreteFunction DiscreteFunctionType;
83 typedef ColCompMatrix<typename AssembledOperatorType::MatrixType::MatrixBaseType,long int> CCSMatrixType;
84 typedef typename DiscreteFunctionType::DofType DofType;
85 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
87 using BaseType :: parameter_;
92 UMFPACKInverseOperator(
const SolverParameter ¶meter = SolverParameter(Parameter::container()) )
93 : BaseType(parameter),
95 UMF_Symbolic( nullptr ),
96 UMF_Numeric( nullptr )
98 Caller::defaults(UMF_Control);
99 UMF_Control[UMFPACK_PRL] = 4;
103 ~UMFPACKInverseOperator()
108 void bind(
const OperatorType& op )
112 BaseType::bind( op );
113 assert( assembledOperator_ );
123 template<
typename... A>
124 void prepare(A... )
const
126 if(assembledOperator_ && !ccsmat_)
128 ccsmat_ = std::make_unique<CCSMatrixType>(assembledOperator_->exportMatrix());
134 virtual void finalize()
138 getCCSMatrix().free();
140 Caller::free_symbolic(&UMF_Symbolic); UMF_Symbolic =
nullptr;
141 Caller::free_numeric(&UMF_Numeric); UMF_Numeric =
nullptr;
151 int apply(
const DofType* arg, DofType* dest)
const
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() )
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;
169 const_cast<ThisType*
>(
this)->finalize();
180 int apply(
const SolverDiscreteFunctionType& arg, SolverDiscreteFunctionType& dest)
const
182 return apply(arg.leakPointer(), dest.leakPointer());
191 template<
typename... DFs>
192 int apply(
const TupleDiscreteFunction<DFs...>& arg,TupleDiscreteFunction<DFs...>& dest)
const
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());
201 apply(vecArg.data(),vecDest.data());
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++));});
208 void printTexInfo(std::ostream& out)
const
210 out<<
"Solver: UMFPACK direct solver"<<std::endl;
213 void setMaxIterations (
int ) {}
218 CCSMatrixType& getCCSMatrix()
const
225 void printDecompositionInfo()
const
227 Caller::report_info(UMF_Control,UMF_Decomposition_Info);
230 UMFPACKInverseOperator(
const UMFPACKInverseOperator &other)
233 UMF_Symbolic(other.UMF_Symbolic),
234 UMF_Numeric(other.UMF_Numeric)
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];
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];
248 typedef typename Dune::UMFPackMethodChooser<DofType> Caller;
251 void decompose()
const
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() )
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;
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76