Discontinuous Galerkin Library
#include "dg/algorithm.h"
Loading...
Searching...
No Matches
multiply.h
Go to the documentation of this file.
1#pragma once
2
3#include "operator.h"
4#include "dg/functors.h"
5#include "dg/blas1.h"
6#include "tensor.h"
7
13namespace dg
14{
16// nvcc does not like local classes so we need to define these globally:
17// \f$ y_i \leftarrow \lambda T_{ij} x_i + \mu y_i\f$
18struct TensorMultiply2d
19{
20 template<class VL, class V0, class V1, class V2, class VM, class V3, class V4>
21 DG_DEVICE
22 void operator() ( VL lambda, V0 t00, V0 t01, V0 t10, V0 t11,
23 V1 in0, V2 in1, VM mu, V3& out0, V4& out1) const
24 {
25 auto tmp0 = DG_FMA(t00,in0 , t01*in1);
26 auto tmp1 = DG_FMA(t10,in0 , t11*in1);
27 auto temp = out1*mu;
28 out1 = DG_FMA( lambda, tmp1, temp);
29 temp = out0*mu;
30 out0 = DG_FMA( lambda, tmp0, temp);
31 }
32};
33// \f$ y_i \leftarrow \lambda T_{ij} x_i + \mu y_i\f$
34struct TensorMultiply3d
35{
36 template<class VL, class V0, class V1, class V2, class V3, class VM, class V4, class V5, class V6>
37 DG_DEVICE
38 void operator() ( VL lambda,
39 V0 t00, V0 t01, V0 t02,
40 V0 t10, V0 t11, V0 t12,
41 V0 t20, V0 t21, V0 t22,
42 V1 in0, V2 in1, V3 in2,
43 VM mu,
44 V4& out0, V5& out1, V6& out2) const
45 {
46 auto tmp0 = DG_FMA( t00,in0 , (DG_FMA( t01,in1 , t02*in2)));
47 auto tmp1 = DG_FMA( t10,in0 , (DG_FMA( t11,in1 , t12*in2)));
48 auto tmp2 = DG_FMA( t20,in0 , (DG_FMA( t21,in1 , t22*in2)));
49 auto temp = out2*mu;
50 out2 = DG_FMA( lambda, tmp2, temp);
51 temp = out1*mu;
52 out1 = DG_FMA( lambda, tmp1, temp);
53 temp = out0*mu;
54 out0 = DG_FMA( lambda, tmp0, temp);
55 }
56};
57// \f$ y_i \leftarrow \lambda T^{-1}_{ij} x_i + \mu y_i\f$
58struct InverseTensorMultiply2d
59{
60 template<class VL, class V0, class V1, class V2, class VM, class V3, class V4>
61 DG_DEVICE
62 void operator() ( VL lambda, V0 t00, V0 t01, V0 t10, V0 t11,
63 V1 in0, V2 in1, VM mu, V3& out0, V4& out1) const
64 {
65 auto dett = DG_FMA( t00,t11 , (-t10*t01));
66 auto tmp0 = DG_FMA( in0,t11 , (-in1*t01));
67 auto tmp1 = DG_FMA( t00,in1 , (-t10*in0));
68 auto temp = out1*mu;
69 out1 = DG_FMA( lambda, tmp1/dett, temp);
70 temp = out0*mu;
71 out0 = DG_FMA( lambda, tmp0/dett, temp);
72 }
73};
74// \f$ y_i \leftarrow \lambda T^{-1}_{ij} x_i + \mu y_i\f$
75struct InverseTensorMultiply3d
76{
77 template<class VL, class V0, class V1, class V2, class V3, class VM, class V4, class V5, class V6>
78 DG_DEVICE
79 void operator() ( VL lambda,
80 V0 t00, V0 t01, V0 t02,
81 V0 t10, V0 t11, V0 t12,
82 V0 t20, V0 t21, V0 t22,
83 V1 in0, V2 in1, V3 in2,
84 VM mu,
85 V4& out0, V5& out1, V6& out2) const
86 {
87 auto dett = t00*DG_FMA(t11, t22, (-t12*t21))
88 -t01*DG_FMA(t10, t22, (-t20*t12))
89 +t02*DG_FMA(t10, t21, (-t20*t11));
90
91 auto tmp0 = in0*DG_FMA(t11, t22, (-t12*t21))
92 -t01*DG_FMA(in1, t22, (-in2*t12))
93 +t02*DG_FMA(in1, t21, (-in2*t11));
94 auto tmp1 = t00*DG_FMA(in1, t22, (-t12*in2))
95 -in0*DG_FMA(t10, t22, (-t20*t12))
96 +t02*DG_FMA(t10, in2, (-t20*in1));
97 auto tmp2 = t00*DG_FMA(t11, in2, (-in1*t21))
98 -t01*DG_FMA(t10, in2, (-t20*in1))
99 +in0*DG_FMA(t10, t21, (-t20*t11));
100 auto temp = out2*mu;
101 out2 = DG_FMA( lambda, tmp2/dett, temp);
102 temp = out1*mu;
103 out1 = DG_FMA( lambda, tmp1/dett, temp);
104 temp = out0*mu;
105 out0 = DG_FMA( lambda, tmp0/dett, temp);
106 }
107};
108//\f$ y = t_{00} t_{11} - t_{10}t_{01} \f$
109struct TensorDeterminant2d
110{
111 template<class value_type>
112 DG_DEVICE
113 value_type operator() ( value_type t00, value_type t01,
114 value_type t10, value_type t11) const
115 {
116 return DG_FMA( t00,t11 , (-t10*t01));
117 }
118};
119//\f$ y = t_{00} t_{11}t_{22} + t_{01}t_{12}t_{20} + t_{02}t_{10}t_{21} - t_{02}t_{11}t_{20} - t_{01}t_{10}t_{22} - t_{00}t_{12}t_{21} \f$
120struct TensorDeterminant3d
121{
122 template<class value_type>
123 DG_DEVICE
124 value_type operator() ( value_type t00, value_type t01, value_type t02,
125 value_type t10, value_type t11, value_type t12,
126 value_type t20, value_type t21, value_type t22) const
127 {
128 return t00* DG_FMA( t11,t22 , (-t21*t12))
129 -t01* DG_FMA( t10,t22 , (-t20*t12))
130 +t02* DG_FMA( t10,t21 , (-t20*t11));
131 }
132};
133
134// \f$ y = \lambda\mu v_i T_{ij} w_j \f$
135struct TensorDot2d
136{
137 template<class VL, class V0, class V1, class V2, class VM, class V3, class V4>
138 DG_DEVICE
139 auto operator() (
140 VL lambda, V0 v0, V1 v1,
141 V2 t00, V2 t01,
142 V2 t10, V2 t11,
143 VM mu, V3 w0, V4 w1
144 ) const
145 {
146 auto tmp0 = DG_FMA(t00,w0 , t01*w1);
147 auto tmp1 = DG_FMA(t10,w0 , t11*w1);
148 return lambda*mu*DG_FMA(v0,tmp0 , v1*tmp1);
149 }
150};
151// \f$ y = \lambda \mu v_i T_{ij} w_j \f$
152struct TensorDot3d
153{
154 template<class VL, class V0, class V1, class V2, class V3, class VM, class V4, class V5, class V6>
155 DG_DEVICE
156 auto operator() (
157 VL lambda,
158 V0 v0, V1 v1, V2 v2,
159 V3 t00, V3 t01, V3 t02,
160 V3 t10, V3 t11, V3 t12,
161 V3 t20, V3 t21, V3 t22,
162 VM mu,
163 V4 w0, V5 w1, V6 w2) const
164 {
165 auto tmp0 = DG_FMA( t00,w0 , (DG_FMA( t01,w1 , t02*w2)));
166 auto tmp1 = DG_FMA( t10,w0 , (DG_FMA( t11,w1 , t12*w2)));
167 auto tmp2 = DG_FMA( t20,w0 , (DG_FMA( t21,w1 , t22*w2)));
168 return lambda*mu*DG_FMA(v0,tmp0 , DG_FMA(v1,tmp1 , v2*tmp2));
169 }
170};
172
178namespace tensor
179{
182
191template<class ContainerType0, class ContainerType1>
192void scal( SparseTensor<ContainerType0>& t, const ContainerType1& mu)
193{
194 unsigned size=t.values().size();
195 for( unsigned i=0; i<size; i++)
196 dg::blas1::pointwiseDot( mu, t.values()[i], t.values()[i]);
197}
198
214template<class ContainerTypeL, class ContainerType0, class ContainerType1, class ContainerType2, class ContainerTypeM, class ContainerType3, class ContainerType4>
215void multiply2d( const ContainerTypeL& lambda, const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, const ContainerTypeM& mu, ContainerType3& out0, ContainerType4& out1)
216{
217 dg::blas1::subroutine( dg::TensorMultiply2d(), lambda,
218 t.value(0,0), t.value(0,1),
219 t.value(1,0), t.value(1,1),
220 in0, in1, mu, out0, out1);
221}
222
239template<class ContainerTypeL, class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerTypeM, class ContainerType4, class ContainerType5, class ContainerType6>
240void multiply3d( const ContainerTypeL& lambda, const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, const ContainerType3& in2, const ContainerTypeM& mu, ContainerType4& out0, ContainerType5& out1, ContainerType6& out2)
241{
242 dg::blas1::subroutine(dg::TensorMultiply3d(),
243 lambda, t.value(0,0), t.value(0,1), t.value(0,2),
244 t.value(1,0), t.value(1,1), t.value(1,2),
245 t.value(2,0), t.value(2,1), t.value(2,2),
246 in0, in1, in2,
247 mu, out0, out1, out2);
248}
249
265template<class ContainerTypeL, class ContainerType0, class ContainerType1, class ContainerType2, class ContainerTypeM, class ContainerType3, class ContainerType4>
266void inv_multiply2d( const ContainerTypeL& lambda, const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, const ContainerTypeM& mu, ContainerType3& out0, ContainerType4& out1)
267{
268 dg::blas1::subroutine( dg::InverseTensorMultiply2d(), lambda,
269 t.value(0,0), t.value(0,1),
270 t.value(1,0), t.value(1,1),
271 in0, in1, mu, out0, out1);
272}
273
291template<class ContainerTypeL, class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerTypeM, class ContainerType4, class ContainerType5, class ContainerType6>
292void inv_multiply3d( const ContainerTypeL& lambda, const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, const ContainerType3& in2, const ContainerTypeM& mu, ContainerType4& out0, ContainerType5& out1, ContainerType6& out2)
293{
294 dg::blas1::subroutine( dg::InverseTensorMultiply3d(),
295 lambda, t.value(0,0), t.value(0,1), t.value(0,2),
296 t.value(1,0), t.value(1,1), t.value(1,2),
297 t.value(2,0), t.value(2,1), t.value(2,2),
298 in0, in1, in2,
299 mu, out0, out1, out2);
300}
301
311template<class ContainerType>
312ContainerType determinant2d( const SparseTensor<ContainerType>& t)
313{
314 ContainerType det = t.value(0,0);
315 dg::blas1::evaluate( det, dg::equals(), dg::TensorDeterminant2d(),
316 t.value(0,0), t.value(0,1),
317 t.value(1,0), t.value(1,1));
318 return det;
319}
320
331template<class ContainerType>
332ContainerType determinant( const SparseTensor<ContainerType>& t)
333{
334 ContainerType det = t.value(0,0);
335 dg::blas1::evaluate( det, dg::equals(), TensorDeterminant3d(),
336 t.value(0,0), t.value(0,1), t.value(0,2),
337 t.value(1,0), t.value(1,1), t.value(1,2),
338 t.value(2,0), t.value(2,1), t.value(2,2));
339 return det;
340}
341
361template<class ContainerType>
362ContainerType volume2d( const SparseTensor<ContainerType>& t)
363{
364 ContainerType vol=determinant2d(t);
365 dg::blas1::transform(vol, vol, dg::InvSqrt<get_value_type<ContainerType>>());
366 return vol;
367}
368
388template<class ContainerType>
389ContainerType volume( const SparseTensor<ContainerType>& t)
390{
391 ContainerType vol=determinant(t);
392 dg::blas1::transform(vol, vol, dg::InvSqrt<get_value_type<ContainerType>>());
393 return vol;
394}
395
396//For convenience
410template<class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerType4>
411void multiply2d( const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, ContainerType3& out0, ContainerType4& out1)
412{
413 multiply2d( 1, t, in0, in1, 0., out0, out1);
414}
415
430template<class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerType4, class ContainerType5, class ContainerType6>
431void multiply3d( const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, const ContainerType3& in2, ContainerType4& out0, ContainerType5& out1, ContainerType6& out2)
432{
433 multiply3d( 1., t, in0, in1, in2, 0., out0, out1, out2);
434}
435
448template<class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerType4>
449void inv_multiply2d( const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, ContainerType3& out0, ContainerType4& out1)
450{
451 inv_multiply2d( 1., t, in0, in1, out0, out1);
452}
453
469template<class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerType4, class ContainerType5, class ContainerType6>
470void inv_multiply3d( const SparseTensor<ContainerType0>& t, const ContainerType1& in0, const ContainerType2& in1, const ContainerType3& in2, ContainerType4& out0, ContainerType5& out1, ContainerType6& out2)
471{
472 inv_multiply3d( 1., t, in0, in1, in2, 0., out0, out1, out2);
473}
474
492template<class ContainerTypeL, class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerTypeM, class ContainerType4, class ContainerType5, class value_type0, class value_type1>
493void scalar_product2d(
494 value_type0 alpha,
495 const ContainerTypeL& lambda,
496 const ContainerType0& v0,
497 const ContainerType1& v1,
498 const SparseTensor<ContainerType2>& t,
499 const ContainerTypeM& mu,
500 const ContainerType3& w0,
501 const ContainerType4& w1,
502 value_type1 beta,
503 ContainerType5& y)
504{
505 dg::blas1::evaluate( y,
506 dg::Axpby( alpha, beta),
507 dg::TensorDot2d(),
508 lambda, v0, v1,
509 t.value(0,0), t.value(0,1),
510 t.value(1,0), t.value(1,1),
511 mu, w0, w1);
512}
513
532template<class ContainerTypeL, class ContainerType0, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerTypeM, class ContainerType4, class ContainerType5, class ContainerType6, class ContainerType7, class value_type0, class value_type1>
533void scalar_product3d(
534 value_type0 alpha,
535 const ContainerTypeL& lambda,
536 const ContainerType0& v0,
537 const ContainerType1& v1,
538 const ContainerType2& v2,
539 const SparseTensor<ContainerType3>& t,
540 const ContainerTypeM& mu,
541 const ContainerType4& w0,
542 const ContainerType5& w1,
543 const ContainerType6& w2,
544 value_type1 beta,
545 ContainerType7& y)
546{
547 dg::blas1::evaluate( y,
548 dg::Axpby( alpha, beta),
549 dg::TensorDot3d(),
550 lambda,
551 v0, v1, v2,
552 t.value(0,0), t.value(0,1), t.value(0,2),
553 t.value(1,0), t.value(1,1), t.value(1,2),
554 t.value(2,0), t.value(2,1), t.value(2,2),
555 mu,
556 w0, w1, w2);
557}
558
560}//namespace tensor
561}//namespace dg
dg::blas1::pointwiseDot
void pointwiseDot(value_type alpha, const ContainerType1 &x1, const ContainerType2 &x2, value_type1 beta, ContainerType &y)
Definition blas1.h:406
dg::blas1::transform
void transform(const ContainerType1 &x, ContainerType &y, UnaryOp op)
Definition blas1.h:585
dg::blas1::subroutine
void subroutine(Subroutine f, ContainerType &&x, ContainerTypes &&... xs)
; Customizable and generic blas1 function
Definition blas1.h:677
dg::blas1::evaluate
void evaluate(ContainerType &y, BinarySubroutine f, Functor g, const ContainerType0 &x0, const ContainerTypes &...xs)
Definition blas1.h:612
dg::coo2d::y
@ y
y direction
dg::get_value_type
typename TensorTraits< std::decay_t< Vector > >::value_type get_value_type
Definition tensor_traits.h:45
dg::tensor::inv_multiply3d
void inv_multiply3d(const ContainerTypeL &lambda, const SparseTensor< ContainerType0 > &t, const ContainerType1 &in0, const ContainerType2 &in1, const ContainerType3 &in2, const ContainerTypeM &mu, ContainerType4 &out0, ContainerType5 &out1, ContainerType6 &out2)
i
Definition multiply.h:292
dg::tensor::scalar_product2d
void scalar_product2d(value_type0 alpha, const ContainerTypeL &lambda, const ContainerType0 &v0, const ContainerType1 &v1, const SparseTensor< ContainerType2 > &t, const ContainerTypeM &mu, const ContainerType3 &w0, const ContainerType4 &w1, value_type1 beta, ContainerType5 &y)
Definition multiply.h:493
dg::tensor::determinant2d
ContainerType determinant2d(const SparseTensor< ContainerType > &t)
Definition multiply.h:312
dg::tensor::inv_multiply2d
void inv_multiply2d(const ContainerTypeL &lambda, const SparseTensor< ContainerType0 > &t, const ContainerType1 &in0, const ContainerType2 &in1, const ContainerTypeM &mu, ContainerType3 &out0, ContainerType4 &out1)
Definition multiply.h:266
dg::tensor::scalar_product3d
void scalar_product3d(value_type0 alpha, const ContainerTypeL &lambda, const ContainerType0 &v0, const ContainerType1 &v1, const ContainerType2 &v2, const SparseTensor< ContainerType3 > &t, const ContainerTypeM &mu, const ContainerType4 &w0, const ContainerType5 &w1, const ContainerType6 &w2, value_type1 beta, ContainerType7 &y)
Definition multiply.h:533
dg::tensor::multiply2d
void multiply2d(const ContainerTypeL &lambda, const SparseTensor< ContainerType0 > &t, const ContainerType1 &in0, const ContainerType2 &in1, const ContainerTypeM &mu, ContainerType3 &out0, ContainerType4 &out1)
Definition multiply.h:215
dg::tensor::determinant
ContainerType determinant(const SparseTensor< ContainerType > &t)
Definition multiply.h:332
dg::tensor::volume
ContainerType volume(const SparseTensor< ContainerType > &t)
Definition multiply.h:389
dg::tensor::multiply3d
void multiply3d(const ContainerTypeL &lambda, const SparseTensor< ContainerType0 > &t, const ContainerType1 &in0, const ContainerType2 &in1, const ContainerType3 &in2, const ContainerTypeM &mu, ContainerType4 &out0, ContainerType5 &out1, ContainerType6 &out2)
Definition multiply.h:240
dg::tensor::scal
void scal(SparseTensor< ContainerType0 > &t, const ContainerType1 &mu)
Definition multiply.h:192
dg::tensor::volume2d
ContainerType volume2d(const SparseTensor< ContainerType > &t)
Definition multiply.h:362
DG_DEVICE
#define DG_DEVICE
Expands to __host__ __device__ if compiled with nvcc else is empty.
Definition dg_doc.h:378
dg
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
dg::InvSqrt
Definition functors.h:83
dg::SparseTensor
Class for 2x2 and 3x3 matrices sharing elements.
Definition tensor.h:51
dg::SparseTensor::values
std::vector< container > & values()
Return write access to the values array.
Definition tensor.h:151
dg::SparseTensor::value
const container & value(size_t i, size_t j) const
Read access the underlying container.
Definition tensor.h:141
dg::equals
Definition subroutines.h:22