16 #ifndef dealii_differentiation_sd_symengine_optimizer_h
17 #define dealii_differentiation_sd_symengine_optimizer_h
21 #ifdef DEAL_II_WITH_SYMENGINE
25 # include <symengine/basic.h>
26 # include <symengine/dict.h>
27 # include <symengine/symengine_exception.h>
28 # include <symengine/symengine_rcp.h>
31 # include <symengine/lambda_double.h>
32 # include <symengine/visitor.h>
33 # ifdef HAVE_SYMENGINE_LLVM
34 # include <symengine/llvm_double.h>
48 # include <boost/serialization/split_member.hpp>
49 # include <boost/type_traits.hpp>
54 # include <type_traits>
76 "SymEngine has not been built with LLVM support.");
83 "The SymEngine LLVM optimizer does not (yet) support the "
84 "selected return type.");
90 template <
typename ReturnType>
121 template <
class StreamType>
183 static_cast<unsigned int>(f2));
215 static_cast<unsigned int>(f2));
261 const bool use_agg_opt =
263 const int opt_level = (use_agg_opt ? 3 : 2);
273 template <
class StreamType>
277 s <<
" OptimizationFlags|";
300 template <
typename ReturnType,
typename T =
void>
313 template <
typename ReturnType,
typename T =
void>
317 # ifdef HAVE_SYMENGINE_LLVM
327 template <
typename ReturnType,
typename T =
void>
328 struct LLVMOptimizer;
347 template <
typename ReturnType,
typename Optimizer,
typename T =
void>
359 template <
typename ReturnType_,
typename T =
void>
360 struct SupportedOptimizerTypeTraits
362 static const bool is_supported =
false;
364 using ReturnType = void;
370 template <
typename ReturnType_>
371 struct SupportedOptimizerTypeTraits<
373 std::enable_if_t<std::is_arithmetic<ReturnType_>::value>>
375 static const bool is_supported =
true;
378 typename std::conditional<std::is_same<ReturnType_, float>::value,
386 template <
typename ReturnType_>
387 struct SupportedOptimizerTypeTraits<
390 boost::is_complex<ReturnType_>::value &&
391 std::is_arithmetic<typename ReturnType_::value_type>::value>>
393 static const bool is_supported =
true;
395 using ReturnType =
typename std::conditional<
396 std::is_same<ReturnType_, std::complex<float>>::value,
398 std::complex<double>>::type;
403 template <
typename ReturnType_>
404 struct DictionaryOptimizer<ReturnType_,
405 std::enable_if_t<SupportedOptimizerTypeTraits<
406 ReturnType_>::is_supported>>
409 typename SupportedOptimizerTypeTraits<ReturnType_>::ReturnType;
411 internal::DictionarySubstitutionVisitor<ReturnType, SD::Expression>;
424 const SymEngine::vec_basic & independent_symbols,
425 const SymEngine::vec_basic & dependent_functions,
429 optimizer.init(independent_symbols,
440 template <
class Archive>
442 save(Archive & archive,
443 const unsigned int version,
446 optimizer.save(archive, version);
455 template <
class Archive>
457 load(Archive & archive,
458 const unsigned int version,
460 const SymEngine::vec_basic & ,
461 const SymEngine::vec_basic & ,
464 optimizer.load(archive, version);
484 template <
typename Stream>
486 print(Stream & stream,
488 const bool print_independent_symbols =
false,
489 const bool print_dependent_functions =
false,
490 const bool print_cse_reductions =
true)
492 optimizer.print(stream,
493 print_independent_symbols,
494 print_dependent_functions,
495 print_cse_reductions);
501 template <
typename ReturnType_>
502 struct LambdaOptimizer<ReturnType_,
503 std::enable_if_t<SupportedOptimizerTypeTraits<
504 ReturnType_>::is_supported>>
507 typename std::conditional<!boost::is_complex<ReturnType_>::value,
509 std::complex<double>>::type;
511 !boost::is_complex<ReturnType_>::value,
512 SymEngine::LambdaRealDoubleVisitor,
513 SymEngine::LambdaComplexDoubleVisitor>::type;
526 const SymEngine::vec_basic & independent_symbols,
527 const SymEngine::vec_basic & dependent_functions,
531 optimizer.init(independent_symbols,
542 template <
class Archive>
554 template <
class Archive>
559 const SymEngine::vec_basic & independent_symbols,
560 const SymEngine::vec_basic & dependent_functions,
563 initialize(optimizer,
586 template <
typename StreamType>
600 # ifdef HAVE_SYMENGINE_LLVM
601 template <
typename ReturnType_>
602 struct LLVMOptimizer<
604 std::enable_if_t<std::is_arithmetic<ReturnType_>::value>>
607 typename std::conditional<std::is_same<ReturnType_, float>::value,
611 typename std::conditional<std::is_same<ReturnType_, float>::value,
612 SymEngine::LLVMFloatVisitor,
613 SymEngine::LLVMDoubleVisitor>::type;
619 static const bool supported_by_LLVM =
true;
632 const SymEngine::vec_basic & independent_symbols,
633 const SymEngine::vec_basic & dependent_functions,
638 optimizer.init(independent_symbols,
650 template <
class Archive>
652 save(Archive &archive,
656 const std::string llvm_compiled_function = optimizer.dumps();
657 archive & llvm_compiled_function;
666 template <
class Archive>
668 load(Archive &archive,
671 const SymEngine::vec_basic & ,
672 const SymEngine::vec_basic & ,
675 std::string llvm_compiled_function;
676 archive & llvm_compiled_function;
677 optimizer.loads(llvm_compiled_function);
697 template <
typename StreamType>
715 template <
typename ReturnType_>
716 struct LLVMOptimizer<
719 boost::is_complex<ReturnType_>::value &&
720 std::is_arithmetic<typename ReturnType_::value_type>::value>>
724 using ReturnType =
typename LambdaOptimizer<ReturnType_>::ReturnType;
732 static const bool supported_by_LLVM =
false;
745 const SymEngine::vec_basic & ,
746 const SymEngine::vec_basic & ,
758 template <
class Archive>
773 template <
class Archive>
778 const SymEngine::vec_basic & ,
779 const SymEngine::vec_basic & ,
802 template <
typename StreamType>
819 template <
typename ReturnType,
typename Optimizer>
820 struct OptimizerHelper<
824 std::is_same<ReturnType, typename Optimizer::ReturnType>::value>>
836 const SymEngine::vec_basic & independent_symbols,
837 const SymEngine::vec_basic & dependent_functions,
845 Optimizer::initialize(*optimizer,
868 std::vector<ReturnType> & output_values,
869 const std::vector<ReturnType> & substitution_values)
872 optimizer->call(output_values.data(), substitution_values.data());
881 template <
class Archive>
883 save(Archive & archive,
884 const unsigned int version,
892 Optimizer::save(archive, version, *optimizer);
901 template <
class Archive>
903 load(Archive & archive,
904 const unsigned int version,
906 const SymEngine::vec_basic & independent_symbols,
907 const SymEngine::vec_basic & dependent_functions,
915 Optimizer::load(archive,
940 template <
typename Stream>
942 print(Stream & stream,
944 const bool print_independent_symbols =
false,
945 const bool print_dependent_functions =
false,
946 const bool print_cse_reductions =
true)
953 Optimizer::print(stream,
955 print_independent_symbols,
956 print_dependent_functions,
957 print_cse_reductions);
961 template <
typename ReturnType,
typename Optimizer>
962 struct OptimizerHelper<
966 !std::is_same<ReturnType, typename Optimizer::ReturnType>::value>>
978 const SymEngine::vec_basic & independent_symbols,
979 const SymEngine::vec_basic & dependent_functions,
985 optimizer->init(independent_symbols,
1007 std::vector<ReturnType> & output_values,
1008 const std::vector<ReturnType> & substitution_values)
1014 std::vector<typename Optimizer::ReturnType> int_outputs(
1015 output_values.size());
1016 std::vector<typename Optimizer::ReturnType> int_inputs(
1017 substitution_values.size());
1020 substitution_values.end(),
1021 int_inputs.begin());
1022 optimizer->call(int_outputs.data(), int_inputs.data());
1025 output_values.begin());
1034 template <
class Archive>
1036 save(Archive & archive,
1037 const unsigned int version,
1041 Optimizer::save(archive, version, *optimizer);
1050 template <
class Archive>
1052 load(Archive & archive,
1053 const unsigned int version,
1055 const SymEngine::vec_basic & independent_symbols,
1056 const SymEngine::vec_basic & dependent_functions,
1064 Optimizer::load(archive,
1067 independent_symbols,
1068 dependent_functions,
1069 optimization_flags);
1089 template <
typename Stream>
1091 print(Stream & stream,
1093 const bool print_cse_reductions =
true,
1094 const bool print_independent_symbols =
false,
1095 const bool print_dependent_functions =
false)
1099 optimizer->print(stream,
1100 print_independent_symbols,
1101 print_dependent_functions,
1102 print_cse_reductions);
1133 template <
typename NumberType,
1136 template <
int,
int,
typename>
1138 TensorType<rank, dim, NumberType>
1140 const TensorType<rank, dim, Expression> &symbol_tensor,
1141 const std::vector<NumberType> & cached_evaluation,
1144 TensorType<rank, dim, NumberType> out;
1145 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
1148 out.unrolled_to_component_indices(i));
1150 optimizer.
extract(symbol_tensor[indices], cached_evaluation);
1178 template <
typename NumberType,
int dim>
1182 const std::vector<NumberType> & cached_evaluation,
1186 for (
unsigned int i = 0;
1187 i < SymmetricTensor<2, dim>::n_independent_components;
1189 for (
unsigned int j = 0;
1190 j < SymmetricTensor<2, dim>::n_independent_components;
1194 make_rank_4_tensor_indices<dim>(i, j);
1196 optimizer.
extract(symbol_tensor[indices], cached_evaluation);
1219 template <
typename NumberType,
typename T>
1245 template <
typename NumberType,
typename T>
1274 template <
typename NumberType,
typename T,
typename... Args>
1278 const Args &...other_functions)
1298 template <
int,
int,
typename>
1302 const TensorType<rank, dim, Expression> &symbol_tensor)
1305 out.reserve(symbol_tensor.n_independent_components);
1306 for (
unsigned int i = 0; i < symbol_tensor.n_independent_components;
1310 symbol_tensor.unrolled_to_component_indices(i));
1311 out.push_back(symbol_tensor[indices].get_RCP());
1333 for (
unsigned int i = 0;
1334 i < SymmetricTensor<2, dim>::n_independent_components;
1336 for (
unsigned int j = 0;
1337 j < SymmetricTensor<2, dim>::n_independent_components;
1341 make_rank_4_tensor_indices<dim>(i, j);
1342 out.push_back(symbol_tensor[indices].get_RCP());
1441 template <
typename ReturnType>
1518 template <
typename Stream>
1520 print(Stream &stream,
const bool print_cse =
false)
const;
1530 template <
class Archive>
1532 save(Archive &archive,
const unsigned int version)
const;
1547 template <
class Archive>
1549 load(Archive &archive,
const unsigned int version);
1573 template <
class Archive>
1579 BOOST_SERIALIZATION_SPLIT_MEMBER()
1662 template <
int rank,
int dim>
1670 template <
int rank,
int dim>
1698 template <
typename T>
1715 template <
typename T,
typename... Args>
1848 const std::vector<ReturnType> &
values)
const;
1861 substitute(
const SymEngine::vec_basic & symbols,
1862 const std::vector<ReturnType> &
values)
const;
1901 const std::vector<ReturnType> &
1922 std::vector<ReturnType>
1923 evaluate(
const std::vector<Expression> &funcs)
const;
1933 template <
int rank,
int dim>
1946 template <
int rank,
int dim>
1960 const std::vector<ReturnType> &cached_evaluation)
const;
1970 std::vector<ReturnType>
1971 extract(
const std::vector<Expression> &funcs,
1972 const std::vector<ReturnType> &cached_evaluation)
const;
1982 template <
int rank,
int dim>
1985 const std::vector<ReturnType> & cached_evaluation)
const;
1995 template <
int rank,
int dim>
1998 const std::vector<ReturnType> &cached_evaluation)
const;
2046 const SymEngine::RCP<const SymEngine::Basic> &
function)
const;
2143 substitute(
const std::vector<ReturnType> &substitution_values)
const;
2154 template <
typename ReturnType>
2155 template <
typename Stream>
2161 stream <<
"Method? " << optimization_method() <<
'\n';
2162 stream <<
"Flags: " << optimization_flags() <<
'\n';
2163 stream <<
"Optimized? " << (optimized() ?
"Yes" :
"No") <<
'\n';
2164 stream <<
"Values substituted? " << values_substituted() <<
"\n\n";
2167 stream <<
"Symbols (" << n_independent_variables()
2168 <<
" independent variables):" <<
'\n';
2170 for (SD::types::substitution_map::const_iterator it =
2171 independent_variables_symbols.begin();
2172 it != independent_variables_symbols.end();
2175 stream << cntr <<
": " << it->first <<
'\n';
2177 stream <<
'\n' << std::flush;
2180 stream <<
"Functions (" << n_dependent_variables()
2181 <<
" dependent variables):" <<
'\n';
2183 for (
typename SD::types::symbol_vector::const_iterator it =
2184 dependent_variables_functions.begin();
2185 it != dependent_variables_functions.end();
2188 stream << cntr <<
": " << (*it) <<
'\n';
2190 stream <<
'\n' << std::flush;
2196 const bool print_cse_reductions =
true;
2197 const bool print_independent_symbols =
false;
2198 const bool print_dependent_functions =
false;
2202 Assert(
dynamic_cast<typename internal::DictionaryOptimizer<
2204 ExcMessage(
"Cannot cast optimizer to Dictionary type."));
2206 internal::OptimizerHelper<
2208 internal::DictionaryOptimizer<ReturnType>>::
2210 dynamic_cast<typename internal::DictionaryOptimizer<
2212 print_independent_symbols,
2213 print_dependent_functions,
2214 print_cse_reductions);
2216 stream <<
'\n' << std::flush;
2220 Assert(
dynamic_cast<typename internal::LambdaOptimizer<
2222 ExcMessage(
"Cannot cast optimizer to Lambda type."));
2224 internal::OptimizerHelper<ReturnType,
2225 internal::LambdaOptimizer<ReturnType>>::
2227 dynamic_cast<typename internal::LambdaOptimizer<
2229 print_independent_symbols,
2230 print_dependent_functions,
2231 print_cse_reductions);
2233 # ifdef HAVE_SYMENGINE_LLVM
2236 Assert(
dynamic_cast<typename internal::LLVMOptimizer<
2238 ExcMessage(
"Cannot cast optimizer to LLVM type."));
2240 internal::OptimizerHelper<ReturnType,
2241 internal::LLVMOptimizer<ReturnType>>::
2243 dynamic_cast<typename internal::LLVMOptimizer<
2245 print_independent_symbols,
2246 print_dependent_functions,
2247 print_cse_reductions);
2256 if (values_substituted())
2258 stream <<
"Evaluated functions:" <<
'\n';
2259 stream << std::flush;
2261 for (
typename std::vector<ReturnType>::const_iterator it =
2262 dependent_variables_output.begin();
2263 it != dependent_variables_output.end();
2266 stream << cntr <<
": " << (*it) <<
'\n';
2268 stream <<
'\n' << std::flush;
2274 template <
typename ReturnType>
2275 template <
class Archive>
2278 const unsigned int version)
const
2283 static_cast<typename std::underlying_type<OptimizerType>::type
>(
2289 static_cast<typename std::underlying_type<OptimizationFlags>::type
>(
2296 ar &independent_variables_symbols;
2297 ar &dependent_variables_functions;
2299 ar &dependent_variables_output;
2300 ar &map_dep_expr_vec_entry;
2301 ar &ready_for_value_extraction;
2304 has_been_serialized =
true;
2305 ar &has_been_serialized;
2314 *opt =
dynamic_cast<typename internal::DictionaryOptimizer<
2319 internal::OptimizerHelper<
2321 internal::DictionaryOptimizer<ReturnType>>::save(ar, version, opt);
2324 *opt =
dynamic_cast<typename internal::LambdaOptimizer<
2329 internal::OptimizerHelper<
2331 internal::LambdaOptimizer<ReturnType>>::save(ar, version, opt);
2333 # ifdef HAVE_SYMENGINE_LLVM
2335 *opt =
dynamic_cast<typename internal::LLVMOptimizer<
2340 internal::OptimizerHelper<
2342 internal::LLVMOptimizer<ReturnType>>::save(ar, version, opt);
2353 template <
typename ReturnType>
2354 template <
class Archive>
2366 typename std::underlying_type<OptimizerType>::type m;
2371 typename std::underlying_type<OptimizationFlags>::type f;
2378 ar &independent_variables_symbols;
2379 ar &dependent_variables_functions;
2381 ar &dependent_variables_output;
2382 ar &map_dep_expr_vec_entry;
2383 ar &ready_for_value_extraction;
2385 ar &has_been_serialized;
2392 create_optimizer(optimizer);
2402 *opt =
dynamic_cast<typename internal::DictionaryOptimizer<
2407 internal::OptimizerHelper<ReturnType,
2408 internal::DictionaryOptimizer<ReturnType>>::
2415 dependent_variables_functions),
2416 optimization_flags());
2419 *opt =
dynamic_cast<typename internal::LambdaOptimizer<
2424 internal::OptimizerHelper<ReturnType,
2425 internal::LambdaOptimizer<ReturnType>>::
2432 dependent_variables_functions),
2433 optimization_flags());
2435 # ifdef HAVE_SYMENGINE_LLVM
2437 *opt =
dynamic_cast<typename internal::LLVMOptimizer<
2442 internal::OptimizerHelper<ReturnType,
2443 internal::LLVMOptimizer<ReturnType>>::
2450 dependent_variables_functions),
2451 optimization_flags());
2462 template <
typename ReturnType>
2463 template <
int rank,
int dim>
2468 Assert(optimized() ==
false,
2470 "Cannot register functions once the optimizer is finalised."));
2472 register_vector_functions(
2478 template <
typename ReturnType>
2479 template <
int rank,
int dim>
2484 Assert(optimized() ==
false,
2486 "Cannot register functions once the optimizer is finalised."));
2488 register_vector_functions(
2494 template <
typename ReturnType>
2495 template <
typename T,
typename... Args>
2499 const Args &...other_functions)
2507 template <
typename ReturnType>
2508 template <
typename T>
2518 template <
typename ReturnType>
2519 template <
int rank,
int dim>
2523 const std::vector<ReturnType> & cached_evaluation)
const
2532 template <
typename ReturnType>
2533 template <
int rank,
int dim>
2539 values_substituted() ==
true,
2541 "The optimizer is not configured to perform evaluation. "
2542 "This action can only performed after substitute() has been called."));
2544 return extract(funcs, dependent_variables_output);
2549 template <
typename ReturnType>
2550 template <
int rank,
int dim>
2554 const std::vector<ReturnType> & cached_evaluation)
const
2563 template <
typename ReturnType>
2564 template <
int rank,
int dim>
2570 values_substituted() ==
true,
2572 "The optimizer is not configured to perform evaluation. "
2573 "This action can only performed after substitute() has been called."));
2575 return extract(funcs, dependent_variables_output);
bool use_symbolic_CSE() const
types::substitution_map independent_variables_symbols
types::symbol_vector dependent_variables_functions
void substitute(const types::substitution_map &substitution_map) const
Tensor< rank, dim, ReturnType > evaluate(const Tensor< rank, dim, Expression > &funcs) const
enum OptimizationFlags flags
void register_functions(const T &functions, const Args &...other_functions)
void register_scalar_function(const SD::Expression &function)
const types::symbol_vector & get_dependent_functions() const
void create_optimizer(std::unique_ptr< SymEngine::Visitor > &optimizer)
bool ready_for_value_extraction
Tensor< rank, dim, ReturnType > extract(const Tensor< rank, dim, Expression > &funcs, const std::vector< ReturnType > &cached_evaluation) const
void print(Stream &stream, const bool print_cse=false) const
enum OptimizerType optimization_method() const
void copy_from(const BatchOptimizer &other)
void register_function(const Tensor< rank, dim, Expression > &function_tensor)
void set_optimization_method(const enum OptimizerType &optimization_method, const enum OptimizationFlags &optimization_flags=OptimizationFlags::optimize_all)
void save(Archive &archive, const unsigned int version) const
enum OptimizationFlags optimization_flags() const
void register_functions(const types::symbol_vector &functions)
BatchOptimizer(BatchOptimizer &&)=default
SymmetricTensor< rank, dim, ReturnType > evaluate(const SymmetricTensor< rank, dim, Expression > &funcs) const
std::vector< ReturnType > dependent_variables_output
~BatchOptimizer()=default
enum OptimizerType method
std::size_t n_dependent_variables() const
void serialize(Archive &archive, const unsigned int version)
void register_symbols(const types::substitution_map &substitution_map)
const std::vector< ReturnType > & evaluate() const
std::map< SD::Expression, std::size_t, SD::types::internal::ExpressionKeyLess > map_dependent_expression_to_vector_entry_t
SymmetricTensor< rank, dim, ReturnType > extract(const SymmetricTensor< rank, dim, Expression > &funcs, const std::vector< ReturnType > &cached_evaluation) const
void register_function(const Expression &function)
std::unique_ptr< SymEngine::Visitor > optimizer
ReturnType extract(const Expression &func, const std::vector< ReturnType > &cached_evaluation) const
void load(Archive &archive, const unsigned int version)
void register_functions(const std::vector< T > &functions)
bool is_valid_nonunique_dependent_variable(const SD::Expression &function) const
void register_vector_functions(const types::symbol_vector &functions)
std::size_t n_independent_variables() const
types::symbol_vector get_independent_symbols() const
void register_function(const SymmetricTensor< rank, dim, Expression > &function_tensor)
map_dependent_expression_to_vector_entry_t map_dep_expr_vec_entry
bool values_substituted() const
static constexpr unsigned int n_independent_components
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcSymEngineLLVMNotAvailable()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcSymEngineLLVMReturnTypeNotSupported()
#define AssertThrow(cond, exc)
SD::types::symbol_vector extract_symbols(const SD::types::substitution_map &substitution_values)
SymEngine::vec_basic convert_expression_vector_to_basic_vector(const SD::types::symbol_vector &symbol_vector)
TensorType< rank, dim, NumberType > tensor_evaluate_optimized(const TensorType< rank, dim, Expression > &symbol_tensor, const std::vector< NumberType > &cached_evaluation, const BatchOptimizer< NumberType > &optimizer)
bool use_symbolic_CSE(const enum OptimizationFlags &flags)
types::symbol_vector unroll_to_expression_vector(const TensorType< rank, dim, Expression > &symbol_tensor)
int get_LLVM_optimization_level(const enum OptimizationFlags &flags)
void register_functions(BatchOptimizer< NumberType > &optimizer, const T &function)
std::vector< SD::Expression > symbol_vector
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
OptimizationFlags & operator|=(OptimizationFlags &f1, const OptimizationFlags f2)
Expression operator|(const Expression &lhs, const Expression &rhs)
Expression operator&(const Expression &lhs, const Expression &rhs)
OptimizationFlags & operator&=(OptimizationFlags &f1, const OptimizationFlags f2)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
std::ostream & operator<<(std::ostream &stream, const Expression &expression)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)