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>(
260 Tensor_Dim2) +
">: 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>
284 template <
int Tensor_Dim1,
int Tensor_Dim2>
297 template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
300 template <
int S,
class T,
class L,
class A>
302 static inline auto get(ublas::matrix<T, L, A> &data) {
304 if (data.size1() != 6)
306 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
307 "of rows should be 6 but is " +
308 boost::lexical_cast<std::string>(data.size1()));
311 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
316 template <
int S,
class T,
class L,
class A>
318 static inline auto get(ublas::matrix<T, L, A> &data) {
320 if (data.size1() != 3)
322 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
323 "of rows should be 3 but is " +
324 boost::lexical_cast<std::string>(data.size1()));
327 &data(0, 0), &data(1, 0), &data(2, 0));
334 template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
339 template <
int Tensor_Dim,
int S = 1>
345 template <
int Tensor_Dim01,
int Tensor_Dim23,
int S,
class T,
class L,
class A>
348 template <
int S,
class T,
class A>
350 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
352 if (data.size1() != 1)
354 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
355 "of rows should be 1 but is " +
356 boost::lexical_cast<std::string>(data.size1()));
362 template <
int S,
class T,
class A>
364 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
366 if (data.size1() != 9) {
368 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
369 "of rows should be 9 but is " +
370 boost::lexical_cast<std::string>(data.size1()));
374 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
375 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
379 template <
int S,
class T,
class A>
381 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
383 if (data.size1() != 36) {
384 cerr << data.size1() << endl;
386 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
387 "of rows should be 36 but is " +
388 boost::lexical_cast<std::string>(data.size1()));
392 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
393 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
394 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
395 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
396 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
397 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
398 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
415 template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1,
class T,
class L,
419 static_assert(!std::is_same<T, T>::value,
420 "Such getFTensor4DdgFromMat specialisation is not implemented");
423 template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1>
430 template <
int Tensor_Dim01,
int Tensor_Dim2,
int S,
class T,
class L,
class A>
433 template <
int S,
class T,
class A>
435 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
437 if (data.size1() != 1)
439 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
440 "of rows should be 1 but is " +
441 boost::lexical_cast<std::string>(data.size1()));
447 template <
int S,
class T,
class A>
449 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
451 if (data.size1() != 6) {
453 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
454 "of rows should be 6 but is " +
455 boost::lexical_cast<std::string>(data.size1()));
459 &data(0, 0), &data(1, 0), &data(2, 0),
460 &data(3, 0), &data(4, 0), &data(5, 0)};
464 template <
int S,
class T,
class A>
466 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
468 if (data.size1() != 18) {
469 cerr << data.size1() << endl;
471 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
472 "of rows should be 18 but is " +
473 boost::lexical_cast<std::string>(data.size1()));
477 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
478 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
479 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
480 &data(15, 0), &data(16, 0), &data(17, 0)};
496 template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1,
class T,
class L,
500 static_assert(!std::is_same<T, T>::value,
501 "Such getFTensor3DgFromMat specialisation is not implemented");
504 template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1>
510 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
511 int S,
class T,
class L,
class A>
514 template <
int S,
class T,
class A>
516 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
518 if (data.size1() != 1)
520 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
521 "of rows should be 1 but is " +
522 boost::lexical_cast<std::string>(data.size1()));
529 template <
int S,
class T,
class A>
531 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
533 if (data.size1() != 16) {
535 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
536 "of rows should be 16 but is " +
537 boost::lexical_cast<std::string>(data.size1()));
541 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
542 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
543 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
544 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
548 template <
int S,
class T,
class A>
550 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
552 if (data.size1() != 81) {
553 cerr << data.size1() << endl;
555 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
556 "of rows should be 81 but is " +
557 boost::lexical_cast<std::string>(data.size1()));
561 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
562 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
563 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
564 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
565 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
566 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
567 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
568 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
569 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
570 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
571 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
572 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
573 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
574 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
575 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
576 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
595 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
596 int S = 1,
class T,
class L,
class A>
598 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
600 static_assert(!std::is_same<T, T>::value,
601 "Such getFTensor4FromMat specialisation is not implemented");
604 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
612 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S,
class T,
616 template <
int S,
class T,
class A>
618 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
620 if (data.size1() != 1)
622 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
623 "of rows should be 1 but is " +
624 boost::lexical_cast<std::string>(data.size1()));
631 template <
int S,
class T,
class A>
633 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
635 if (data.size1() != 8)
637 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
638 "of rows should be 8 but is " +
639 boost::lexical_cast<std::string>(data.size1()));
642 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
643 &data(5, 0), &data(6, 0), &data(7, 0)
649 template <
int S,
class T,
class A>
651 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
653 if (data.size1() != 12)
655 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
656 "of rows should be 12 but is " +
657 boost::lexical_cast<std::string>(data.size1()));
660 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
661 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
662 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
666 template <
int S,
class T,
class A>
668 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
670 if (data.size1() != 12)
672 "getFTensor3FromMat<2, 2, 3>: wrong size of data matrix, number "
673 "of rows should be 12 but is " +
674 boost::lexical_cast<std::string>(data.size1()));
677 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
678 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
679 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
683 template <
int S,
class T,
class A>
685 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
687 if (data.size1() != 27)
689 "getFTensor3FromMat<3, 3, 3>: wrong size of data matrix, number "
690 "of rows should be 27 but is " +
691 boost::lexical_cast<std::string>(data.size1()));
694 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
695 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
696 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
697 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
698 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
699 &data(25, 0), &data(26, 0)};
703 template <
int S,
class T,
class A>
705 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
707 if (data.size1() != 54)
709 "getFTensor3FromMat<6, 3, 3>: wrong size of data matrix, number "
710 "of rows should be 54 but is " +
711 boost::lexical_cast<std::string>(data.size1()));
713 std::array<double *, 54> ptrs;
714 for (
auto i = 0;
i != 54; ++
i)
715 ptrs[
i] = &data(
i, 0);
720 template <
int S,
class T,
class A>
722 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
724 if (data.size1() != 54)
726 "getFTensor3FromMat<3, 3, 6>: wrong size of data matrix, number "
727 "of rows should be 54 but is " +
728 boost::lexical_cast<std::string>(data.size1()));
730 std::array<double *, 54> ptrs;
731 for (
auto i = 0;
i != 54; ++
i)
732 ptrs[
i] = &data(
i, 0);
751 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1,
class T,
754 Tensor_Dim1, Tensor_Dim2>
756 static_assert(!std::is_same<T, T>::value,
757 "Such getFTensor3FromMat specialisation is not implemented");
760 template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1>
771 inline static auto get(T *ptr) {
779 inline static auto get(T *ptr) {
787 inline static auto get(T *ptr) {
789 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
800 template <
int DIM,
int S = DIM>
807 template <
int DIM,
int S = DIM>
814 template <
int DIM,
int S = DIM>
824 inline static auto get(T *ptr) {
826 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
832 inline static auto get(T *ptr) {
834 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
835 &ptr[8], &ptr[9], &ptr[10], &ptr[11], &ptr[12], &ptr[13], &ptr[14],
836 &ptr[15], &ptr[16], &ptr[17], &ptr[18], &ptr[19], &ptr[20], &ptr[21],
837 &ptr[22], &ptr[23], &ptr[24], &ptr[25], &ptr[26], &ptr[27], &ptr[28],
838 &ptr[29], &ptr[30], &ptr[31], &ptr[32], &ptr[33], &ptr[34], &ptr[35]);
844 inline static auto get(T *ptr) {
846 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
853 inline static auto get(T *ptr) {
861 inline static auto get(T *ptr) {
869 inline static auto get(T *ptr) {
876 inline static auto get(T *ptr) {
888 template <
int DIM1,
int DIM2,
int S = DIM1 * DIM2>
900 template <
int DIM1,
int DIM2,
int S = DIM1 * DIM2>
912 template <
int DIM1,
int DIM2>
916 "Such getFTensor2HVecFromPtr is not implemented");
949 template <
int DIM1,
int DIM2,
int DIM3>
954 "Such getFTensor3FromPtr is not implemented");
961 &ptr[0], &ptr[1], &ptr[2],
963 &ptr[3], &ptr[4], &ptr[5],
965 &ptr[6], &ptr[7], &ptr[8],
967 &ptr[9], &ptr[10], &ptr[11]
976 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
977 &ptr[8], &ptr[9], &ptr[10], &ptr[11], &ptr[12], &ptr[13], &ptr[14],
978 &ptr[15], &ptr[16], &ptr[17], &ptr[18], &ptr[19], &ptr[20], &ptr[21],
979 &ptr[22], &ptr[23], &ptr[24], &ptr[25], &ptr[26]);
993 static_assert(DIM,
"Such getFTensor2SymmetricFromPtr is not implemented");
1000 &ptr[0], &ptr[1], &ptr[2],
1011 &ptr[0], &ptr[1], &ptr[2]);
1027 static_assert(DIM,
"Such getFTensor2SymmetricFromPtr is not implemented");
1034 &ptr[0], &ptr[1], &ptr[2],
1045 &ptr[0], &ptr[1], &ptr[2]);
1061 "Such getFTensor2SymmetricLowerFromPtr is not implemented");
1079 &ptr[0], &ptr[1], &ptr[3]);
1086 template <
typename V>
static inline auto get(V &data) {
1087 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1094 template <
typename V>
static inline auto get(V &data) {
1095 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1102 template <
typename V>
static inline auto get(V &data) {
1103 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1111 template <
typename V>
static inline auto get(V &data) {
1112 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1114 &data[2], &data[3]};
1120 template <
typename V>
static inline auto get(V &data) {
1121 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1123 &data[0], &data[1], &data[2], &data[3], &data[4], &data[5]};
1129 template <
typename V>
static inline auto get(V &data) {
1130 using T =
typename std::remove_reference<decltype(data[0])>
::type;
1132 &data[0], &data[1], &data[2], &data[3], &data[4],
1133 &data[5], &data[6], &data[7], &data[8]};
1151 template <
int DIM,
int S = 0>
1158 template <
int DIM,
int S>
1173 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1180 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0), &data(rr + 3, 0)};
1187 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0),
1188 &data(rr + 3, 0), &data(rr + 4, 0), &data(rr + 5, 0),
1189 &data(rr + 6, 0), &data(rr + 7, 0), &data(rr + 8, 0)};
1202 template <
int DIM,
int S>
1205 static_assert(DIM != DIM,
"not implemented");
1220 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1234 template <
int DIM1,
int DIM2,
int S,
class T,
class L,
class A>
1237 template <
int DIM1,
int DIM2,
class T,
class L,
class A>
1240 template <
int S,
class T,
class L,
class A>
1243 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1246 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1248 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1252 template <
int S,
class T,
class L,
class A>
1255 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1258 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1259 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1260 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
1264 template <
class T,
class L,
class A>
1267 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1268 const size_t cc,
const int ss = 0) {
1270 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1272 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
1276 template <
class T,
class L,
class A>
1279 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1280 const size_t cc,
const int ss = 0) {
1282 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1283 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1284 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
1289 template <
int DIM1,
int DIM2,
int S>
1296 template <
int DIM1,
int DIM2>
1305 template <
int S,
typename T,
typename L,
typename A>
1308 const size_t rr,
const size_t cc = 0) {
1312 template <
int S,
typename T,
typename L,
typename A>
1315 const size_t rr,
const size_t cc = 0) {
1321 template <
int DIM1,
int DIM2,
int S>
1339 const size_t M = mat.size1();
1340 const size_t N = mat.size2();
1344 "The input matrix for inverse computation is not square %d != %d",
1347 int *ipv =
new int[
N];
1349 double *work =
new double[lwork];
1354 "lapack error info = %d", info);
1358 "lapack error info = %d", info);
1375 const size_t M = mat.size1();
1376 const size_t N = mat.size2();
1378 if (
M == 0 ||
M !=
N)
1380 "The input matrix for inverse computation is not square %d != %d",
1387 int *ipiv =
new int[
M];
1389 &*
f.data().begin(), nrhs);
1392 SETERRQ1(PETSC_COMM_SELF, 1,
"error lapack solve dgesv info = %d", info);
1410 auto mat_copy = mat;
1428 const size_t M = mat.size1();
1429 const size_t N = mat.size2();
1431 if (
M == 0 ||
M !=
N)
1434 "The input matrix for eigen value computation is not square %d != %d",
1436 if (eig.size() !=
M)
1437 eig.resize(
M,
false);
1442 const int size = (
M + 2) *
M;
1444 double *work =
new double[size];
1446 if (
lapack_dsyev(
'V',
'U',
n, &*eigen_vec.data().begin(), lda,
1447 &*eig.data().begin(), work, lwork) > 0)
1449 "The algorithm failed to compute eigenvalues.");
1462 template <
int DIM,
typename T1,
typename T2>
1469 const int lda = DIM;
1470 const int lwork = (DIM + 2) * DIM;
1471 std::array<
double, (DIM + 2) * DIM> work;
1476 "The algorithm failed to compute eigenvalues.");
1489 template <
int DIM,
typename T1,
typename T2,
typename T3>
1495 for (
int ii = 0; ii != DIM; ii++)
1496 for (
int jj = 0; jj != DIM; jj++)
1497 eigen_vec(ii, jj) = mat(ii, jj);
1512 return t(0, 0) *
t(1, 1) *
t(2, 2) +
t(1, 0) *
t(2, 1) *
t(0, 2) +
1513 t(2, 0) *
t(0, 1) *
t(1, 2) -
t(0, 0) *
t(2, 1) *
t(1, 2) -
1514 t(2, 0) *
t(1, 1) *
t(0, 2) -
t(1, 0) *
t(0, 1) *
t(2, 2);
1525 return t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0);
1541 template <
typename T,
int DIM>
1549 template <
typename T,
int DIM>
1558 template <
int Tensor_Dim,
class T,
class L,
class A>
1560 ublas::vector<T, A> &det_data,
1561 ublas::matrix<T, L, A> &inv_jac_data) {
1564 "Specialization for this template not yet implemented");
1570 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
1577 template <
class T1,
class T2>
1588 template <
class T1,
class T2>
1599 template <
class T1,
class T2,
class T3>
1602 const auto inv_det = 1. / det;
1603 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1604 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1605 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1606 inv_t(1, 0) = (
t(1, 2) *
t(2, 0) -
t(1, 0) *
t(2, 2)) * inv_det;
1607 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1608 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1609 inv_t(2, 0) = (
t(1, 0) *
t(2, 1) -
t(1, 1) *
t(2, 0)) * inv_det;
1610 inv_t(2, 1) = (
t(0, 1) *
t(2, 0) -
t(0, 0) *
t(2, 1)) * inv_det;
1611 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1619 template <
class T1,
class T2,
class T3>
1622 const auto inv_det = 1. / det;
1623 inv_t(0, 0) =
t(1, 1) * inv_det;
1624 inv_t(0, 1) = -
t(0, 1) * inv_det;
1625 inv_t(1, 0) = -
t(1, 0) * inv_det;
1626 inv_t(1, 1) =
t(0, 0) * inv_det;
1642 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) / det;
1643 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) / det;
1644 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) / det;
1645 inv_t(1, 0) = (
t(1, 2) *
t(2, 0) -
t(1, 0) *
t(2, 2)) / det;
1646 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) / det;
1647 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) / det;
1648 inv_t(2, 0) = (
t(1, 0) *
t(2, 1) -
t(1, 1) *
t(2, 0)) / det;
1649 inv_t(2, 1) = (
t(0, 1) *
t(2, 0) -
t(0, 0) *
t(2, 1)) / det;
1650 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) / det;
1662 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>,
double,
1667 const auto inv_det = 1. / det;
1668 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1669 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1670 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1671 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1672 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1673 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1685 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>,
adouble,
1690 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) / det;
1691 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) / det;
1692 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) / det;
1693 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) / det;
1694 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) / det;
1695 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) / det;
1708 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>,
double,
1713 const auto inv_det = 1. / det;
1714 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1715 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1716 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1717 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1718 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1719 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1723 template <
typename T1,
typename T2,
typename T3,
int DIM>
1726 template <
typename T1,
typename T2,
typename T3>
1733 template <
typename T1,
typename T2,
typename T3>
1740 template <
typename T1,
typename T2,
typename T3,
int DIM>
1749 template <
typename T1,
typename T2,
typename T3,
int DIM>
1755 DIM>::invert(
t, det, inv_t);
1763 template <
typename Iterator>
1765 return (*it)->getEnt();
1787 template <
typename Extractor,
typename Iterator>
1789 Iterator end_iter) {
1790 moab::Range::iterator hint =
r.begin();
1791 while (begin_iter != end_iter) {
1793 auto bi = Extractor::extract(begin_iter);
1794 Iterator pj = begin_iter;
1795 while (pj != end_iter && (bi +
j) == Extractor::extract(pj)) {
1799 hint =
r.insert(hint, bi, bi + (
j - 1));
1823 template <
typename MI,
typename MO = Modify_change_nothing>
1827 for (
auto it = mi.begin(); it != mi.end(); ++it) {
1828 if (!
const_cast<MI &
>(mi).modify(it, mo))
1830 "Houston we have a problem");
1858 return boost::make_shared<TempMeshset>(moab);
1915 template <
typename I>
1922 if (tester(first)) {
1924 auto second = first;
1927 while (second != s) {
1934 CHKERR inserter(first, second);
1960 template <
typename Dest = void,
typename... Arg>
1962 if constexpr (std::is_same<void, Dest>::value)
1963 return std::array<std::common_type_t<std::decay_t<Arg>...>,
sizeof...(Arg)>{
1964 {std::forward<Arg>(arg)...}};
1966 return std::array<Dest,
sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
1971 #endif //__TEMPLATES_HPP__