Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81190935
test_tensor_algebra.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Sep 4, 16:20
Size
8 KB
Mime Type
text/x-c
Expires
Fri, Sep 6, 16:20 (2 d)
Engine
blob
Format
Raw Data
Handle
20528798
Attached To
rMUSPECTRE µSpectre
test_tensor_algebra.cc
View Options
/**
* @file test_tensor_algebra.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 05 Nov 2017
*
* @brief Tests for the tensor algebra functions
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <iomanip>
#include <unsupported/Eigen/CXX11/Tensor>
#include "common/tensor_algebra.hh"
#include "tests.hh"
#include "tests/test_goodies.hh"
namespace
muSpectre
{
BOOST_AUTO_TEST_SUITE
(
tensor_algebra
)
auto
TerrNorm
=
[](
auto
&&
t
){
return
Eigen
::
Tensor
<
Real
,
0
>
(
t
.
abs
().
sum
())();
};
/* ---------------------------------------------------------------------- */
BOOST_AUTO_TEST_CASE
(
tensor_outer_product_test
)
{
constexpr
Dim_t
dim
{
2
};
Eigen
::
TensorFixedSize
<
Real
,
Eigen
::
Sizes
<
dim
,
dim
>>
A
,
B
;
// use prime numbers so that every multiple is uniquely identifiable
A
.
setValues
({{
1
,
2
},
{
3
,
7
}});
B
.
setValues
({{
11
,
13
},
{
17
,
19
}});
Eigen
::
TensorFixedSize
<
Real
,
Eigen
::
Sizes
<
dim
,
dim
,
dim
,
dim
>>
Res1
,
Res2
,
Res3
;
for
(
Dim_t
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
Dim_t
j
=
0
;
j
<
dim
;
++
j
)
{
for
(
Dim_t
k
=
0
;
k
<
dim
;
++
k
)
{
for
(
Dim_t
l
=
0
;
l
<
dim
;
++
l
)
{
Res1
(
i
,
j
,
k
,
l
)
=
A
(
i
,
j
)
*
B
(
k
,
l
);
Res2
(
i
,
j
,
k
,
l
)
=
A
(
i
,
k
)
*
B
(
j
,
l
);
Res3
(
i
,
j
,
k
,
l
)
=
A
(
i
,
l
)
*
B
(
j
,
k
);
}
}
}
}
Real
error
=
TerrNorm
(
Res1
-
Tensors
::
outer
<
dim
>
(
A
,
B
));
BOOST_CHECK_LT
(
error
,
tol
);
error
=
TerrNorm
(
Res2
-
Tensors
::
outer_under
<
dim
>
(
A
,
B
));
BOOST_CHECK_LT
(
error
,
tol
);
error
=
TerrNorm
(
Res3
-
Tensors
::
outer_over
<
dim
>
(
A
,
B
));
if
(
error
>
tol
)
{
std
::
cout
<<
"reference:"
<<
std
::
endl
<<
Res3
<<
std
::
endl
;
std
::
cout
<<
"result:"
<<
std
::
endl
<<
Tensors
::
outer_over
<
dim
>
(
A
,
B
)
<<
std
::
endl
;
std
::
cout
<<
"A:"
<<
std
::
endl
<<
A
<<
std
::
endl
;
std
::
cout
<<
"B"
<<
std
::
endl
<<
B
<<
std
::
endl
;
decltype
(
Res3
)
tmp
=
Tensors
::
outer_over
<
dim
>
(
A
,
B
);
for
(
Dim_t
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
Dim_t
j
=
0
;
j
<
dim
;
++
j
)
{
for
(
Dim_t
k
=
0
;
k
<
dim
;
++
k
)
{
for
(
Dim_t
l
=
0
;
l
<
dim
;
++
l
)
{
std
::
cout
<<
"for ("
<<
i
<<
", "
<<
j
<<
", "
<<
k
<<
", "
<<
l
<<
"), ref: "
<<
std
::
setw
(
3
)
<<
Res3
(
i
,
j
,
k
,
l
)
<<
", res: "
<<
std
::
setw
(
3
)
<<
tmp
(
i
,
j
,
k
,
l
)
<<
std
::
endl
;
}
}
}
}
}
BOOST_CHECK_LT
(
error
,
tol
);
error
=
TerrNorm
(
Res3
-
Tensors
::
outer
<
dim
>
(
A
,
B
));
BOOST_CHECK_GT
(
error
,
tol
);
};
BOOST_FIXTURE_TEST_CASE_TEMPLATE
(
outer_products
,
Fix
,
testGoodies
::
dimlist
,
Fix
)
{
constexpr
auto
dim
{
Fix
::
dim
};
using
T2
=
Tensors
::
Tens2_t
<
dim
>
;
using
M2
=
Matrices
::
Tens2_t
<
dim
>
;
using
Map2
=
Eigen
::
Map
<
const
M2
>
;
using
T4
=
Tensors
::
Tens4_t
<
dim
>
;
using
M4
=
Matrices
::
Tens4_t
<
dim
>
;
using
Map4
=
Eigen
::
Map
<
const
M4
>
;
T2
A
,
B
;
T4
RT
;
A
.
setRandom
();
B
.
setRandom
();
Map2
Amap
(
A
.
data
());
Map2
Bmap
(
B
.
data
());
M2
C
,
D
;
M4
RM
;
C
=
Amap
;
D
=
Bmap
;
auto
error
=
[](
const
T4
&
A
,
const
M4
&
B
)
{
return
(
B
-
Map4
(
A
.
data
())).
norm
();};
// Check outer product
RT
=
Tensors
::
outer
<
dim
>
(
A
,
B
);
RM
=
Matrices
::
outer
(
C
,
D
);
BOOST_CHECK_LT
(
error
(
RT
,
RM
),
tol
);
// Check outer_under product
RT
=
Tensors
::
outer_under
<
dim
>
(
A
,
B
);
RM
=
Matrices
::
outer_under
(
C
,
D
);
BOOST_CHECK_LT
(
error
(
RT
,
RM
),
tol
);
// Check outer_over product
RT
=
Tensors
::
outer_over
<
dim
>
(
A
,
B
);
RM
=
Matrices
::
outer_over
(
C
,
D
);
BOOST_CHECK_LT
(
error
(
RT
,
RM
),
tol
);
}
BOOST_AUTO_TEST_CASE
(
tensor_multiplication
)
{
constexpr
Dim_t
dim
{
2
};
using
Strain_t
=
Eigen
::
TensorFixedSize
<
Real
,
Eigen
::
Sizes
<
dim
,
dim
>>
;
Strain_t
A
,
B
;
A
.
setValues
({{
1
,
2
},
{
3
,
7
}});
B
.
setValues
({{
11
,
13
},
{
17
,
19
}});
Strain_t
FF1
=
A
*
B
;
// element-wise multiplication
std
::
array
<
Eigen
::
IndexPair
<
int
>
,
1
>
prod_dims
{
Eigen
::
IndexPair
<
int
>
{
1
,
0
}};
Strain_t
FF2
=
A
.
contract
(
B
,
prod_dims
);
// correct option 1
Strain_t
FF3
;
using
Mat_t
=
Eigen
::
Map
<
Eigen
::
Matrix
<
Real
,
dim
,
dim
>>
;
// following only works for evaluated tensors (which already have data())
Mat_t
(
FF3
.
data
())
=
Mat_t
(
A
.
data
())
*
Mat_t
(
B
.
data
());
Strain_t
ref
;
ref
.
setZero
();
for
(
Dim_t
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
Dim_t
j
=
0
;
j
<
dim
;
++
j
)
{
for
(
Dim_t
a
=
0
;
a
<
dim
;
++
a
)
{
ref
(
i
,
j
)
+=
A
(
i
,
a
)
*
B
(
a
,
j
);
}
}
}
using
Strain_tw
=
Eigen
::
TensorFixedSize
<
Real
,
Eigen
::
Sizes
<
dim
+
1
,
dim
+
1
>>
;
Strain_tw
C
;
C
.
setConstant
(
100
);
//static_assert(!std::is_convertible<Strain_t, Strain_tw>::value,
// "Tensors not size-protected");
if
(
std
::
is_convertible
<
Strain_t
,
Strain_tw
>::
value
)
{
//std::cout << "this is not good, should I abandon Tensors?";
}
// this test seems useless. I use to detect if Eigen changed the
// default tensor product
Real
error
=
TerrNorm
(
FF1
-
ref
);
if
(
error
<
tol
)
{
std
::
cout
<<
"A ="
<<
std
::
endl
<<
A
<<
std
::
endl
;
std
::
cout
<<
"B ="
<<
std
::
endl
<<
B
<<
std
::
endl
;
std
::
cout
<<
"FF1 ="
<<
std
::
endl
<<
FF1
<<
std
::
endl
;
std
::
cout
<<
"ref ="
<<
std
::
endl
<<
ref
<<
std
::
endl
;
}
BOOST_CHECK_GT
(
error
,
tol
);
error
=
TerrNorm
(
FF2
-
ref
);
if
(
error
>
tol
)
{
std
::
cout
<<
"A ="
<<
std
::
endl
<<
A
<<
std
::
endl
;
std
::
cout
<<
"B ="
<<
std
::
endl
<<
B
<<
std
::
endl
;
std
::
cout
<<
"FF2 ="
<<
std
::
endl
<<
FF2
<<
std
::
endl
;
std
::
cout
<<
"ref ="
<<
std
::
endl
<<
ref
<<
std
::
endl
;
}
BOOST_CHECK_LT
(
error
,
tol
);
error
=
TerrNorm
(
FF3
-
ref
);
if
(
error
>
tol
)
{
std
::
cout
<<
"A ="
<<
std
::
endl
<<
A
<<
std
::
endl
;
std
::
cout
<<
"B ="
<<
std
::
endl
<<
B
<<
std
::
endl
;
std
::
cout
<<
"FF3 ="
<<
std
::
endl
<<
FF3
<<
std
::
endl
;
std
::
cout
<<
"ref ="
<<
std
::
endl
<<
ref
<<
std
::
endl
;
}
BOOST_CHECK_LT
(
error
,
tol
);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE
(
test_tensmult
,
Fix
,
testGoodies
::
dimlist
,
Fix
)
{
using
namespace
Matrices
;
constexpr
Dim_t
dim
{
Fix
::
dim
};
using
T4
=
Tens4_t
<
dim
>
;
using
T2
=
Tens2_t
<
dim
>
;
using
V2
=
Eigen
::
Matrix
<
Real
,
dim
*
dim
,
1
>
;
T4
C
;
C
.
setRandom
();
T2
E
;
E
.
setRandom
();
Eigen
::
Map
<
const
V2
>
Ev
(
E
.
data
());
T2
R
=
tensmult
(
C
,
E
);
auto
error
=
(
Eigen
::
Map
<
const
V2
>
(
R
.
data
())
-
C
*
Ev
).
norm
();
BOOST_CHECK_LT
(
error
,
tol
);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE
(
test_tracer
,
Fix
,
testGoodies
::
dimlist
,
Fix
)
{
using
namespace
Matrices
;
constexpr
Dim_t
dim
{
Fix
::
dim
};
using
T2
=
Tens2_t
<
dim
>
;
auto
tracer
=
Itrac
<
dim
>
();
T2
F
;
F
.
setRandom
();
auto
Ftrac
=
tensmult
(
tracer
,
F
);
auto
error
=
(
Ftrac
-
F
.
trace
()
*
F
.
Identity
()).
norm
();
BOOST_CHECK_LT
(
error
,
tol
);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE
(
test_identity
,
Fix
,
testGoodies
::
dimlist
,
Fix
)
{
using
namespace
Matrices
;
constexpr
Dim_t
dim
{
Fix
::
dim
};
using
T2
=
Tens2_t
<
dim
>
;
auto
ident
=
Iiden
<
dim
>
();
T2
F
;
F
.
setRandom
();
auto
Fiden
=
tensmult
(
ident
,
F
);
auto
error
=
(
Fiden
-
F
).
norm
();
BOOST_CHECK_LT
(
error
,
tol
);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE
(
test_transposer
,
Fix
,
testGoodies
::
dimlist
,
Fix
)
{
using
namespace
Matrices
;
constexpr
Dim_t
dim
{
Fix
::
dim
};
using
T2
=
Tens2_t
<
dim
>
;
auto
trnst
=
Itrns
<
dim
>
();
T2
F
;
F
.
setRandom
();
auto
Ftrns
=
tensmult
(
trnst
,
F
);
auto
error
=
(
Ftrns
-
F
.
transpose
()).
norm
();
BOOST_CHECK_LT
(
error
,
tol
);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE
(
test_symmetriser
,
Fix
,
testGoodies
::
dimlist
,
Fix
)
{
using
namespace
Matrices
;
constexpr
Dim_t
dim
{
Fix
::
dim
};
using
T2
=
Tens2_t
<
dim
>
;
auto
symmt
=
Isymm
<
dim
>
();
T2
F
;
F
.
setRandom
();
auto
Fsymm
=
tensmult
(
symmt
,
F
);
auto
error
=
(
Fsymm
-
.5
*
(
F
+
F
.
transpose
())).
norm
();
BOOST_CHECK_LT
(
error
,
tol
);
}
BOOST_AUTO_TEST_SUITE_END
();
}
// muSpectre
Event Timeline
Log In to Comment