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> 11 #include <dune/fem/function/adaptivefunction.hh> 12 #include <dune/fem/function/blockvectorfunction.hh> 13 #include <dune/fem/function/tuplediscretefunction.hh> 14 #include <dune/fem/io/parameter.hh> 15 #include <dune/fem/operator/common/operator.hh> 16 #include <dune/fem/operator/linear/spoperator.hh> 17 #include <dune/fem/operator/matrix/colcompspmatrix.hh> 18 #include <dune/fem/operator/linear/spoperator.hh> 19 #include <dune/fem/solver/inverseoperatorinterface.hh> 20 21 #if HAVE_SUITESPARSE_UMFPACK 22 #include <dune/fem/misc/umfpack.hh> 23 24 namespace Dune 25 { 26 namespace Fem 27 { 28 29 /** @addtogroup DirectSolver 30 * 31 * In this section implementations of direct solvers 32 * for solving linear systems of the from 33 * \f$A x = b\f$, where \f$A\f$ is a Mapping or 34 * Operator and \f$x\f$ and \f$b\f$ are discrete functions 35 * (see DiscreteFunctionInterface) can be found. 36 */ 37 38 template<class DiscreteFunction, class Matrix> 39 class UMFPACKInverseOperator; 40 41 template<class DiscreteFunction, class Matrix> 42 struct UMFPACKInverseOperatorTraits 43 { 44 typedef DiscreteFunction DiscreteFunctionType; 45 typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType; 46 47 typedef Dune::Fem::Operator< DiscreteFunction, DiscreteFunction > OperatorType; 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 56 /** \class UMFPACKInverseOperator 57 * \ingroup DirectSolver 58 * \brief The %UMFPack direct sparse solver 59 * %UMFPack will always go double precision and supports complex numbers. 60 * Details on UMFPack can be found on http://www.cise.ufl.edu/research/sparse/umfpack/ 61 * \note This will only work if dune-fem has been configured to use UMFPACK 62 */ 63 template<class DiscreteFunction, 64 class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > > 65 class UMFPACKInverseOperator : 66 public InverseOperatorInterface< UMFPACKInverseOperatorTraits< DiscreteFunction, Matrix > > 67 { 68 public: 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 89 /** \brief Constructor. 90 * \param[in] parameter parameters for the solver 91 */ UMFPACKInverseOperator(const SolverParameter & parameter=SolverParameter (Parameter::container ()))92 UMFPACKInverseOperator(const SolverParameter ¶meter = 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. ~UMFPACKInverseOperator()103 ~UMFPACKInverseOperator() 104 { 105 unbind(); 106 } 107 bind(const OperatorType & op)108 void bind( const OperatorType& op ) 109 { 110 // clear old storage 111 finalize(); 112 BaseType::bind( op ); 113 assert( assembledOperator_ ); 114 } 115 unbind()116 void unbind () 117 { 118 finalize(); 119 BaseType::unbind(); 120 } 121 122 // \brief Decompose matrix. 123 template<typename... A> prepare(A...) const124 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. finalize()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 145 /** \brief Solve the system. 146 * \param[in] arg right hand side 147 * \param[out] dest solution 148 * \warning You have to decompose the matrix before calling the apply (using the method prepare) 149 * and you have free the decompistion when is not needed anymore (using the method finalize). 150 */ apply(const DofType * arg,DofType * dest) const151 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 174 /** \brief Solve the system. 175 * \param[in] arg right hand side 176 * \param[out] dest solution 177 * \warning You have to decompose the matrix before calling the apply (using the method prepare) 178 * and you have free the decompistion when is not needed anymore (using the method finalize). 179 */ apply(const SolverDiscreteFunctionType & arg,SolverDiscreteFunctionType & dest) const180 int apply(const SolverDiscreteFunctionType& arg, SolverDiscreteFunctionType& dest) const 181 { 182 return apply(arg.leakPointer(), dest.leakPointer()); 183 } 184 185 /** \brief Solve the system. 186 * \param[in] arg right hand side 187 * \param[out] dest solution 188 * \warning You have to decompose the matrix before calling the apply (using the method prepare) 189 * and you have free the decompistion when is not needed anymore (using the method finalize). 190 */ 191 template<typename... DFs> apply(const TupleDiscreteFunction<DFs...> & arg,TupleDiscreteFunction<DFs...> & dest) const192 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 printTexInfo(std::ostream & out) const208 void printTexInfo(std::ostream& out) const 209 { 210 out<<"Solver: UMFPACK direct solver"<<std::endl; 211 } 212 setMaxIterations(int)213 void setMaxIterations ( int ) {} 214 215 /** \brief Get CCS matrix of the operator to solve. 216 * \warning It is up to the user to preserve consistency. 217 */ getCCSMatrix() const218 CCSMatrixType& getCCSMatrix() const 219 { 220 assert( ccsmat_ ); 221 return *ccsmat_; 222 } 223 224 // \brief Print some statistics about the UMFPACK decomposition. printDecompositionInfo() const225 void printDecompositionInfo() const 226 { 227 Caller::report_info(UMF_Control,UMF_Decomposition_Info); 228 } 229 UMFPACKInverseOperator(const UMFPACKInverseOperator & other)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 240 protected: 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. decompose() const251 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 280