5 #ifndef __TEMPLATES_HPP__
6 #define __TEMPLATES_HPP__
10 template <
typename T>
using ShardVec = boost::shared_ptr<std::vector<T>>;
34 ublas::shallow_array_adaptor<T>(
n, ptr));
56 template <
typename T1>
60 n,
m, ublas::shallow_array_adaptor<T>(
n *
m, ptr));
71 template <
class KeyExtractor1,
class KeyExtractor2>
struct KeyFromKey {
76 const KeyExtractor2 &key2_ = KeyExtractor2())
88 template <
typename id_type>
struct LtBit {
89 inline bool operator()(
const id_type &valueA,
const id_type &valueB)
const {
90 return valueA.to_ulong() < valueB.to_ulong();
94 template <
typename id_type>
struct EqBit {
95 inline bool operator()(
const id_type &valueA,
const id_type &valueB)
const {
96 return valueA.to_ulong() == valueB.to_ulong();
102 return value.to_ulong();
106 template <
class X>
inline std::string
toString(X x) {
107 std::ostringstream buffer;
113 static inline auto get(ublas::vector<T, A> &data) {
134 template <
int S = 1,
class T,
class A>
139 template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
142 template <
int S,
class T,
class A>
144 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
146 if (data.size1() != 3)
148 "getFTensor1FromMat<3>: wrong size of data matrix, number of "
149 "rows should be 3 but is " +
150 boost::lexical_cast<std::string>(data.size1()));
153 &data(0, 0), &data(1, 0), &data(2, 0));
157 template <
int S,
class T,
class A>
159 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
161 if (data.size1() != 4)
163 "getFTensor1FromMat<4>: wrong size of data matrix, number of "
164 "rows should be 4 but is " +
165 boost::lexical_cast<std::string>(data.size1()));
168 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
172 template <
int S,
class T,
class A>
174 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
176 if (data.size1() != 6)
178 "getFTensor1FromMat<6>: wrong size of data matrix, number of "
179 "rows should be 6 but is " +
180 boost::lexical_cast<std::string>(data.size1()));
183 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
188 template <
int S,
class T,
class A>
190 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
192 if (data.size1() != 9)
194 "getFTensor1FromMat<9>: wrong size of data matrix, number of "
195 "rows should be 9 but is " +
196 boost::lexical_cast<std::string>(data.size1()));
199 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
200 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
204 template <
int S,
class T,
class A>
206 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
208 if (data.size1() != 2)
210 "getFTensor1FromMat<2>: wrong size of data matrix, number of "
211 "rows should be 2 but is " +
212 boost::lexical_cast<std::string>(data.size1()));
219 template <
int S,
class T,
class A>
221 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
223 if (data.size1() != 1)
225 "getFTensor1FromMat<1>: wrong size of data matrix, number of "
226 "rows should be 1 but is " +
227 boost::lexical_cast<std::string>(data.size1()));
236 template <
int Tensor_Dim,
int S = 1,
class T,
class L,
class A>
239 static_assert(!std::is_same<T, T>::value,
"not implemented");
245 template <
int Tensor_Dim,
int S = 1>
251 template <
int Tensor_Dim1,
int Tensor_Dim2,
int S,
class T,
class L,
class A>
253 static inline auto get(ublas::matrix<T, L, A> &data) {
255 if (data.size1() != Tensor_Dim1 * Tensor_Dim2) {
257 "getFTensor2FromMat<" +
258 boost::lexical_cast<std::string>(Tensor_Dim1) +
"," +
259 boost::lexical_cast<std::string>(Tensor_Dim2) +
260 ">: wrong size of rows in data matrix, should be " +
261 boost::lexical_cast<std::string>(Tensor_Dim1 * Tensor_Dim2) +
262 " but is " + boost::lexical_cast<std::string>(data.size1()));
265 std::array<double *, Tensor_Dim1 * Tensor_Dim2> ptrs;
266 for (
auto i = 0;
i != Tensor_Dim1 * Tensor_Dim2; ++
i)
267 ptrs[
i] = &data(
i, 0);
276 template <
int Tensor_Dim1,
int Tensor_Dim2>
283 template <
int Tensor_Dim1,
int Tensor_Dim2>
296 template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
299 template <
int S,
class T,
class L,
class A>
301 static inline auto get(ublas::matrix<T, L, A> &data) {
303 if (data.size1() != 6)
305 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
306 "of rows should be 6 but is " +
307 boost::lexical_cast<std::string>(data.size1()));
310 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
315 template <
int S,
class T,
class L,
class A>
317 static inline auto get(ublas::matrix<T, L, A> &data) {
319 if (data.size1() != 3)
321 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
322 "of rows should be 3 but is " +
323 boost::lexical_cast<std::string>(data.size1()));
326 &data(0, 0), &data(1, 0), &data(2, 0));
333 template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
338 template <
int Tensor_Dim,
int S = 1>
344 template <
int Tensor_Dim01,
int Tensor_Dim23,
int S,
class T,
class L,
class A>
347 template <
int S,
class T,
class A>
349 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
351 if (data.size1() != 1)
353 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
354 "of rows should be 1 but is " +
355 boost::lexical_cast<std::string>(data.size1()));
361 template <
int S,
class T,
class A>
363 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
365 if (data.size1() != 9) {
367 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
368 "of rows should be 9 but is " +
369 boost::lexical_cast<std::string>(data.size1()));
373 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
374 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
378 template <
int S,
class T,
class A>
380 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
382 if (data.size1() != 36) {
383 cerr << data.size1() << endl;
385 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
386 "of rows should be 36 but is " +
387 boost::lexical_cast<std::string>(data.size1()));
391 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
392 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
393 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
394 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
395 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
396 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
397 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
414 template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1,
class T,
class L,
418 static_assert(!std::is_same<T, T>::value,
419 "Such getFTensor4DdgFromMat specialisation is not implemented");
422 template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1>
429 template <
int Tensor_Dim01,
int Tensor_Dim2,
int S,
class T,
class L,
class A>
432 template <
int S,
class T,
class A>
434 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
436 if (data.size1() != 1)
438 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
439 "of rows should be 1 but is " +
440 boost::lexical_cast<std::string>(data.size1()));
446 template <
int S,
class T,
class A>
448 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
450 if (data.size1() != 6) {
452 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
453 "of rows should be 6 but is " +
454 boost::lexical_cast<std::string>(data.size1()));
458 &data(0, 0), &data(1, 0), &data(2, 0),
459 &data(3, 0), &data(4, 0), &data(5, 0)};
463 template <
int S,
class T,
class A>
465 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
467 if (data.size1() != 18) {
468 cerr << data.size1() << endl;
470 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
471 "of rows should be 18 but is " +
472 boost::lexical_cast<std::string>(data.size1()));
476 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
477 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
478 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
479 &data(15, 0), &data(16, 0), &data(17, 0)};
495 template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1,
class T,
class L,
499 static_assert(!std::is_same<T, T>::value,
500 "Such getFTensor3DgFromMat specialisation is not implemented");
503 template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1>
509 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
510 int S,
class T,
class L,
class A>
513 template <
int S,
class T,
class A>
515 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
517 if (data.size1() != 1)
519 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
520 "of rows should be 1 but is " +
521 boost::lexical_cast<std::string>(data.size1()));
528 template <
int S,
class T,
class A>
530 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
532 if (data.size1() != 16) {
534 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
535 "of rows should be 16 but is " +
536 boost::lexical_cast<std::string>(data.size1()));
540 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
541 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
542 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
543 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
547 template <
int S,
class T,
class A>
549 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
551 if (data.size1() != 81) {
552 cerr << data.size1() << endl;
554 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
555 "of rows should be 81 but is " +
556 boost::lexical_cast<std::string>(data.size1()));
560 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
561 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
562 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
563 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
564 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
565 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
566 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
567 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
568 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
569 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
570 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
571 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
572 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
573 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
574 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
575 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
594 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
595 int S = 1,
class T,
class L,
class A>
597 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
599 static_assert(!std::is_same<T, T>::value,
600 "Such getFTensor4FromMat specialisation is not implemented");
603 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
611 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S,
class T,
615 template <
int S,
class T,
class A>
617 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
619 if (data.size1() != 1)
621 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
622 "of rows should be 1 but is " +
623 boost::lexical_cast<std::string>(data.size1()));
630 template <
int S,
class T,
class A>
632 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
634 if (data.size1() != 8)
636 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
637 "of rows should be 8 but is " +
638 boost::lexical_cast<std::string>(data.size1()));
641 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
642 &data(5, 0), &data(6, 0), &data(7, 0)
648 template <
int S,
class T,
class A>
650 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
652 if (data.size1() != 12)
654 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
655 "of rows should be 12 but is " +
656 boost::lexical_cast<std::string>(data.size1()));
659 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
660 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
661 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
665 template <
int S,
class T,
class A>
667 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
669 if (data.size1() != 12)
671 "getFTensor3FromMat<2, 2, 3>: wrong size of data matrix, number "
672 "of rows should be 12 but is " +
673 boost::lexical_cast<std::string>(data.size1()));
676 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
677 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
678 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
682 template <
int S,
class T,
class A>
684 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
686 if (data.size1() != 27)
688 "getFTensor3FromMat<3, 3, 3>: wrong size of data matrix, number "
689 "of rows should be 27 but is " +
690 boost::lexical_cast<std::string>(data.size1()));
693 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
694 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
695 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
696 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
697 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
698 &data(25, 0), &data(26, 0)};
702 template <
int S,
class T,
class A>
704 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
706 if (data.size1() != 54)
708 "getFTensor3FromMat<6, 3, 3>: wrong size of data matrix, number "
709 "of rows should be 54 but is " +
710 boost::lexical_cast<std::string>(data.size1()));
712 std::array<double *, 54> ptrs;
713 for (
auto i = 0;
i != 54; ++
i)
714 ptrs[
i] = &data(
i, 0);
719 template <
int S,
class T,
class A>
721 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
723 if (data.size1() != 54)
725 "getFTensor3FromMat<3, 3, 6>: wrong size of data matrix, number "
726 "of rows should be 54 but is " +
727 boost::lexical_cast<std::string>(data.size1()));
729 std::array<double *, 54> ptrs;
730 for (
auto i = 0;
i != 54; ++
i)
731 ptrs[
i] = &data(
i, 0);
750 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1,
class T,
753 Tensor_Dim1, Tensor_Dim2>
755 static_assert(!std::is_same<T, T>::value,
756 "Such getFTensor3FromMat specialisation is not implemented");
759 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1>
770 inline static auto get(T *ptr) {
777 inline static auto get(T *ptr) {
785 inline static auto get(T *ptr) {
793 inline static auto get(T *ptr) {
801 inline static auto get(T *ptr) {
803 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5);
814 template <
int DIM,
int S = DIM>
821 template <
int DIM,
int S = DIM>
828 template <
int DIM,
int S = DIM>
838 inline static auto get(T *ptr) {
840 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5);
846 inline static auto get(T *ptr) {
848 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
849 ptr + 8, ptr + 9, ptr + 10, ptr + 11, ptr + 12, ptr + 13, ptr + 14,
850 ptr + 15, ptr + 16, ptr + 17, ptr + 18, ptr + 19, ptr + 20, ptr + 21,
851 ptr + 22, ptr + 23, ptr + 24, ptr + 25, ptr + 26, ptr + 27, ptr + 28,
852 ptr + 29, ptr + 30, ptr + 31, ptr + 32, ptr + 33, ptr + 34, ptr + 35);
858 inline static auto get(T *ptr) {
860 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
867 inline static auto get(T *ptr) {
875 inline static auto get(T *ptr) {
883 inline static auto get(T *ptr) {
890 inline static auto get(T *ptr) {
902 template <
int DIM1,
int DIM2,
int S = DIM1 * DIM2>
914 template <
int DIM1,
int DIM2,
int S = DIM1 * DIM2>
926 template <
int DIM1,
int DIM2>
930 "Such getFTensor2HVecFromPtr is not implemented");
963 template <
int DIM1,
int DIM2,
int DIM3>
968 "Such getFTensor3FromPtr is not implemented");
975 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
976 ptr + 8, ptr + 9, ptr + 10, ptr + 11);
983 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
984 ptr + 8, ptr + 9, ptr + 10, ptr + 11, ptr + 12, ptr + 13, ptr + 14,
985 ptr + 15, ptr + 16, ptr + 17, ptr + 18, ptr + 19, ptr + 20, ptr + 21,
986 ptr + 22, ptr + 23, ptr + 24, ptr + 25, ptr + 26);
1000 static_assert(DIM,
"Such getFTensor2SymmetricFromPtr is not implemented");
1007 ptr + 0, ptr + 1, ptr + 2,
1018 &ptr[0], &ptr[1], &ptr[2]);
1034 static_assert(DIM,
"Such getFTensor2SymmetricFromPtr is not implemented");
1041 ptr + 0, ptr + 1, ptr + 2,
1052 ptr + 0, ptr + 1, ptr + 2);
1068 "Such getFTensor2SymmetricLowerFromPtr is not implemented");
1086 ptr + 0, ptr + 1, ptr + 3);
1093 template <
typename V>
static inline auto get(V &data) {
1094 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1101 template <
typename V>
static inline auto get(V &data) {
1102 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1109 template <
typename V>
static inline auto get(V &data) {
1110 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1118 template <
typename V>
static inline auto get(V &data) {
1119 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1121 &data[2], &data[3]};
1127 template <
typename V>
static inline auto get(V &data) {
1128 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1130 &data[0], &data[1], &data[2], &data[3], &data[4], &data[5]};
1136 template <
typename V>
static inline auto get(V &data) {
1137 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1139 &data[0], &data[1], &data[2], &data[3], &data[4],
1140 &data[5], &data[6], &data[7], &data[8]};
1158 template <
int DIM,
int S = 0>
1165 template <
int DIM,
int S>
1180 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1187 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0), &data(rr + 3, 0)};
1194 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0),
1195 &data(rr + 3, 0), &data(rr + 4, 0), &data(rr + 5, 0),
1196 &data(rr + 6, 0), &data(rr + 7, 0), &data(rr + 8, 0)};
1209 template <
int DIM,
int S>
1212 static_assert(DIM != DIM,
"not implemented");
1227 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1241 template <
int DIM1,
int DIM2,
int S,
class T,
class L,
class A>
1244 template <
int DIM1,
int DIM2,
class T,
class L,
class A>
1247 template <
int S,
class T,
class L,
class A>
1250 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1253 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1255 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1259 template <
int S,
class T,
class L,
class A>
1262 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1265 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1266 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1267 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
1271 template <
class T,
class L,
class A>
1274 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1275 const size_t cc,
const int ss = 0) {
1277 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1279 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
1283 template <
class T,
class L,
class A>
1286 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1287 const size_t cc,
const int ss = 0) {
1289 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1290 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1291 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
1296 template <
int DIM1,
int DIM2,
int S>
1303 template <
int DIM1,
int DIM2>
1312 template <
int S,
typename T,
typename L,
typename A>
1315 const size_t rr,
const size_t cc = 0) {
1319 template <
int S,
typename T,
typename L,
typename A>
1322 const size_t rr,
const size_t cc = 0) {
1328 template <
int DIM1,
int DIM2,
int S>
1346 const size_t M = mat.size1();
1347 const size_t N = mat.size2();
1351 "The input matrix for inverse computation is not square %d != %d",
1354 int *ipv =
new int[
N];
1356 double *work =
new double[lwork];
1361 "lapack error info = %d", info);
1365 "lapack error info = %d", info);
1382 const size_t M = mat.size1();
1383 const size_t N = mat.size2();
1385 if (
M == 0 ||
M !=
N)
1387 "The input matrix for inverse computation is not square %d != %d",
1394 int *ipiv =
new int[
M];
1396 &*
f.data().begin(), nrhs);
1399 SETERRQ1(PETSC_COMM_SELF, 1,
"error lapack solve dgesv info = %d", info);
1417 auto mat_copy = mat;
1435 const size_t M = mat.size1();
1436 const size_t N = mat.size2();
1438 if (
M == 0 ||
M !=
N)
1441 "The input matrix for eigen value computation is not square %d != %d",
1443 if (eig.size() !=
M)
1444 eig.resize(
M,
false);
1449 const int size = (
M + 2) *
M;
1451 double *work =
new double[size];
1453 if (
lapack_dsyev(
'V',
'U',
n, &*eigen_vec.data().begin(), lda,
1454 &*eig.data().begin(), work, lwork) > 0)
1456 "The algorithm failed to compute eigenvalues.");
1469 template <
int DIM,
typename T1,
typename T2>
1476 const int lda = DIM;
1477 const int lwork = (DIM + 2) * DIM;
1478 std::array<
double, (DIM + 2) * DIM> work;
1483 "The algorithm failed to compute eigenvalues.");
1496 template <
int DIM,
typename T1,
typename T2,
typename T3>
1502 for (
int ii = 0; ii != DIM; ii++)
1503 for (
int jj = 0; jj != DIM; jj++)
1504 eigen_vec(ii, jj) = mat(ii, jj);
1519 return t(0, 0) *
t(1, 1) *
t(2, 2) +
t(1, 0) *
t(2, 1) *
t(0, 2) +
1520 t(2, 0) *
t(0, 1) *
t(1, 2) -
t(0, 0) *
t(2, 1) *
t(1, 2) -
1521 t(2, 0) *
t(1, 1) *
t(0, 2) -
t(1, 0) *
t(0, 1) *
t(2, 2);
1532 return t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0);
1548 template <
typename T,
int DIM>
1556 template <
typename T,
int DIM>
1565 template <
int Tensor_Dim,
class T,
class L,
class A>
1567 ublas::vector<T, A> &det_data,
1568 ublas::matrix<T, L, A> &inv_jac_data) {
1571 "Specialization for this template not yet implemented");
1577 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
1584 template <
class T1,
class T2>
1595 template <
class T1,
class T2>
1606 template <
class T1,
class T2,
class T3>
1609 const auto inv_det = 1. / det;
1610 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1611 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1612 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1613 inv_t(1, 0) = (
t(1, 2) *
t(2, 0) -
t(1, 0) *
t(2, 2)) * inv_det;
1614 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1615 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1616 inv_t(2, 0) = (
t(1, 0) *
t(2, 1) -
t(1, 1) *
t(2, 0)) * inv_det;
1617 inv_t(2, 1) = (
t(0, 1) *
t(2, 0) -
t(0, 0) *
t(2, 1)) * inv_det;
1618 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1626 template <
class T1,
class T2,
class T3>
1629 const auto inv_det = 1. / det;
1630 inv_t(0, 0) =
t(1, 1) * inv_det;
1631 inv_t(0, 1) = -
t(0, 1) * inv_det;
1632 inv_t(1, 0) = -
t(1, 0) * inv_det;
1633 inv_t(1, 1) =
t(0, 0) * inv_det;
1649 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) / det;
1650 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) / det;
1651 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) / det;
1652 inv_t(1, 0) = (
t(1, 2) *
t(2, 0) -
t(1, 0) *
t(2, 2)) / det;
1653 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) / det;
1654 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) / det;
1655 inv_t(2, 0) = (
t(1, 0) *
t(2, 1) -
t(1, 1) *
t(2, 0)) / det;
1656 inv_t(2, 1) = (
t(0, 1) *
t(2, 0) -
t(0, 0) *
t(2, 1)) / det;
1657 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) / det;
1669 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>,
double,
1674 const auto inv_det = 1. / det;
1675 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1676 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1677 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1678 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1679 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1680 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1692 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>,
adouble,
1697 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) / det;
1698 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) / det;
1699 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) / det;
1700 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) / det;
1701 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) / det;
1702 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) / det;
1715 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>,
double,
1720 const auto inv_det = 1. / det;
1721 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1722 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1723 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1724 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1725 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1726 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1730 template <
typename T1,
typename T2,
typename T3,
int DIM>
1733 template <
typename T1,
typename T2,
typename T3>
1740 template <
typename T1,
typename T2,
typename T3>
1747 template <
typename T1,
typename T2,
typename T3,
int DIM>
1756 template <
typename T1,
typename T2,
typename T3,
int DIM>
1762 DIM>::invert(
t, det, inv_t);
1770 template <
typename Iterator>
1772 return (*it)->getEnt();
1794 template <
typename Extractor,
typename Iterator>
1796 Iterator end_iter) {
1797 moab::Range::iterator hint =
r.begin();
1798 while (begin_iter != end_iter) {
1800 auto bi = Extractor::extract(begin_iter);
1801 Iterator pj = begin_iter;
1802 while (pj != end_iter && (bi +
j) == Extractor::extract(pj)) {
1806 hint =
r.insert(hint, bi, bi + (
j - 1));
1830 template <
typename MI,
typename MO = Modify_change_nothing>
1834 for (
auto it = mi.begin(); it != mi.end(); ++it) {
1835 if (!
const_cast<MI &
>(mi).modify(it, mo))
1837 "Houston we have a problem");
1865 return boost::make_shared<TempMeshset>(moab);
1922 template <
typename I>
1929 if (tester(first)) {
1931 auto second = first;
1934 while (second != s) {
1941 CHKERR inserter(first, second);
1967 template <
typename Dest = void,
typename... Arg>
1969 if constexpr (std::is_same<void, Dest>::value)
1970 return std::array<std::common_type_t<std::decay_t<Arg>...>,
sizeof...(Arg)>{
1971 {std::forward<Arg>(arg)...}};
1973 return std::array<Dest,
sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
1985 #define FTENSOR_INDEX(DIM, I) I##_FTIndex<DIM> I;
1989 #endif //__TEMPLATES_HPP__