Reference documentation for deal.II version GIT relicensing-245-g36f19064f7 2024-03-29 07:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
symengine_optimizer.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_differentiation_sd_symengine_optimizer_h
16#define dealii_differentiation_sd_symengine_optimizer_h
17
18#include <deal.II/base/config.h>
19
20#ifdef DEAL_II_WITH_SYMENGINE
21
22// Low level
23# include <symengine/basic.h>
24# include <symengine/dict.h>
25# include <symengine/symengine_exception.h>
26# include <symengine/symengine_rcp.h>
27
28// Optimization
29# include <symengine/lambda_double.h>
30# include <symengine/visitor.h>
31# ifdef HAVE_SYMENGINE_LLVM
32# include <symengine/llvm_double.h>
33# endif
34
38
44
45# include <boost/serialization/split_member.hpp>
46# include <boost/type_traits.hpp>
47
48# include <algorithm>
49# include <map>
50# include <memory>
51# include <type_traits>
52# include <utility>
53# include <vector>
54
55
57
58
59namespace Differentiation
60{
61 namespace SD
62 {
73 "SymEngine has not been built with LLVM support.");
74
80 "The SymEngine LLVM optimizer does not (yet) support the "
81 "selected return type.");
82
86 // Forward declarations
87 template <typename ReturnType>
88 class BatchOptimizer;
89
90
96 enum class OptimizerType
97 {
106 lambda,
111 llvm
112 };
113
114
118 template <typename StreamType>
119 inline StreamType &
120 operator<<(StreamType &s, OptimizerType o)
121 {
123 s << "dictionary";
124 else if (o == OptimizerType::lambda)
125 s << "lambda";
126 else if (o == OptimizerType::llvm)
127 s << "llvm";
128 else
129 {
130 Assert(false, ExcMessage("Unknown optimization method."));
131 }
132
133 return s;
134 }
135
136
142 enum class OptimizationFlags : unsigned char
143 {
151 optimize_cse = 0x0001,
156 optimize_aggressive = 0x0002,
161 };
162
163
172 // This operator exists since if it did not then the result of the bit-or
173 // <tt>operator |</tt> would be an integer which would in turn trigger a
174 // compiler warning when we tried to assign it to an object of type
175 // OptimizationFlags.
176 inline OptimizationFlags
178 {
179 return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) |
180 static_cast<unsigned int>(f2));
181 }
182
183
188 inline OptimizationFlags &
190 {
191 f1 = f1 | f2;
192 return f1;
193 }
194
195
204 // This operator exists since if it did not then the result of the bit-or
205 // <tt>operator |</tt> would be an integer which would in turn trigger a
206 // compiler warning when we tried to assign it to an object of type
207 // OptimizationFlags.
208 inline OptimizationFlags
210 {
211 return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) &
212 static_cast<unsigned int>(f2));
213 }
214
215
220 inline OptimizationFlags &
222 {
223 f1 = f1 & f2;
224 return f1;
225 }
226
227
228 namespace internal
229 {
234 inline bool
236 {
237 return static_cast<int>(flags & OptimizationFlags::optimize_cse);
238 }
239
244 inline int
246 {
247 // With the LLVM compiler there exists the opportunity to tune
248 // the level of optimizations performed during compilation.
249 // By default SymEngine sets this at "opt_level=2", which one
250 // presumes targets -O2. Here we are a bit more specific about
251 // want we want it to do:
252 // - Normal compilation: -02 (default settings)
253 // - Aggressive mode: -03 (the whole lot!)
254 // In theory we could also target
255 // - Debug mode: -O0 (no optimizations)
256 // but this doesn't make much sense since SymEngine is a
257 // tested external library.
258 const bool use_agg_opt =
259 static_cast<int>(flags & OptimizationFlags::optimize_aggressive);
260 const int opt_level = (use_agg_opt ? 3 : 2);
261 return opt_level;
262 }
263 } // namespace internal
264
265
270 template <typename StreamType>
271 inline StreamType &
272 operator<<(StreamType &s, OptimizationFlags o)
273 {
274 s << " OptimizationFlags|";
275 if (static_cast<unsigned int>(o & OptimizationFlags::optimize_cse))
276 s << "cse|";
277
278 // LLVM optimization level
279 s << "-O" + std::to_string(internal::get_LLVM_optimization_level(o)) +
280 "|";
281
282 return s;
283 }
284
285
286 namespace internal
287 {
297 template <typename ReturnType, typename T = void>
299
300
310 template <typename ReturnType, typename T = void>
312
313
314# ifdef HAVE_SYMENGINE_LLVM
324 template <typename ReturnType, typename T = void>
325 struct LLVMOptimizer;
326# endif // HAVE_SYMENGINE_LLVM
327
328
344 template <typename ReturnType, typename Optimizer, typename T = void>
346
347
348# ifndef DOXYGEN
349
350
351 /* ----------- Specializations for the Optimizers ----------- */
352
353
354 // A helper struct to type trait detection for the optimizers that
355 // will be defined next.
356 template <typename ReturnType_, typename T = void>
357 struct SupportedOptimizerTypeTraits
358 {
359 static const bool is_supported = false;
360
361 using ReturnType = void;
362 };
363
364
365
366 // Specialization for arithmetic types
367 template <typename ReturnType_>
368 struct SupportedOptimizerTypeTraits<
369 ReturnType_,
370 std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
371 {
372 static const bool is_supported = true;
373
374 using ReturnType = typename std::
375 conditional<std::is_same_v<ReturnType_, float>, float, double>::type;
376 };
377
378
379
380 // Specialization for complex arithmetic types
381 template <typename ReturnType_>
382 struct SupportedOptimizerTypeTraits<
383 ReturnType_,
384 std::enable_if_t<
385 boost::is_complex<ReturnType_>::value &&
386 std::is_arithmetic_v<typename ReturnType_::value_type>>>
387 {
388 static const bool is_supported = true;
389
390 using ReturnType =
391 std::conditional_t<std::is_same_v<ReturnType_, std::complex<float>>,
392 std::complex<float>,
393 std::complex<double>>;
394 };
395
396
397
398 template <typename ReturnType_>
399 struct DictionaryOptimizer<ReturnType_,
400 std::enable_if_t<SupportedOptimizerTypeTraits<
401 ReturnType_>::is_supported>>
402 {
403 using ReturnType =
404 typename SupportedOptimizerTypeTraits<ReturnType_>::ReturnType;
405 using OptimizerType =
406 internal::DictionarySubstitutionVisitor<ReturnType, SD::Expression>;
407
408
417 static void
418 initialize(OptimizerType &optimizer,
419 const SymEngine::vec_basic &independent_symbols,
420 const SymEngine::vec_basic &dependent_functions,
421 const enum OptimizationFlags &optimization_flags)
422 {
423 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
424 optimizer.init(independent_symbols,
425 dependent_functions,
426 use_symbolic_cse);
427 }
428
429
430
435 template <class Archive>
436 static void
437 save(Archive &archive,
438 const unsigned int version,
439 OptimizerType &optimizer)
440 {
441 optimizer.save(archive, version);
442 }
443
444
445
450 template <class Archive>
451 static void
452 load(Archive &archive,
453 const unsigned int version,
454 OptimizerType &optimizer,
455 const SymEngine::vec_basic & /*independent_symbols*/,
456 const SymEngine::vec_basic & /*dependent_functions*/,
457 const enum OptimizationFlags & /*optimization_flags*/)
458 {
459 optimizer.load(archive, version);
460 }
461
462
463
479 template <typename Stream>
480 static void
481 print(Stream &stream,
482 const OptimizerType &optimizer,
483 const bool print_independent_symbols = false,
484 const bool print_dependent_functions = false,
485 const bool print_cse_reductions = true)
486 {
487 optimizer.print(stream,
488 print_independent_symbols,
489 print_dependent_functions,
490 print_cse_reductions);
491 }
492 };
493
494
495
496 template <typename ReturnType_>
497 struct LambdaOptimizer<ReturnType_,
498 std::enable_if_t<SupportedOptimizerTypeTraits<
499 ReturnType_>::is_supported>>
500 {
501 using ReturnType =
502 std::conditional_t<!boost::is_complex<ReturnType_>::value,
503 double,
504 std::complex<double>>;
505 using OptimizerType =
506 std::conditional_t<!boost::is_complex<ReturnType_>::value,
507 SymEngine::LambdaRealDoubleVisitor,
508 SymEngine::LambdaComplexDoubleVisitor>;
509
510
519 static void
520 initialize(OptimizerType &optimizer,
521 const SymEngine::vec_basic &independent_symbols,
522 const SymEngine::vec_basic &dependent_functions,
523 const enum OptimizationFlags &optimization_flags)
524 {
525 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
526 optimizer.init(independent_symbols,
527 dependent_functions,
528 use_symbolic_cse);
529 }
530
531
532
537 template <class Archive>
538 static void
539 save(Archive & /*archive*/,
540 const unsigned int /*version*/,
541 OptimizerType & /*optimizer*/)
542 {}
543
544
549 template <class Archive>
550 static void
551 load(Archive & /*archive*/,
552 const unsigned int /*version*/,
553 OptimizerType &optimizer,
554 const SymEngine::vec_basic &independent_symbols,
555 const SymEngine::vec_basic &dependent_functions,
556 const enum OptimizationFlags &optimization_flags)
557 {
558 initialize(optimizer,
559 independent_symbols,
560 dependent_functions,
561 optimization_flags);
562 }
563
564
565
581 template <typename StreamType>
582 static void
583 print(StreamType & /*stream*/,
584 const OptimizerType & /*optimizer*/,
585 const bool /*print_independent_symbols*/ = false,
586 const bool /*print_dependent_functions*/ = false,
587 const bool /*print_cse_reductions*/ = true)
588 {
589 // No built-in print function
590 }
591 };
592
593
594
595# ifdef HAVE_SYMENGINE_LLVM
596 template <typename ReturnType_>
597 struct LLVMOptimizer<ReturnType_,
598 std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
599 {
600 using ReturnType = typename std::
601 conditional<std::is_same_v<ReturnType_, float>, float, double>::type;
602 using OptimizerType =
603 std::conditional_t<std::is_same_v<ReturnType_, float>,
604 SymEngine::LLVMFloatVisitor,
605 SymEngine::LLVMDoubleVisitor>;
606
611 static const bool supported_by_LLVM = true;
612
613
622 static void
623 initialize(OptimizerType &optimizer,
624 const SymEngine::vec_basic &independent_symbols,
625 const SymEngine::vec_basic &dependent_functions,
626 const enum OptimizationFlags &optimization_flags)
627 {
628 const int opt_level = get_LLVM_optimization_level(optimization_flags);
629 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
630 optimizer.init(independent_symbols,
631 dependent_functions,
632 use_symbolic_cse,
633 opt_level);
634 }
635
636
637
642 template <class Archive>
643 static void
644 save(Archive &archive,
645 const unsigned int /*version*/,
646 OptimizerType &optimizer)
647 {
648 const std::string llvm_compiled_function = optimizer.dumps();
649 archive &llvm_compiled_function;
650 }
651
652
653
658 template <class Archive>
659 static void
660 load(Archive &archive,
661 const unsigned int /*version*/,
662 OptimizerType &optimizer,
663 const SymEngine::vec_basic & /*independent_symbols*/,
664 const SymEngine::vec_basic & /*dependent_functions*/,
665 const enum OptimizationFlags & /*optimization_flags*/)
666 {
667 std::string llvm_compiled_function;
668 archive &llvm_compiled_function;
669 optimizer.loads(llvm_compiled_function);
670 }
671
672
673
689 template <typename StreamType>
690 static void
691 print(StreamType & /*stream*/,
692 const OptimizerType & /*optimizer*/,
693 const bool /*print_independent_symbols*/ = false,
694 const bool /*print_dependent_functions*/ = false,
695 const bool /*print_cse_reductions*/ = true)
696 {
697 // No built-in print function
698 }
699 };
700
701
702 // There is no LLVM optimizer built with complex number support.
703 // So we fall back to the LambdaDouble case as a type (required
704 // at compile time), but offer no implementation. We expect that
705 // the calling class does not create this type: This can be done by
706 // checking the `supported_by_LLVM` flag.
707 template <typename ReturnType_>
708 struct LLVMOptimizer<
709 ReturnType_,
710 std::enable_if_t<
711 boost::is_complex<ReturnType_>::value &&
712 std::is_arithmetic_v<typename ReturnType_::value_type>>>
713 {
714 // Since there is no working implementation, these are dummy types
715 // that help with templating in the calling function.
716 using ReturnType = typename LambdaOptimizer<ReturnType_>::ReturnType;
717 using OptimizerType =
718 typename LambdaOptimizer<ReturnType_>::OptimizerType;
719
724 static const bool supported_by_LLVM = false;
725
726
735 static void
736 initialize(OptimizerType & /*optimizer*/,
737 const SymEngine::vec_basic & /*independent_symbols*/,
738 const SymEngine::vec_basic & /*dependent_functions*/,
739 const enum OptimizationFlags & /*optimization_flags*/)
740 {
742 }
743
744
745
750 template <class Archive>
751 static void
752 save(Archive & /*archive*/,
753 const unsigned int /*version*/,
754 OptimizerType & /*optimizer*/)
755 {
757 }
758
759
760
765 template <class Archive>
766 static void
767 load(Archive & /*archive*/,
768 const unsigned int /*version*/,
769 OptimizerType & /*optimizer*/,
770 const SymEngine::vec_basic & /*independent_symbols*/,
771 const SymEngine::vec_basic & /*dependent_functions*/,
772 const enum OptimizationFlags & /*optimization_flags*/)
773 {
775 }
776
777
778
794 template <typename StreamType>
795 static void
796 print(StreamType & /*stream*/,
797 const OptimizerType & /*optimizer*/,
798 const bool /*print_independent_symbols*/ = false,
799 const bool /*print_dependent_functions*/ = false,
800 const bool /*print_cse_reductions*/ = true)
801 {
803 }
804 };
805# endif // HAVE_SYMENGINE_LLVM
806
807
808 /* ----------- Specializations for OptimizerHelper ----------- */
809
810
811 template <typename ReturnType, typename Optimizer>
812 struct OptimizerHelper<
813 ReturnType,
814 Optimizer,
815 std::enable_if_t<
816 std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
817 {
826 static void
827 initialize(typename Optimizer::OptimizerType *optimizer,
828 const SymEngine::vec_basic &independent_symbols,
829 const SymEngine::vec_basic &dependent_functions,
830 const enum OptimizationFlags &optimization_flags)
831 {
832 Assert(optimizer, ExcNotInitialized());
833
834 // Some optimizers don't have the same interface for
835 // initialization, we filter them out through the specializations
836 // of the Optimizer class
837 Optimizer::initialize(*optimizer,
838 independent_symbols,
839 dependent_functions,
840 optimization_flags);
841 }
842
843
844
858 static void
859 substitute(typename Optimizer::OptimizerType *optimizer,
860 std::vector<ReturnType> &output_values,
861 const std::vector<ReturnType> &substitution_values)
862 {
863 Assert(optimizer, ExcNotInitialized());
864 optimizer->call(output_values.data(), substitution_values.data());
865 }
866
867
868
873 template <class Archive>
874 static void
875 save(Archive &archive,
876 const unsigned int version,
877 typename Optimizer::OptimizerType *optimizer)
878 {
879 Assert(optimizer, ExcNotInitialized());
880
881 // Some optimizers don't have the same interface for
882 // serialization, we filter them out through the specializations
883 // of the Optimizer class
884 Optimizer::save(archive, version, *optimizer);
885 }
886
887
888
893 template <class Archive>
894 static void
895 load(Archive &archive,
896 const unsigned int version,
897 typename Optimizer::OptimizerType *optimizer,
898 const SymEngine::vec_basic &independent_symbols,
899 const SymEngine::vec_basic &dependent_functions,
900 const enum OptimizationFlags &optimization_flags)
901 {
902 Assert(optimizer, ExcNotInitialized());
903
904 // Some optimizers don't have the same interface for
905 // serialization, we filter them out through the specializations
906 // of the Optimizer class
907 Optimizer::load(archive,
908 version,
909 *optimizer,
910 independent_symbols,
911 dependent_functions,
912 optimization_flags);
913 }
914
915
916
932 template <typename Stream>
933 static void
934 print(Stream &stream,
935 typename Optimizer::OptimizerType *optimizer,
936 const bool print_independent_symbols = false,
937 const bool print_dependent_functions = false,
938 const bool print_cse_reductions = true)
939 {
940 Assert(optimizer, ExcNotInitialized());
941
942 // Some optimizers don't have a print function, so
943 // we filter them out through the specializations of
944 // the Optimizer class
945 Optimizer::print(stream,
946 *optimizer,
947 print_independent_symbols,
948 print_dependent_functions,
949 print_cse_reductions);
950 }
951 };
952
953 template <typename ReturnType, typename Optimizer>
954 struct OptimizerHelper<
955 ReturnType,
956 Optimizer,
957 std::enable_if_t<
958 !std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
959 {
968 static void
969 initialize(typename Optimizer::OptimizerType *optimizer,
970 const SymEngine::vec_basic &independent_symbols,
971 const SymEngine::vec_basic &dependent_functions,
972 const enum OptimizationFlags &optimization_flags)
973 {
974 Assert(optimizer, ExcNotInitialized());
975
976 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
977 optimizer->init(independent_symbols,
978 dependent_functions,
979 use_symbolic_cse);
980 }
981
982
983
997 static void
998 substitute(typename Optimizer::OptimizerType *optimizer,
999 std::vector<ReturnType> &output_values,
1000 const std::vector<ReturnType> &substitution_values)
1001 {
1002 Assert(optimizer, ExcNotInitialized());
1003
1004 // Intermediate values to accommodate the difference in
1005 // value types.
1006 std::vector<typename Optimizer::ReturnType> int_outputs(
1007 output_values.size());
1008 std::vector<typename Optimizer::ReturnType> int_inputs(
1009 substitution_values.size());
1010
1011 std::copy(substitution_values.begin(),
1012 substitution_values.end(),
1013 int_inputs.begin());
1014 optimizer->call(int_outputs.data(), int_inputs.data());
1015 std::copy(int_outputs.begin(),
1016 int_outputs.end(),
1017 output_values.begin());
1018 }
1019
1020
1021
1026 template <class Archive>
1027 static void
1028 save(Archive &archive,
1029 const unsigned int version,
1030 typename Optimizer::OptimizerType *optimizer)
1031 {
1032 Assert(optimizer, ExcNotInitialized());
1033 Optimizer::save(archive, version, *optimizer);
1034 }
1035
1036
1037
1042 template <class Archive>
1043 static void
1044 load(Archive &archive,
1045 const unsigned int version,
1046 typename Optimizer::OptimizerType *optimizer,
1047 const SymEngine::vec_basic &independent_symbols,
1048 const SymEngine::vec_basic &dependent_functions,
1049 const enum OptimizationFlags &optimization_flags)
1050 {
1051 Assert(optimizer, ExcNotInitialized());
1052
1053 // Some optimizers don't have the same interface for
1054 // serialization, we filter them out through the specializations
1055 // of the Optimizer class
1056 Optimizer::load(archive,
1057 version,
1058 *optimizer,
1059 independent_symbols,
1060 dependent_functions,
1061 optimization_flags);
1062 }
1063
1064
1065
1081 template <typename Stream>
1082 static void
1083 print(Stream &stream,
1084 typename Optimizer::OptimizerType *optimizer,
1085 const bool print_cse_reductions = true,
1086 const bool print_independent_symbols = false,
1087 const bool print_dependent_functions = false)
1088 {
1089 Assert(optimizer, ExcNotInitialized());
1090
1091 optimizer->print(stream,
1092 print_independent_symbols,
1093 print_dependent_functions,
1094 print_cse_reductions);
1095 }
1096 };
1097
1098# endif // DOXYGEN
1099
1100
1101 /* -------------------- Utility functions ---------------------- */
1102
1103
1125 template <typename NumberType,
1126 int rank,
1127 int dim,
1128 template <int, int, typename>
1129 class TensorType>
1130 TensorType<rank, dim, NumberType>
1132 const TensorType<rank, dim, Expression> &symbol_tensor,
1133 const std::vector<NumberType> &cached_evaluation,
1134 const BatchOptimizer<NumberType> &optimizer)
1135 {
1136 TensorType<rank, dim, NumberType> out;
1137 for (unsigned int i = 0; i < out.n_independent_components; ++i)
1138 {
1139 const TableIndices<rank> indices(
1140 out.unrolled_to_component_indices(i));
1141 out[indices] =
1142 optimizer.extract(symbol_tensor[indices], cached_evaluation);
1143 }
1144 return out;
1145 }
1146
1147
1170 template <typename NumberType, int dim>
1173 const SymmetricTensor<4, dim, Expression> &symbol_tensor,
1174 const std::vector<NumberType> &cached_evaluation,
1175 const BatchOptimizer<NumberType> &optimizer)
1176 {
1178 for (unsigned int i = 0;
1179 i < SymmetricTensor<2, dim>::n_independent_components;
1180 ++i)
1181 for (unsigned int j = 0;
1182 j < SymmetricTensor<2, dim>::n_independent_components;
1183 ++j)
1184 {
1185 const TableIndices<4> indices =
1186 make_rank_4_tensor_indices<dim>(i, j);
1187 out[indices] =
1188 optimizer.extract(symbol_tensor[indices], cached_evaluation);
1189 }
1190 return out;
1191 }
1192
1193
1211 template <typename NumberType, typename T>
1212 void
1214 const T &function)
1215 {
1216 optimizer.register_function(function);
1217 }
1218
1219
1237 template <typename NumberType, typename T>
1238 void
1240 const std::vector<T> &functions)
1241 {
1242 for (const auto &function : functions)
1243 register_functions(optimizer, function);
1244 }
1245
1246
1266 template <typename NumberType, typename T, typename... Args>
1267 void
1269 const T &function,
1270 const Args &...other_functions)
1271 {
1272 register_functions(optimizer, function);
1273 register_functions(optimizer, other_functions...);
1274 }
1275
1276
1288 template <int rank,
1289 int dim,
1290 template <int, int, typename>
1291 class TensorType>
1294 const TensorType<rank, dim, Expression> &symbol_tensor)
1295 {
1297 out.reserve(symbol_tensor.n_independent_components);
1298 for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
1299 ++i)
1300 {
1301 const TableIndices<rank> indices(
1302 symbol_tensor.unrolled_to_component_indices(i));
1303 out.push_back(symbol_tensor[indices].get_RCP());
1304 }
1305 return out;
1306 }
1307
1308
1318 template <int dim>
1321 const SymmetricTensor<4, dim, Expression> &symbol_tensor)
1322 {
1324 out.reserve(symbol_tensor.n_independent_components);
1325 for (unsigned int i = 0;
1326 i < SymmetricTensor<2, dim>::n_independent_components;
1327 ++i)
1328 for (unsigned int j = 0;
1329 j < SymmetricTensor<2, dim>::n_independent_components;
1330 ++j)
1331 {
1332 const TableIndices<4> indices =
1333 make_rank_4_tensor_indices<dim>(i, j);
1334 out.push_back(symbol_tensor[indices].get_RCP());
1335 }
1336 return out;
1337 }
1338
1339 } // namespace internal
1340
1341
1342
1433 template <typename ReturnType>
1435 {
1436 public:
1445
1466
1476 BatchOptimizer(const BatchOptimizer &other);
1477
1481 BatchOptimizer(BatchOptimizer &&) noexcept = default;
1482
1486 ~BatchOptimizer() = default;
1487
1498 void
1499 copy_from(const BatchOptimizer &other);
1500
1510 template <typename Stream>
1511 void
1512 print(Stream &stream, const bool print_cse = false) const;
1513
1522 template <class Archive>
1523 void
1524 save(Archive &archive, const unsigned int version) const;
1525
1539 template <class Archive>
1540 void
1541 load(Archive &archive, const unsigned int version);
1542
1543# ifdef DOXYGEN
1565 template <class Archive>
1566 void
1567 serialize(Archive &archive, const unsigned int version);
1568# else
1569 // This macro defines the serialize() method that is compatible with
1570 // the templated save() and load() method that have been implemented.
1571 BOOST_SERIALIZATION_SPLIT_MEMBER()
1572# endif
1573
1584 void
1585 register_symbols(const types::substitution_map &substitution_map);
1586
1592 void
1593 register_symbols(const SymEngine::map_basic_basic &substitution_map);
1594
1605 void
1606 register_symbols(const types::symbol_vector &symbols);
1607
1618 void
1619 register_symbols(const SymEngine::vec_basic &symbols);
1620
1627
1633 std::size_t
1635
1647 void
1648 register_function(const Expression &function);
1649
1654 template <int rank, int dim>
1655 void
1657
1662 template <int rank, int dim>
1663 void
1665 const SymmetricTensor<rank, dim, Expression> &function_tensor);
1666
1671 void
1672 register_functions(const types::symbol_vector &functions);
1673
1678 void
1679 register_functions(const SymEngine::vec_basic &functions);
1680
1690 template <typename T>
1691 void
1692 register_functions(const std::vector<T> &functions);
1693
1707 template <typename T, typename... Args>
1708 void
1709 register_functions(const T &functions, const Args &...other_functions);
1710
1715 const types::symbol_vector &
1717
1724 std::size_t
1725 n_dependent_variables() const;
1726
1745 void
1749
1754 enum OptimizerType
1755 optimization_method() const;
1756
1762 optimization_flags() const;
1763
1769 bool
1770 use_symbolic_CSE() const;
1771
1787 void
1788 optimize();
1789
1794 bool
1795 optimized() const;
1796
1813 void
1814 substitute(const types::substitution_map &substitution_map) const;
1815
1825 void
1826 substitute(const SymEngine::map_basic_basic &substitution_map) const;
1827
1838 void
1839 substitute(const types::symbol_vector &symbols,
1840 const std::vector<ReturnType> &values) const;
1841
1852 void
1853 substitute(const SymEngine::vec_basic &symbols,
1854 const std::vector<ReturnType> &values) const;
1855
1861 bool
1862 values_substituted() const;
1863
1893 const std::vector<ReturnType> &
1894 evaluate() const;
1895
1903 ReturnType
1904 evaluate(const Expression &func) const;
1905
1914 std::vector<ReturnType>
1915 evaluate(const std::vector<Expression> &funcs) const;
1916
1925 template <int rank, int dim>
1928
1929
1938 template <int rank, int dim>
1941
1942
1950 ReturnType
1951 extract(const Expression &func,
1952 const std::vector<ReturnType> &cached_evaluation) const;
1953
1954
1962 std::vector<ReturnType>
1963 extract(const std::vector<Expression> &funcs,
1964 const std::vector<ReturnType> &cached_evaluation) const;
1965
1966
1974 template <int rank, int dim>
1977 const std::vector<ReturnType> &cached_evaluation) const;
1978
1979
1987 template <int rank, int dim>
1990 const std::vector<ReturnType> &cached_evaluation) const;
1991
1994 private:
1999
2005
2016
2023
2028 bool
2030 const SD::Expression &function) const;
2031
2036 bool
2038 const SymEngine::RCP<const SymEngine::Basic> &function) const;
2039
2054 mutable std::vector<ReturnType> dependent_variables_output;
2055
2065 std::map<SD::Expression,
2066 std::size_t,
2068
2074
2081 mutable std::unique_ptr<SymEngine::Visitor> optimizer;
2082
2092
2098
2102 void
2104
2109 void
2111
2115 void
2116 create_optimizer(std::unique_ptr<SymEngine::Visitor> &optimizer);
2117
2134 void
2135 substitute(const std::vector<ReturnType> &substitution_values) const;
2136 };
2137
2138
2139
2140 /* -------------------- inline and template functions ------------------ */
2141
2142
2143# ifndef DOXYGEN
2144
2145
2146 template <typename ReturnType>
2147 template <typename Stream>
2148 void
2149 BatchOptimizer<ReturnType>::print(Stream &stream,
2150 const bool /*print_cse*/) const
2151 {
2152 // Settings
2153 stream << "Method? " << optimization_method() << '\n';
2154 stream << "Flags: " << optimization_flags() << '\n';
2155 stream << "Optimized? " << (optimized() ? "Yes" : "No") << '\n';
2156 stream << "Values substituted? " << values_substituted() << "\n\n";
2157
2158 // Independent variables
2159 stream << "Symbols (" << n_independent_variables()
2160 << " independent variables):" << '\n';
2161 int cntr = 0;
2162 for (SD::types::substitution_map::const_iterator it =
2163 independent_variables_symbols.begin();
2164 it != independent_variables_symbols.end();
2165 ++it, ++cntr)
2166 {
2167 stream << cntr << ": " << it->first << '\n';
2168 }
2169 stream << '\n' << std::flush;
2170
2171 // Dependent functions
2172 stream << "Functions (" << n_dependent_variables()
2173 << " dependent variables):" << '\n';
2174 cntr = 0;
2175 for (typename SD::types::symbol_vector::const_iterator it =
2176 dependent_variables_functions.begin();
2177 it != dependent_variables_functions.end();
2178 ++it, ++cntr)
2179 {
2180 stream << cntr << ": " << (*it) << '\n';
2181 }
2182 stream << '\n' << std::flush;
2183
2184 // Common subexpression
2185 if (optimized() == true && use_symbolic_CSE() == true)
2186 {
2187 Assert(optimizer, ExcNotInitialized());
2188 const bool print_cse_reductions = true;
2189 const bool print_independent_symbols = false;
2190 const bool print_dependent_functions = false;
2191
2192 if (optimization_method() == OptimizerType::dictionary)
2193 {
2194 Assert(dynamic_cast<typename internal::DictionaryOptimizer<
2195 ReturnType>::OptimizerType *>(optimizer.get()),
2196 ExcMessage("Cannot cast optimizer to Dictionary type."));
2197
2198 internal::OptimizerHelper<
2199 ReturnType,
2200 internal::DictionaryOptimizer<ReturnType>>::
2201 print(stream,
2202 dynamic_cast<typename internal::DictionaryOptimizer<
2203 ReturnType>::OptimizerType *>(optimizer.get()),
2204 print_independent_symbols,
2205 print_dependent_functions,
2206 print_cse_reductions);
2207
2208 stream << '\n' << std::flush;
2209 }
2210 else if (optimization_method() == OptimizerType::lambda)
2211 {
2212 Assert(dynamic_cast<typename internal::LambdaOptimizer<
2213 ReturnType>::OptimizerType *>(optimizer.get()),
2214 ExcMessage("Cannot cast optimizer to Lambda type."));
2215
2216 internal::OptimizerHelper<ReturnType,
2217 internal::LambdaOptimizer<ReturnType>>::
2218 print(stream,
2219 dynamic_cast<typename internal::LambdaOptimizer<
2220 ReturnType>::OptimizerType *>(optimizer.get()),
2221 print_independent_symbols,
2222 print_dependent_functions,
2223 print_cse_reductions);
2224 }
2225# ifdef HAVE_SYMENGINE_LLVM
2226 else if (optimization_method() == OptimizerType::llvm)
2227 {
2228 Assert(dynamic_cast<typename internal::LLVMOptimizer<
2229 ReturnType>::OptimizerType *>(optimizer.get()),
2230 ExcMessage("Cannot cast optimizer to LLVM type."));
2231
2232 internal::OptimizerHelper<ReturnType,
2233 internal::LLVMOptimizer<ReturnType>>::
2234 print(stream,
2235 dynamic_cast<typename internal::LLVMOptimizer<
2236 ReturnType>::OptimizerType *>(optimizer.get()),
2237 print_independent_symbols,
2238 print_dependent_functions,
2239 print_cse_reductions);
2240 }
2241# endif // HAVE_SYMENGINE_LLVM
2242 else
2243 {
2244 AssertThrow(false, ExcMessage("Unknown optimizer type."));
2245 }
2246 }
2247
2248 if (values_substituted())
2249 {
2250 stream << "Evaluated functions:" << '\n';
2251 stream << std::flush;
2252 cntr = 0;
2253 for (typename std::vector<ReturnType>::const_iterator it =
2254 dependent_variables_output.begin();
2255 it != dependent_variables_output.end();
2256 ++it, ++cntr)
2257 {
2258 stream << cntr << ": " << (*it) << '\n';
2259 }
2260 stream << '\n' << std::flush;
2261 }
2262 }
2263
2264
2265
2266 template <typename ReturnType>
2267 template <class Archive>
2268 void
2270 const unsigned int version) const
2271 {
2272 // Serialize enum classes...
2273 {
2274 const auto m =
2275 static_cast<typename std::underlying_type<OptimizerType>::type>(
2276 method);
2277 ar &m;
2278 }
2279 {
2280 const auto f =
2281 static_cast<typename std::underlying_type<OptimizationFlags>::type>(
2282 flags);
2283 ar &f;
2284 }
2285
2286 // Important: Independent variables must always be
2287 // serialized before the dependent variables.
2288 ar &independent_variables_symbols;
2289 ar &dependent_variables_functions;
2290
2291 ar &dependent_variables_output;
2292 ar &map_dep_expr_vec_entry;
2293 ar &ready_for_value_extraction;
2294
2295 // Mark that we've saved this class at some point.
2296 has_been_serialized = true;
2297 ar &has_been_serialized;
2298
2299 // When we serialize the optimizer itself, we have to (unfortunately)
2300 // provide it with sufficient information to rebuild itself from scratch.
2301 // This is because only two of the three optimization classes support
2302 // real serialization (i.e. have save/load capability).
2303 const SD::types::symbol_vector symbol_vec =
2304 Utilities::extract_symbols(independent_variables_symbols);
2305 if (typename internal::DictionaryOptimizer<ReturnType>::OptimizerType
2306 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2307 ReturnType>::OptimizerType *>(optimizer.get()))
2308 {
2309 Assert(optimization_method() == OptimizerType::dictionary,
2311 internal::OptimizerHelper<
2312 ReturnType,
2313 internal::DictionaryOptimizer<ReturnType>>::save(ar, version, opt);
2314 }
2315 else if (typename internal::LambdaOptimizer<ReturnType>::OptimizerType
2316 *opt = dynamic_cast<typename internal::LambdaOptimizer<
2317 ReturnType>::OptimizerType *>(optimizer.get()))
2318 {
2319 Assert(optimization_method() == OptimizerType::lambda,
2321 internal::OptimizerHelper<
2322 ReturnType,
2323 internal::LambdaOptimizer<ReturnType>>::save(ar, version, opt);
2324 }
2325# ifdef HAVE_SYMENGINE_LLVM
2326 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
2327 *opt = dynamic_cast<typename internal::LLVMOptimizer<
2328 ReturnType>::OptimizerType *>(optimizer.get()))
2329 {
2330 Assert(optimization_method() == OptimizerType::llvm,
2332 internal::OptimizerHelper<
2333 ReturnType,
2334 internal::LLVMOptimizer<ReturnType>>::save(ar, version, opt);
2335 }
2336# endif
2337 else
2338 {
2339 AssertThrow(false, ExcMessage("Unknown optimizer type."));
2340 }
2341 }
2342
2343
2344
2345 template <typename ReturnType>
2346 template <class Archive>
2347 void
2348 BatchOptimizer<ReturnType>::load(Archive &ar, const unsigned int version)
2349 {
2350 Assert(independent_variables_symbols.empty(), ExcInternalError());
2351 Assert(dependent_variables_functions.empty(), ExcInternalError());
2352 Assert(dependent_variables_output.empty(), ExcInternalError());
2353 Assert(map_dep_expr_vec_entry.empty(), ExcInternalError());
2354 Assert(ready_for_value_extraction == false, ExcInternalError());
2355
2356 // Deserialize enum classes...
2357 {
2358 typename std::underlying_type<OptimizerType>::type m;
2359 ar &m;
2360 method = static_cast<OptimizerType>(m);
2361 }
2362 {
2363 typename std::underlying_type<OptimizationFlags>::type f;
2364 ar &f;
2365 flags = static_cast<OptimizationFlags>(f);
2366 }
2367
2368 // Important: Independent variables must always be
2369 // deserialized before the dependent variables.
2370 ar &independent_variables_symbols;
2371 ar &dependent_variables_functions;
2372
2373 ar &dependent_variables_output;
2374 ar &map_dep_expr_vec_entry;
2375 ar &ready_for_value_extraction;
2376
2377 ar &has_been_serialized;
2378
2379 // If we're reading in data, then create the optimizer
2380 // and then deserialize it.
2381 Assert(!optimizer, ExcInternalError());
2382
2383 // Create and configure the optimizer
2384 create_optimizer(optimizer);
2385 Assert(optimizer, ExcNotInitialized());
2386
2387 // When we deserialize the optimizer itself, we have to (unfortunately)
2388 // provide it with sufficient information to rebuild itself from scratch.
2389 // This is because only two of the three optimization classes support
2390 // real serialization (i.e. have save/load capability).
2391 const SD::types::symbol_vector symbol_vec =
2392 Utilities::extract_symbols(independent_variables_symbols);
2393 if (typename internal::DictionaryOptimizer<ReturnType>::OptimizerType
2394 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2395 ReturnType>::OptimizerType *>(optimizer.get()))
2396 {
2397 Assert(optimization_method() == OptimizerType::dictionary,
2399 internal::OptimizerHelper<ReturnType,
2400 internal::DictionaryOptimizer<ReturnType>>::
2401 load(ar,
2402 version,
2403 opt,
2405 symbol_vec),
2407 dependent_variables_functions),
2408 optimization_flags());
2409 }
2410 else if (typename internal::LambdaOptimizer<ReturnType>::OptimizerType
2411 *opt = dynamic_cast<typename internal::LambdaOptimizer<
2412 ReturnType>::OptimizerType *>(optimizer.get()))
2413 {
2414 Assert(optimization_method() == OptimizerType::lambda,
2416 internal::OptimizerHelper<ReturnType,
2417 internal::LambdaOptimizer<ReturnType>>::
2418 load(ar,
2419 version,
2420 opt,
2422 symbol_vec),
2424 dependent_variables_functions),
2425 optimization_flags());
2426 }
2427# ifdef HAVE_SYMENGINE_LLVM
2428 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
2429 *opt = dynamic_cast<typename internal::LLVMOptimizer<
2430 ReturnType>::OptimizerType *>(optimizer.get()))
2431 {
2432 Assert(optimization_method() == OptimizerType::llvm,
2434 internal::OptimizerHelper<ReturnType,
2435 internal::LLVMOptimizer<ReturnType>>::
2436 load(ar,
2437 version,
2438 opt,
2440 symbol_vec),
2442 dependent_variables_functions),
2443 optimization_flags());
2444 }
2445# endif
2446 else
2447 {
2448 AssertThrow(false, ExcMessage("Unknown optimizer type."));
2449 }
2450 }
2451
2452
2453
2454 template <typename ReturnType>
2455 template <int rank, int dim>
2456 void
2458 const Tensor<rank, dim, Expression> &function_tensor)
2459 {
2460 Assert(optimized() == false,
2461 ExcMessage(
2462 "Cannot register functions once the optimizer is finalised."));
2463
2464 register_vector_functions(
2465 internal::unroll_to_expression_vector(function_tensor));
2466 }
2467
2468
2469
2470 template <typename ReturnType>
2471 template <int rank, int dim>
2472 void
2474 const SymmetricTensor<rank, dim, Expression> &function_tensor)
2475 {
2476 Assert(optimized() == false,
2477 ExcMessage(
2478 "Cannot register functions once the optimizer is finalised."));
2479
2480 register_vector_functions(
2481 internal::unroll_to_expression_vector(function_tensor));
2482 }
2483
2484
2485
2486 template <typename ReturnType>
2487 template <typename T, typename... Args>
2488 void
2490 const T &functions,
2491 const Args &...other_functions)
2492 {
2493 internal::register_functions(*this, functions);
2494 internal::register_functions(*this, other_functions...);
2495 }
2496
2497
2498
2499 template <typename ReturnType>
2500 template <typename T>
2501 void
2503 const std::vector<T> &functions)
2504 {
2505 internal::register_functions(*this, functions);
2506 }
2507
2508
2509
2510 template <typename ReturnType>
2511 template <int rank, int dim>
2514 const Tensor<rank, dim, Expression> &funcs,
2515 const std::vector<ReturnType> &cached_evaluation) const
2516 {
2518 cached_evaluation,
2519 *this);
2520 }
2521
2522
2523
2524 template <typename ReturnType>
2525 template <int rank, int dim>
2528 const Tensor<rank, dim, Expression> &funcs) const
2529 {
2530 Assert(
2531 values_substituted() == true,
2532 ExcMessage(
2533 "The optimizer is not configured to perform evaluation. "
2534 "This action can only performed after substitute() has been called."));
2535
2536 return extract(funcs, dependent_variables_output);
2537 }
2538
2539
2540
2541 template <typename ReturnType>
2542 template <int rank, int dim>
2546 const std::vector<ReturnType> &cached_evaluation) const
2547 {
2549 cached_evaluation,
2550 *this);
2551 }
2552
2553
2554
2555 template <typename ReturnType>
2556 template <int rank, int dim>
2559 const SymmetricTensor<rank, dim, Expression> &funcs) const
2560 {
2561 Assert(
2562 values_substituted() == true,
2563 ExcMessage(
2564 "The optimizer is not configured to perform evaluation. "
2565 "This action can only performed after substitute() has been called."));
2566
2567 return extract(funcs, dependent_variables_output);
2568 }
2569
2570# endif // DOXYGEN
2571
2572 } // namespace SD
2573} // namespace Differentiation
2574
2575
2577
2578#endif // DEAL_II_WITH_SYMENGINE
2579
2580#endif
SymmetricTensor< rank, dim, ReturnType > extract(const SymmetricTensor< rank, dim, Expression > &funcs, const std::vector< ReturnType > &cached_evaluation) const
types::substitution_map independent_variables_symbols
void substitute(const types::substitution_map &substitution_map) const
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)
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)
SymmetricTensor< rank, dim, ReturnType > evaluate(const SymmetricTensor< rank, dim, Expression > &funcs) const
void save(Archive &archive, const unsigned int version) const
enum OptimizationFlags optimization_flags() const
void register_functions(const types::symbol_vector &functions)
std::vector< ReturnType > dependent_variables_output
Tensor< rank, dim, ReturnType > extract(const Tensor< rank, dim, Expression > &funcs, const std::vector< ReturnType > &cached_evaluation) 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
void register_function(const Expression &function)
Tensor< rank, dim, ReturnType > evaluate(const Tensor< rank, dim, Expression > &funcs) const
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)
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
static constexpr unsigned int n_independent_components
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcSymEngineLLVMReturnTypeNotSupported()
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcSymEngineLLVMNotAvailable()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#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)
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
std::ostream & operator<<(std::ostream &stream, const Expression &expression)
OptimizationFlags & operator&=(OptimizationFlags &f1, const OptimizationFlags f2)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
STL namespace.