MFEM v4.9.0
Finite element discretization library
Loading...
Searching...
No Matches
ttensor.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_TEMPLATE_TENSOR
13#define MFEM_TEMPLATE_TENSOR
14
15#include "../config/tconfig.hpp"
16#include "../linalg/simd.hpp"
17#include "../general/tassign.hpp"
18#include "tlayout.hpp"
19#include "tmatrix.hpp"
20
21// Templated tensor implementation (up to order 4)
22
23namespace mfem
24{
25
26// Element-wise tensor operations
27
28namespace internal
29{
30
31template <int Rank>
32struct TensorOps;
33
34template <>
35struct TensorOps<1> // rank = 1
36{
37 // Assign: A {=,+=,*=} scalar_value
38 template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
39 typename scalar_t>
40 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
41 const scalar_t value)
42 {
43 MFEM_STATIC_ASSERT(A_layout_t::rank == 1, "invalid rank");
44 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
45 {
46 mfem::Assign<Op>(A_data[A_layout.ind(i1)], value);
47 }
48 }
49
50 // Assign: A {=,+=,*=} scalar_value, host+device version
51 template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
52 typename scalar_t>
53 MFEM_HOST_DEVICE
54 static void AssignHD(const A_layout_t &A_layout, A_data_t &A_data,
55 const scalar_t value)
56 {
57 MFEM_STATIC_ASSERT(A_layout_t::rank == 1, "invalid rank");
58 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
59 {
60 mfem::AssignHD<Op>(A_data[A_layout.ind(i1)], value);
61 }
62 }
63
64 // Assign: A {=,+=,*=} B
65 template <AssignOp::Type Op,
66 typename A_layout_t, typename A_data_t,
67 typename B_layout_t, typename B_data_t>
68 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
69 const B_layout_t &B_layout, const B_data_t &B_data)
70 {
71 MFEM_STATIC_ASSERT(A_layout_t::rank == 1 && B_layout_t::rank == 1,
72 "invalid ranks");
73 MFEM_STATIC_ASSERT(A_layout_t::dim_1 == B_layout_t::dim_1,
74 "invalid dimensions");
75 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
76 {
77 mfem::Assign<Op>(A_data[A_layout.ind(i1)], B_data[B_layout.ind(i1)]);
78 }
79 }
80};
81
82template <>
83struct TensorOps<2> // rank = 2
84{
85 // Assign: A {=,+=,*=} scalar_value
86 template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
87 typename scalar_t>
88 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
89 scalar_t value)
90 {
91 MFEM_STATIC_ASSERT(A_layout_t::rank == 2, "invalid rank");
92 for (int i2 = 0; i2 < A_layout_t::dim_2; i2++)
93 {
94 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
95 {
96 mfem::Assign<Op>(A_data[A_layout.ind(i1,i2)], value);
97 }
98 }
99 }
100
101 // Assign: A {=,+=,*=} scalar_value, host+device version
102 template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
103 typename scalar_t>
104 MFEM_HOST_DEVICE
105 static void AssignHD(const A_layout_t &A_layout, A_data_t &A_data,
106 scalar_t value)
107 {
108 MFEM_STATIC_ASSERT(A_layout_t::rank == 2, "invalid rank");
109 for (int i2 = 0; i2 < A_layout_t::dim_2; i2++)
110 {
111 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
112 {
113 mfem::AssignHD<Op>(A_data[A_layout.ind(i1,i2)], value);
114 }
115 }
116 }
117
118 // Assign: A {=,+=,*=} B
119 template <AssignOp::Type Op,
120 typename A_layout_t, typename A_data_t,
121 typename B_layout_t, typename B_data_t>
122 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
123 const B_layout_t &B_layout, const B_data_t &B_data)
124 {
125 MFEM_STATIC_ASSERT(A_layout_t::rank == 2 && B_layout_t::rank == 2,
126 "invalid ranks");
127 MFEM_STATIC_ASSERT(A_layout_t::dim_1 == B_layout_t::dim_1 &&
128 A_layout_t::dim_2 == B_layout_t::dim_2,
129 "invalid dimensions");
130 for (int i2 = 0; i2 < A_layout_t::dim_2; i2++)
131 {
132 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
133 {
134 mfem::Assign<Op>(A_data[A_layout.ind(i1,i2)],
135 B_data[B_layout.ind(i1,i2)]);
136 }
137 }
138 }
139};
140
141template <>
142struct TensorOps<3> // rank = 3
143{
144 // Assign: A {=,+=,*=} scalar_value
145 template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
146 typename scalar_t>
147 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
148 scalar_t value)
149 {
150 MFEM_STATIC_ASSERT(A_layout_t::rank == 3, "invalid rank");
151 for (int i3 = 0; i3 < A_layout_t::dim_3; i3++)
152 {
153 for (int i2 = 0; i2 < A_layout_t::dim_2; i2++)
154 {
155 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
156 {
157 mfem::Assign<Op>(A_data[A_layout.ind(i1,i2,i3)], value);
158 }
159 }
160 }
161 }
162
163 // Assign: A {=,+=,*=} B
164 template <AssignOp::Type Op,
165 typename A_layout_t, typename A_data_t,
166 typename B_layout_t, typename B_data_t>
167 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
168 const B_layout_t &B_layout, const B_data_t &B_data)
169 {
170 MFEM_STATIC_ASSERT(A_layout_t::rank == 3 && B_layout_t::rank == 3,
171 "invalid ranks");
172 MFEM_STATIC_ASSERT(A_layout_t::dim_1 == B_layout_t::dim_1 &&
173 A_layout_t::dim_2 == B_layout_t::dim_2 &&
174 A_layout_t::dim_3 == B_layout_t::dim_3,
175 "invalid dimensions");
176 for (int i3 = 0; i3 < A_layout_t::dim_3; i3++)
177 {
178 for (int i2 = 0; i2 < A_layout_t::dim_2; i2++)
179 {
180 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
181 {
182 mfem::Assign<Op>(A_data[A_layout.ind(i1,i2,i3)],
183 B_data[B_layout.ind(i1,i2,i3)]);
184 }
185 }
186 }
187 }
188};
189
190template <>
191struct TensorOps<4> // rank = 4
192{
193 // Assign: A {=,+=,*=} scalar_value
194 template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
195 typename scalar_t>
196 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
197 scalar_t value)
198 {
199 MFEM_STATIC_ASSERT(A_layout_t::rank == 4, "invalid rank");
200 for (int i4 = 0; i4 < A_layout_t::dim_4; i4++)
201 {
202 for (int i3 = 0; i3 < A_layout_t::dim_3; i3++)
203 {
204 for (int i2 = 0; i2 < A_layout_t::dim_2; i2++)
205 {
206 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
207 {
208 mfem::Assign<Op>(A_data[A_layout.ind(i1,i2,i3,i4)], value);
209 }
210 }
211 }
212 }
213 }
214
215 // Assign: A {=,+=,*=} B
216 template <AssignOp::Type Op,
217 typename A_layout_t, typename A_data_t,
218 typename B_layout_t, typename B_data_t>
219 static void Assign(const A_layout_t &A_layout, A_data_t &A_data,
220 const B_layout_t &B_layout, const B_data_t &B_data)
221 {
222 MFEM_STATIC_ASSERT(A_layout_t::rank == 4 && B_layout_t::rank == 4,
223 "invalid ranks");
224 MFEM_STATIC_ASSERT(A_layout_t::dim_1 == B_layout_t::dim_1 &&
225 A_layout_t::dim_2 == B_layout_t::dim_2 &&
226 A_layout_t::dim_3 == B_layout_t::dim_3 &&
227 A_layout_t::dim_4 == B_layout_t::dim_4,
228 "invalid dimensions");
229 for (int i4 = 0; i4 < A_layout_t::dim_4; i4++)
230 {
231 for (int i3 = 0; i3 < A_layout_t::dim_3; i3++)
232 {
233 for (int i2 = 0; i2 < A_layout_t::dim_2; i2++)
234 {
235 for (int i1 = 0; i1 < A_layout_t::dim_1; i1++)
236 {
237 mfem::Assign<Op>(A_data[A_layout.ind(i1,i2,i3,i4)],
238 B_data[B_layout.ind(i1,i2,i3,i4)]);
239 }
240 }
241 }
242 }
243 }
244};
245
246} // namespace mfem::internal
247
248// Tensor or sub-tensor assign function: A {=,+=,*=} scalar_value.
249template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
250 typename scalar_t>
251inline void TAssign(const A_layout_t &A_layout, A_data_t &A_data,
252 const scalar_t value)
253{
254 internal::TensorOps<A_layout_t::rank>::
255 template Assign<Op>(A_layout, A_data, value);
256}
257
258// Tensor or sub-tensor assign function: A {=,+=,*=} scalar_value.
259// Host+device version.
260template <AssignOp::Type Op, typename A_layout_t, typename A_data_t,
261 typename scalar_t>
262MFEM_HOST_DEVICE
263inline void TAssignHD(const A_layout_t &A_layout, A_data_t &A_data,
264 const scalar_t value)
265{
266 internal::TensorOps<A_layout_t::rank>::
267 template AssignHD<Op>(A_layout, A_data, value);
268}
269
270// Tensor assign function: A {=,+=,*=} B that allows different input and output
271// layouts. With suitable layouts this function can be used to permute
272// (transpose) tensors, extract sub-tensors, etc.
273template <AssignOp::Type Op,
274 typename A_layout_t, typename A_data_t,
275 typename B_layout_t, typename B_data_t>
276inline void TAssign(const A_layout_t &A_layout, A_data_t &A_data,
277 const B_layout_t &B_layout, const B_data_t &B_data)
278{
279 internal::TensorOps<A_layout_t::rank>::
280 template Assign<Op>(A_layout, A_data, B_layout, B_data);
281}
282
283// classes TVector, TMatrix, TTensor3, TTensor4
284
285template <int S, typename data_t = double, bool align = false>
286struct TVector
287{
288public:
289 static const int size = S;
290 static const int aligned_size = align ? MFEM_ALIGN_SIZE(S,data_t) : size;
291 typedef data_t data_type;
292 data_t data[aligned_size>0?aligned_size:1];
293
294 typedef StridedLayout1D<S,1> layout_type;
295 static const layout_type layout;
296
297 data_t &operator[](int i) { return data[i]; }
298 const data_t &operator[](int i) const { return data[i]; }
299
300 template <AssignOp::Type Op>
301 void Assign(const data_t d)
302 {
303 TAssign<Op>(layout, data, d);
304 }
305
306 template <AssignOp::Type Op, typename src_data_t>
307 void Assign(const src_data_t &src)
308 {
309 TAssign<Op>(layout, data, layout, src);
310 }
311
312 template <AssignOp::Type Op, typename dest_data_t>
313 void AssignTo(dest_data_t &dest)
314 {
315 TAssign<Op>(layout, dest, layout, data);
316 }
317
318 void Set(const data_t d)
319 {
320 Assign<AssignOp::Set>(d);
321 }
322
323 template <typename src_data_t>
324 void Set(const src_data_t &src)
325 {
326 Assign<AssignOp::Set>(src);
327 }
328
329 template <typename dest_data_t>
330 void Assemble(dest_data_t &dest) const
331 {
332 AssignTo<AssignOp::Add>(dest);
333 }
334
335 void Scale(const data_t scale)
336 {
337 Assign<AssignOp::Mult>(scale);
338 }
339};
340
341template <int S, typename data_t, bool align>
342const typename TVector<S,data_t,align>::layout_type
343TVector<S,data_t,align>::layout = layout_type();
344
345
346template <int N1, int N2, typename data_t = double, bool align = false>
347struct TMatrix : public TVector<N1*N2,data_t,align>
348{
349 typedef TVector<N1*N2,data_t,align> base_class;
350 using base_class::size;
351 using base_class::data;
352
353 typedef ColumnMajorLayout2D<N1,N2> layout_type;
354 static const layout_type layout;
355 static inline int ind(int i1, int i2) { return layout.ind(i1,i2); }
356
357 data_t &operator()(int i, int j) { return data[ind(i,j)]; }
358 const data_t &operator()(int i, int j) const { return data[ind(i,j)]; }
359
360 inline data_t Det() const
361 {
362 return TDet<data_t>(layout, data);
363 }
364
365 inline void Adjugate(TMatrix<N1,N2,data_t> &adj) const
366 {
367 TAdjugate<data_t>(layout, data, layout, adj.data);
368 }
369
370 // Compute the adjugate and the determinant of a (small) matrix.
371 inline data_t AdjDet(TMatrix<N2,N1,data_t> &adj) const
372 {
373 return TAdjDet<data_t>(layout, data, layout, adj.data);
374 }
375};
376
377template <int N1, int N2, typename data_t, bool align>
378const typename TMatrix<N1,N2,data_t,align>::layout_type
379TMatrix<N1,N2,data_t,align>::layout = layout_type();
380
381
382template <int N1, int N2, int N3, typename data_t = double, bool align = false>
383struct TTensor3 : TVector<N1*N2*N3,data_t,align>
384{
385 typedef TVector<N1*N2*N3,data_t,align> base_class;
386 using base_class::size;
387 using base_class::data;
388
389 typedef ColumnMajorLayout3D<N1,N2,N3> layout_type;
390 static const layout_type layout;
391 static inline int ind(int i1, int i2, int i3)
392 { return layout.ind(i1,i2,i3); }
393
394 data_t &operator()(int i, int j, int k) { return data[ind(i,j,k)]; }
395 const data_t &operator()(int i, int j, int k) const
396 { return data[ind(i,j,k)]; }
397};
398
399template <int N1, int N2, int N3, typename data_t, bool align>
400const typename TTensor3<N1,N2,N3,data_t,align>::layout_type
401TTensor3<N1,N2,N3,data_t,align>::layout = layout_type();
402
403template <int N1, int N2, int N3, int N4, typename data_t = double,
404 bool align = false>
405struct TTensor4 : TVector<N1*N2*N3*N4,data_t,align>
406{
407 typedef TVector<N1*N2*N3*N4,data_t,align> base_class;
408 using base_class::size;
409 using base_class::data;
410
411 typedef ColumnMajorLayout4D<N1,N2,N3,N4> layout_type;
412 static const layout_type layout;
413 static inline int ind(int i1, int i2, int i3, int i4)
414 { return layout.ind(i1,i2,i3,i4); }
415
416 data_t &operator()(int i, int j, int k, int l)
417 { return data[ind(i,j,k,l)]; }
418 const data_t &operator()(int i, int j, int k, int l) const
419 { return data[ind(i,j,k,l)]; }
420};
421
422template <int N1, int N2, int N3, int N4, typename data_t, bool align>
423const typename TTensor4<N1,N2,N3,N4,data_t,align>::layout_type
424TTensor4<N1,N2,N3,N4,data_t,align>::layout = layout_type();
425
426
427// Tensor products
428
429// C_{i,j,k} {=|+=} \sum_s A_{s,j} B_{i,s,k}
430template <bool Add,
431 typename A_layout_t, typename A_data_t,
432 typename B_layout_t, typename B_data_t,
433 typename C_layout_t, typename C_data_t>
434MFEM_ALWAYS_INLINE inline
435void Mult_1_2(const A_layout_t &A_layout, const A_data_t &A_data,
436 const B_layout_t &B_layout, const B_data_t &B_data,
437 const C_layout_t &C_layout, C_data_t &C_data)
438{
439 MFEM_STATIC_ASSERT(A_layout_t::rank == 2 && B_layout_t::rank == 3 &&
440 C_layout_t::rank == 3, "invalid ranks");
441 const int B3 = B_layout_t::dim_3;
442 const int C3 = C_layout_t::dim_3;
443 MFEM_STATIC_ASSERT(B3 == C3, "invalid dimensions");
444 for (int k = 0; k < B3; k++)
445 {
446 Mult_AB<Add>(B_layout.ind3(k), B_data,
447 A_layout, A_data,
448 C_layout.ind3(k), C_data);
449 }
450}
451
452// C_{i,j,k} {=|+=} \sum_s A_{i,s} B_{s,j,k}
453template <bool Add,
454 typename A_layout_t, typename A_data_t,
455 typename B_layout_t, typename B_data_t,
456 typename C_layout_t, typename C_data_t>
457MFEM_ALWAYS_INLINE inline
458void Mult_2_1(const A_layout_t &A_layout, const A_data_t &A_data,
459 const B_layout_t &B_layout, const B_data_t &B_data,
460 const C_layout_t &C_layout, C_data_t &C_data)
461{
462 MFEM_STATIC_ASSERT(A_layout_t::rank == 2 && B_layout_t::rank == 3 &&
463 C_layout_t::rank == 3, "invalid ranks");
464 Mult_AB<Add>(A_layout, A_data,
465 B_layout.merge_23(), B_data,
466 C_layout.merge_23(), C_data);
467}
468
469// C_{i,k,j,l} {=|+=} \sum_s A_{s,i} A_{s,j} B_{k,s,l}
470template <bool Add,
471 typename A_layout_t, typename A_data_t,
472 typename B_layout_t, typename B_data_t,
473 typename C_layout_t, typename C_data_t>
474MFEM_ALWAYS_INLINE inline
475void TensorAssemble(const A_layout_t &A_layout, const A_data_t &A_data,
476 const B_layout_t &B_layout, const B_data_t &B_data,
477 const C_layout_t &C_layout, C_data_t &C_data)
478{
479 MFEM_STATIC_ASSERT(A_layout_t::rank == 2 && B_layout_t::rank == 3 &&
480 C_layout_t::rank == 4, "invalid ranks");
481 const int A1 = A_layout_t::dim_1;
482 const int A2 = A_layout_t::dim_2;
483 const int B1 = B_layout_t::dim_1;
484 const int B2 = B_layout_t::dim_2;
485 const int B3 = B_layout_t::dim_3;
486 const int C1 = C_layout_t::dim_1;
487 const int C2 = C_layout_t::dim_2;
488 const int C3 = C_layout_t::dim_3;
489 const int C4 = C_layout_t::dim_4;
490 MFEM_STATIC_ASSERT(A1 == B2 && A2 == C1 && A2 == C3 && B1 == C2 && B3 == C4,
491 "invalid dimensions");
492
493#if 1
494 // Impl == 3
495 MFEM_FLOPS_ADD(3*A1*A2*A2*B1*B3);
496 if (!Add) { TAssign<AssignOp::Set>(C_layout, C_data, 0.0); }
497 for (int j = 0; j < A2; j++)
498 {
499 for (int i = 0; i < A2; i++)
500 {
501 for (int l = 0; l < B3; l++)
502 {
503 for (int k = 0; k < B1; k++)
504 {
505 for (int s = 0; s < A1; s++)
506 {
507 // C(i,k,j,l) += A(s,i) * A(s,j) * B(k,s,l);
508 C_data[C_layout.ind(i,k,j,l)] +=
509 A_data[A_layout.ind(s,i)] *
510 A_data[A_layout.ind(s,j)] *
511 B_data[B_layout.ind(k,s,l)];
512 }
513 }
514 }
515 }
516 }
517#else
518 // Impl == 1
519 if (!Add) { TAssign<AssignOp::Set>(C_layout, C_data, 0.0); }
520 for (int s = 0; s < A1; s++)
521 {
522 for (int i = 0; i < A2; i++)
523 {
524 for (int k = 0; k < B1; k++)
525 {
526 for (int j = 0; j < A2; j++)
527 {
528 for (int l = 0; l < B3; l++)
529 {
530 // C(i,k,j,l) += A(s,i) * A(s,j) * B(k,s,l);
531 C_data[C_layout.ind(i,k,j,l)] +=
532 A_data[A_layout.ind(s,i)] *
533 A_data[A_layout.ind(s,j)] *
534 B_data[B_layout.ind(k,s,l)];
535 }
536 }
537 }
538 }
539 }
540#endif
541}
542
543// D_{i,k,j,l} {=|+=} \sum_s A_{i,s} B_{s,j} C_{k,s,l}
544template <bool Add,
545 typename A_layout_t, typename A_data_t,
546 typename B_layout_t, typename B_data_t,
547 typename C_layout_t, typename C_data_t,
548 typename D_layout_t, typename D_data_t>
549MFEM_ALWAYS_INLINE inline
550void TensorAssemble(const A_layout_t &A_layout, const A_data_t &A_data,
551 const B_layout_t &B_layout, const B_data_t &B_data,
552 const C_layout_t &C_layout, const C_data_t &C_data,
553 const D_layout_t &D_layout, D_data_t &D_data)
554{
555 MFEM_STATIC_ASSERT(A_layout_t::rank == 2 && B_layout_t::rank == 2 &&
556 C_layout_t::rank == 3 && D_layout_t::rank == 4,
557 "invalid ranks");
558 const int A1 = A_layout_t::dim_1;
559 const int A2 = A_layout_t::dim_2;
560 const int B1 = B_layout_t::dim_1;
561 const int B2 = B_layout_t::dim_2;
562 const int C1 = C_layout_t::dim_1;
563 const int C2 = C_layout_t::dim_2;
564 const int C3 = C_layout_t::dim_3;
565 const int D1 = D_layout_t::dim_1;
566 const int D2 = D_layout_t::dim_2;
567 const int D3 = D_layout_t::dim_3;
568 const int D4 = D_layout_t::dim_4;
569 MFEM_STATIC_ASSERT(A2 == B1 && A2 == C2 && A1 == D1 && B2 == D3 &&
570 C1 == D2 && C3 == D4, "invalid dimensions");
571
572#if 0
573 TTensor4<A1,C1,A2,C3> H;
574 // H_{i,k,s,l} = A_{i,s} C_{k,s,l}
575 for (int l = 0; l < C3; l++)
576 {
577 for (int s = 0; s < B1; s++)
578 {
579 for (int k = 0; k < C1; k++)
580 {
581 for (int i = 0; i < A1; i++)
582 {
583 H(i,k,s,l) = A_data[A_layout.ind(i,s)]*
584 C_data[C_layout.ind(k,s,l)];
585 }
586 }
587 }
588 }
589 // D_{(i,k),j,l} = \sum_s B_{s,j} H_{(i,k),s,l}
590 Mult_1_2<Add>(B_layout, B_data, H.layout.merge_12(), H,
591 D_layout.merge_12(), D_data);
592#elif 1
593 MFEM_FLOPS_ADD(A1*B1*C1*C3); // computation of H(l)
594 for (int l = 0; l < C3; l++)
595 {
596 TTensor3<A1,C1,A2,typename C_data_t::data_type> H;
597 // H(l)_{i,k,s} = A_{i,s} C_{k,s,l}
598 for (int s = 0; s < B1; s++)
599 {
600 for (int k = 0; k < C1; k++)
601 {
602 for (int i = 0; i < A1; i++)
603 {
604 H(i,k,s) = A_data[A_layout.ind(i,s)]*
605 C_data[C_layout.ind(k,s,l)];
606 }
607 }
608 }
609 // D_{(i,k),j,l} = \sum_s H(l)_{(i,k),s} B_{s,j}
610 Mult_AB<Add>(H.layout.merge_12(), H, B_layout, B_data,
611 D_layout.merge_12().ind3(l), D_data);
612 }
613#else
614 TTensor4<B1,C1,B2,C3> F;
615 for (int l = 0; l < C3; l++)
616 {
617 for (int j = 0; j < B2; j++)
618 {
619 for (int k = 0; k < C1; k++)
620 {
621 for (int s = 0; s < B1; s++)
622 {
623 F(s,k,j,l) = B_data[B_layout.ind(s,j)]*
624 C_data[C_layout.ind(k,s,l)];
625 }
626 }
627 }
628 }
629 Mult_AB<Add>(A_layout, A_data, F.layout.merge_34().merge_23(), F,
630 D_layout.merge_34().merge_23(), D_data);
631#endif
632}
633
634
635// C_{i,j,k,l} {=|+=} A_{i,j,k} B_{j,l}
636template <AssignOp::Type Op,
637 typename A_layout_t, typename A_data_t,
638 typename B_layout_t, typename B_data_t,
639 typename C_layout_t, typename C_data_t>
640MFEM_ALWAYS_INLINE inline
641void TensorProduct(const A_layout_t &a, const A_data_t &A,
642 const B_layout_t &b, const B_data_t &B,
643 const C_layout_t &c, C_data_t &C)
644{
645 const int A1 = A_layout_t::dim_1;
646 const int A2 = A_layout_t::dim_2;
647 const int A3 = A_layout_t::dim_3;
648 const int B1 = B_layout_t::dim_1;
649 const int B2 = B_layout_t::dim_2;
650 const int C1 = C_layout_t::dim_1;
651 const int C2 = C_layout_t::dim_2;
652 const int C3 = C_layout_t::dim_3;
653 const int C4 = C_layout_t::dim_4;
654 MFEM_STATIC_ASSERT(A1 == C1 && A2 == B1 && A2 == C2 && A3 == C3 && B2 == C4,
655 "invalid dimensions");
656
657 MFEM_FLOPS_ADD(A1*A2*A3*B2);
658 for (int l = 0; l < B2; l++)
659 {
660 for (int k = 0; k < A3; k++)
661 {
662 for (int j = 0; j < A2; j++)
663 {
664 for (int i = 0; i < A1; i++)
665 {
666 mfem::Assign<Op>(C[c.ind(i,j,k,l)],
667 A[a.ind(i,j,k)]*B[b.ind(j,l)]);
668 }
669 }
670 }
671 }
672}
673
674} // namespace mfem
675
676#endif // MFEM_TEMPLATE_TENSOR
b
real_t b
Definition lissajous.cpp:42
a
real_t a
Definition lissajous.cpp:41
mfem::Mult_1_2
MFEM_ALWAYS_INLINE void Mult_1_2(const A_layout_t &A_layout, const A_data_t &A_data, const B_layout_t &B_layout, const B_data_t &B_data, const C_layout_t &C_layout, C_data_t &C_data)
Definition ttensor.hpp:435
mfem::TAssignHD
MFEM_HOST_DEVICE void TAssignHD(const A_layout_t &A_layout, A_data_t &A_data, const scalar_t value)
Definition ttensor.hpp:263
mfem::Assign
lvalue_t & Assign(lvalue_t &a, const rvalue_t &b)
Definition tassign.hpp:137
mfem::TDet
scalar_t TDet(const layout_t &a, const data_t &A)
Definition tmatrix.hpp:572
mfem::Mult_AB
MFEM_ALWAYS_INLINE void Mult_AB(const A_layout_t &A_layout, const A_data_t &A_data, const B_layout_t &B_layout, const B_data_t &B_data, const C_layout_t &C_layout, C_data_t &C_data)
Definition tmatrix.hpp:229
mfem::Mult_2_1
MFEM_ALWAYS_INLINE void Mult_2_1(const A_layout_t &A_layout, const A_data_t &A_data, const B_layout_t &B_layout, const B_data_t &B_data, const C_layout_t &C_layout, C_data_t &C_data)
Definition ttensor.hpp:458
mfem::TensorProduct
MFEM_ALWAYS_INLINE void TensorProduct(const A_layout_t &a, const A_data_t &A, const B_layout_t &b, const B_data_t &B, const C_layout_t &c, C_data_t &C)
Definition ttensor.hpp:641
mfem::TAdjugate
void TAdjugate(const A_layout_t &a, const A_data_t &A, const B_layout_t &b, B_data_t &B)
Definition tmatrix.hpp:619
mfem::TAssign
void TAssign(const A_layout_t &A_layout, A_data_t &A_data, const scalar_t value)
Definition ttensor.hpp:251
mfem::TAdjDet
scalar_t TAdjDet(const A_layout_t &a, const A_data_t &A, const B_layout_t &b, B_data_t &B)
Definition tmatrix.hpp:648
mfem::TensorAssemble
MFEM_ALWAYS_INLINE void TensorAssemble(const A_layout_t &A_layout, const A_data_t &A_data, const B_layout_t &B_layout, const B_data_t &B_data, const C_layout_t &C_layout, C_data_t &C_data)
Definition ttensor.hpp:475
mfem::Add
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
Definition densemat.cpp:2339
mfem::AssignHD
MFEM_HOST_DEVICE lvalue_t & AssignHD(lvalue_t &a, const rvalue_t &b)
Definition tassign.hpp:144
mfem::StridedLayout2D::ind
static MFEM_HOST_DEVICE int ind(int i1, int i2)
Definition tlayout.hpp:106
mfem::StridedLayout3D::merge_12
static StridedLayout2D< N1 *N2, S1, N3, S3 > merge_12()
Definition tlayout.hpp:267
mfem::StridedLayout3D::ind
static int ind(int i1, int i2, int i3)
Definition tlayout.hpp:248
mfem::StridedLayout4D::ind
static int ind(int i1, int i2, int i3, int i4)
Definition tlayout.hpp:395
mfem::StridedLayout4D::merge_34
static StridedLayout3D< N1, S1, N2, S2, N3 *N4, S3 > merge_34()
Definition tlayout.hpp:419
mfem::StridedLayout4D::merge_12
static StridedLayout3D< N1 *N2, S1, N3, S3, N4, S4 > merge_12()
Definition tlayout.hpp:412
mfem::TMatrix::layout
static const layout_type layout
Definition ttensor.hpp:354
mfem::TMatrix::operator()
data_t & operator()(int i, int j)
Definition ttensor.hpp:357
mfem::TMatrix::data
data_t data[aligned_size >0?aligned_size:1]
Definition ttensor.hpp:292
mfem::TMatrix::operator()
const data_t & operator()(int i, int j) const
Definition ttensor.hpp:358
mfem::TMatrix::ind
static int ind(int i1, int i2)
Definition ttensor.hpp:355
mfem::TMatrix::layout_type
ColumnMajorLayout2D< N1, N2 > layout_type
Definition ttensor.hpp:353
mfem::TMatrix::Det
data_t Det() const
Definition ttensor.hpp:360
mfem::TMatrix::base_class
TVector< N1 *N2, data_t, align > base_class
Definition ttensor.hpp:349
mfem::TMatrix::Adjugate
void Adjugate(TMatrix< N1, N2, data_t > &adj) const
Definition ttensor.hpp:365
mfem::TMatrix::AdjDet
data_t AdjDet(TMatrix< N2, N1, data_t > &adj) const
Definition ttensor.hpp:371
mfem::TTensor3::layout_type
ColumnMajorLayout3D< N1, N2, N3 > layout_type
Definition ttensor.hpp:389
mfem::TTensor3::data
data_t data[aligned_size >0?aligned_size:1]
Definition ttensor.hpp:292
mfem::TTensor3::layout
static const layout_type layout
Definition ttensor.hpp:390
mfem::TTensor3::operator()
data_t & operator()(int i, int j, int k)
Definition ttensor.hpp:394
mfem::TTensor3::operator()
const data_t & operator()(int i, int j, int k) const
Definition ttensor.hpp:395
mfem::TTensor3::base_class
TVector< N1 *N2 *N3, data_t, align > base_class
Definition ttensor.hpp:385
mfem::TTensor3::ind
static int ind(int i1, int i2, int i3)
Definition ttensor.hpp:391
mfem::TTensor4::operator()
const data_t & operator()(int i, int j, int k, int l) const
Definition ttensor.hpp:418
mfem::TTensor4::data
data_t data[aligned_size >0?aligned_size:1]
Definition ttensor.hpp:292
mfem::TTensor4::operator()
data_t & operator()(int i, int j, int k, int l)
Definition ttensor.hpp:416
mfem::TTensor4::base_class
TVector< N1 *N2 *N3 *N4, data_t, align > base_class
Definition ttensor.hpp:407
mfem::TTensor4::layout_type
ColumnMajorLayout4D< N1, N2, N3, N4 > layout_type
Definition ttensor.hpp:411
mfem::TTensor4::ind
static int ind(int i1, int i2, int i3, int i4)
Definition ttensor.hpp:413
mfem::TTensor4::layout
static const layout_type layout
Definition ttensor.hpp:412
mfem::TVector::aligned_size
static const int aligned_size
Definition ttensor.hpp:290
mfem::TVector::Scale
void Scale(const data_t scale)
Definition ttensor.hpp:335
mfem::TVector::layout
static const layout_type layout
Definition ttensor.hpp:295
mfem::TVector::operator[]
const data_t & operator[](int i) const
Definition ttensor.hpp:298
mfem::TVector::data
data_t data[aligned_size >0?aligned_size:1]
Definition ttensor.hpp:292
mfem::TVector::operator[]
data_t & operator[](int i)
Definition ttensor.hpp:297
mfem::TVector::Assign
void Assign(const src_data_t &src)
Definition ttensor.hpp:307
mfem::TVector::AssignTo
void AssignTo(dest_data_t &dest)
Definition ttensor.hpp:313
mfem::TVector::Assemble
void Assemble(dest_data_t &dest) const
Definition ttensor.hpp:330
mfem::TVector::size
static const int size
Definition ttensor.hpp:289
mfem::TVector::Set
void Set(const data_t d)
Definition ttensor.hpp:318
mfem::TVector::layout_type
StridedLayout1D< S, 1 > layout_type
Definition ttensor.hpp:294
mfem::TVector::Assign
void Assign(const data_t d)
Definition ttensor.hpp:301
mfem::TVector::data_type
data_t data_type
Definition ttensor.hpp:291
mfem::TVector::Set
void Set(const src_data_t &src)
Definition ttensor.hpp:324
Generated on Thu Dec 11 2025 11:52:55 for MFEM by doxygen 1.11.0