Reference documentation for deal.II version GIT a189bc2bdf 2022-12-07 02:45: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 - 2022 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 
24 // Low level
25 # include <symengine/basic.h>
26 # include <symengine/dict.h>
27 # include <symengine/symengine_exception.h>
28 # include <symengine/symengine_rcp.h>
29 
30 // Optimization
31 # include <symengine/lambda_double.h>
32 # include <symengine/visitor.h>
33 # ifdef HAVE_SYMENGINE_LLVM
34 # include <symengine/llvm_double.h>
35 # endif
37 
38 # include <deal.II/base/exceptions.h>
39 # include <deal.II/base/logstream.h>
40 # include <deal.II/base/utilities.h>
41 
47 
48 # include <boost/serialization/split_member.hpp>
49 # include <boost/type_traits.hpp>
50 
51 # include <algorithm>
52 # include <map>
53 # include <memory>
54 # include <type_traits>
55 # include <utility>
56 # include <vector>
57 
58 
60 
61 
62 namespace Differentiation
63 {
64  namespace SD
65  {
76  "SymEngine has not been built with LLVM support.");
77 
83  "The SymEngine LLVM optimizer does not (yet) support the "
84  "selected return type.");
85 
89  // Forward declarations
90  template <typename ReturnType>
91  class BatchOptimizer;
92 
93 
99  enum class OptimizerType
100  {
104  dictionary,
109  lambda,
114  llvm
115  };
116 
117 
121  template <class StreamType>
122  inline StreamType &
123  operator<<(StreamType &s, OptimizerType o)
124  {
125  if (o == OptimizerType::dictionary)
126  s << "dictionary";
127  else if (o == OptimizerType::lambda)
128  s << "lambda";
129  else if (o == OptimizerType::llvm)
130  s << "llvm";
131  else
132  {
133  Assert(false, ExcMessage("Unknown optimization method."));
134  }
135 
136  return s;
137  }
138 
139 
145  enum class OptimizationFlags : unsigned char
146  {
150  optimize_default = 0,
154  optimize_cse = 0x0001,
159  optimize_aggressive = 0x0002,
164  };
165 
166 
175  // This operator exists since if it did not then the result of the bit-or
176  // <tt>operator |</tt> would be an integer which would in turn trigger a
177  // compiler warning when we tried to assign it to an object of type
178  // OptimizationFlags.
179  inline OptimizationFlags
181  {
182  return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) |
183  static_cast<unsigned int>(f2));
184  }
185 
186 
191  inline OptimizationFlags &
193  {
194  f1 = f1 | f2;
195  return f1;
196  }
197 
198 
207  // This operator exists since if it did not then the result of the bit-or
208  // <tt>operator |</tt> would be an integer which would in turn trigger a
209  // compiler warning when we tried to assign it to an object of type
210  // OptimizationFlags.
211  inline OptimizationFlags
213  {
214  return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) &
215  static_cast<unsigned int>(f2));
216  }
217 
218 
223  inline OptimizationFlags &
225  {
226  f1 = f1 & f2;
227  return f1;
228  }
229 
230 
231  namespace internal
232  {
237  inline bool
239  {
240  return static_cast<int>(flags & OptimizationFlags::optimize_cse);
241  }
242 
247  inline int
249  {
250  // With the LLVM compiler there exists the opportunity to tune
251  // the level of optimizations performed during compilation.
252  // By default SymEngine sets this at "opt_level=2", which one
253  // presumes targets -O2. Here we are a bit more specific about
254  // want we want it to do:
255  // - Normal compilation: -02 (default settings)
256  // - Aggressive mode: -03 (the whole lot!)
257  // In theory we could also target
258  // - Debug mode: -O0 (no optimizations)
259  // but this doesn't make much sense since SymEngine is a
260  // tested external library.
261  const bool use_agg_opt =
262  static_cast<int>(flags & OptimizationFlags::optimize_aggressive);
263  const int opt_level = (use_agg_opt ? 3 : 2);
264  return opt_level;
265  }
266  } // namespace internal
267 
268 
273  template <class StreamType>
274  inline StreamType &
275  operator<<(StreamType &s, OptimizationFlags o)
276  {
277  s << " OptimizationFlags|";
278  if (static_cast<unsigned int>(o & OptimizationFlags::optimize_cse))
279  s << "cse|";
280 
281  // LLVM optimization level
283  "|";
284 
285  return s;
286  }
287 
288 
289  namespace internal
290  {
300  template <typename ReturnType, typename T = void>
302 
303 
313  template <typename ReturnType, typename T = void>
315 
316 
317 # ifdef HAVE_SYMENGINE_LLVM
327  template <typename ReturnType, typename T = void>
328  struct LLVMOptimizer;
329 # endif // HAVE_SYMENGINE_LLVM
330 
331 
347  template <typename ReturnType, typename Optimizer, typename T = void>
349 
350 
351 # ifndef DOXYGEN
352 
353 
354  /* ----------- Specializations for the Optimizers ----------- */
355 
356 
357  // A helper struct to type trait detection for the optimizers that
358  // will be defined next.
359  template <typename ReturnType_, typename T = void>
360  struct SupportedOptimizerTypeTraits
361  {
362  static const bool is_supported = false;
363 
364  using ReturnType = void;
365  };
366 
367 
368 
369  // Specialization for arithmetic types
370  template <typename ReturnType_>
371  struct SupportedOptimizerTypeTraits<
372  ReturnType_,
373  std::enable_if_t<std::is_arithmetic<ReturnType_>::value>>
374  {
375  static const bool is_supported = true;
376 
377  using ReturnType =
378  typename std::conditional<std::is_same<ReturnType_, float>::value,
379  float,
380  double>::type;
381  };
382 
383 
384 
385  // Specialization for complex arithmetic types
386  template <typename ReturnType_>
387  struct SupportedOptimizerTypeTraits<
388  ReturnType_,
389  std::enable_if_t<
390  boost::is_complex<ReturnType_>::value &&
391  std::is_arithmetic<typename ReturnType_::value_type>::value>>
392  {
393  static const bool is_supported = true;
394 
395  using ReturnType = typename std::conditional<
396  std::is_same<ReturnType_, std::complex<float>>::value,
397  std::complex<float>,
398  std::complex<double>>::type;
399  };
400 
401 
402 
403  template <typename ReturnType_>
404  struct DictionaryOptimizer<ReturnType_,
405  std::enable_if_t<SupportedOptimizerTypeTraits<
406  ReturnType_>::is_supported>>
407  {
408  using ReturnType =
409  typename SupportedOptimizerTypeTraits<ReturnType_>::ReturnType;
410  using OptimizerType =
411  internal::DictionarySubstitutionVisitor<ReturnType, SD::Expression>;
412 
413 
422  static void
423  initialize(OptimizerType & optimizer,
424  const SymEngine::vec_basic & independent_symbols,
425  const SymEngine::vec_basic & dependent_functions,
426  const enum OptimizationFlags &optimization_flags)
427  {
428  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
429  optimizer.init(independent_symbols,
430  dependent_functions,
431  use_symbolic_cse);
432  }
433 
434 
435 
440  template <class Archive>
441  static void
442  save(Archive & archive,
443  const unsigned int version,
444  OptimizerType & optimizer)
445  {
446  optimizer.save(archive, version);
447  }
448 
449 
450 
455  template <class Archive>
456  static void
457  load(Archive & archive,
458  const unsigned int version,
459  OptimizerType & optimizer,
460  const SymEngine::vec_basic & /*independent_symbols*/,
461  const SymEngine::vec_basic & /*dependent_functions*/,
462  const enum OptimizationFlags & /*optimization_flags*/)
463  {
464  optimizer.load(archive, version);
465  }
466 
467 
468 
484  template <typename Stream>
485  static void
486  print(Stream & stream,
487  const OptimizerType &optimizer,
488  const bool print_independent_symbols = false,
489  const bool print_dependent_functions = false,
490  const bool print_cse_reductions = true)
491  {
492  optimizer.print(stream,
493  print_independent_symbols,
494  print_dependent_functions,
495  print_cse_reductions);
496  }
497  };
498 
499 
500 
501  template <typename ReturnType_>
502  struct LambdaOptimizer<ReturnType_,
503  std::enable_if_t<SupportedOptimizerTypeTraits<
504  ReturnType_>::is_supported>>
505  {
506  using ReturnType =
507  typename std::conditional<!boost::is_complex<ReturnType_>::value,
508  double,
509  std::complex<double>>::type;
510  using OptimizerType = typename std::conditional<
511  !boost::is_complex<ReturnType_>::value,
512  SymEngine::LambdaRealDoubleVisitor,
513  SymEngine::LambdaComplexDoubleVisitor>::type;
514 
515 
524  static void
525  initialize(OptimizerType & optimizer,
526  const SymEngine::vec_basic & independent_symbols,
527  const SymEngine::vec_basic & dependent_functions,
528  const enum OptimizationFlags &optimization_flags)
529  {
530  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
531  optimizer.init(independent_symbols,
532  dependent_functions,
533  use_symbolic_cse);
534  }
535 
536 
537 
542  template <class Archive>
543  static void
544  save(Archive & /*archive*/,
545  const unsigned int /*version*/,
546  OptimizerType & /*optimizer*/)
547  {}
548 
549 
554  template <class Archive>
555  static void
556  load(Archive & /*archive*/,
557  const unsigned int /*version*/,
558  OptimizerType & optimizer,
559  const SymEngine::vec_basic & independent_symbols,
560  const SymEngine::vec_basic & dependent_functions,
561  const enum OptimizationFlags &optimization_flags)
562  {
563  initialize(optimizer,
564  independent_symbols,
565  dependent_functions,
566  optimization_flags);
567  }
568 
569 
570 
586  template <typename StreamType>
587  static void
588  print(StreamType & /*stream*/,
589  const OptimizerType & /*optimizer*/,
590  const bool /*print_independent_symbols*/ = false,
591  const bool /*print_dependent_functions*/ = false,
592  const bool /*print_cse_reductions*/ = true)
593  {
594  // No built-in print function
595  }
596  };
597 
598 
599 
600 # ifdef HAVE_SYMENGINE_LLVM
601  template <typename ReturnType_>
602  struct LLVMOptimizer<
603  ReturnType_,
604  std::enable_if_t<std::is_arithmetic<ReturnType_>::value>>
605  {
606  using ReturnType =
607  typename std::conditional<std::is_same<ReturnType_, float>::value,
608  float,
609  double>::type;
610  using OptimizerType =
611  typename std::conditional<std::is_same<ReturnType_, float>::value,
612  SymEngine::LLVMFloatVisitor,
613  SymEngine::LLVMDoubleVisitor>::type;
614 
619  static const bool supported_by_LLVM = true;
620 
621 
630  static void
631  initialize(OptimizerType & optimizer,
632  const SymEngine::vec_basic & independent_symbols,
633  const SymEngine::vec_basic & dependent_functions,
634  const enum OptimizationFlags &optimization_flags)
635  {
636  const int opt_level = get_LLVM_optimization_level(optimization_flags);
637  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
638  optimizer.init(independent_symbols,
639  dependent_functions,
640  use_symbolic_cse,
641  opt_level);
642  }
643 
644 
645 
650  template <class Archive>
651  static void
652  save(Archive &archive,
653  const unsigned int /*version*/,
654  OptimizerType &optimizer)
655  {
656  const std::string llvm_compiled_function = optimizer.dumps();
657  archive & llvm_compiled_function;
658  }
659 
660 
661 
666  template <class Archive>
667  static void
668  load(Archive &archive,
669  const unsigned int /*version*/,
670  OptimizerType &optimizer,
671  const SymEngine::vec_basic & /*independent_symbols*/,
672  const SymEngine::vec_basic & /*dependent_functions*/,
673  const enum OptimizationFlags & /*optimization_flags*/)
674  {
675  std::string llvm_compiled_function;
676  archive & llvm_compiled_function;
677  optimizer.loads(llvm_compiled_function);
678  }
679 
680 
681 
697  template <typename StreamType>
698  static void
699  print(StreamType & /*stream*/,
700  const OptimizerType & /*optimizer*/,
701  const bool /*print_independent_symbols*/ = false,
702  const bool /*print_dependent_functions*/ = false,
703  const bool /*print_cse_reductions*/ = true)
704  {
705  // No built-in print function
706  }
707  };
708 
709 
710  // There is no LLVM optimizer built with complex number support.
711  // So we fall back to the LambdaDouble case as a type (required
712  // at compile time), but offer no implementation. We expect that
713  // the calling class does not create this type: This can be done by
714  // checking the `supported_by_LLVM` flag.
715  template <typename ReturnType_>
716  struct LLVMOptimizer<
717  ReturnType_,
718  std::enable_if_t<
719  boost::is_complex<ReturnType_>::value &&
720  std::is_arithmetic<typename ReturnType_::value_type>::value>>
721  {
722  // Since there is no working implementation, these are dummy types
723  // that help with templating in the calling function.
724  using ReturnType = typename LambdaOptimizer<ReturnType_>::ReturnType;
725  using OptimizerType =
727 
732  static const bool supported_by_LLVM = false;
733 
734 
743  static void
744  initialize(OptimizerType & /*optimizer*/,
745  const SymEngine::vec_basic & /*independent_symbols*/,
746  const SymEngine::vec_basic & /*dependent_functions*/,
747  const enum OptimizationFlags & /*optimization_flags*/)
748  {
749  AssertThrow(false, ExcNotImplemented());
750  }
751 
752 
753 
758  template <class Archive>
759  static void
760  save(Archive & /*archive*/,
761  const unsigned int /*version*/,
762  OptimizerType & /*optimizer*/)
763  {
764  AssertThrow(false, ExcNotImplemented());
765  }
766 
767 
768 
773  template <class Archive>
774  static void
775  load(Archive & /*archive*/,
776  const unsigned int /*version*/,
777  OptimizerType & /*optimizer*/,
778  const SymEngine::vec_basic & /*independent_symbols*/,
779  const SymEngine::vec_basic & /*dependent_functions*/,
780  const enum OptimizationFlags & /*optimization_flags*/)
781  {
782  AssertThrow(false, ExcNotImplemented());
783  }
784 
785 
786 
802  template <typename StreamType>
803  static void
804  print(StreamType & /*stream*/,
805  const OptimizerType & /*optimizer*/,
806  const bool /*print_independent_symbols*/ = false,
807  const bool /*print_dependent_functions*/ = false,
808  const bool /*print_cse_reductions*/ = true)
809  {
810  AssertThrow(false, ExcNotImplemented());
811  }
812  };
813 # endif // HAVE_SYMENGINE_LLVM
814 
815 
816  /* ----------- Specializations for OptimizerHelper ----------- */
817 
818 
819  template <typename ReturnType, typename Optimizer>
820  struct OptimizerHelper<
821  ReturnType,
822  Optimizer,
823  std::enable_if_t<
824  std::is_same<ReturnType, typename Optimizer::ReturnType>::value>>
825  {
834  static void
835  initialize(typename Optimizer::OptimizerType *optimizer,
836  const SymEngine::vec_basic & independent_symbols,
837  const SymEngine::vec_basic & dependent_functions,
838  const enum OptimizationFlags & optimization_flags)
839  {
840  Assert(optimizer, ExcNotInitialized());
841 
842  // Some optimizers don't have the same interface for
843  // initialization, we filter them out through the specializations
844  // of the Optimizer class
845  Optimizer::initialize(*optimizer,
846  independent_symbols,
847  dependent_functions,
848  optimization_flags);
849  }
850 
851 
852 
866  static void
867  substitute(typename Optimizer::OptimizerType *optimizer,
868  std::vector<ReturnType> & output_values,
869  const std::vector<ReturnType> & substitution_values)
870  {
871  Assert(optimizer, ExcNotInitialized());
872  optimizer->call(output_values.data(), substitution_values.data());
873  }
874 
875 
876 
881  template <class Archive>
882  static void
883  save(Archive & archive,
884  const unsigned int version,
885  typename Optimizer::OptimizerType *optimizer)
886  {
887  Assert(optimizer, ExcNotInitialized());
888 
889  // Some optimizers don't have the same interface for
890  // serialization, we filter them out through the specializations
891  // of the Optimizer class
892  Optimizer::save(archive, version, *optimizer);
893  }
894 
895 
896 
901  template <class Archive>
902  static void
903  load(Archive & archive,
904  const unsigned int version,
905  typename Optimizer::OptimizerType *optimizer,
906  const SymEngine::vec_basic & independent_symbols,
907  const SymEngine::vec_basic & dependent_functions,
908  const enum OptimizationFlags & optimization_flags)
909  {
910  Assert(optimizer, ExcNotInitialized());
911 
912  // Some optimizers don't have the same interface for
913  // serialization, we filter them out through the specializations
914  // of the Optimizer class
915  Optimizer::load(archive,
916  version,
917  *optimizer,
918  independent_symbols,
919  dependent_functions,
920  optimization_flags);
921  }
922 
923 
924 
940  template <typename Stream>
941  static void
942  print(Stream & stream,
943  typename Optimizer::OptimizerType *optimizer,
944  const bool print_independent_symbols = false,
945  const bool print_dependent_functions = false,
946  const bool print_cse_reductions = true)
947  {
948  Assert(optimizer, ExcNotInitialized());
949 
950  // Some optimizers don't have a print function, so
951  // we filter them out through the specializations of
952  // the Optimizer class
953  Optimizer::print(stream,
954  *optimizer,
955  print_independent_symbols,
956  print_dependent_functions,
957  print_cse_reductions);
958  }
959  };
960 
961  template <typename ReturnType, typename Optimizer>
962  struct OptimizerHelper<
963  ReturnType,
964  Optimizer,
965  std::enable_if_t<
966  !std::is_same<ReturnType, typename Optimizer::ReturnType>::value>>
967  {
976  static void
977  initialize(typename Optimizer::OptimizerType *optimizer,
978  const SymEngine::vec_basic & independent_symbols,
979  const SymEngine::vec_basic & dependent_functions,
980  const enum OptimizationFlags & optimization_flags)
981  {
982  Assert(optimizer, ExcNotInitialized());
983 
984  const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
985  optimizer->init(independent_symbols,
986  dependent_functions,
987  use_symbolic_cse);
988  }
989 
990 
991 
1005  static void
1006  substitute(typename Optimizer::OptimizerType *optimizer,
1007  std::vector<ReturnType> & output_values,
1008  const std::vector<ReturnType> & substitution_values)
1009  {
1010  Assert(optimizer, ExcNotInitialized());
1011 
1012  // Intermediate values to accommodate the difference in
1013  // value types.
1014  std::vector<typename Optimizer::ReturnType> int_outputs(
1015  output_values.size());
1016  std::vector<typename Optimizer::ReturnType> int_inputs(
1017  substitution_values.size());
1018 
1019  std::copy(substitution_values.begin(),
1020  substitution_values.end(),
1021  int_inputs.begin());
1022  optimizer->call(int_outputs.data(), int_inputs.data());
1023  std::copy(int_outputs.begin(),
1024  int_outputs.end(),
1025  output_values.begin());
1026  }
1027 
1028 
1029 
1034  template <class Archive>
1035  static void
1036  save(Archive & archive,
1037  const unsigned int version,
1038  typename Optimizer::OptimizerType *optimizer)
1039  {
1040  Assert(optimizer, ExcNotInitialized());
1041  Optimizer::save(archive, version, *optimizer);
1042  }
1043 
1044 
1045 
1050  template <class Archive>
1051  static void
1052  load(Archive & archive,
1053  const unsigned int version,
1054  typename Optimizer::OptimizerType *optimizer,
1055  const SymEngine::vec_basic & independent_symbols,
1056  const SymEngine::vec_basic & dependent_functions,
1057  const enum OptimizationFlags & optimization_flags)
1058  {
1059  Assert(optimizer, ExcNotInitialized());
1060 
1061  // Some optimizers don't have the same interface for
1062  // serialization, we filter them out through the specializations
1063  // of the Optimizer class
1064  Optimizer::load(archive,
1065  version,
1066  *optimizer,
1067  independent_symbols,
1068  dependent_functions,
1069  optimization_flags);
1070  }
1071 
1072 
1073 
1089  template <typename Stream>
1090  static void
1091  print(Stream & stream,
1092  typename Optimizer::OptimizerType *optimizer,
1093  const bool print_cse_reductions = true,
1094  const bool print_independent_symbols = false,
1095  const bool print_dependent_functions = false)
1096  {
1097  Assert(optimizer, ExcNotInitialized());
1098 
1099  optimizer->print(stream,
1100  print_independent_symbols,
1101  print_dependent_functions,
1102  print_cse_reductions);
1103  }
1104  };
1105 
1106 # endif // DOXYGEN
1107 
1108 
1109  /* -------------------- Utility functions ---------------------- */
1110 
1111 
1133  template <typename NumberType,
1134  int rank,
1135  int dim,
1136  template <int, int, typename>
1137  class TensorType>
1138  TensorType<rank, dim, NumberType>
1140  const TensorType<rank, dim, Expression> &symbol_tensor,
1141  const std::vector<NumberType> & cached_evaluation,
1142  const BatchOptimizer<NumberType> & optimizer)
1143  {
1144  TensorType<rank, dim, NumberType> out;
1145  for (unsigned int i = 0; i < out.n_independent_components; ++i)
1146  {
1147  const TableIndices<rank> indices(
1148  out.unrolled_to_component_indices(i));
1149  out[indices] =
1150  optimizer.extract(symbol_tensor[indices], cached_evaluation);
1151  }
1152  return out;
1153  }
1154 
1155 
1178  template <typename NumberType, int dim>
1181  const SymmetricTensor<4, dim, Expression> &symbol_tensor,
1182  const std::vector<NumberType> & cached_evaluation,
1183  const BatchOptimizer<NumberType> & optimizer)
1184  {
1186  for (unsigned int i = 0;
1187  i < SymmetricTensor<2, dim>::n_independent_components;
1188  ++i)
1189  for (unsigned int j = 0;
1190  j < SymmetricTensor<2, dim>::n_independent_components;
1191  ++j)
1192  {
1193  const TableIndices<4> indices =
1194  make_rank_4_tensor_indices<dim>(i, j);
1195  out[indices] =
1196  optimizer.extract(symbol_tensor[indices], cached_evaluation);
1197  }
1198  return out;
1199  }
1200 
1201 
1219  template <typename NumberType, typename T>
1220  void
1222  const T & function)
1223  {
1224  optimizer.register_function(function);
1225  }
1226 
1227 
1245  template <typename NumberType, typename T>
1246  void
1248  const std::vector<T> & functions)
1249  {
1250  for (const auto &function : functions)
1251  register_functions(optimizer, function);
1252  }
1253 
1254 
1274  template <typename NumberType, typename T, typename... Args>
1275  void
1277  const T & function,
1278  const Args &...other_functions)
1279  {
1280  register_functions(optimizer, function);
1281  register_functions(optimizer, other_functions...);
1282  }
1283 
1284 
1296  template <int rank,
1297  int dim,
1298  template <int, int, typename>
1299  class TensorType>
1302  const TensorType<rank, dim, Expression> &symbol_tensor)
1303  {
1305  out.reserve(symbol_tensor.n_independent_components);
1306  for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
1307  ++i)
1308  {
1309  const TableIndices<rank> indices(
1310  symbol_tensor.unrolled_to_component_indices(i));
1311  out.push_back(symbol_tensor[indices].get_RCP());
1312  }
1313  return out;
1314  }
1315 
1316 
1326  template <int dim>
1329  const SymmetricTensor<4, dim, Expression> &symbol_tensor)
1330  {
1332  out.reserve(symbol_tensor.n_independent_components);
1333  for (unsigned int i = 0;
1334  i < SymmetricTensor<2, dim>::n_independent_components;
1335  ++i)
1336  for (unsigned int j = 0;
1337  j < SymmetricTensor<2, dim>::n_independent_components;
1338  ++j)
1339  {
1340  const TableIndices<4> indices =
1341  make_rank_4_tensor_indices<dim>(i, j);
1342  out.push_back(symbol_tensor[indices].get_RCP());
1343  }
1344  return out;
1345  }
1346 
1347  } // namespace internal
1348 
1349 
1350 
1441  template <typename ReturnType>
1443  {
1444  public:
1452  BatchOptimizer();
1453 
1474 
1484  BatchOptimizer(const BatchOptimizer &other);
1485 
1490 
1494  ~BatchOptimizer() = default;
1495 
1506  void
1507  copy_from(const BatchOptimizer &other);
1508 
1518  template <typename Stream>
1519  void
1520  print(Stream &stream, const bool print_cse = false) const;
1521 
1530  template <class Archive>
1531  void
1532  save(Archive &archive, const unsigned int version) const;
1533 
1547  template <class Archive>
1548  void
1549  load(Archive &archive, const unsigned int version);
1550 
1551 # ifdef DOXYGEN
1573  template <class Archive>
1574  void
1575  serialize(Archive &archive, const unsigned int version);
1576 # else
1577  // This macro defines the serialize() method that is compatible with
1578  // the templated save() and load() method that have been implemented.
1579  BOOST_SERIALIZATION_SPLIT_MEMBER()
1580 # endif
1581 
1592  void
1594 
1600  void
1601  register_symbols(const SymEngine::map_basic_basic &substitution_map);
1602 
1613  void
1614  register_symbols(const types::symbol_vector &symbols);
1615 
1626  void
1627  register_symbols(const SymEngine::vec_basic &symbols);
1628 
1634  get_independent_symbols() const;
1635 
1641  std::size_t
1642  n_independent_variables() const;
1643 
1655  void
1656  register_function(const Expression &function);
1657 
1662  template <int rank, int dim>
1663  void
1665 
1670  template <int rank, int dim>
1671  void
1673  const SymmetricTensor<rank, dim, Expression> &function_tensor);
1674 
1679  void
1681 
1686  void
1687  register_functions(const SymEngine::vec_basic &functions);
1688 
1698  template <typename T>
1699  void
1700  register_functions(const std::vector<T> &functions);
1701 
1715  template <typename T, typename... Args>
1716  void
1717  register_functions(const T &functions, const Args &...other_functions);
1718 
1723  const types::symbol_vector &
1724  get_dependent_functions() const;
1725 
1732  std::size_t
1733  n_dependent_variables() const;
1734 
1753  void
1757 
1762  enum OptimizerType
1763  optimization_method() const;
1764 
1769  enum OptimizationFlags
1770  optimization_flags() const;
1771 
1777  bool
1778  use_symbolic_CSE() const;
1779 
1795  void
1796  optimize();
1797 
1802  bool
1803  optimized() const;
1804 
1821  void
1823 
1833  void
1834  substitute(const SymEngine::map_basic_basic &substitution_map) const;
1835 
1846  void
1847  substitute(const types::symbol_vector & symbols,
1848  const std::vector<ReturnType> &values) const;
1849 
1860  void
1861  substitute(const SymEngine::vec_basic & symbols,
1862  const std::vector<ReturnType> &values) const;
1863 
1869  bool
1870  values_substituted() const;
1871 
1901  const std::vector<ReturnType> &
1902  evaluate() const;
1903 
1911  ReturnType
1912  evaluate(const Expression &func) const;
1913 
1922  std::vector<ReturnType>
1923  evaluate(const std::vector<Expression> &funcs) const;
1924 
1933  template <int rank, int dim>
1936 
1937 
1946  template <int rank, int dim>
1949 
1950 
1958  ReturnType
1959  extract(const Expression & func,
1960  const std::vector<ReturnType> &cached_evaluation) const;
1961 
1962 
1970  std::vector<ReturnType>
1971  extract(const std::vector<Expression> &funcs,
1972  const std::vector<ReturnType> &cached_evaluation) const;
1973 
1974 
1982  template <int rank, int dim>
1985  const std::vector<ReturnType> & cached_evaluation) const;
1986 
1987 
1995  template <int rank, int dim>
1998  const std::vector<ReturnType> &cached_evaluation) const;
1999 
2002  private:
2006  enum OptimizerType method;
2007 
2012  enum OptimizationFlags flags;
2013 
2024 
2031 
2036  bool
2038  const SD::Expression &function) const;
2039 
2044  bool
2046  const SymEngine::RCP<const SymEngine::Basic> &function) const;
2047 
2062  mutable std::vector<ReturnType> dependent_variables_output;
2063 
2073  std::map<SD::Expression,
2074  std::size_t,
2076 
2082 
2089  mutable std::unique_ptr<SymEngine::Visitor> optimizer;
2090 
2100 
2105  mutable bool has_been_serialized;
2106 
2110  void
2111  register_scalar_function(const SD::Expression &function);
2112 
2117  void
2119 
2123  void
2124  create_optimizer(std::unique_ptr<SymEngine::Visitor> &optimizer);
2125 
2142  void
2143  substitute(const std::vector<ReturnType> &substitution_values) const;
2144  };
2145 
2146 
2147 
2148  /* -------------------- inline and template functions ------------------ */
2149 
2150 
2151 # ifndef DOXYGEN
2152 
2153 
2154  template <typename ReturnType>
2155  template <typename Stream>
2156  void
2157  BatchOptimizer<ReturnType>::print(Stream &stream,
2158  const bool /*print_cse*/) const
2159  {
2160  // Settings
2161  stream << "Method? " << optimization_method() << '\n';
2162  stream << "Flags: " << optimization_flags() << '\n';
2163  stream << "Optimized? " << (optimized() ? "Yes" : "No") << '\n';
2164  stream << "Values substituted? " << values_substituted() << "\n\n";
2165 
2166  // Independent variables
2167  stream << "Symbols (" << n_independent_variables()
2168  << " independent variables):" << '\n';
2169  int cntr = 0;
2170  for (SD::types::substitution_map::const_iterator it =
2171  independent_variables_symbols.begin();
2172  it != independent_variables_symbols.end();
2173  ++it, ++cntr)
2174  {
2175  stream << cntr << ": " << it->first << '\n';
2176  }
2177  stream << '\n' << std::flush;
2178 
2179  // Dependent functions
2180  stream << "Functions (" << n_dependent_variables()
2181  << " dependent variables):" << '\n';
2182  cntr = 0;
2183  for (typename SD::types::symbol_vector::const_iterator it =
2184  dependent_variables_functions.begin();
2185  it != dependent_variables_functions.end();
2186  ++it, ++cntr)
2187  {
2188  stream << cntr << ": " << (*it) << '\n';
2189  }
2190  stream << '\n' << std::flush;
2191 
2192  // Common subexpression
2193  if (optimized() == true && use_symbolic_CSE() == true)
2194  {
2195  Assert(optimizer, ExcNotInitialized());
2196  const bool print_cse_reductions = true;
2197  const bool print_independent_symbols = false;
2198  const bool print_dependent_functions = false;
2199 
2200  if (optimization_method() == OptimizerType::dictionary)
2201  {
2202  Assert(dynamic_cast<typename internal::DictionaryOptimizer<
2203  ReturnType>::OptimizerType *>(optimizer.get()),
2204  ExcMessage("Cannot cast optimizer to Dictionary type."));
2205 
2206  internal::OptimizerHelper<
2207  ReturnType,
2208  internal::DictionaryOptimizer<ReturnType>>::
2209  print(stream,
2210  dynamic_cast<typename internal::DictionaryOptimizer<
2211  ReturnType>::OptimizerType *>(optimizer.get()),
2212  print_independent_symbols,
2213  print_dependent_functions,
2214  print_cse_reductions);
2215 
2216  stream << '\n' << std::flush;
2217  }
2218  else if (optimization_method() == OptimizerType::lambda)
2219  {
2220  Assert(dynamic_cast<typename internal::LambdaOptimizer<
2221  ReturnType>::OptimizerType *>(optimizer.get()),
2222  ExcMessage("Cannot cast optimizer to Lambda type."));
2223 
2224  internal::OptimizerHelper<ReturnType,
2225  internal::LambdaOptimizer<ReturnType>>::
2226  print(stream,
2227  dynamic_cast<typename internal::LambdaOptimizer<
2228  ReturnType>::OptimizerType *>(optimizer.get()),
2229  print_independent_symbols,
2230  print_dependent_functions,
2231  print_cse_reductions);
2232  }
2233 # ifdef HAVE_SYMENGINE_LLVM
2234  else if (optimization_method() == OptimizerType::llvm)
2235  {
2236  Assert(dynamic_cast<typename internal::LLVMOptimizer<
2237  ReturnType>::OptimizerType *>(optimizer.get()),
2238  ExcMessage("Cannot cast optimizer to LLVM type."));
2239 
2240  internal::OptimizerHelper<ReturnType,
2241  internal::LLVMOptimizer<ReturnType>>::
2242  print(stream,
2243  dynamic_cast<typename internal::LLVMOptimizer<
2244  ReturnType>::OptimizerType *>(optimizer.get()),
2245  print_independent_symbols,
2246  print_dependent_functions,
2247  print_cse_reductions);
2248  }
2249 # endif // HAVE_SYMENGINE_LLVM
2250  else
2251  {
2252  AssertThrow(false, ExcMessage("Unknown optimizer type."));
2253  }
2254  }
2255 
2256  if (values_substituted())
2257  {
2258  stream << "Evaluated functions:" << '\n';
2259  stream << std::flush;
2260  cntr = 0;
2261  for (typename std::vector<ReturnType>::const_iterator it =
2262  dependent_variables_output.begin();
2263  it != dependent_variables_output.end();
2264  ++it, ++cntr)
2265  {
2266  stream << cntr << ": " << (*it) << '\n';
2267  }
2268  stream << '\n' << std::flush;
2269  }
2270  }
2271 
2272 
2273 
2274  template <typename ReturnType>
2275  template <class Archive>
2276  void
2277  BatchOptimizer<ReturnType>::save(Archive & ar,
2278  const unsigned int version) const
2279  {
2280  // Serialize enum classes...
2281  {
2282  const auto m =
2283  static_cast<typename std::underlying_type<OptimizerType>::type>(
2284  method);
2285  ar &m;
2286  }
2287  {
2288  const auto f =
2289  static_cast<typename std::underlying_type<OptimizationFlags>::type>(
2290  flags);
2291  ar &f;
2292  }
2293 
2294  // Important: Independent variables must always be
2295  // serialized before the dependent variables.
2296  ar &independent_variables_symbols;
2297  ar &dependent_variables_functions;
2298 
2299  ar &dependent_variables_output;
2300  ar &map_dep_expr_vec_entry;
2301  ar &ready_for_value_extraction;
2302 
2303  // Mark that we've saved this class at some point.
2304  has_been_serialized = true;
2305  ar &has_been_serialized;
2306 
2307  // When we serialize the optimizer itself, we have to (unfortunately)
2308  // provide it with sufficient information to rebuild itself from scratch.
2309  // This is because only two of the three optimization classes support
2310  // real serialization (i.e. have save/load capability).
2311  const SD::types::symbol_vector symbol_vec =
2312  Utilities::extract_symbols(independent_variables_symbols);
2314  *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2315  ReturnType>::OptimizerType *>(optimizer.get()))
2316  {
2317  Assert(optimization_method() == OptimizerType::dictionary,
2318  ExcInternalError());
2319  internal::OptimizerHelper<
2320  ReturnType,
2321  internal::DictionaryOptimizer<ReturnType>>::save(ar, version, opt);
2322  }
2324  *opt = dynamic_cast<typename internal::LambdaOptimizer<
2325  ReturnType>::OptimizerType *>(optimizer.get()))
2326  {
2327  Assert(optimization_method() == OptimizerType::lambda,
2328  ExcInternalError());
2329  internal::OptimizerHelper<
2330  ReturnType,
2331  internal::LambdaOptimizer<ReturnType>>::save(ar, version, opt);
2332  }
2333 # ifdef HAVE_SYMENGINE_LLVM
2335  *opt = dynamic_cast<typename internal::LLVMOptimizer<
2336  ReturnType>::OptimizerType *>(optimizer.get()))
2337  {
2338  Assert(optimization_method() == OptimizerType::llvm,
2339  ExcInternalError());
2340  internal::OptimizerHelper<
2341  ReturnType,
2342  internal::LLVMOptimizer<ReturnType>>::save(ar, version, opt);
2343  }
2344 # endif
2345  else
2346  {
2347  AssertThrow(false, ExcMessage("Unknown optimizer type."));
2348  }
2349  }
2350 
2351 
2352 
2353  template <typename ReturnType>
2354  template <class Archive>
2355  void
2356  BatchOptimizer<ReturnType>::load(Archive &ar, const unsigned int version)
2357  {
2358  Assert(independent_variables_symbols.empty(), ExcInternalError());
2359  Assert(dependent_variables_functions.empty(), ExcInternalError());
2360  Assert(dependent_variables_output.empty(), ExcInternalError());
2361  Assert(map_dep_expr_vec_entry.empty(), ExcInternalError());
2362  Assert(ready_for_value_extraction == false, ExcInternalError());
2363 
2364  // Deserialize enum classes...
2365  {
2366  typename std::underlying_type<OptimizerType>::type m;
2367  ar & m;
2368  method = static_cast<OptimizerType>(m);
2369  }
2370  {
2371  typename std::underlying_type<OptimizationFlags>::type f;
2372  ar & f;
2373  flags = static_cast<OptimizationFlags>(f);
2374  }
2375 
2376  // Important: Independent variables must always be
2377  // deserialized before the dependent variables.
2378  ar &independent_variables_symbols;
2379  ar &dependent_variables_functions;
2380 
2381  ar &dependent_variables_output;
2382  ar &map_dep_expr_vec_entry;
2383  ar &ready_for_value_extraction;
2384 
2385  ar &has_been_serialized;
2386 
2387  // If we're reading in data, then create the optimizer
2388  // and then deserialize it.
2389  Assert(!optimizer, ExcInternalError());
2390 
2391  // Create and configure the optimizer
2392  create_optimizer(optimizer);
2393  Assert(optimizer, ExcNotInitialized());
2394 
2395  // When we deserialize the optimizer itself, we have to (unfortunately)
2396  // provide it with sufficient information to rebuild itself from scratch.
2397  // This is because only two of the three optimization classes support
2398  // real serialization (i.e. have save/load capability).
2399  const SD::types::symbol_vector symbol_vec =
2400  Utilities::extract_symbols(independent_variables_symbols);
2402  *opt = dynamic_cast<typename internal::DictionaryOptimizer<
2403  ReturnType>::OptimizerType *>(optimizer.get()))
2404  {
2405  Assert(optimization_method() == OptimizerType::dictionary,
2406  ExcInternalError());
2407  internal::OptimizerHelper<ReturnType,
2408  internal::DictionaryOptimizer<ReturnType>>::
2409  load(ar,
2410  version,
2411  opt,
2413  symbol_vec),
2415  dependent_variables_functions),
2416  optimization_flags());
2417  }
2419  *opt = dynamic_cast<typename internal::LambdaOptimizer<
2420  ReturnType>::OptimizerType *>(optimizer.get()))
2421  {
2422  Assert(optimization_method() == OptimizerType::lambda,
2423  ExcInternalError());
2424  internal::OptimizerHelper<ReturnType,
2425  internal::LambdaOptimizer<ReturnType>>::
2426  load(ar,
2427  version,
2428  opt,
2430  symbol_vec),
2432  dependent_variables_functions),
2433  optimization_flags());
2434  }
2435 # ifdef HAVE_SYMENGINE_LLVM
2437  *opt = dynamic_cast<typename internal::LLVMOptimizer<
2438  ReturnType>::OptimizerType *>(optimizer.get()))
2439  {
2440  Assert(optimization_method() == OptimizerType::llvm,
2441  ExcInternalError());
2442  internal::OptimizerHelper<ReturnType,
2443  internal::LLVMOptimizer<ReturnType>>::
2444  load(ar,
2445  version,
2446  opt,
2448  symbol_vec),
2450  dependent_variables_functions),
2451  optimization_flags());
2452  }
2453 # endif
2454  else
2455  {
2456  AssertThrow(false, ExcMessage("Unknown optimizer type."));
2457  }
2458  }
2459 
2460 
2461 
2462  template <typename ReturnType>
2463  template <int rank, int dim>
2464  void
2466  const Tensor<rank, dim, Expression> &function_tensor)
2467  {
2468  Assert(optimized() == false,
2469  ExcMessage(
2470  "Cannot register functions once the optimizer is finalised."));
2471 
2472  register_vector_functions(
2473  internal::unroll_to_expression_vector(function_tensor));
2474  }
2475 
2476 
2477 
2478  template <typename ReturnType>
2479  template <int rank, int dim>
2480  void
2482  const SymmetricTensor<rank, dim, Expression> &function_tensor)
2483  {
2484  Assert(optimized() == false,
2485  ExcMessage(
2486  "Cannot register functions once the optimizer is finalised."));
2487 
2488  register_vector_functions(
2489  internal::unroll_to_expression_vector(function_tensor));
2490  }
2491 
2492 
2493 
2494  template <typename ReturnType>
2495  template <typename T, typename... Args>
2496  void
2498  const T &functions,
2499  const Args &...other_functions)
2500  {
2502  internal::register_functions(*this, other_functions...);
2503  }
2504 
2505 
2506 
2507  template <typename ReturnType>
2508  template <typename T>
2509  void
2511  const std::vector<T> &functions)
2512  {
2514  }
2515 
2516 
2517 
2518  template <typename ReturnType>
2519  template <int rank, int dim>
2522  const Tensor<rank, dim, Expression> &funcs,
2523  const std::vector<ReturnType> & cached_evaluation) const
2524  {
2526  cached_evaluation,
2527  *this);
2528  }
2529 
2530 
2531 
2532  template <typename ReturnType>
2533  template <int rank, int dim>
2536  const Tensor<rank, dim, Expression> &funcs) const
2537  {
2538  Assert(
2539  values_substituted() == true,
2540  ExcMessage(
2541  "The optimizer is not configured to perform evaluation. "
2542  "This action can only performed after substitute() has been called."));
2543 
2544  return extract(funcs, dependent_variables_output);
2545  }
2546 
2547 
2548 
2549  template <typename ReturnType>
2550  template <int rank, int dim>
2554  const std::vector<ReturnType> & cached_evaluation) const
2555  {
2557  cached_evaluation,
2558  *this);
2559  }
2560 
2561 
2562 
2563  template <typename ReturnType>
2564  template <int rank, int dim>
2567  const SymmetricTensor<rank, dim, Expression> &funcs) const
2568  {
2569  Assert(
2570  values_substituted() == true,
2571  ExcMessage(
2572  "The optimizer is not configured to perform evaluation. "
2573  "This action can only performed after substitute() has been called."));
2574 
2575  return extract(funcs, dependent_variables_output);
2576  }
2577 
2578 # endif // DOXYGEN
2579 
2580  } // namespace SD
2581 } // namespace Differentiation
2582 
2583 
2585 
2586 #endif // DEAL_II_WITH_SYMENGINE
2587 
2588 #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)
BatchOptimizer(BatchOptimizer &&)=default
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
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:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:512
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcSymEngineLLVMNotAvailable()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcSymEngineLLVMReturnTypeNotSupported()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
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:2394
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)