Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19: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 - 2024 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
37
43
44# include <boost/serialization/split_member.hpp>
45# include <boost/type_traits.hpp>
46
47# include <algorithm>
48# include <map>
49# include <memory>
50# include <type_traits>
51# include <utility>
52# include <vector>
53
54
56
57
58namespace Differentiation
59{
60 namespace SD
61 {
72 "SymEngine has not been built with LLVM support.");
73
79 "The SymEngine LLVM optimizer does not (yet) support the "
80 "selected return type.");
81
85 // Forward declarations
86 template <typename ReturnType>
87 class BatchOptimizer;
88
89
95 enum class OptimizerType
96 {
105 lambda,
110 llvm
111 };
112
113
117 template <typename StreamType>
118 inline StreamType &
119 operator<<(StreamType &s, OptimizerType o)
120 {
122 s << "dictionary";
123 else if (o == OptimizerType::lambda)
124 s << "lambda";
125 else if (o == OptimizerType::llvm)
126 s << "llvm";
127 else
128 {
129 Assert(false, ExcMessage("Unknown optimization method."));
130 }
131
132 return s;
133 }
134
135
141 enum class OptimizationFlags : unsigned char
142 {
150 optimize_cse = 0x0001,
155 optimize_aggressive = 0x0002,
160 };
161
162
171 // This operator exists since if it did not then the result of the bit-or
172 // <tt>operator |</tt> would be an integer which would in turn trigger a
173 // compiler warning when we tried to assign it to an object of type
174 // OptimizationFlags.
175 inline OptimizationFlags
177 {
178 return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) |
179 static_cast<unsigned int>(f2));
180 }
181
182
187 inline OptimizationFlags &
189 {
190 f1 = f1 | f2;
191 return f1;
192 }
193
194
203 // This operator exists since if it did not then the result of the bit-or
204 // <tt>operator |</tt> would be an integer which would in turn trigger a
205 // compiler warning when we tried to assign it to an object of type
206 // OptimizationFlags.
207 inline OptimizationFlags
209 {
210 return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) &
211 static_cast<unsigned int>(f2));
212 }
213
214
219 inline OptimizationFlags &
221 {
222 f1 = f1 & f2;
223 return f1;
224 }
225
226
227 namespace internal
228 {
233 inline bool
235 {
236 return static_cast<int>(flags & OptimizationFlags::optimize_cse);
237 }
238
243 inline int
245 {
246 // With the LLVM compiler there exists the opportunity to tune
247 // the level of optimizations performed during compilation.
248 // By default SymEngine sets this at "opt_level=2", which one
249 // presumes targets -O2. Here we are a bit more specific about
250 // want we want it to do:
251 // - Normal compilation: -02 (default settings)
252 // - Aggressive mode: -03 (the whole lot!)
253 // In theory we could also target
254 // - Debug mode: -O0 (no optimizations)
255 // but this doesn't make much sense since SymEngine is a
256 // tested external library.
257 const bool use_agg_opt =
258 static_cast<int>(flags & OptimizationFlags::optimize_aggressive);
259 const int opt_level = (use_agg_opt ? 3 : 2);
260 return opt_level;
261 }
262 } // namespace internal
263
264
269 template <typename StreamType>
270 inline StreamType &
271 operator<<(StreamType &s, OptimizationFlags o)
272 {
273 s << " OptimizationFlags|";
274 if (static_cast<unsigned int>(o & OptimizationFlags::optimize_cse))
275 s << "cse|";
276
277 // LLVM optimization level
278 s << "-O" + std::to_string(internal::get_LLVM_optimization_level(o)) +
279 "|";
280
281 return s;
282 }
283
284
285 namespace internal
286 {
296 template <typename ReturnType, typename T = void>
298
299
309 template <typename ReturnType, typename T = void>
311
312
313# ifdef HAVE_SYMENGINE_LLVM
323 template <typename ReturnType, typename T = void>
324 struct LLVMOptimizer;
325# endif // HAVE_SYMENGINE_LLVM
326
327
343 template <typename ReturnType, typename Optimizer, typename T = void>
345
346
347# ifndef DOXYGEN
348
349
350 /* ----------- Specializations for the Optimizers ----------- */
351
352
353 // A helper struct to type trait detection for the optimizers that
354 // will be defined next.
355 template <typename ReturnType_, typename T = void>
356 struct SupportedOptimizerTypeTraits
357 {
358 static const bool is_supported = false;
359
360 using ReturnType = void;
361 };
362
363
364
365 // Specialization for arithmetic types
366 template <typename ReturnType_>
367 struct SupportedOptimizerTypeTraits<
368 ReturnType_,
369 std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
370 {
371 static const bool is_supported = true;
372
373 using ReturnType =
374 std::conditional_t<std::is_same_v<ReturnType_, float>, float, double>;
375 };
376
377
378
379 // Specialization for complex arithmetic types
380 template <typename ReturnType_>
381 struct SupportedOptimizerTypeTraits<
382 ReturnType_,
383 std::enable_if_t<
384 boost::is_complex<ReturnType_>::value &&
385 std::is_arithmetic_v<typename ReturnType_::value_type>>>
386 {
387 static const bool is_supported = true;
388
389 using ReturnType =
390 std::conditional_t<std::is_same_v<ReturnType_, std::complex<float>>,
391 std::complex<float>,
392 std::complex<double>>;
393 };
394
395
396
397 template <typename ReturnType_>
398 struct DictionaryOptimizer<ReturnType_,
399 std::enable_if_t<SupportedOptimizerTypeTraits<
400 ReturnType_>::is_supported>>
401 {
402 using ReturnType =
403 typename SupportedOptimizerTypeTraits<ReturnType_>::ReturnType;
404 using OptimizerType =
405 internal::DictionarySubstitutionVisitor<ReturnType, SD::Expression>;
406
407
416 static void
417 initialize(OptimizerType &optimizer,
418 const SymEngine::vec_basic &independent_symbols,
419 const SymEngine::vec_basic &dependent_functions,
420 const enum OptimizationFlags &optimization_flags)
421 {
422 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
423 optimizer.init(independent_symbols,
424 dependent_functions,
425 use_symbolic_cse);
426 }
427
428
429
434 template <class Archive>
435 static void
436 save(Archive &archive,
437 const unsigned int version,
438 OptimizerType &optimizer)
439 {
440 optimizer.save(archive, version);
441 }
442
443
444
449 template <class Archive>
450 static void
451 load(Archive &archive,
452 const unsigned int version,
453 OptimizerType &optimizer,
454 const SymEngine::vec_basic & /*independent_symbols*/,
455 const SymEngine::vec_basic & /*dependent_functions*/,
456 const enum OptimizationFlags & /*optimization_flags*/)
457 {
458 optimizer.load(archive, version);
459 }
460
461
462
478 template <typename Stream>
479 static void
480 print(Stream &stream,
481 const OptimizerType &optimizer,
482 const bool print_independent_symbols = false,
483 const bool print_dependent_functions = false,
484 const bool print_cse_reductions = true)
485 {
486 optimizer.print(stream,
487 print_independent_symbols,
488 print_dependent_functions,
489 print_cse_reductions);
490 }
491 };
492
493
494
495 template <typename ReturnType_>
496 struct LambdaOptimizer<ReturnType_,
497 std::enable_if_t<SupportedOptimizerTypeTraits<
498 ReturnType_>::is_supported>>
499 {
500 using ReturnType =
501 std::conditional_t<!boost::is_complex<ReturnType_>::value,
502 double,
503 std::complex<double>>;
504 using OptimizerType =
505 std::conditional_t<!boost::is_complex<ReturnType_>::value,
506 SymEngine::LambdaRealDoubleVisitor,
507 SymEngine::LambdaComplexDoubleVisitor>;
508
509
518 static void
519 initialize(OptimizerType &optimizer,
520 const SymEngine::vec_basic &independent_symbols,
521 const SymEngine::vec_basic &dependent_functions,
522 const enum OptimizationFlags &optimization_flags)
523 {
524 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
525 optimizer.init(independent_symbols,
526 dependent_functions,
527 use_symbolic_cse);
528 }
529
530
531
536 template <class Archive>
537 static void
538 save(Archive & /*archive*/,
539 const unsigned int /*version*/,
540 OptimizerType & /*optimizer*/)
541 {}
542
543
548 template <class Archive>
549 static void
550 load(Archive & /*archive*/,
551 const unsigned int /*version*/,
552 OptimizerType &optimizer,
553 const SymEngine::vec_basic &independent_symbols,
554 const SymEngine::vec_basic &dependent_functions,
555 const enum OptimizationFlags &optimization_flags)
556 {
557 initialize(optimizer,
558 independent_symbols,
559 dependent_functions,
560 optimization_flags);
561 }
562
563
564
580 template <typename StreamType>
581 static void
582 print(StreamType & /*stream*/,
583 const OptimizerType & /*optimizer*/,
584 const bool /*print_independent_symbols*/ = false,
585 const bool /*print_dependent_functions*/ = false,
586 const bool /*print_cse_reductions*/ = true)
587 {
588 // No built-in print function
589 }
590 };
591
592
593
594# ifdef HAVE_SYMENGINE_LLVM
595 template <typename ReturnType_>
596 struct LLVMOptimizer<ReturnType_,
597 std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
598 {
599 using ReturnType =
600 std::conditional_t<std::is_same_v<ReturnType_, float>, float, double>;
601 using OptimizerType =
602 std::conditional_t<std::is_same_v<ReturnType_, float>,
603 SymEngine::LLVMFloatVisitor,
604 SymEngine::LLVMDoubleVisitor>;
605
610 static const bool supported_by_LLVM = true;
611
612
621 static void
622 initialize(OptimizerType &optimizer,
623 const SymEngine::vec_basic &independent_symbols,
624 const SymEngine::vec_basic &dependent_functions,
625 const enum OptimizationFlags &optimization_flags)
626 {
627 const int opt_level = get_LLVM_optimization_level(optimization_flags);
628 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
629 optimizer.init(independent_symbols,
630 dependent_functions,
631 use_symbolic_cse,
632 opt_level);
633 }
634
635
636
641 template <class Archive>
642 static void
643 save(Archive &archive,
644 const unsigned int /*version*/,
645 OptimizerType &optimizer)
646 {
647 const std::string llvm_compiled_function = optimizer.dumps();
648 archive &llvm_compiled_function;
649 }
650
651
652
657 template <class Archive>
658 static void
659 load(Archive &archive,
660 const unsigned int /*version*/,
661 OptimizerType &optimizer,
662 const SymEngine::vec_basic & /*independent_symbols*/,
663 const SymEngine::vec_basic & /*dependent_functions*/,
664 const enum OptimizationFlags & /*optimization_flags*/)
665 {
666 std::string llvm_compiled_function;
667 archive &llvm_compiled_function;
668 optimizer.loads(llvm_compiled_function);
669 }
670
671
672
688 template <typename StreamType>
689 static void
690 print(StreamType & /*stream*/,
691 const OptimizerType & /*optimizer*/,
692 const bool /*print_independent_symbols*/ = false,
693 const bool /*print_dependent_functions*/ = false,
694 const bool /*print_cse_reductions*/ = true)
695 {
696 // No built-in print function
697 }
698 };
699
700
701 // There is no LLVM optimizer built with complex number support.
702 // So we fall back to the LambdaDouble case as a type (required
703 // at compile time), but offer no implementation. We expect that
704 // the calling class does not create this type: This can be done by
705 // checking the `supported_by_LLVM` flag.
706 template <typename ReturnType_>
707 struct LLVMOptimizer<
708 ReturnType_,
709 std::enable_if_t<
710 boost::is_complex<ReturnType_>::value &&
711 std::is_arithmetic_v<typename ReturnType_::value_type>>>
712 {
713 // Since there is no working implementation, these are dummy types
714 // that help with templating in the calling function.
715 using ReturnType = typename LambdaOptimizer<ReturnType_>::ReturnType;
716 using OptimizerType =
717 typename LambdaOptimizer<ReturnType_>::OptimizerType;
718
723 static const bool supported_by_LLVM = false;
724
725
734 static void
735 initialize(OptimizerType & /*optimizer*/,
736 const SymEngine::vec_basic & /*independent_symbols*/,
737 const SymEngine::vec_basic & /*dependent_functions*/,
738 const enum OptimizationFlags & /*optimization_flags*/)
739 {
741 }
742
743
744
749 template <class Archive>
750 static void
751 save(Archive & /*archive*/,
752 const unsigned int /*version*/,
753 OptimizerType & /*optimizer*/)
754 {
756 }
757
758
759
764 template <class Archive>
765 static void
766 load(Archive & /*archive*/,
767 const unsigned int /*version*/,
768 OptimizerType & /*optimizer*/,
769 const SymEngine::vec_basic & /*independent_symbols*/,
770 const SymEngine::vec_basic & /*dependent_functions*/,
771 const enum OptimizationFlags & /*optimization_flags*/)
772 {
774 }
775
776
777
793 template <typename StreamType>
794 static void
795 print(StreamType & /*stream*/,
796 const OptimizerType & /*optimizer*/,
797 const bool /*print_independent_symbols*/ = false,
798 const bool /*print_dependent_functions*/ = false,
799 const bool /*print_cse_reductions*/ = true)
800 {
802 }
803 };
804# endif // HAVE_SYMENGINE_LLVM
805
806
807 /* ----------- Specializations for OptimizerHelper ----------- */
808
809
810 template <typename ReturnType, typename Optimizer>
811 struct OptimizerHelper<
812 ReturnType,
813 Optimizer,
814 std::enable_if_t<
815 std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
816 {
825 static void
826 initialize(typename Optimizer::OptimizerType *optimizer,
827 const SymEngine::vec_basic &independent_symbols,
828 const SymEngine::vec_basic &dependent_functions,
829 const enum OptimizationFlags &optimization_flags)
830 {
831 Assert(optimizer, ExcNotInitialized());
832
833 // Some optimizers don't have the same interface for
834 // initialization, we filter them out through the specializations
835 // of the Optimizer class
836 Optimizer::initialize(*optimizer,
837 independent_symbols,
838 dependent_functions,
839 optimization_flags);
840 }
841
842
843
857 static void
858 substitute(typename Optimizer::OptimizerType *optimizer,
859 std::vector<ReturnType> &output_values,
860 const std::vector<ReturnType> &substitution_values)
861 {
862 Assert(optimizer, ExcNotInitialized());
863 optimizer->call(output_values.data(), substitution_values.data());
864 }
865
866
867
872 template <class Archive>
873 static void
874 save(Archive &archive,
875 const unsigned int version,
876 typename Optimizer::OptimizerType *optimizer)
877 {
878 Assert(optimizer, ExcNotInitialized());
879
880 // Some optimizers don't have the same interface for
881 // serialization, we filter them out through the specializations
882 // of the Optimizer class
883 Optimizer::save(archive, version, *optimizer);
884 }
885
886
887
892 template <class Archive>
893 static void
894 load(Archive &archive,
895 const unsigned int version,
896 typename Optimizer::OptimizerType *optimizer,
897 const SymEngine::vec_basic &independent_symbols,
898 const SymEngine::vec_basic &dependent_functions,
899 const enum OptimizationFlags &optimization_flags)
900 {
901 Assert(optimizer, ExcNotInitialized());
902
903 // Some optimizers don't have the same interface for
904 // serialization, we filter them out through the specializations
905 // of the Optimizer class
906 Optimizer::load(archive,
907 version,
908 *optimizer,
909 independent_symbols,
910 dependent_functions,
911 optimization_flags);
912 }
913
914
915
931 template <typename Stream>
932 static void
933 print(Stream &stream,
934 typename Optimizer::OptimizerType *optimizer,
935 const bool print_independent_symbols = false,
936 const bool print_dependent_functions = false,
937 const bool print_cse_reductions = true)
938 {
939 Assert(optimizer, ExcNotInitialized());
940
941 // Some optimizers don't have a print function, so
942 // we filter them out through the specializations of
943 // the Optimizer class
944 Optimizer::print(stream,
945 *optimizer,
946 print_independent_symbols,
947 print_dependent_functions,
948 print_cse_reductions);
949 }
950 };
951
952 template <typename ReturnType, typename Optimizer>
953 struct OptimizerHelper<
954 ReturnType,
955 Optimizer,
956 std::enable_if_t<
957 !std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
958 {
967 static void
968 initialize(typename Optimizer::OptimizerType *optimizer,
969 const SymEngine::vec_basic &independent_symbols,
970 const SymEngine::vec_basic &dependent_functions,
971 const enum OptimizationFlags &optimization_flags)
972 {
973 Assert(optimizer, ExcNotInitialized());
974
975 const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
976 optimizer->init(independent_symbols,
977 dependent_functions,
978 use_symbolic_cse);
979 }
980
981
982
996 static void
997 substitute(typename Optimizer::OptimizerType *optimizer,
998 std::vector<ReturnType> &output_values,
999 const std::vector<ReturnType> &substitution_values)
1000 {
1001 Assert(optimizer, ExcNotInitialized());
1002
1003 // Intermediate values to accommodate the difference in
1004 // value types.
1005 std::vector<typename Optimizer::ReturnType> int_outputs(
1006 output_values.size());
1007 std::vector<typename Optimizer::ReturnType> int_inputs(
1008 substitution_values.size());
1009
1010 std::copy(substitution_values.begin(),
1011 substitution_values.end(),
1012 int_inputs.begin());
1013 optimizer->call(int_outputs.data(), int_inputs.data());
1014 std::copy(int_outputs.begin(),
1015 int_outputs.end(),
1016 output_values.begin());
1017 }
1018
1019
1020
1025 template <class Archive>
1026 static void
1027 save(Archive &archive,
1028 const unsigned int version,
1029 typename Optimizer::OptimizerType *optimizer)
1030 {
1031 Assert(optimizer, ExcNotInitialized());
1032 Optimizer::save(archive, version, *optimizer);
1033 }
1034
1035
1036
1041 template <class Archive>
1042 static void
1043 load(Archive &archive,
1044 const unsigned int version,
1045 typename Optimizer::OptimizerType *optimizer,
1046 const SymEngine::vec_basic &independent_symbols,
1047 const SymEngine::vec_basic &dependent_functions,
1048 const enum OptimizationFlags &optimization_flags)
1049 {
1050 Assert(optimizer, ExcNotInitialized());
1051
1052 // Some optimizers don't have the same interface for
1053 // serialization, we filter them out through the specializations
1054 // of the Optimizer class
1055 Optimizer::load(archive,
1056 version,
1057 *optimizer,
1058 independent_symbols,
1059 dependent_functions,
1060 optimization_flags);
1061 }
1062
1063
1064
1080 template <typename Stream>
1081 static void
1082 print(Stream &stream,
1083 typename Optimizer::OptimizerType *optimizer,
1084 const bool print_cse_reductions = true,
1085 const bool print_independent_symbols = false,
1086 const bool print_dependent_functions = false)
1087 {
1088 Assert(optimizer, ExcNotInitialized());
1089
1090 optimizer->print(stream,
1091 print_independent_symbols,
1092 print_dependent_functions,
1093 print_cse_reductions);
1094 }
1095 };
1096
1097# endif // DOXYGEN
1098
1099
1100 /* -------------------- Utility functions ---------------------- */
1101
1102
1124 template <typename NumberType,
1125 int rank,
1126 int dim,
1127 template <int, int, typename>
1128 class TensorType>
1129 TensorType<rank, dim, NumberType>
1131 const TensorType<rank, dim, Expression> &symbol_tensor,
1132 const std::vector<NumberType> &cached_evaluation,
1133 const BatchOptimizer<NumberType> &optimizer)
1134 {
1135 TensorType<rank, dim, NumberType> out;
1136 for (unsigned int i = 0; i < out.n_independent_components; ++i)
1137 {
1138 const TableIndices<rank> indices(
1139 out.unrolled_to_component_indices(i));
1140 out[indices] =
1141 optimizer.extract(symbol_tensor[indices], cached_evaluation);
1142 }
1143 return out;
1144 }
1145
1146
1169 template <typename NumberType, int dim>
1172 const SymmetricTensor<4, dim, Expression> &symbol_tensor,
1173 const std::vector<NumberType> &cached_evaluation,
1174 const BatchOptimizer<NumberType> &optimizer)
1175 {
1177 for (unsigned int i = 0;
1178 i < SymmetricTensor<2, dim>::n_independent_components;
1179 ++i)
1180 for (unsigned int j = 0;
1181 j < SymmetricTensor<2, dim>::n_independent_components;
1182 ++j)
1183 {
1184 const TableIndices<4> indices =
1185 make_rank_4_tensor_indices<dim>(i, j);
1186 out[indices] =
1187 optimizer.extract(symbol_tensor[indices], cached_evaluation);
1188 }
1189 return out;
1190 }
1191
1192
1210 template <typename NumberType, typename T>
1211 void
1213 const T &function)
1214 {
1215 optimizer.register_function(function);
1216 }
1217
1218
1236 template <typename NumberType, typename T>
1237 void
1239 const std::vector<T> &functions)
1240 {
1241 for (const auto &function : functions)
1242 register_functions(optimizer, function);
1243 }
1244
1245
1265 template <typename NumberType, typename T, typename... Args>
1266 void
1268 const T &function,
1269 const Args &...other_functions)
1270 {
1271 register_functions(optimizer, function);
1272 register_functions(optimizer, other_functions...);
1273 }
1274
1275
1287 template <int rank,
1288 int dim,
1289 template <int, int, typename>
1290 class TensorType>
1293 const TensorType<rank, dim, Expression> &symbol_tensor)
1294 {
1296 out.reserve(symbol_tensor.n_independent_components);
1297 for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
1298 ++i)
1299 {
1300 const TableIndices<rank> indices(
1301 symbol_tensor.unrolled_to_component_indices(i));
1302 out.push_back(symbol_tensor[indices].get_RCP());
1303 }
1304 return out;
1305 }
1306
1307
1317 template <int dim>
1320 const SymmetricTensor<4, dim, Expression> &symbol_tensor)
1321 {
1323 out.reserve(symbol_tensor.n_independent_components);
1324 for (unsigned int i = 0;
1325 i < SymmetricTensor<2, dim>::n_independent_components;
1326 ++i)
1327 for (unsigned int j = 0;
1328 j < SymmetricTensor<2, dim>::n_independent_components;
1329 ++j)
1330 {
1331 const TableIndices<4> indices =
1332 make_rank_4_tensor_indices<dim>(i, j);
1333 out.push_back(symbol_tensor[indices].get_RCP());
1334 }
1335 return out;
1336 }
1337
1338 } // namespace internal
1339
1340
1341
1432 template <typename ReturnType>
1434 {
1435 public:
1444
1465
1475 BatchOptimizer(const BatchOptimizer &other);
1476
1480 BatchOptimizer(BatchOptimizer &&) noexcept = default;
1481
1485 ~BatchOptimizer() = default;
1486
1497 void
1498 copy_from(const BatchOptimizer &other);
1499
1509 template <typename Stream>
1510 void
1511 print(Stream &stream, const bool print_cse = false) const;
1512
1521 template <class Archive>
1522 void
1523 save(Archive &archive, const unsigned int version) const;
1524
1538 template <class Archive>
1539 void
1540 load(Archive &archive, const unsigned int version);
1541
1542# ifdef DOXYGEN
1564 template <class Archive>
1565 void
1566 serialize(Archive &archive, const unsigned int version);
1567# else
1568 // This macro defines the serialize() method that is compatible with
1569 // the templated save() and load() method that have been implemented.
1570 BOOST_SERIALIZATION_SPLIT_MEMBER()
1571# endif
1572
1583 void
1584 register_symbols(const types::substitution_map &substitution_map);
1585
1591 void
1592 register_symbols(const SymEngine::map_basic_basic &substitution_map);
1593
1604 void
1605 register_symbols(const types::symbol_vector &symbols);
1606
1617 void
1618 register_symbols(const SymEngine::vec_basic &symbols);
1619
1626
1632 std::size_t
1634
1646 void
1647 register_function(const Expression &function);
1648
1653 template <int rank, int dim>
1654 void
1656
1661 template <int rank, int dim>
1662 void
1664 const SymmetricTensor<rank, dim, Expression> &function_tensor);
1665
1670 void
1671 register_functions(const types::symbol_vector &functions);
1672
1677 void
1678 register_functions(const SymEngine::vec_basic &functions);
1679
1689 template <typename T>
1690 void
1691 register_functions(const std::vector<T> &functions);
1692
1706 template <typename T, typename... Args>
1707 void
1708 register_functions(const T &functions, const Args &...other_functions);
1709
1714 const types::symbol_vector &
1716
1723 std::size_t
1724 n_dependent_variables() const;
1725
1744 void
1748
1753 enum OptimizerType
1754 optimization_method() const;
1755
1761 optimization_flags() const;
1762
1768 bool
1769 use_symbolic_CSE() const;
1770
1786 void
1787 optimize();
1788
1793 bool
1794 optimized() const;
1795
1812 void
1813 substitute(const types::substitution_map &substitution_map) const;
1814
1824 void
1825 substitute(const SymEngine::map_basic_basic &substitution_map) const;
1826
1837 void
1838 substitute(const types::symbol_vector &symbols,
1839 const std::vector<ReturnType> &values) const;
1840
1851 void
1852 substitute(const SymEngine::vec_basic &symbols,
1853 const std::vector<ReturnType> &values) const;
1854
1860 bool
1861 values_substituted() const;
1862
1892 const std::vector<ReturnType> &
1893 evaluate() const;
1894
1902 ReturnType
1903 evaluate(const Expression &func) const;
1904
1913 std::vector<ReturnType>
1914 evaluate(const std::vector<Expression> &funcs) const;
1915
1924 template <int rank, int dim>
1927
1928
1937 template <int rank, int dim>
1940
1941
1949 ReturnType
1950 extract(const Expression &func,
1951 const std::vector<ReturnType> &cached_evaluation) const;
1952
1953
1961 std::vector<ReturnType>
1962 extract(const std::vector<Expression> &funcs,
1963 const std::vector<ReturnType> &cached_evaluation) const;
1964
1965
1973 template <int rank, int dim>
1976 const std::vector<ReturnType> &cached_evaluation) const;
1977
1978
1986 template <int rank, int dim>
1989 const std::vector<ReturnType> &cached_evaluation) const;
1990
1993 private:
1998
2004
2015
2022
2027 bool
2029 const SD::Expression &function) const;
2030
2035 bool
2037 const SymEngine::RCP<const SymEngine::Basic> &function) const;
2038
2053 mutable std::vector<ReturnType> dependent_variables_output;
2054
2064 std::map<SD::Expression,
2065 std::size_t,
2067
2073
2080 mutable std::unique_ptr<SymEngine::Visitor> optimizer;
2081
2091
2097
2101 void
2103
2108 void
2110
2114 void
2115 create_optimizer(std::unique_ptr<SymEngine::Visitor> &optimizer);
2116
2133 void
2134 substitute(const std::vector<ReturnType> &substitution_values) const;
2135 };
2136
2137
2138
2139 /* -------------------- inline and template functions ------------------ */
2140
2141
2142# ifndef DOXYGEN
2143
2144
2145 template <typename ReturnType>
2146 template <typename Stream>
2147 void
2148 BatchOptimizer<ReturnType>::print(Stream &stream,
2149 const bool /*print_cse*/) const
2150 {
2151 // Settings
2152 stream << "Method? " << optimization_method() << '\n';
2153 stream << "Flags: " << optimization_flags() << '\n';
2154 stream << "Optimized? " << (optimized() ? "Yes" : "No") << '\n';
2155 stream << "Values substituted? " << values_substituted() << "\n\n";
2156
2157 // Independent variables
2158 stream << "Symbols (" << n_independent_variables()
2159 << " independent variables):" << '\n';
2160 int cntr = 0;
2161 for (SD::types::substitution_map::const_iterator it =
2162 independent_variables_symbols.begin();
2163 it != independent_variables_symbols.end();
2164 ++it, ++cntr)
2165 {
2166 stream << cntr << ": " << it->first << '\n';
2167 }
2168 stream << '\n' << std::flush;
2169
2170 // Dependent functions
2171 stream << "Functions (" << n_dependent_variables()
2172 << " dependent variables):" << '\n';
2173 cntr = 0;
2174 for (typename SD::types::symbol_vector::const_iterator it =
2175 dependent_variables_functions.begin();
2176 it != dependent_variables_functions.end();
2177 ++it, ++cntr)
2178 {
2179 stream << cntr << ": " << (*it) << '\n';
2180 }
2181 stream << '\n' << std::flush;
2182
2183 // Common subexpression
2184 if (optimized() == true && use_symbolic_CSE() == true)
2185 {
2186 Assert(optimizer, ExcNotInitialized());
2187 const bool print_cse_reductions = true;
2188 const bool print_independent_symbols = false;
2189 const bool print_dependent_functions = false;
2190
2191 if (optimization_method() == OptimizerType::dictionary)
2192 {
2193 Assert(dynamic_cast<typename internal::DictionaryOptimizer<
2194 ReturnType>::OptimizerType *>(optimizer.get()),
2195 ExcMessage("Cannot cast optimizer to Dictionary type."));
2196
2197 internal::OptimizerHelper<
2198 ReturnType,
2199 internal::DictionaryOptimizer<ReturnType>>::
2200 print(stream,
2201 dynamic_cast<typename internal::DictionaryOptimizer<
2202 ReturnType>::OptimizerType *>(optimizer.get()),
2203 print_independent_symbols,
2204 print_dependent_functions,
2205 print_cse_reductions);
2206
2207 stream << '\n' << std::flush;
2208 }
2209 else if (optimization_method() == OptimizerType::lambda)
2210 {
2211 Assert(dynamic_cast<typename internal::LambdaOptimizer<
2212 ReturnType>::OptimizerType *>(optimizer.get()),
2213 ExcMessage("Cannot cast optimizer to Lambda type."));
2214
2215 internal::OptimizerHelper<ReturnType,
2216 internal::LambdaOptimizer<ReturnType>>::
2217 print(stream,
2218 dynamic_cast<typename internal::LambdaOptimizer<
2219 ReturnType>::OptimizerType *>(optimizer.get()),
2220 print_independent_symbols,
2221 print_dependent_functions,
2222 print_cse_reductions);
2223 }
2224# ifdef HAVE_SYMENGINE_LLVM
2225 else if (optimization_method() == OptimizerType::llvm)
2226 {
2227 Assert(dynamic_cast<typename internal::LLVMOptimizer<
2228 ReturnType>::OptimizerType *>(optimizer.get()),
2229 ExcMessage("Cannot cast optimizer to LLVM type."));
2230
2231 internal::OptimizerHelper<ReturnType,
2232 internal::LLVMOptimizer<ReturnType>>::
2233 print(stream,
2234 dynamic_cast<typename internal::LLVMOptimizer<
2235 ReturnType>::OptimizerType *>(optimizer.get()),
2236 print_independent_symbols,
2237 print_dependent_functions,
2238 print_cse_reductions);
2239 }
2240# endif // HAVE_SYMENGINE_LLVM
2241 else
2242 {
2243 AssertThrow(false, ExcMessage("Unknown optimizer type."));
2244 }
2245 }
2246
2247 if (values_substituted())
2248 {
2249 stream << "Evaluated functions:" << '\n';
2250 stream << std::flush;
2251 cntr = 0;
2252 for (typename std::vector<ReturnType>::const_iterator it =
2253 dependent_variables_output.begin();
2254 it != dependent_variables_output.end();
2255 ++it, ++cntr)
2256 {
2257 stream << cntr << ": " << (*it) << '\n';
2258 }
2259 stream << '\n' << std::flush;
2260 }
2261 }
2262
2263
2264
2265 template <typename ReturnType>
2266 template <class Archive>
2267 void
2269 const unsigned int version) const
2270 {
2271 // Serialize enum classes...
2272 {
2273 const auto m =
2274 static_cast<std::underlying_type_t<OptimizerType>>(method);
2275 ar &m;
2276 }
2277 {
2278 const auto f =
2279 static_cast<std::underlying_type_t<OptimizationFlags>>(flags);
2280 ar &f;
2281 }
2282
2283 // Important: Independent variables must always be
2284 // serialized before the dependent variables.
2285 ar &independent_variables_symbols;
2286 ar &dependent_variables_functions;
2287
2288 ar &dependent_variables_output;
2289 ar &map_dep_expr_vec_entry;
2290 ar &ready_for_value_extraction;
2291
2292 // Mark that we've saved this class at some point.
2293 has_been_serialized = true;
2294 ar &has_been_serialized;
2295
2296 // When we serialize the optimizer itself, we have to (unfortunately)
2297 // provide it with sufficient information to rebuild itself from scratch.
2298 // This is because only two of the three optimization classes support
2299 // real serialization (i.e. have save/load capability).
2300 const SD::types::symbol_vector symbol_vec =
2301 Utilities::extract_symbols(independent_variables_symbols);
2302 if (typename internal::DictionaryOptimizer<ReturnType>::OptimizerType
2303 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2304 ReturnType>::OptimizerType *>(optimizer.get()))
2305 {
2306 Assert(optimization_method() == OptimizerType::dictionary,
2308 internal::OptimizerHelper<
2309 ReturnType,
2310 internal::DictionaryOptimizer<ReturnType>>::save(ar, version, opt);
2311 }
2312 else if (typename internal::LambdaOptimizer<ReturnType>::OptimizerType
2313 *opt = dynamic_cast<typename internal::LambdaOptimizer<
2314 ReturnType>::OptimizerType *>(optimizer.get()))
2315 {
2316 Assert(optimization_method() == OptimizerType::lambda,
2318 internal::OptimizerHelper<
2319 ReturnType,
2320 internal::LambdaOptimizer<ReturnType>>::save(ar, version, opt);
2321 }
2322# ifdef HAVE_SYMENGINE_LLVM
2323 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
2324 *opt = dynamic_cast<typename internal::LLVMOptimizer<
2325 ReturnType>::OptimizerType *>(optimizer.get()))
2326 {
2327 Assert(optimization_method() == OptimizerType::llvm,
2329 internal::OptimizerHelper<
2330 ReturnType,
2331 internal::LLVMOptimizer<ReturnType>>::save(ar, version, opt);
2332 }
2333# endif
2334 else
2335 {
2336 AssertThrow(false, ExcMessage("Unknown optimizer type."));
2337 }
2338 }
2339
2340
2341
2342 template <typename ReturnType>
2343 template <class Archive>
2344 void
2345 BatchOptimizer<ReturnType>::load(Archive &ar, const unsigned int version)
2346 {
2347 Assert(independent_variables_symbols.empty(), ExcInternalError());
2348 Assert(dependent_variables_functions.empty(), ExcInternalError());
2349 Assert(dependent_variables_output.empty(), ExcInternalError());
2350 Assert(map_dep_expr_vec_entry.empty(), ExcInternalError());
2351 Assert(ready_for_value_extraction == false, ExcInternalError());
2352
2353 // Deserialize enum classes...
2354 {
2355 std::underlying_type_t<OptimizerType> m;
2356 ar &m;
2357 method = static_cast<OptimizerType>(m);
2358 }
2359 {
2360 std::underlying_type_t<OptimizationFlags> f;
2361 ar &f;
2362 flags = static_cast<OptimizationFlags>(f);
2363 }
2364
2365 // Important: Independent variables must always be
2366 // deserialized before the dependent variables.
2367 ar &independent_variables_symbols;
2368 ar &dependent_variables_functions;
2369
2370 ar &dependent_variables_output;
2371 ar &map_dep_expr_vec_entry;
2372 ar &ready_for_value_extraction;
2373
2374 ar &has_been_serialized;
2375
2376 // If we're reading in data, then create the optimizer
2377 // and then deserialize it.
2378 Assert(!optimizer, ExcInternalError());
2379
2380 // Create and configure the optimizer
2381 create_optimizer(optimizer);
2382 Assert(optimizer, ExcNotInitialized());
2383
2384 // When we deserialize the optimizer itself, we have to (unfortunately)
2385 // provide it with sufficient information to rebuild itself from scratch.
2386 // This is because only two of the three optimization classes support
2387 // real serialization (i.e. have save/load capability).
2388 const SD::types::symbol_vector symbol_vec =
2389 Utilities::extract_symbols(independent_variables_symbols);
2390 if (typename internal::DictionaryOptimizer<ReturnType>::OptimizerType
2391 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2392 ReturnType>::OptimizerType *>(optimizer.get()))
2393 {
2394 Assert(optimization_method() == OptimizerType::dictionary,
2396 internal::OptimizerHelper<ReturnType,
2397 internal::DictionaryOptimizer<ReturnType>>::
2398 load(ar,
2399 version,
2400 opt,
2402 symbol_vec),
2404 dependent_variables_functions),
2405 optimization_flags());
2406 }
2407 else if (typename internal::LambdaOptimizer<ReturnType>::OptimizerType
2408 *opt = dynamic_cast<typename internal::LambdaOptimizer<
2409 ReturnType>::OptimizerType *>(optimizer.get()))
2410 {
2411 Assert(optimization_method() == OptimizerType::lambda,
2413 internal::OptimizerHelper<ReturnType,
2414 internal::LambdaOptimizer<ReturnType>>::
2415 load(ar,
2416 version,
2417 opt,
2419 symbol_vec),
2421 dependent_variables_functions),
2422 optimization_flags());
2423 }
2424# ifdef HAVE_SYMENGINE_LLVM
2425 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
2426 *opt = dynamic_cast<typename internal::LLVMOptimizer<
2427 ReturnType>::OptimizerType *>(optimizer.get()))
2428 {
2429 Assert(optimization_method() == OptimizerType::llvm,
2431 internal::OptimizerHelper<ReturnType,
2432 internal::LLVMOptimizer<ReturnType>>::
2433 load(ar,
2434 version,
2435 opt,
2437 symbol_vec),
2439 dependent_variables_functions),
2440 optimization_flags());
2441 }
2442# endif
2443 else
2444 {
2445 AssertThrow(false, ExcMessage("Unknown optimizer type."));
2446 }
2447 }
2448
2449
2450
2451 template <typename ReturnType>
2452 template <int rank, int dim>
2453 void
2455 const Tensor<rank, dim, Expression> &function_tensor)
2456 {
2457 Assert(optimized() == false,
2458 ExcMessage(
2459 "Cannot register functions once the optimizer is finalised."));
2460
2461 register_vector_functions(
2462 internal::unroll_to_expression_vector(function_tensor));
2463 }
2464
2465
2466
2467 template <typename ReturnType>
2468 template <int rank, int dim>
2469 void
2471 const SymmetricTensor<rank, dim, Expression> &function_tensor)
2472 {
2473 Assert(optimized() == false,
2474 ExcMessage(
2475 "Cannot register functions once the optimizer is finalised."));
2476
2477 register_vector_functions(
2478 internal::unroll_to_expression_vector(function_tensor));
2479 }
2480
2481
2482
2483 template <typename ReturnType>
2484 template <typename T, typename... Args>
2485 void
2487 const T &functions,
2488 const Args &...other_functions)
2489 {
2490 internal::register_functions(*this, functions);
2491 internal::register_functions(*this, other_functions...);
2492 }
2493
2494
2495
2496 template <typename ReturnType>
2497 template <typename T>
2498 void
2500 const std::vector<T> &functions)
2501 {
2502 internal::register_functions(*this, functions);
2503 }
2504
2505
2506
2507 template <typename ReturnType>
2508 template <int rank, int dim>
2511 const Tensor<rank, dim, Expression> &funcs,
2512 const std::vector<ReturnType> &cached_evaluation) const
2513 {
2515 cached_evaluation,
2516 *this);
2517 }
2518
2519
2520
2521 template <typename ReturnType>
2522 template <int rank, int dim>
2525 const Tensor<rank, dim, Expression> &funcs) const
2526 {
2527 Assert(
2528 values_substituted() == true,
2529 ExcMessage(
2530 "The optimizer is not configured to perform evaluation. "
2531 "This action can only performed after substitute() has been called."));
2532
2533 return extract(funcs, dependent_variables_output);
2534 }
2535
2536
2537
2538 template <typename ReturnType>
2539 template <int rank, int dim>
2543 const std::vector<ReturnType> &cached_evaluation) const
2544 {
2546 cached_evaluation,
2547 *this);
2548 }
2549
2550
2551
2552 template <typename ReturnType>
2553 template <int rank, int dim>
2556 const SymmetricTensor<rank, dim, Expression> &funcs) const
2557 {
2558 Assert(
2559 values_substituted() == true,
2560 ExcMessage(
2561 "The optimizer is not configured to perform evaluation. "
2562 "This action can only performed after substitute() has been called."));
2563
2564 return extract(funcs, dependent_variables_output);
2565 }
2566
2567# endif // DOXYGEN
2568
2569 } // namespace SD
2570} // namespace Differentiation
2571
2572
2574
2575#endif // DEAL_II_WITH_SYMENGINE
2576
2577#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.