Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F94146067
aka_tensor.hh
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, Dec 4, 06:33
Size
15 KB
Mime Type
text/x-c++
Expires
Fri, Dec 6, 06:33 (1 d, 12 h)
Engine
blob
Format
Raw Data
Handle
22532685
Attached To
rAKA akantu
aka_tensor.hh
View Options
/**
* Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#ifndef AKA_TENSOR_HH_
#define AKA_TENSOR_HH_
namespace
akantu
{
/* -------------------------------------------------------------------------- */
/* TensorBase */
/* -------------------------------------------------------------------------- */
template
<
typename
T
,
Int
ndim
>
class
TensorBase
:
public
TensorTrait
<
ndim
>
{
using
RetType
=
TensorBase
<
T
,
ndim
>
;
static_assert
(
ndim
>
2
,
"TensorBase cannot by used for dimensions < 3"
);
protected
:
using
idx_t
=
Idx
;
template
<
typename
...
Args
>
using
valid_args_t
=
typename
std
::
enable_if
<
aka
::
conjunction
<
aka
::
disjunction
<
std
::
is_integral
<
Args
>
,
std
::
is_enum
<
Args
>>
...
>::
value
and
ndim
==
sizeof
...(
Args
),
int
>::
type
;
public
:
using
proxy
=
TensorBase
<
T
,
ndim
>
;
using
size_type
=
Idx
;
template
<
Int
_ndim
=
ndim
,
std
::
enable_if_t
<
_ndim
==
1
or
_ndim
==
2
,
int
>
=
0
>
TensorBase
()
{
n
.
fill
(
0
);
}
TensorBase
()
{
n
.
fill
(
0
);
}
template
<
typename
...
Args
,
valid_args_t
<
Args
...
>
=
0
>
constexpr
TensorBase
(
Args
...
args
)
:
n
{
idx_t
(
args
)...},
_size
(
detail
::
product_all
(
args
...))
{}
constexpr
TensorBase
(
const
TensorBase
&
other
)
:
n
(
other
.
n
),
_size
(
other
.
_size
),
values
(
other
.
values
)
{}
constexpr
TensorBase
(
TensorBase
&&
other
)
noexcept
:
n
(
std
::
move
(
other
.
n
)),
_size
(
std
::
exchange
(
other
.
_size
,
0
)),
values
(
std
::
exchange
(
other
.
values
,
nullptr
))
{}
protected
:
template
<
typename
Array
,
idx_t
...
I
>
constexpr
auto
check_indices
(
const
Array
&
idx
,
std
::
integer_sequence
<
idx_t
,
I
...
>
/* for_template_deduction */
)
const
{
bool
result
=
true
;
(
void
)
std
::
initializer_list
<
int
>
{(
result
&=
idx
[
I
]
<
n
[
I
],
0
)...};
return
result
;
}
template
<
typename
...
Args
>
constexpr
auto
compute_index
(
Args
...
args
)
const
{
std
::
array
<
idx_t
,
sizeof
...(
Args
)
>
idx
{
idx_t
(
args
)...};
assert
(
check_indices
(
idx
,
std
::
make_integer_sequence
<
idx_t
,
sizeof
...(
Args
)
>
{})
&&
"The indexes are out of bound"
);
idx_t
index
=
0
,
i
=
(
sizeof
...(
Args
))
-
1
;
for
(;
i
>
0
;
i
--
)
{
index
+=
idx
[
i
];
if
(
i
>
0
)
{
index
*=
n
[
i
-
1
];
}
}
return
index
+
idx
[
0
];
}
template
<
typename
S
,
int
...
I
>
constexpr
auto
get_slice
(
idx_t
s
,
std
::
index_sequence
<
I
...
>
)
{
return
S
(
this
->
values
+
s
*
detail
::
product_all
(
n
[
I
]...),
n
[
I
]...);
}
template
<
typename
S
,
std
::
size_t
...
I
>
constexpr
auto
get_slice
(
idx_t
s
,
std
::
index_sequence
<
I
...
>
)
const
{
return
S
(
this
->
values
+
s
*
detail
::
product_all
(
n
[
I
]...),
n
[
I
]...);
}
public
:
template
<
typename
...
Args
,
valid_args_t
<
Args
...
>
=
0
>
inline
auto
operator
()(
Args
...
args
)
->
T
&
{
return
*
(
this
->
values
+
compute_index
(
std
::
move
(
args
)...));
}
template
<
typename
...
Args
,
valid_args_t
<
Args
...
>
=
0
>
inline
auto
operator
()(
Args
...
args
)
const
->
const
T
&
{
return
*
(
this
->
values
+
compute_index
(
std
::
move
(
args
)...));
}
template
<
class
R
=
T
,
idx_t
_ndim
=
ndim
,
std
::
enable_if_t
<
(
_ndim
>
3
)
and
not
std
::
is_const
<
R
>::
value
>
*
=
nullptr
>
inline
auto
operator
()(
idx_t
s
)
{
return
get_slice
<
TensorProxy
<
T
,
ndim
-
1
>>
(
s
,
std
::
make_index_sequence
<
ndim
-
1
>
());
}
template
<
idx_t
_ndim
=
ndim
,
std
::
enable_if_t
<
(
_ndim
>
3
)
>
*
=
nullptr
>
inline
auto
operator
()(
idx_t
s
)
const
{
return
get_slice
<
TensorProxy
<
T
,
ndim
-
1
>>
(
s
,
std
::
make_index_sequence
<
ndim
-
1
>
());
}
template
<
class
R
=
T
,
idx_t
_ndim
=
ndim
,
std
::
enable_if_t
<
(
_ndim
==
3
)
and
not
std
::
is_const
<
R
>::
value
>
*
=
nullptr
>
inline
auto
operator
()(
idx_t
s
)
{
return
get_slice
<
Eigen
::
Map
<
Eigen
::
Matrix
<
T
,
Eigen
::
Dynamic
,
Eigen
::
Dynamic
>>>
(
s
,
std
::
make_index_sequence
<
ndim
-
1
>
());
}
template
<
idx_t
_ndim
=
ndim
,
std
::
enable_if_t
<
_ndim
==
3
>
*
=
nullptr
>
inline
auto
operator
()(
idx_t
s
)
const
{
return
get_slice
<
Eigen
::
Map
<
const
Eigen
::
Matrix
<
std
::
remove_const_t
<
T
>
,
Eigen
::
Dynamic
,
Eigen
::
Dynamic
>>>
(
s
,
std
::
make_index_sequence
<
ndim
-
1
>
());
}
protected
:
template
<
class
Operator
>
auto
transform
(
Operator
&&
op
)
->
RetType
&
{
std
::
transform
(
this
->
values
,
this
->
values
+
this
->
_size
,
this
->
values
,
std
::
forward
<
Operator
>
(
op
));
return
*
(
static_cast
<
RetType
*>
(
this
));
}
template
<
class
Other
,
class
Operator
>
auto
transform
(
Other
&&
other
,
Operator
&&
op
)
->
RetType
&
{
AKANTU_DEBUG_ASSERT
(
_size
==
other
.
size
(),
"The two tensors do not have the same size "
<<
this
->
_size
<<
" != "
<<
other
.
_size
);
std
::
transform
(
this
->
values
,
this
->
values
+
this
->
_size
,
other
.
values
,
this
->
values
,
std
::
forward
<
Operator
>
(
op
));
return
*
(
static_cast
<
RetType
*>
(
this
));
}
template
<
class
Operator
>
auto
accumulate
(
T
init
,
Operator
&&
op
)
->
T
{
return
std
::
accumulate
(
this
->
values
,
this
->
values
+
this
->
_size
,
std
::
move
(
init
),
std
::
forward
<
Operator
>
(
op
));
}
template
<
class
Other
,
class
Init
,
class
Accumulate
,
class
Operator
>
auto
transform_reduce
(
Other
&&
other
,
T
init
,
Accumulate
&&
acc
,
Operator
&&
op
)
->
T
{
return
std
::
inner_product
(
this
->
values
,
this
->
values
+
this
->
_size
,
other
.
data
(),
std
::
move
(
init
),
std
::
forward
<
Accumulate
>
(
acc
),
std
::
forward
<
Operator
>
(
op
));
}
// element wise arithmetic operators -----------------------------------------
public
:
inline
decltype
(
auto
)
operator
+=
(
const
TensorBase
&
other
)
{
return
transform
(
other
,
[](
auto
&&
a
,
auto
&&
b
)
{
return
a
+
b
;
});
}
/* ------------------------------------------------------------------------ */
inline
auto
operator
-=
(
const
TensorBase
&
other
)
->
TensorBase
&
{
return
transform
(
other
,
[](
auto
&&
a
,
auto
&&
b
)
{
return
a
-
b
;
});
}
/* ------------------------------------------------------------------------ */
inline
auto
operator
+=
(
const
T
&
x
)
->
TensorBase
&
{
return
transform
([
&
x
](
auto
&&
a
)
{
return
a
+
x
;
});
}
/* ------------------------------------------------------------------------ */
inline
auto
operator
-=
(
const
T
&
x
)
->
TensorBase
&
{
return
transform
([
&
x
](
auto
&&
a
)
{
return
a
-
x
;
});
}
/* ------------------------------------------------------------------------ */
inline
auto
operator
*=
(
const
T
&
x
)
->
TensorBase
&
{
return
transform
([
&
x
](
auto
&&
a
)
{
return
a
*
x
;
});
}
/* ---------------------------------------------------------------------- */
inline
auto
operator
/=
(
const
T
&
x
)
->
TensorBase
&
{
return
transform
([
&
x
](
auto
&&
a
)
{
return
a
/
x
;
});
}
/// Y = \alpha X + Y
inline
auto
aXplusY
(
const
TensorBase
&
other
,
const
T
alpha
=
1.
)
->
TensorBase
&
{
return
transform
(
other
,
[
&
alpha
](
auto
&&
a
,
auto
&&
b
)
{
return
alpha
*
a
+
b
;
});
}
/* ------------------------------------------------------------------------ */
auto
data
()
->
T
*
{
return
values
;
}
auto
data
()
const
->
const
T
*
{
return
values
;
}
// clang-format off
[[
deprecated
(
"use data instead to be stl compatible"
)]]
auto
storage
()
->
T
*
{
return
values
;
}
[[
deprecated
(
"use data instead to be stl compatible"
)]]
auto
storage
()
const
->
const
T
*
{
return
values
;
}
// clang-format on
auto
size
()
const
{
return
_size
;
}
auto
size
(
idx_t
i
)
const
{
AKANTU_DEBUG_ASSERT
(
i
<
ndim
,
"This tensor has only "
<<
ndim
<<
" dimensions, not "
<<
(
i
+
1
));
return
n
[
i
];
};
inline
void
set
(
const
T
&
t
)
{
std
::
fill_n
(
values
,
_size
,
t
);
};
inline
void
zero
()
{
set
(
T
());
};
public
:
/// "Entrywise" norm norm<L_p> @f[ \|\boldsymbol{T}\|_p = \left(
/// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}}
/// @f]
template
<
Int
norm_type
,
std
::
enable_if_t
<
norm_type
==
Eigen
::
Infinity
>
*
=
nullptr
>
auto
lpNorm
()
const
->
T
{
return
accumulate
(
T
(),
[](
auto
&&
init
,
auto
&&
a
)
{
return
init
+
std
::
abs
(
a
);
});
}
template
<
Int
norm_type
,
std
::
enable_if_t
<
norm_type
==
1
>
*
=
nullptr
>
auto
lpNorm
()
const
->
T
{
return
accumulate
(
T
(),
[](
auto
&&
init
,
auto
&&
a
)
{
return
std
::
max
(
init
,
std
::
abs
(
a
));
});
}
template
<
Int
norm_type
,
std
::
enable_if_t
<
norm_type
==
2
>
*
=
nullptr
>
auto
norm
()
const
->
T
{
return
std
::
sqrt
(
accumulate
(
T
(),
[](
auto
&&
init
,
auto
&&
a
)
{
return
init
+
a
*
a
;
}));
}
template
<
Int
norm_type
,
std
::
enable_if_t
<
(
norm_type
>
2
)
>
*
=
nullptr
>
auto
norm
()
const
->
T
{
return
std
::
pow
(
accumulate
(
T
(),
[](
auto
&&
init
,
auto
&&
a
)
{
return
init
+
std
::
pow
(
a
,
norm_type
);
}),
1.
/
norm_type
);
}
auto
norm
()
const
->
T
{
return
lpNorm
<
2
>
();
}
protected
:
template
<
Int
N
,
typename
...
Args
,
std
::
enable_if_t
<
(
sizeof
...(
Args
)
==
ndim
),
int
>
=
0
>
void
serialize
(
std
::
ostream
&
stream
,
Args
...
args
)
const
{
stream
<<
this
->
operator
()(
std
::
move
(
args
)...);
}
template
<
Int
N
,
typename
...
Args
,
std
::
enable_if_t
<
(
sizeof
...(
Args
)
<
ndim
),
int
>
=
0
>
void
serialize
(
std
::
ostream
&
stream
,
Args
...
args
)
const
{
stream
<<
"["
;
for
(
idx_t
i
=
0
;
i
<
n
[
N
];
++
i
)
{
if
(
i
!=
0
)
{
stream
<<
","
;
}
serialize
<
N
+
1
>
(
stream
,
std
::
move
(
args
)...,
i
);
}
stream
<<
"]"
;
}
public
:
void
printself
(
std
::
ostream
&
stream
)
const
{
serialize
<
0
>
(
stream
);
};
protected
:
template
<
std
::
size_t
...
I
>
constexpr
decltype
(
auto
)
begin
(
std
::
index_sequence
<
I
...
>
)
{
return
view_iterator
<
ViewIteratorHelper_t
<
sizeof
...(
I
),
T
>>
(
values
,
n
[
I
]...);
}
template
<
std
::
size_t
...
I
>
constexpr
decltype
(
auto
)
end
(
std
::
index_sequence
<
I
...
>
)
{
return
view_iterator
<
ViewIteratorHelper_t
<
sizeof
...(
I
),
T
>>
(
values
+
_size
,
n
[
I
]...);
}
template
<
std
::
size_t
...
I
>
constexpr
decltype
(
auto
)
begin
(
std
::
index_sequence
<
I
...
>
)
const
{
return
const_view_iterator
<
ViewIteratorHelper_t
<
sizeof
...(
I
),
T
>>
(
values
,
n
[
I
]...);
}
template
<
std
::
size_t
...
I
>
constexpr
decltype
(
auto
)
end
(
std
::
index_sequence
<
I
...
>
)
const
{
return
const_view_iterator
<
ViewIteratorHelper_t
<
sizeof
...(
I
),
T
>>
(
values
+
_size
,
n
[
I
]...);
}
public
:
decltype
(
auto
)
begin
()
{
return
begin
(
std
::
make_index_sequence
<
ndim
-
1
>
{});
}
decltype
(
auto
)
end
()
{
return
end
(
std
::
make_index_sequence
<
ndim
-
1
>
{});
}
decltype
(
auto
)
begin
()
const
{
return
begin
(
std
::
make_index_sequence
<
ndim
-
1
>
{});
}
decltype
(
auto
)
end
()
const
{
return
end
(
std
::
make_index_sequence
<
ndim
-
1
>
{});
}
protected
:
// size per dimension
std
::
array
<
idx_t
,
ndim
>
n
;
// total storage size
idx_t
_size
{
0
};
// actual data location
T
*
values
{
nullptr
};
};
/* -------------------------------------------------------------------------- */
/* TensorProxy */
/* -------------------------------------------------------------------------- */
template
<
typename
T
,
Int
ndim
>
class
TensorProxy
:
public
TensorBase
<
T
,
ndim
>
{
private
:
using
parent
=
TensorBase
<
T
,
ndim
>
;
public
:
// proxy constructor
template
<
typename
...
Args
>
constexpr
TensorProxy
(
T
*
data
=
reinterpret_cast
<
T
*>
(
0xdeadbeef
),
Args
...
args
)
:
parent
(
args
...)
{
this
->
values
=
data
;
}
constexpr
TensorProxy
(
const
TensorProxy
<
T
,
ndim
>
&
other
)
:
parent
(
other
)
{
this
->
values
=
other
.
values
;
}
constexpr
TensorProxy
(
const
Tensor
<
T
,
ndim
>
&
other
)
:
parent
(
other
)
{
this
->
values
=
other
.
values
;
}
// move constructors ---------------------------------------------------------
// proxy -> proxy
TensorProxy
(
TensorProxy
&&
other
)
noexcept
:
parent
(
other
)
{}
auto
operator
=
(
const
TensorBase
<
T
,
ndim
>
&
other
)
->
TensorProxy
&
{
AKANTU_DEBUG_ASSERT
(
other
.
size
()
==
this
->
size
(),
"You are trying to copy too a tensors proxy with the wrong size "
<<
this
->
_size
<<
" != "
<<
other
.
_size
);
static_assert
(
std
::
is_trivially_copyable
<
T
>
{},
"Cannot copy a tensor on non trivial types"
);
std
::
copy
(
other
.
values
,
other
.
values
+
this
->
_size
,
this
->
values
);
return
*
this
;
}
};
/* -------------------------------------------------------------------------- */
/* Tensor */
/* -------------------------------------------------------------------------- */
template
<
typename
T
,
Int
ndim
>
class
Tensor
:
public
TensorBase
<
T
,
ndim
>
{
private
:
using
parent
=
TensorBase
<
T
,
ndim
>
;
public
:
template
<
typename
...
Args
>
constexpr
Tensor
(
Args
...
args
)
:
parent
(
args
...)
{
static_assert
(
std
::
is_trivially_constructible
<
T
>
{},
"Cannot create a tensor on non trivially constructible types"
);
this
->
values
=
new
T
[
this
->
_size
];
}
/* ------------------------------------------------------------------------ */
virtual
~
Tensor
()
{
delete
[]
this
->
values
;
}
// copy constructors ---------------------------------------------------------
constexpr
Tensor
(
const
Tensor
&
other
)
:
parent
(
other
)
{
this
->
values
=
new
T
[
this
->
_size
];
std
::
copy
(
other
.
values
,
other
.
values
+
this
->
_size
,
this
->
values
);
}
constexpr
explicit
Tensor
(
const
TensorProxy
<
T
,
ndim
>
&
other
)
:
parent
(
other
)
{
// static_assert(false, "Copying data are you sure");
this
->
values
=
new
T
[
this
->
_size
];
std
::
copy
(
other
.
values
,
other
.
values
+
this
->
_size
,
this
->
values
);
}
// move constructors ---------------------------------------------------------
// proxy -> proxy, non proxy -> non proxy
Tensor
(
Tensor
&&
other
)
noexcept
:
parent
(
other
)
{}
// copy operator -------------------------------------------------------------
/// operator= copy-and-swap
auto
operator
=
(
const
TensorBase
<
T
,
ndim
>
&
other
)
->
Tensor
&
{
if
(
&
other
==
this
)
return
*
this
;
std
::
cout
<<
"Warning: operator= delete data"
<<
std
::
endl
;
delete
[]
this
->
values
;
this
->
n
=
other
.
n
;
this
->
_size
=
other
.
_size
;
static_assert
(
std
::
is_trivially_constructible
<
T
>
{},
"Cannot create a tensor on non trivially constructible types"
);
this
->
values
=
new
T
[
this
->
_size
];
static_assert
(
std
::
is_trivially_copyable
<
T
>
{},
"Cannot copy a tensor on non trivial types"
);
std
::
copy
(
other
.
values
,
other
.
values
+
this
->
_size
,
this
->
values
);
return
*
this
;
}
};
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
using
Tensor3
=
Tensor
<
T
,
3
>
;
template
<
typename
T
>
using
Tensor3Proxy
=
TensorProxy
<
T
,
3
>
;
template
<
typename
T
>
using
Tensor3Base
=
TensorBase
<
T
,
3
>
;
}
// namespace akantu
#endif
// AKA_TENSOR_HH_
Event Timeline
Log In to Comment