Reference documentation for deal.II version GIT c14369f203 2023-10-01 07:40: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\}}\)
symengine_optimizer.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_differentiation_sd_symengine_optimizer_h
17 #define dealii_differentiation_sd_symengine_optimizer_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_WITH_SYMENGINE
22 
23 // Low level
24 # include <symengine/basic.h>
25 # include <symengine/dict.h>
26 # include <symengine/symengine_exception.h>
27 # include <symengine/symengine_rcp.h>
28 
29 // Optimization
30 # include <symengine/lambda_double.h>
31 # include <symengine/visitor.h>
32 # ifdef HAVE_SYMENGINE_LLVM
33 # include <symengine/llvm_double.h>
34 # endif
35 
36 # include <deal.II/base/exceptions.h>
37 # include <deal.II/base/logstream.h>
38 # include <deal.II/base/utilities.h>
39 
45 
46 # include <boost/serialization/split_member.hpp>
47 # include <boost/type_traits.hpp>
48 
49 # include <algorithm>
50 # include <map>
51 # include <memory>
52 # include <type_traits>
53 # include <utility>
54 # include <vector>
55 
56 
58 
59 
60 namespace Differentiation
61 {
62  namespace SD
63  {
74  "SymEngine has not been built with LLVM support.");
75 
81  "The SymEngine LLVM optimizer does not (yet) support the "
82  "selected return type.");
83 
87  // Forward declarations
88  template <typename ReturnType>
89  class BatchOptimizer;
90 
91 
97  enum class OptimizerType
98  {
102  dictionary,
107  lambda,
112  llvm
113  };
114 
115 
119  template <typename StreamType>
120  inline StreamType &
121  operator<<(StreamType &s, OptimizerType o)
122  {
123  if (o == OptimizerType::dictionary)
124  s << "dictionary";
125  else if (o == OptimizerType::lambda)
126  s << "lambda";
127  else if (o == OptimizerType::llvm)
128  s << "llvm";
129  else
130  {
131  Assert(false, ExcMessage("Unknown optimization method."));
132  }
133 
134  return s;
135  }
136 
137 
143  enum class OptimizationFlags : unsigned char
144  {
148  optimize_default = 0,
152  optimize_cse = 0x0001,
157  optimize_aggressive = 0x0002,
162  };
163 
164 
173  // This operator exists since if it did not then the result of the bit-or
174  // <tt>operator |</tt> would be an integer which would in turn trigger a
175  // compiler warning when we tried to assign it to an object of type
176  // OptimizationFlags.
177  inline OptimizationFlags
179  {
180  return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) |
181  static_cast<unsigned int>(f2));
182  }
183 
184 
189  inline OptimizationFlags &
191  {
192  f1 = f1 | f2;
193  return f1;
194  }
195 
196 
205  // This operator exists since if it did not then the result of the bit-or
206  // <tt>operator |</tt> would be an integer which would in turn trigger a
207  // compiler warning when we tried to assign it to an object of type
208  // OptimizationFlags.
209  inline OptimizationFlags
211  {
212  return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) &
213  static_cast<unsigned int>(f2));
214  }
215 
216 
221  inline OptimizationFlags &
223  {
224  f1 = f1 & f2;
225  return f1;
226  }
227 
228 
229  namespace internal
230  {
235  inline bool
237  {
238  return static_cast<int>(flags & OptimizationFlags::optimize_cse);
239  }
240 
245  inline int
247  {
248  // With the LLVM compiler there exists the opportunity to tune
249  // the level of optimizations performed during compilation.
250  // By default SymEngine sets this at "opt_level=2", which one
251  // presumes targets -O2. Here we are a bit more specific about
252  // want we want it to do:
253  // - Normal compilation: -02 (default settings)
254  // - Aggressive mode: -03 (the whole lot!)
255  // In theory we could also target
256  // - Debug mode: -O0 (no optimizations)
257  // but this doesn't make much sense since SymEngine is a
258  // tested external library.
259  const bool use_agg_opt =
260  static_cast<int>(flags & OptimizationFlags::optimize_aggressive);
261  const int opt_level = (use_agg_opt ? 3 : 2);
262  return opt_level;
263  }
264  } // namespace internal
265 
266 
271  template <typename StreamType>
272  inline StreamType &
273  operator<<(StreamType &s, OptimizationFlags o)
274  {
275  s << " OptimizationFlags|";
276  if (static_cast<unsigned int>(o & OptimizationFlags::optimize_cse))
277  s << "cse|";
278 
279  // LLVM optimization level
281  "|";
282 
283  return s;
284  }
285 
286 
287  namespace internal
288  {
298  template <typename ReturnType, typename T = void>
300 
301 
311  template <typename ReturnType, typename T = void>
313 
314 
315 # ifdef HAVE_SYMENGINE_LLVM
325  template <typename ReturnType, typename T = void>
326  struct LLVMOptimizer;
327 # endif // HAVE_SYMENGINE_LLVM
328 
329 
345  template <typename ReturnType, typename Optimizer, typename T = void>
347 
348 
349 # ifndef DOXYGEN
350 
351 
352  /* ----------- Specializations for the Optimizers ----------- */
353 
354 
355  // A helper struct to type trait detection for the optimizers that
356  // will be defined next.
357  template <typename ReturnType_, typename T = void>
358  struct SupportedOptimizerTypeTraits
359  {
360  static const bool is_supported = false;
361 
362  using ReturnType = void;
363  };
364 
365 
366 
367  // Specialization for arithmetic types
368  template <typename ReturnType_>
369  struct SupportedOptimizerTypeTraits<
370  ReturnType_,
371  std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
372  {
373  static const bool is_supported = true;
374 
375  using ReturnType = typename std::
376  conditional<std::is_same_v<ReturnType_, float>, float, double>::type;
377  };
378 
379 
380 
381  // Specialization for complex arithmetic types
382  template <typename ReturnType_>
383  struct SupportedOptimizerTypeTraits<
384  ReturnType_,
385  std::enable_if_t<
386  boost::is_complex<ReturnType_>::value &&
387  std::is_arithmetic_v<typename ReturnType_::value_type>>>
388  {
389  static const bool is_supported = true;
390 
391  using ReturnType = typename std::conditional<
392  std::is_same_v<ReturnType_, std::complex<float>>,
393  std::complex<float>,
394  std::complex<double>>::type;
395  };
396 
397 
398 
399  template <typename ReturnType_>
400  struct DictionaryOptimizer<ReturnType_,
401  std::enable_if_t<SupportedOptimizerTypeTraits<
402  ReturnType_>::is_supported>>
403  {
404  using ReturnType =
405  typename SupportedOptimizerTypeTraits<ReturnType_>::ReturnType;
406  using OptimizerType =
407  internal::DictionarySubstitutionVisitor<ReturnType, SD::Expression>;
408 
409 
418  static void
419  initialize(OptimizerType &optimizer,
420  const SymEngine::vec_basic &independent_symbols,
421  const SymEngine::vec_basic &dependent_functions,
422  const enum OptimizationFlags &optimization_flags)
423  {
424  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
425  optimizer.init(independent_symbols,
426  dependent_functions,
427  use_symbolic_cse);
428  }
429 
430 
431 
436  template <class Archive>
437  static void
438  save(Archive &archive,
439  const unsigned int version,
440  OptimizerType &optimizer)
441  {
442  optimizer.save(archive, version);
443  }
444 
445 
446 
451  template <class Archive>
452  static void
453  load(Archive &archive,
454  const unsigned int version,
455  OptimizerType &optimizer,
456  const SymEngine::vec_basic & /*independent_symbols*/,
457  const SymEngine::vec_basic & /*dependent_functions*/,
458  const enum OptimizationFlags & /*optimization_flags*/)
459  {
460  optimizer.load(archive, version);
461  }
462 
463 
464 
480  template <typename Stream>
481  static void
482  print(Stream &stream,
483  const OptimizerType &optimizer,
484  const bool print_independent_symbols = false,
485  const bool print_dependent_functions = false,
486  const bool print_cse_reductions = true)
487  {
488  optimizer.print(stream,
489  print_independent_symbols,
490  print_dependent_functions,
491  print_cse_reductions);
492  }
493  };
494 
495 
496 
497  template <typename ReturnType_>
498  struct LambdaOptimizer<ReturnType_,
499  std::enable_if_t<SupportedOptimizerTypeTraits<
500  ReturnType_>::is_supported>>
501  {
502  using ReturnType =
503  typename std::conditional<!boost::is_complex<ReturnType_>::value,
504  double,
505  std::complex<double>>::type;
506  using OptimizerType = typename std::conditional<
507  !boost::is_complex<ReturnType_>::value,
508  SymEngine::LambdaRealDoubleVisitor,
509  SymEngine::LambdaComplexDoubleVisitor>::type;
510 
511 
520  static void
521  initialize(OptimizerType &optimizer,
522  const SymEngine::vec_basic &independent_symbols,
523  const SymEngine::vec_basic &dependent_functions,
524  const enum OptimizationFlags &optimization_flags)
525  {
526  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
527  optimizer.init(independent_symbols,
528  dependent_functions,
529  use_symbolic_cse);
530  }
531 
532 
533 
538  template <class Archive>
539  static void
540  save(Archive & /*archive*/,
541  const unsigned int /*version*/,
542  OptimizerType & /*optimizer*/)
543  {}
544 
545 
550  template <class Archive>
551  static void
552  load(Archive & /*archive*/,
553  const unsigned int /*version*/,
554  OptimizerType &optimizer,
555  const SymEngine::vec_basic &independent_symbols,
556  const SymEngine::vec_basic &dependent_functions,
557  const enum OptimizationFlags &optimization_flags)
558  {
559  initialize(optimizer,
560  independent_symbols,
561  dependent_functions,
562  optimization_flags);
563  }
564 
565 
566 
582  template <typename StreamType>
583  static void
584  print(StreamType & /*stream*/,
585  const OptimizerType & /*optimizer*/,
586  const bool /*print_independent_symbols*/ = false,
587  const bool /*print_dependent_functions*/ = false,
588  const bool /*print_cse_reductions*/ = true)
589  {
590  // No built-in print function
591  }
592  };
593 
594 
595 
596 # ifdef HAVE_SYMENGINE_LLVM
597  template <typename ReturnType_>
598  struct LLVMOptimizer<ReturnType_,
599  std::enable_if_t<std::is_arithmetic_v<ReturnType_>>>
600  {
601  using ReturnType = typename std::
602  conditional<std::is_same_v<ReturnType_, float>, float, double>::type;
603  using OptimizerType =
604  typename std::conditional<std::is_same_v<ReturnType_, float>,
605  SymEngine::LLVMFloatVisitor,
606  SymEngine::LLVMDoubleVisitor>::type;
607 
612  static const bool supported_by_LLVM = true;
613 
614 
623  static void
624  initialize(OptimizerType &optimizer,
625  const SymEngine::vec_basic &independent_symbols,
626  const SymEngine::vec_basic &dependent_functions,
627  const enum OptimizationFlags &optimization_flags)
628  {
629  const int opt_level = get_LLVM_optimization_level(optimization_flags);
630  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
631  optimizer.init(independent_symbols,
632  dependent_functions,
633  use_symbolic_cse,
634  opt_level);
635  }
636 
637 
638 
643  template <class Archive>
644  static void
645  save(Archive &archive,
646  const unsigned int /*version*/,
647  OptimizerType &optimizer)
648  {
649  const std::string llvm_compiled_function = optimizer.dumps();
650  archive &llvm_compiled_function;
651  }
652 
653 
654 
659  template <class Archive>
660  static void
661  load(Archive &archive,
662  const unsigned int /*version*/,
663  OptimizerType &optimizer,
664  const SymEngine::vec_basic & /*independent_symbols*/,
665  const SymEngine::vec_basic & /*dependent_functions*/,
666  const enum OptimizationFlags & /*optimization_flags*/)
667  {
668  std::string llvm_compiled_function;
669  archive &llvm_compiled_function;
670  optimizer.loads(llvm_compiled_function);
671  }
672 
673 
674 
690  template <typename StreamType>
691  static void
692  print(StreamType & /*stream*/,
693  const OptimizerType & /*optimizer*/,
694  const bool /*print_independent_symbols*/ = false,
695  const bool /*print_dependent_functions*/ = false,
696  const bool /*print_cse_reductions*/ = true)
697  {
698  // No built-in print function
699  }
700  };
701 
702 
703  // There is no LLVM optimizer built with complex number support.
704  // So we fall back to the LambdaDouble case as a type (required
705  // at compile time), but offer no implementation. We expect that
706  // the calling class does not create this type: This can be done by
707  // checking the `supported_by_LLVM` flag.
708  template <typename ReturnType_>
709  struct LLVMOptimizer<
710  ReturnType_,
711  std::enable_if_t<
712  boost::is_complex<ReturnType_>::value &&
713  std::is_arithmetic_v<typename ReturnType_::value_type>>>
714  {
715  // Since there is no working implementation, these are dummy types
716  // that help with templating in the calling function.
717  using ReturnType = typename LambdaOptimizer<ReturnType_>::ReturnType;
718  using OptimizerType =
720 
725  static const bool supported_by_LLVM = false;
726 
727 
736  static void
737  initialize(OptimizerType & /*optimizer*/,
738  const SymEngine::vec_basic & /*independent_symbols*/,
739  const SymEngine::vec_basic & /*dependent_functions*/,
740  const enum OptimizationFlags & /*optimization_flags*/)
741  {
742  AssertThrow(false, ExcNotImplemented());
743  }
744 
745 
746 
751  template <class Archive>
752  static void
753  save(Archive & /*archive*/,
754  const unsigned int /*version*/,
755  OptimizerType & /*optimizer*/)
756  {
757  AssertThrow(false, ExcNotImplemented());
758  }
759 
760 
761 
766  template <class Archive>
767  static void
768  load(Archive & /*archive*/,
769  const unsigned int /*version*/,
770  OptimizerType & /*optimizer*/,
771  const SymEngine::vec_basic & /*independent_symbols*/,
772  const SymEngine::vec_basic & /*dependent_functions*/,
773  const enum OptimizationFlags & /*optimization_flags*/)
774  {
775  AssertThrow(false, ExcNotImplemented());
776  }
777 
778 
779 
795  template <typename StreamType>
796  static void
797  print(StreamType & /*stream*/,
798  const OptimizerType & /*optimizer*/,
799  const bool /*print_independent_symbols*/ = false,
800  const bool /*print_dependent_functions*/ = false,
801  const bool /*print_cse_reductions*/ = true)
802  {
803  AssertThrow(false, ExcNotImplemented());
804  }
805  };
806 # endif // HAVE_SYMENGINE_LLVM
807 
808 
809  /* ----------- Specializations for OptimizerHelper ----------- */
810 
811 
812  template <typename ReturnType, typename Optimizer>
813  struct OptimizerHelper<
814  ReturnType,
815  Optimizer,
816  std::enable_if_t<
817  std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
818  {
827  static void
828  initialize(typename Optimizer::OptimizerType *optimizer,
829  const SymEngine::vec_basic &independent_symbols,
830  const SymEngine::vec_basic &dependent_functions,
831  const enum OptimizationFlags &optimization_flags)
832  {
833  Assert(optimizer, ExcNotInitialized());
834 
835  // Some optimizers don't have the same interface for
836  // initialization, we filter them out through the specializations
837  // of the Optimizer class
838  Optimizer::initialize(*optimizer,
839  independent_symbols,
840  dependent_functions,
841  optimization_flags);
842  }
843 
844 
845 
859  static void
860  substitute(typename Optimizer::OptimizerType *optimizer,
861  std::vector<ReturnType> &output_values,
862  const std::vector<ReturnType> &substitution_values)
863  {
864  Assert(optimizer, ExcNotInitialized());
865  optimizer->call(output_values.data(), substitution_values.data());
866  }
867 
868 
869 
874  template <class Archive>
875  static void
876  save(Archive &archive,
877  const unsigned int version,
878  typename Optimizer::OptimizerType *optimizer)
879  {
880  Assert(optimizer, ExcNotInitialized());
881 
882  // Some optimizers don't have the same interface for
883  // serialization, we filter them out through the specializations
884  // of the Optimizer class
885  Optimizer::save(archive, version, *optimizer);
886  }
887 
888 
889 
894  template <class Archive>
895  static void
896  load(Archive &archive,
897  const unsigned int version,
898  typename Optimizer::OptimizerType *optimizer,
899  const SymEngine::vec_basic &independent_symbols,
900  const SymEngine::vec_basic &dependent_functions,
901  const enum OptimizationFlags &optimization_flags)
902  {
903  Assert(optimizer, ExcNotInitialized());
904 
905  // Some optimizers don't have the same interface for
906  // serialization, we filter them out through the specializations
907  // of the Optimizer class
908  Optimizer::load(archive,
909  version,
910  *optimizer,
911  independent_symbols,
912  dependent_functions,
913  optimization_flags);
914  }
915 
916 
917 
933  template <typename Stream>
934  static void
935  print(Stream &stream,
936  typename Optimizer::OptimizerType *optimizer,
937  const bool print_independent_symbols = false,
938  const bool print_dependent_functions = false,
939  const bool print_cse_reductions = true)
940  {
941  Assert(optimizer, ExcNotInitialized());
942 
943  // Some optimizers don't have a print function, so
944  // we filter them out through the specializations of
945  // the Optimizer class
946  Optimizer::print(stream,
947  *optimizer,
948  print_independent_symbols,
949  print_dependent_functions,
950  print_cse_reductions);
951  }
952  };
953 
954  template <typename ReturnType, typename Optimizer>
955  struct OptimizerHelper<
956  ReturnType,
957  Optimizer,
958  std::enable_if_t<
959  !std::is_same_v<ReturnType, typename Optimizer::ReturnType>>>
960  {
969  static void
970  initialize(typename Optimizer::OptimizerType *optimizer,
971  const SymEngine::vec_basic &independent_symbols,
972  const SymEngine::vec_basic &dependent_functions,
973  const enum OptimizationFlags &optimization_flags)
974  {
975  Assert(optimizer, ExcNotInitialized());
976 
977  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
978  optimizer->init(independent_symbols,
979  dependent_functions,
980  use_symbolic_cse);
981  }
982 
983 
984 
998  static void
999  substitute(typename Optimizer::OptimizerType *optimizer,
1000  std::vector<ReturnType> &output_values,
1001  const std::vector<ReturnType> &substitution_values)
1002  {
1003  Assert(optimizer, ExcNotInitialized());
1004 
1005  // Intermediate values to accommodate the difference in
1006  // value types.
1007  std::vector<typename Optimizer::ReturnType> int_outputs(
1008  output_values.size());
1009  std::vector<typename Optimizer::ReturnType> int_inputs(
1010  substitution_values.size());
1011 
1012  std::copy(substitution_values.begin(),
1013  substitution_values.end(),
1014  int_inputs.begin());
1015  optimizer->call(int_outputs.data(), int_inputs.data());
1016  std::copy(int_outputs.begin(),
1017  int_outputs.end(),
1018  output_values.begin());
1019  }
1020 
1021 
1022 
1027  template <class Archive>
1028  static void
1029  save(Archive &archive,
1030  const unsigned int version,
1031  typename Optimizer::OptimizerType *optimizer)
1032  {
1033  Assert(optimizer, ExcNotInitialized());
1034  Optimizer::save(archive, version, *optimizer);
1035  }
1036 
1037 
1038 
1043  template <class Archive>
1044  static void
1045  load(Archive &archive,
1046  const unsigned int version,
1047  typename Optimizer::OptimizerType *optimizer,
1048  const SymEngine::vec_basic &independent_symbols,
1049  const SymEngine::vec_basic &dependent_functions,
1050  const enum OptimizationFlags &optimization_flags)
1051  {
1052  Assert(optimizer, ExcNotInitialized());
1053 
1054  // Some optimizers don't have the same interface for
1055  // serialization, we filter them out through the specializations
1056  // of the Optimizer class
1057  Optimizer::load(archive,
1058  version,
1059  *optimizer,
1060  independent_symbols,
1061  dependent_functions,
1062  optimization_flags);
1063  }
1064 
1065 
1066 
1082  template <typename Stream>
1083  static void
1084  print(Stream &stream,
1085  typename Optimizer::OptimizerType *optimizer,
1086  const bool print_cse_reductions = true,
1087  const bool print_independent_symbols = false,
1088  const bool print_dependent_functions = false)
1089  {
1090  Assert(optimizer, ExcNotInitialized());
1091 
1092  optimizer->print(stream,
1093  print_independent_symbols,
1094  print_dependent_functions,
1095  print_cse_reductions);
1096  }
1097  };
1098 
1099 # endif // DOXYGEN
1100 
1101 
1102  /* -------------------- Utility functions ---------------------- */
1103 
1104 
1126  template <typename NumberType,
1127  int rank,
1128  int dim,
1129  template <int, int, typename>
1130  class TensorType>
1131  TensorType<rank, dim, NumberType>
1133  const TensorType<rank, dim, Expression> &symbol_tensor,
1134  const std::vector<NumberType> &cached_evaluation,
1135  const BatchOptimizer<NumberType> &optimizer)
1136  {
1137  TensorType<rank, dim, NumberType> out;
1138  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1139  {
1140  const TableIndices<rank> indices(
1141  out.unrolled_to_component_indices(i));
1142  out[indices] =
1143  optimizer.extract(symbol_tensor[indices], cached_evaluation);
1144  }
1145  return out;
1146  }
1147 
1148 
1171  template <typename NumberType, int dim>
1174  const SymmetricTensor<4, dim, Expression> &symbol_tensor,
1175  const std::vector<NumberType> &cached_evaluation,
1176  const BatchOptimizer<NumberType> &optimizer)
1177  {
1179  for (unsigned int i = 0;
1180  i < SymmetricTensor<2, dim>::n_independent_components;
1181  ++i)
1182  for (unsigned int j = 0;
1183  j < SymmetricTensor<2, dim>::n_independent_components;
1184  ++j)
1185  {
1186  const TableIndices<4> indices =
1187  make_rank_4_tensor_indices<dim>(i, j);
1188  out[indices] =
1189  optimizer.extract(symbol_tensor[indices], cached_evaluation);
1190  }
1191  return out;
1192  }
1193 
1194 
1212  template <typename NumberType, typename T>
1213  void
1215  const T &function)
1216  {
1217  optimizer.register_function(function);
1218  }
1219 
1220 
1238  template <typename NumberType, typename T>
1239  void
1241  const std::vector<T> &functions)
1242  {
1243  for (const auto &function : functions)
1244  register_functions(optimizer, function);
1245  }
1246 
1247 
1267  template <typename NumberType, typename T, typename... Args>
1268  void
1270  const T &function,
1271  const Args &...other_functions)
1272  {
1273  register_functions(optimizer, function);
1274  register_functions(optimizer, other_functions...);
1275  }
1276 
1277 
1289  template <int rank,
1290  int dim,
1291  template <int, int, typename>
1292  class TensorType>
1295  const TensorType<rank, dim, Expression> &symbol_tensor)
1296  {
1298  out.reserve(symbol_tensor.n_independent_components);
1299  for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
1300  ++i)
1301  {
1302  const TableIndices<rank> indices(
1303  symbol_tensor.unrolled_to_component_indices(i));
1304  out.push_back(symbol_tensor[indices].get_RCP());
1305  }
1306  return out;
1307  }
1308 
1309 
1319  template <int dim>
1322  const SymmetricTensor<4, dim, Expression> &symbol_tensor)
1323  {
1325  out.reserve(symbol_tensor.n_independent_components);
1326  for (unsigned int i = 0;
1327  i < SymmetricTensor<2, dim>::n_independent_components;
1328  ++i)
1329  for (unsigned int j = 0;
1330  j < SymmetricTensor<2, dim>::n_independent_components;
1331  ++j)
1332  {
1333  const TableIndices<4> indices =
1334  make_rank_4_tensor_indices<dim>(i, j);
1335  out.push_back(symbol_tensor[indices].get_RCP());
1336  }
1337  return out;
1338  }
1339 
1340  } // namespace internal
1341 
1342 
1343 
1434  template <typename ReturnType>
1436  {
1437  public:
1445  BatchOptimizer();
1446 
1467 
1477  BatchOptimizer(const BatchOptimizer &other);
1478 
1482  BatchOptimizer(BatchOptimizer &&) noexcept = default;
1483 
1487  ~BatchOptimizer() = default;
1488 
1499  void
1500  copy_from(const BatchOptimizer &other);
1501 
1511  template <typename Stream>
1512  void
1513  print(Stream &stream, const bool print_cse = false) const;
1514 
1523  template <class Archive>
1524  void
1525  save(Archive &archive, const unsigned int version) const;
1526 
1540  template <class Archive>
1541  void
1542  load(Archive &archive, const unsigned int version);
1543 
1544 # ifdef DOXYGEN
1566  template <class Archive>
1567  void
1568  serialize(Archive &archive, const unsigned int version);
1569 # else
1570  // This macro defines the serialize() method that is compatible with
1571  // the templated save() and load() method that have been implemented.
1572  BOOST_SERIALIZATION_SPLIT_MEMBER()
1573 # endif
1574 
1585  void
1587 
1593  void
1594  register_symbols(const SymEngine::map_basic_basic &substitution_map);
1595 
1606  void
1607  register_symbols(const types::symbol_vector &symbols);
1608 
1619  void
1620  register_symbols(const SymEngine::vec_basic &symbols);
1621 
1627  get_independent_symbols() const;
1628 
1634  std::size_t
1635  n_independent_variables() const;
1636 
1648  void
1649  register_function(const Expression &function);
1650 
1655  template <int rank, int dim>
1656  void
1658 
1663  template <int rank, int dim>
1664  void
1666  const SymmetricTensor<rank, dim, Expression> &function_tensor);
1667 
1672  void
1674 
1679  void
1680  register_functions(const SymEngine::vec_basic &functions);
1681 
1691  template <typename T>
1692  void
1693  register_functions(const std::vector<T> &functions);
1694 
1708  template <typename T, typename... Args>
1709  void
1710  register_functions(const T &functions, const Args &...other_functions);
1711 
1716  const types::symbol_vector &
1717  get_dependent_functions() const;
1718 
1725  std::size_t
1726  n_dependent_variables() const;
1727 
1746  void
1750 
1755  enum OptimizerType
1756  optimization_method() const;
1757 
1762  enum OptimizationFlags
1763  optimization_flags() const;
1764 
1770  bool
1771  use_symbolic_CSE() const;
1772 
1788  void
1789  optimize();
1790 
1795  bool
1796  optimized() const;
1797 
1814  void
1816 
1826  void
1827  substitute(const SymEngine::map_basic_basic &substitution_map) const;
1828 
1839  void
1840  substitute(const types::symbol_vector &symbols,
1841  const std::vector<ReturnType> &values) const;
1842 
1853  void
1854  substitute(const SymEngine::vec_basic &symbols,
1855  const std::vector<ReturnType> &values) const;
1856 
1862  bool
1863  values_substituted() const;
1864 
1894  const std::vector<ReturnType> &
1895  evaluate() const;
1896 
1904  ReturnType
1905  evaluate(const Expression &func) const;
1906 
1915  std::vector<ReturnType>
1916  evaluate(const std::vector<Expression> &funcs) const;
1917 
1926  template <int rank, int dim>
1929 
1930 
1939  template <int rank, int dim>
1942 
1943 
1951  ReturnType
1952  extract(const Expression &func,
1953  const std::vector<ReturnType> &cached_evaluation) const;
1954 
1955 
1963  std::vector<ReturnType>
1964  extract(const std::vector<Expression> &funcs,
1965  const std::vector<ReturnType> &cached_evaluation) const;
1966 
1967 
1975  template <int rank, int dim>
1978  const std::vector<ReturnType> &cached_evaluation) const;
1979 
1980 
1988  template <int rank, int dim>
1991  const std::vector<ReturnType> &cached_evaluation) const;
1992 
1995  private:
1999  enum OptimizerType method;
2000 
2005  enum OptimizationFlags flags;
2006 
2017 
2024 
2029  bool
2031  const SD::Expression &function) const;
2032 
2037  bool
2039  const SymEngine::RCP<const SymEngine::Basic> &function) const;
2040 
2055  mutable std::vector<ReturnType> dependent_variables_output;
2056 
2066  std::map<SD::Expression,
2067  std::size_t,
2069 
2075 
2082  mutable std::unique_ptr<SymEngine::Visitor> optimizer;
2083 
2093 
2098  mutable bool has_been_serialized;
2099 
2103  void
2104  register_scalar_function(const SD::Expression &function);
2105 
2110  void
2112 
2116  void
2117  create_optimizer(std::unique_ptr<SymEngine::Visitor> &optimizer);
2118 
2135  void
2136  substitute(const std::vector<ReturnType> &substitution_values) const;
2137  };
2138 
2139 
2140 
2141  /* -------------------- inline and template functions ------------------ */
2142 
2143 
2144 # ifndef DOXYGEN
2145 
2146 
2147  template <typename ReturnType>
2148  template <typename Stream>
2149  void
2150  BatchOptimizer<ReturnType>::print(Stream &stream,
2151  const bool /*print_cse*/) const
2152  {
2153  // Settings
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";
2158 
2159  // Independent variables
2160  stream << "Symbols (" << n_independent_variables()
2161  << " independent variables):" << '\n';
2162  int cntr = 0;
2163  for (SD::types::substitution_map::const_iterator it =
2164  independent_variables_symbols.begin();
2165  it != independent_variables_symbols.end();
2166  ++it, ++cntr)
2167  {
2168  stream << cntr << ": " << it->first << '\n';
2169  }
2170  stream << '\n' << std::flush;
2171 
2172  // Dependent functions
2173  stream << "Functions (" << n_dependent_variables()
2174  << " dependent variables):" << '\n';
2175  cntr = 0;
2176  for (typename SD::types::symbol_vector::const_iterator it =
2177  dependent_variables_functions.begin();
2178  it != dependent_variables_functions.end();
2179  ++it, ++cntr)
2180  {
2181  stream << cntr << ": " << (*it) << '\n';
2182  }
2183  stream << '\n' << std::flush;
2184 
2185  // Common subexpression
2186  if (optimized() == true && use_symbolic_CSE() == true)
2187  {
2188  Assert(optimizer, ExcNotInitialized());
2189  const bool print_cse_reductions = true;
2190  const bool print_independent_symbols = false;
2191  const bool print_dependent_functions = false;
2192 
2193  if (optimization_method() == OptimizerType::dictionary)
2194  {
2195  Assert(dynamic_cast<typename internal::DictionaryOptimizer<
2196  ReturnType>::OptimizerType *>(optimizer.get()),
2197  ExcMessage("Cannot cast optimizer to Dictionary type."));
2198 
2199  internal::OptimizerHelper<
2200  ReturnType,
2201  internal::DictionaryOptimizer<ReturnType>>::
2202  print(stream,
2203  dynamic_cast<typename internal::DictionaryOptimizer<
2204  ReturnType>::OptimizerType *>(optimizer.get()),
2205  print_independent_symbols,
2206  print_dependent_functions,
2207  print_cse_reductions);
2208 
2209  stream << '\n' << std::flush;
2210  }
2211  else if (optimization_method() == OptimizerType::lambda)
2212  {
2213  Assert(dynamic_cast<typename internal::LambdaOptimizer<
2214  ReturnType>::OptimizerType *>(optimizer.get()),
2215  ExcMessage("Cannot cast optimizer to Lambda type."));
2216 
2217  internal::OptimizerHelper<ReturnType,
2218  internal::LambdaOptimizer<ReturnType>>::
2219  print(stream,
2220  dynamic_cast<typename internal::LambdaOptimizer<
2221  ReturnType>::OptimizerType *>(optimizer.get()),
2222  print_independent_symbols,
2223  print_dependent_functions,
2224  print_cse_reductions);
2225  }
2226 # ifdef HAVE_SYMENGINE_LLVM
2227  else if (optimization_method() == OptimizerType::llvm)
2228  {
2229  Assert(dynamic_cast<typename internal::LLVMOptimizer<
2230  ReturnType>::OptimizerType *>(optimizer.get()),
2231  ExcMessage("Cannot cast optimizer to LLVM type."));
2232 
2233  internal::OptimizerHelper<ReturnType,
2234  internal::LLVMOptimizer<ReturnType>>::
2235  print(stream,
2236  dynamic_cast<typename internal::LLVMOptimizer<
2237  ReturnType>::OptimizerType *>(optimizer.get()),
2238  print_independent_symbols,
2239  print_dependent_functions,
2240  print_cse_reductions);
2241  }
2242 # endif // HAVE_SYMENGINE_LLVM
2243  else
2244  {
2245  AssertThrow(false, ExcMessage("Unknown optimizer type."));
2246  }
2247  }
2248 
2249  if (values_substituted())
2250  {
2251  stream << "Evaluated functions:" << '\n';
2252  stream << std::flush;
2253  cntr = 0;
2254  for (typename std::vector<ReturnType>::const_iterator it =
2255  dependent_variables_output.begin();
2256  it != dependent_variables_output.end();
2257  ++it, ++cntr)
2258  {
2259  stream << cntr << ": " << (*it) << '\n';
2260  }
2261  stream << '\n' << std::flush;
2262  }
2263  }
2264 
2265 
2266 
2267  template <typename ReturnType>
2268  template <class Archive>
2269  void
2271  const unsigned int version) const
2272  {
2273  // Serialize enum classes...
2274  {
2275  const auto m =
2276  static_cast<typename std::underlying_type<OptimizerType>::type>(
2277  method);
2278  ar &m;
2279  }
2280  {
2281  const auto f =
2282  static_cast<typename std::underlying_type<OptimizationFlags>::type>(
2283  flags);
2284  ar &f;
2285  }
2286 
2287  // Important: Independent variables must always be
2288  // serialized before the dependent variables.
2289  ar &independent_variables_symbols;
2290  ar &dependent_variables_functions;
2291 
2292  ar &dependent_variables_output;
2293  ar &map_dep_expr_vec_entry;
2294  ar &ready_for_value_extraction;
2295 
2296  // Mark that we've saved this class at some point.
2297  has_been_serialized = true;
2298  ar &has_been_serialized;
2299 
2300  // When we serialize the optimizer itself, we have to (unfortunately)
2301  // provide it with sufficient information to rebuild itself from scratch.
2302  // This is because only two of the three optimization classes support
2303  // real serialization (i.e. have save/load capability).
2304  const SD::types::symbol_vector symbol_vec =
2305  Utilities::extract_symbols(independent_variables_symbols);
2307  *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2308  ReturnType>::OptimizerType *>(optimizer.get()))
2309  {
2310  Assert(optimization_method() == OptimizerType::dictionary,
2311  ExcInternalError());
2312  internal::OptimizerHelper<
2313  ReturnType,
2314  internal::DictionaryOptimizer<ReturnType>>::save(ar, version, opt);
2315  }
2317  *opt = dynamic_cast<typename internal::LambdaOptimizer<
2318  ReturnType>::OptimizerType *>(optimizer.get()))
2319  {
2320  Assert(optimization_method() == OptimizerType::lambda,
2321  ExcInternalError());
2322  internal::OptimizerHelper<
2323  ReturnType,
2324  internal::LambdaOptimizer<ReturnType>>::save(ar, version, opt);
2325  }
2326 # ifdef HAVE_SYMENGINE_LLVM
2328  *opt = dynamic_cast<typename internal::LLVMOptimizer<
2329  ReturnType>::OptimizerType *>(optimizer.get()))
2330  {
2331  Assert(optimization_method() == OptimizerType::llvm,
2332  ExcInternalError());
2333  internal::OptimizerHelper<
2334  ReturnType,
2335  internal::LLVMOptimizer<ReturnType>>::save(ar, version, opt);
2336  }
2337 # endif
2338  else
2339  {
2340  AssertThrow(false, ExcMessage("Unknown optimizer type."));
2341  }
2342  }
2343 
2344 
2345 
2346  template <typename ReturnType>
2347  template <class Archive>
2348  void
2349  BatchOptimizer<ReturnType>::load(Archive &ar, const unsigned int version)
2350  {
2351  Assert(independent_variables_symbols.empty(), ExcInternalError());
2352  Assert(dependent_variables_functions.empty(), ExcInternalError());
2353  Assert(dependent_variables_output.empty(), ExcInternalError());
2354  Assert(map_dep_expr_vec_entry.empty(), ExcInternalError());
2355  Assert(ready_for_value_extraction == false, ExcInternalError());
2356 
2357  // Deserialize enum classes...
2358  {
2359  typename std::underlying_type<OptimizerType>::type m;
2360  ar &m;
2361  method = static_cast<OptimizerType>(m);
2362  }
2363  {
2364  typename std::underlying_type<OptimizationFlags>::type f;
2365  ar &f;
2366  flags = static_cast<OptimizationFlags>(f);
2367  }
2368 
2369  // Important: Independent variables must always be
2370  // deserialized before the dependent variables.
2371  ar &independent_variables_symbols;
2372  ar &dependent_variables_functions;
2373 
2374  ar &dependent_variables_output;
2375  ar &map_dep_expr_vec_entry;
2376  ar &ready_for_value_extraction;
2377 
2378  ar &has_been_serialized;
2379 
2380  // If we're reading in data, then create the optimizer
2381  // and then deserialize it.
2382  Assert(!optimizer, ExcInternalError());
2383 
2384  // Create and configure the optimizer
2385  create_optimizer(optimizer);
2386  Assert(optimizer, ExcNotInitialized());
2387 
2388  // When we deserialize the optimizer itself, we have to (unfortunately)
2389  // provide it with sufficient information to rebuild itself from scratch.
2390  // This is because only two of the three optimization classes support
2391  // real serialization (i.e. have save/load capability).
2392  const SD::types::symbol_vector symbol_vec =
2393  Utilities::extract_symbols(independent_variables_symbols);
2395  *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2396  ReturnType>::OptimizerType *>(optimizer.get()))
2397  {
2398  Assert(optimization_method() == OptimizerType::dictionary,
2399  ExcInternalError());
2400  internal::OptimizerHelper<ReturnType,
2401  internal::DictionaryOptimizer<ReturnType>>::
2402  load(ar,
2403  version,
2404  opt,
2406  symbol_vec),
2408  dependent_variables_functions),
2409  optimization_flags());
2410  }
2412  *opt = dynamic_cast<typename internal::LambdaOptimizer<
2413  ReturnType>::OptimizerType *>(optimizer.get()))
2414  {
2415  Assert(optimization_method() == OptimizerType::lambda,
2416  ExcInternalError());
2417  internal::OptimizerHelper<ReturnType,
2418  internal::LambdaOptimizer<ReturnType>>::
2419  load(ar,
2420  version,
2421  opt,
2423  symbol_vec),
2425  dependent_variables_functions),
2426  optimization_flags());
2427  }
2428 # ifdef HAVE_SYMENGINE_LLVM
2430  *opt = dynamic_cast<typename internal::LLVMOptimizer<
2431  ReturnType>::OptimizerType *>(optimizer.get()))
2432  {
2433  Assert(optimization_method() == OptimizerType::llvm,
2434  ExcInternalError());
2435  internal::OptimizerHelper<ReturnType,
2436  internal::LLVMOptimizer<ReturnType>>::
2437  load(ar,
2438  version,
2439  opt,
2441  symbol_vec),
2443  dependent_variables_functions),
2444  optimization_flags());
2445  }
2446 # endif
2447  else
2448  {
2449  AssertThrow(false, ExcMessage("Unknown optimizer type."));
2450  }
2451  }
2452 
2453 
2454 
2455  template <typename ReturnType>
2456  template <int rank, int dim>
2457  void
2459  const Tensor<rank, dim, Expression> &function_tensor)
2460  {
2461  Assert(optimized() == false,
2462  ExcMessage(
2463  "Cannot register functions once the optimizer is finalised."));
2464 
2465  register_vector_functions(
2466  internal::unroll_to_expression_vector(function_tensor));
2467  }
2468 
2469 
2470 
2471  template <typename ReturnType>
2472  template <int rank, int dim>
2473  void
2475  const SymmetricTensor<rank, dim, Expression> &function_tensor)
2476  {
2477  Assert(optimized() == false,
2478  ExcMessage(
2479  "Cannot register functions once the optimizer is finalised."));
2480 
2481  register_vector_functions(
2482  internal::unroll_to_expression_vector(function_tensor));
2483  }
2484 
2485 
2486 
2487  template <typename ReturnType>
2488  template <typename T, typename... Args>
2489  void
2491  const T &functions,
2492  const Args &...other_functions)
2493  {
2495  internal::register_functions(*this, other_functions...);
2496  }
2497 
2498 
2499 
2500  template <typename ReturnType>
2501  template <typename T>
2502  void
2504  const std::vector<T> &functions)
2505  {
2507  }
2508 
2509 
2510 
2511  template <typename ReturnType>
2512  template <int rank, int dim>
2515  const Tensor<rank, dim, Expression> &funcs,
2516  const std::vector<ReturnType> &cached_evaluation) const
2517  {
2519  cached_evaluation,
2520  *this);
2521  }
2522 
2523 
2524 
2525  template <typename ReturnType>
2526  template <int rank, int dim>
2529  const Tensor<rank, dim, Expression> &funcs) const
2530  {
2531  Assert(
2532  values_substituted() == true,
2533  ExcMessage(
2534  "The optimizer is not configured to perform evaluation. "
2535  "This action can only performed after substitute() has been called."));
2536 
2537  return extract(funcs, dependent_variables_output);
2538  }
2539 
2540 
2541 
2542  template <typename ReturnType>
2543  template <int rank, int dim>
2547  const std::vector<ReturnType> &cached_evaluation) const
2548  {
2550  cached_evaluation,
2551  *this);
2552  }
2553 
2554 
2555 
2556  template <typename ReturnType>
2557  template <int rank, int dim>
2560  const SymmetricTensor<rank, dim, Expression> &funcs) const
2561  {
2562  Assert(
2563  values_substituted() == true,
2564  ExcMessage(
2565  "The optimizer is not configured to perform evaluation. "
2566  "This action can only performed after substitute() has been called."));
2567 
2568  return extract(funcs, dependent_variables_output);
2569  }
2570 
2571 # endif // DOXYGEN
2572 
2573  } // namespace SD
2574 } // namespace Differentiation
2575 
2576 
2578 
2579 #endif // DEAL_II_WITH_SYMENGINE
2580 
2581 #endif
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
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)
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
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)
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
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:490
static ::ExceptionBase & ExcSymEngineLLVMNotAvailable()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcSymEngineLLVMReturnTypeNotSupported()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
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)
static const char T
std::string to_string(const T &t)
Definition: patterns.h:2391
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)