16 #ifndef dealii_differentiation_sd_symengine_optimizer_h
17 #define dealii_differentiation_sd_symengine_optimizer_h
21 #ifdef DEAL_II_WITH_SYMENGINE
24 # include <symengine/basic.h>
25 # include <symengine/dict.h>
26 # include <symengine/symengine_exception.h>
27 # include <symengine/symengine_rcp.h>
30 # include <symengine/lambda_double.h>
31 # include <symengine/visitor.h>
32 # ifdef HAVE_SYMENGINE_LLVM
33 # include <symengine/llvm_double.h>
46 # include <boost/serialization/split_member.hpp>
47 # include <boost/type_traits.hpp>
52 # include <type_traits>
74 "SymEngine has not been built with LLVM support.");
81 "The SymEngine LLVM optimizer does not (yet) support the "
82 "selected return type.");
88 template <
typename ReturnType>
119 template <
typename StreamType>
181 static_cast<unsigned int>(f2));
213 static_cast<unsigned int>(f2));
259 const bool use_agg_opt =
261 const int opt_level = (use_agg_opt ? 3 : 2);
271 template <
typename StreamType>
275 s <<
" OptimizationFlags|";
298 template <
typename ReturnType,
typename T =
void>
311 template <
typename ReturnType,
typename T =
void>
315 # ifdef HAVE_SYMENGINE_LLVM
325 template <
typename ReturnType,
typename T =
void>
326 struct LLVMOptimizer;
345 template <
typename ReturnType,
typename Optimizer,
typename T =
void>
357 template <
typename ReturnType_,
typename T =
void>
358 struct SupportedOptimizerTypeTraits
360 static const bool is_supported =
false;
362 using ReturnType = void;
368 template <
typename ReturnType_>
369 struct SupportedOptimizerTypeTraits<
371 std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
373 static const bool is_supported =
true;
375 using ReturnType =
typename std::
376 conditional<std::is_same_v<ReturnType_, float>, float,
double>::type;
382 template <
typename ReturnType_>
383 struct SupportedOptimizerTypeTraits<
386 boost::is_complex<ReturnType_>::value &&
387 std::is_arithmetic_v<typename ReturnType_::value_type>>>
389 static const bool is_supported =
true;
391 using ReturnType =
typename std::conditional<
392 std::is_same_v<ReturnType_, std::complex<float>>,
394 std::complex<double>>::type;
399 template <
typename ReturnType_>
400 struct DictionaryOptimizer<ReturnType_,
401 std::enable_if_t<SupportedOptimizerTypeTraits<
402 ReturnType_>::is_supported>>
405 typename SupportedOptimizerTypeTraits<ReturnType_>::ReturnType;
407 internal::DictionarySubstitutionVisitor<ReturnType, SD::Expression>;
420 const SymEngine::vec_basic &independent_symbols,
421 const SymEngine::vec_basic &dependent_functions,
425 optimizer.init(independent_symbols,
436 template <
class Archive>
438 save(Archive &archive,
439 const unsigned int version,
442 optimizer.save(archive, version);
451 template <
class Archive>
453 load(Archive &archive,
454 const unsigned int version,
456 const SymEngine::vec_basic & ,
457 const SymEngine::vec_basic & ,
460 optimizer.load(archive, version);
480 template <
typename Stream>
482 print(Stream &stream,
484 const bool print_independent_symbols =
false,
485 const bool print_dependent_functions =
false,
486 const bool print_cse_reductions =
true)
488 optimizer.print(stream,
489 print_independent_symbols,
490 print_dependent_functions,
491 print_cse_reductions);
497 template <
typename ReturnType_>
498 struct LambdaOptimizer<ReturnType_,
499 std::enable_if_t<SupportedOptimizerTypeTraits<
500 ReturnType_>::is_supported>>
503 typename std::conditional<!boost::is_complex<ReturnType_>::value,
505 std::complex<double>>::type;
507 !boost::is_complex<ReturnType_>::value,
508 SymEngine::LambdaRealDoubleVisitor,
509 SymEngine::LambdaComplexDoubleVisitor>::type;
522 const SymEngine::vec_basic &independent_symbols,
523 const SymEngine::vec_basic &dependent_functions,
527 optimizer.init(independent_symbols,
538 template <
class Archive>
550 template <
class Archive>
555 const SymEngine::vec_basic &independent_symbols,
556 const SymEngine::vec_basic &dependent_functions,
559 initialize(optimizer,
582 template <
typename StreamType>
596 # ifdef HAVE_SYMENGINE_LLVM
597 template <
typename ReturnType_>
598 struct LLVMOptimizer<ReturnType_,
599 std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
601 using ReturnType =
typename std::
602 conditional<std::is_same_v<ReturnType_, float>, float,
double>::type;
604 typename std::conditional<std::is_same_v<ReturnType_, float>,
605 SymEngine::LLVMFloatVisitor,
606 SymEngine::LLVMDoubleVisitor>::type;
612 static const bool supported_by_LLVM =
true;
625 const SymEngine::vec_basic &independent_symbols,
626 const SymEngine::vec_basic &dependent_functions,
631 optimizer.init(independent_symbols,
643 template <
class Archive>
645 save(Archive &archive,
649 const std::string llvm_compiled_function = optimizer.dumps();
650 archive &llvm_compiled_function;
659 template <
class Archive>
661 load(Archive &archive,
664 const SymEngine::vec_basic & ,
665 const SymEngine::vec_basic & ,
668 std::string llvm_compiled_function;
669 archive &llvm_compiled_function;
670 optimizer.loads(llvm_compiled_function);
690 template <
typename StreamType>
708 template <
typename ReturnType_>
709 struct LLVMOptimizer<
712 boost::is_complex<ReturnType_>::value &&
713 std::is_arithmetic_v<typename ReturnType_::value_type>>>
717 using ReturnType =
typename LambdaOptimizer<ReturnType_>::ReturnType;
725 static const bool supported_by_LLVM =
false;
738 const SymEngine::vec_basic & ,
739 const SymEngine::vec_basic & ,
751 template <
class Archive>
766 template <
class Archive>
771 const SymEngine::vec_basic & ,
772 const SymEngine::vec_basic & ,
795 template <
typename StreamType>
812 template <
typename ReturnType,
typename Optimizer>
813 struct OptimizerHelper<
817 std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
829 const SymEngine::vec_basic &independent_symbols,
830 const SymEngine::vec_basic &dependent_functions,
838 Optimizer::initialize(*optimizer,
861 std::vector<ReturnType> &output_values,
862 const std::vector<ReturnType> &substitution_values)
865 optimizer->call(output_values.data(), substitution_values.data());
874 template <
class Archive>
876 save(Archive &archive,
877 const unsigned int version,
885 Optimizer::save(archive, version, *optimizer);
894 template <
class Archive>
896 load(Archive &archive,
897 const unsigned int version,
899 const SymEngine::vec_basic &independent_symbols,
900 const SymEngine::vec_basic &dependent_functions,
908 Optimizer::load(archive,
933 template <
typename Stream>
935 print(Stream &stream,
937 const bool print_independent_symbols =
false,
938 const bool print_dependent_functions =
false,
939 const bool print_cse_reductions =
true)
946 Optimizer::print(stream,
948 print_independent_symbols,
949 print_dependent_functions,
950 print_cse_reductions);
954 template <
typename ReturnType,
typename Optimizer>
955 struct OptimizerHelper<
959 !std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
971 const SymEngine::vec_basic &independent_symbols,
972 const SymEngine::vec_basic &dependent_functions,
978 optimizer->init(independent_symbols,
1000 std::vector<ReturnType> &output_values,
1001 const std::vector<ReturnType> &substitution_values)
1007 std::vector<typename Optimizer::ReturnType> int_outputs(
1008 output_values.size());
1009 std::vector<typename Optimizer::ReturnType> int_inputs(
1010 substitution_values.size());
1013 substitution_values.end(),
1014 int_inputs.begin());
1015 optimizer->call(int_outputs.data(), int_inputs.data());
1018 output_values.begin());
1027 template <
class Archive>
1029 save(Archive &archive,
1030 const unsigned int version,
1034 Optimizer::save(archive, version, *optimizer);
1043 template <
class Archive>
1045 load(Archive &archive,
1046 const unsigned int version,
1048 const SymEngine::vec_basic &independent_symbols,
1049 const SymEngine::vec_basic &dependent_functions,
1057 Optimizer::load(archive,
1060 independent_symbols,
1061 dependent_functions,
1062 optimization_flags);
1082 template <
typename Stream>
1084 print(Stream &stream,
1086 const bool print_cse_reductions =
true,
1087 const bool print_independent_symbols =
false,
1088 const bool print_dependent_functions =
false)
1092 optimizer->print(stream,
1093 print_independent_symbols,
1094 print_dependent_functions,
1095 print_cse_reductions);
1126 template <
typename NumberType,
1129 template <
int,
int,
typename>
1131 TensorType<rank, dim, NumberType>
1133 const TensorType<rank, dim, Expression> &symbol_tensor,
1134 const std::vector<NumberType> &cached_evaluation,
1137 TensorType<rank, dim, NumberType> out;
1138 for (
unsigned int i = 0; i < out.n_independent_components; ++i)
1141 out.unrolled_to_component_indices(i));
1143 optimizer.
extract(symbol_tensor[indices], cached_evaluation);
1171 template <
typename NumberType,
int dim>
1175 const std::vector<NumberType> &cached_evaluation,
1179 for (
unsigned int i = 0;
1180 i < SymmetricTensor<2, dim>::n_independent_components;
1182 for (
unsigned int j = 0;
1183 j < SymmetricTensor<2, dim>::n_independent_components;
1187 make_rank_4_tensor_indices<dim>(i, j);
1189 optimizer.
extract(symbol_tensor[indices], cached_evaluation);
1212 template <
typename NumberType,
typename T>
1238 template <
typename NumberType,
typename T>
1267 template <
typename NumberType,
typename T,
typename... Args>
1271 const Args &...other_functions)
1291 template <
int,
int,
typename>
1295 const TensorType<rank, dim, Expression> &symbol_tensor)
1298 out.reserve(symbol_tensor.n_independent_components);
1299 for (
unsigned int i = 0; i < symbol_tensor.n_independent_components;
1303 symbol_tensor.unrolled_to_component_indices(i));
1304 out.push_back(symbol_tensor[indices].get_RCP());
1326 for (
unsigned int i = 0;
1327 i < SymmetricTensor<2, dim>::n_independent_components;
1329 for (
unsigned int j = 0;
1330 j < SymmetricTensor<2, dim>::n_independent_components;
1334 make_rank_4_tensor_indices<dim>(i, j);
1335 out.push_back(symbol_tensor[indices].get_RCP());
1434 template <
typename ReturnType>
1511 template <typename Stream>
1513 print(Stream &stream, const
bool print_cse = false) const;
1523 template <class Archive>
1525 save(Archive &archive, const
unsigned int version) const;
1540 template <class Archive>
1542 load(Archive &archive, const
unsigned int version);
1566 template <
class Archive>
1572 BOOST_SERIALIZATION_SPLIT_MEMBER()
1655 template <
int rank,
int dim>
1663 template <
int rank,
int dim>
1691 template <
typename T>
1708 template <
typename T,
typename... Args>
1841 const std::vector<ReturnType> &
values)
const;
1854 substitute(
const SymEngine::vec_basic &symbols,
1855 const std::vector<ReturnType> &
values)
const;
1894 const std::vector<ReturnType> &
1915 std::vector<ReturnType>
1916 evaluate(
const std::vector<Expression> &funcs)
const;
1926 template <
int rank,
int dim>
1939 template <
int rank,
int dim>
1953 const std::vector<ReturnType> &cached_evaluation)
const;
1963 std::vector<ReturnType>
1964 extract(
const std::vector<Expression> &funcs,
1965 const std::vector<ReturnType> &cached_evaluation)
const;
1975 template <
int rank,
int dim>
1978 const std::vector<ReturnType> &cached_evaluation)
const;
1988 template <
int rank,
int dim>
1991 const std::vector<ReturnType> &cached_evaluation)
const;
2039 const SymEngine::RCP<const SymEngine::Basic> &
function)
const;
2136 substitute(
const std::vector<ReturnType> &substitution_values)
const;
2147 template <
typename ReturnType>
2148 template <
typename Stream>
2154 stream <<
"Method? " << optimization_method() <<
'\n';
2155 stream <<
"Flags: " << optimization_flags() <<
'\n';
2156 stream <<
"Optimized? " << (optimized() ?
"Yes" :
"No") <<
'\n';
2157 stream <<
"Values substituted? " << values_substituted() <<
"\n\n";
2160 stream <<
"Symbols (" << n_independent_variables()
2161 <<
" independent variables):" <<
'\n';
2163 for (SD::types::substitution_map::const_iterator it =
2164 independent_variables_symbols.begin();
2165 it != independent_variables_symbols.end();
2168 stream << cntr <<
": " << it->first <<
'\n';
2170 stream <<
'\n' << std::flush;
2173 stream <<
"Functions (" << n_dependent_variables()
2174 <<
" dependent variables):" <<
'\n';
2176 for (
typename SD::types::symbol_vector::const_iterator it =
2177 dependent_variables_functions.begin();
2178 it != dependent_variables_functions.end();
2181 stream << cntr <<
": " << (*it) <<
'\n';
2183 stream <<
'\n' << std::flush;
2189 const bool print_cse_reductions =
true;
2190 const bool print_independent_symbols =
false;
2191 const bool print_dependent_functions =
false;
2195 Assert(
dynamic_cast<typename internal::DictionaryOptimizer<
2197 ExcMessage(
"Cannot cast optimizer to Dictionary type."));
2199 internal::OptimizerHelper<
2201 internal::DictionaryOptimizer<ReturnType>>::
2203 dynamic_cast<typename internal::DictionaryOptimizer<
2205 print_independent_symbols,
2206 print_dependent_functions,
2207 print_cse_reductions);
2209 stream <<
'\n' << std::flush;
2213 Assert(
dynamic_cast<typename internal::LambdaOptimizer<
2215 ExcMessage(
"Cannot cast optimizer to Lambda type."));
2217 internal::OptimizerHelper<ReturnType,
2218 internal::LambdaOptimizer<ReturnType>>::
2220 dynamic_cast<typename internal::LambdaOptimizer<
2222 print_independent_symbols,
2223 print_dependent_functions,
2224 print_cse_reductions);
2226 # ifdef HAVE_SYMENGINE_LLVM
2229 Assert(
dynamic_cast<typename internal::LLVMOptimizer<
2231 ExcMessage(
"Cannot cast optimizer to LLVM type."));
2233 internal::OptimizerHelper<ReturnType,
2234 internal::LLVMOptimizer<ReturnType>>::
2236 dynamic_cast<typename internal::LLVMOptimizer<
2238 print_independent_symbols,
2239 print_dependent_functions,
2240 print_cse_reductions);
2249 if (values_substituted())
2251 stream <<
"Evaluated functions:" <<
'\n';
2252 stream << std::flush;
2254 for (
typename std::vector<ReturnType>::const_iterator it =
2255 dependent_variables_output.begin();
2256 it != dependent_variables_output.end();
2259 stream << cntr <<
": " << (*it) <<
'\n';
2261 stream <<
'\n' << std::flush;
2267 template <
typename ReturnType>
2268 template <
class Archive>
2271 const unsigned int version)
const
2276 static_cast<typename std::underlying_type<OptimizerType>::type
>(
2282 static_cast<typename std::underlying_type<OptimizationFlags>::type
>(
2289 ar &independent_variables_symbols;
2290 ar &dependent_variables_functions;
2292 ar &dependent_variables_output;
2293 ar &map_dep_expr_vec_entry;
2294 ar &ready_for_value_extraction;
2297 has_been_serialized =
true;
2298 ar &has_been_serialized;
2307 *opt =
dynamic_cast<typename internal::DictionaryOptimizer<
2312 internal::OptimizerHelper<
2314 internal::DictionaryOptimizer<ReturnType>>::save(ar, version, opt);
2317 *opt =
dynamic_cast<typename internal::LambdaOptimizer<
2322 internal::OptimizerHelper<
2324 internal::LambdaOptimizer<ReturnType>>::save(ar, version, opt);
2326 # ifdef HAVE_SYMENGINE_LLVM
2328 *opt =
dynamic_cast<typename internal::LLVMOptimizer<
2333 internal::OptimizerHelper<
2335 internal::LLVMOptimizer<ReturnType>>::save(ar, version, opt);
2346 template <
typename ReturnType>
2347 template <
class Archive>
2359 typename std::underlying_type<OptimizerType>::type m;
2364 typename std::underlying_type<OptimizationFlags>::type f;
2371 ar &independent_variables_symbols;
2372 ar &dependent_variables_functions;
2374 ar &dependent_variables_output;
2375 ar &map_dep_expr_vec_entry;
2376 ar &ready_for_value_extraction;
2378 ar &has_been_serialized;
2385 create_optimizer(optimizer);
2395 *opt =
dynamic_cast<typename internal::DictionaryOptimizer<
2400 internal::OptimizerHelper<ReturnType,
2401 internal::DictionaryOptimizer<ReturnType>>::
2408 dependent_variables_functions),
2409 optimization_flags());
2412 *opt =
dynamic_cast<typename internal::LambdaOptimizer<
2417 internal::OptimizerHelper<ReturnType,
2418 internal::LambdaOptimizer<ReturnType>>::
2425 dependent_variables_functions),
2426 optimization_flags());
2428 # ifdef HAVE_SYMENGINE_LLVM
2430 *opt =
dynamic_cast<typename internal::LLVMOptimizer<
2435 internal::OptimizerHelper<ReturnType,
2436 internal::LLVMOptimizer<ReturnType>>::
2443 dependent_variables_functions),
2444 optimization_flags());
2455 template <
typename ReturnType>
2456 template <
int rank,
int dim>
2461 Assert(optimized() ==
false,
2463 "Cannot register functions once the optimizer is finalised."));
2465 register_vector_functions(
2471 template <
typename ReturnType>
2472 template <
int rank,
int dim>
2477 Assert(optimized() ==
false,
2479 "Cannot register functions once the optimizer is finalised."));
2481 register_vector_functions(
2487 template <
typename ReturnType>
2488 template <
typename T,
typename... Args>
2492 const Args &...other_functions)
2500 template <
typename ReturnType>
2501 template <
typename T>
2511 template <
typename ReturnType>
2512 template <
int rank,
int dim>
2516 const std::vector<ReturnType> &cached_evaluation)
const
2525 template <
typename ReturnType>
2526 template <
int rank,
int dim>
2532 values_substituted() ==
true,
2534 "The optimizer is not configured to perform evaluation. "
2535 "This action can only performed after substitute() has been called."));
2537 return extract(funcs, dependent_variables_output);
2542 template <
typename ReturnType>
2543 template <
int rank,
int dim>
2547 const std::vector<ReturnType> &cached_evaluation)
const
2556 template <
typename ReturnType>
2557 template <
int rank,
int dim>
2563 values_substituted() ==
true,
2565 "The optimizer is not configured to perform evaluation. "
2566 "This action can only performed after substitute() has been called."));
2568 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)
SymmetricTensor< rank, dim, ReturnType > evaluate(const SymmetricTensor< rank, dim, Expression > &funcs) const
std::vector< ReturnType > dependent_variables_output
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
BatchOptimizer(BatchOptimizer &&) noexcept=default
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_NAMESPACE_CLOSE
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)