dune-fem 2.8.0
Loading...
Searching...
No Matches
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>
20
21#if HAVE_SUITESPARSE_UMFPACK
23
24namespace Dune
25{
26namespace Fem
27{
28
38template<class DiscreteFunction, class Matrix>
39class UMFPACKInverseOperator;
40
41template<class DiscreteFunction, class Matrix>
42struct 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
63template<class DiscreteFunction,
64 class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
65class UMFPACKInverseOperator :
66 public InverseOperatorInterface< UMFPACKInverseOperatorTraits< DiscreteFunction, Matrix > >
67{
68public:
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
240protected:
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