Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79505653
aka_types.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
Mon, Aug 26, 19:31
Size
45 KB
Mime Type
text/x-c++
Expires
Wed, Aug 28, 19:31 (2 d)
Engine
blob
Format
Raw Data
Handle
20100653
Attached To
rAKA akantu
aka_types.hh
View Options
/**
* @file aka_types.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 17 2011
* @date last modification: Tue Feb 20 2018
*
* @brief description of the "simple" types
*
* @section LICENSE
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_error.hh"
#include "aka_fwd.hh"
#include "aka_math.hh"
/* -------------------------------------------------------------------------- */
#include <initializer_list>
#include <iomanip>
#include <type_traits>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_TYPES_HH__
#define __AKANTU_AKA_TYPES_HH__
namespace
akantu
{
enum
NormType
{
L_1
=
1
,
L_2
=
2
,
L_inf
=
UInt
(
-
1
)
};
/**
* DimHelper is a class to generalize the setup of a dim array from 3
* values. This gives a common interface in the TensorStorage class
* independently of its derived inheritance (Vector, Matrix, Tensor3)
* @tparam dim
*/
template
<
UInt
dim
>
struct
DimHelper
{
static
inline
void
setDims
(
UInt
m
,
UInt
n
,
UInt
p
,
UInt
dims
[
dim
]);
};
/* -------------------------------------------------------------------------- */
template
<>
struct
DimHelper
<
1
>
{
static
inline
void
setDims
(
UInt
m
,
__attribute__
((
unused
))
UInt
n
,
__attribute__
((
unused
))
UInt
p
,
UInt
dims
[
1
])
{
dims
[
0
]
=
m
;
}
};
/* -------------------------------------------------------------------------- */
template
<>
struct
DimHelper
<
2
>
{
static
inline
void
setDims
(
UInt
m
,
UInt
n
,
__attribute__
((
unused
))
UInt
p
,
UInt
dims
[
2
])
{
dims
[
0
]
=
m
;
dims
[
1
]
=
n
;
}
};
/* -------------------------------------------------------------------------- */
template
<>
struct
DimHelper
<
3
>
{
static
inline
void
setDims
(
UInt
m
,
UInt
n
,
UInt
p
,
UInt
dims
[
3
])
{
dims
[
0
]
=
m
;
dims
[
1
]
=
n
;
dims
[
2
]
=
p
;
}
};
/* -------------------------------------------------------------------------- */
template
<
typename
T
,
UInt
ndim
,
class
RetType
>
class
TensorStorage
;
/* -------------------------------------------------------------------------- */
/* Proxy classes */
/* -------------------------------------------------------------------------- */
namespace
tensors
{
template
<
class
A
,
class
B
>
struct
is_copyable
{
enum
:
bool
{
value
=
false
};
};
template
<
class
A
>
struct
is_copyable
<
A
,
A
>
{
enum
:
bool
{
value
=
true
};
};
template
<
class
A
>
struct
is_copyable
<
A
,
typename
A
::
RetType
>
{
enum
:
bool
{
value
=
true
};
};
template
<
class
A
>
struct
is_copyable
<
A
,
typename
A
::
RetType
::
proxy
>
{
enum
:
bool
{
value
=
true
};
};
}
// namespace tensors
/**
* @class TensorProxy aka_types.hh
* @desc The TensorProxy class is a proxy class to the TensorStorage it handles
* the
* wrapped case. That is to say if an accessor should give access to a Tensor
* wrapped on some data, like the Array<T>::iterator they can return a
* TensorProxy that will be automatically transformed as a TensorStorage wrapped
* on the same data
* @tparam T stored type
* @tparam ndim order of the tensor
* @tparam RetType real derived type
*/
template
<
typename
T
,
UInt
ndim
,
class
_RetType
>
class
TensorProxy
{
protected
:
using
RetTypeProxy
=
typename
_RetType
::
proxy
;
constexpr
TensorProxy
(
T
*
data
,
UInt
m
,
UInt
n
,
UInt
p
)
{
DimHelper
<
ndim
>::
setDims
(
m
,
n
,
p
,
this
->
n
);
this
->
values
=
data
;
}
#ifndef SWIG
template
<
class
Other
,
typename
=
std
::
enable_if_t
<
tensors
::
is_copyable
<
TensorProxy
,
Other
>::
value
>>
explicit
TensorProxy
(
const
Other
&
other
)
{
this
->
values
=
other
.
storage
();
for
(
UInt
i
=
0
;
i
<
ndim
;
++
i
)
this
->
n
[
i
]
=
other
.
size
(
i
);
}
#endif
public
:
using
RetType
=
_RetType
;
UInt
size
(
UInt
i
)
const
{
AKANTU_DEBUG_ASSERT
(
i
<
ndim
,
"This tensor has only "
<<
ndim
<<
" dimensions, not "
<<
(
i
+
1
));
return
n
[
i
];
}
inline
UInt
size
()
const
{
UInt
_size
=
1
;
for
(
UInt
d
=
0
;
d
<
ndim
;
++
d
)
_size
*=
this
->
n
[
d
];
return
_size
;
}
T
*
storage
()
const
{
return
values
;
}
#ifndef SWIG
template
<
class
Other
,
typename
=
std
::
enable_if_t
<
tensors
::
is_copyable
<
TensorProxy
,
Other
>::
value
>>
inline
TensorProxy
&
operator
=
(
const
Other
&
other
)
{
AKANTU_DEBUG_ASSERT
(
other
.
size
()
==
this
->
size
(),
"You are trying to copy two tensors with different sizes"
);
memcpy
(
this
->
values
,
other
.
storage
(),
this
->
size
()
*
sizeof
(
T
));
return
*
this
;
}
#endif
// template <class Other, typename = std::enable_if_t<
// tensors::is_copyable<TensorProxy, Other>::value>>
// inline TensorProxy & operator=(const Other && other) {
// AKANTU_DEBUG_ASSERT(
// other.size() == this->size(),
// "You are trying to copy two tensors with different sizes");
// memcpy(this->values, other.storage(), this->size() * sizeof(T));
// return *this;
// }
template
<
typename
O
>
inline
RetTypeProxy
&
operator
*=
(
const
O
&
o
)
{
RetType
(
*
this
)
*=
o
;
return
static_cast
<
RetTypeProxy
&>
(
*
this
);
}
template
<
typename
O
>
inline
RetTypeProxy
&
operator
/=
(
const
O
&
o
)
{
RetType
(
*
this
)
/=
o
;
return
static_cast
<
RetTypeProxy
&>
(
*
this
);
}
protected
:
T
*
values
;
UInt
n
[
ndim
];
};
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
class
VectorProxy
:
public
TensorProxy
<
T
,
1
,
Vector
<
T
>>
{
using
parent
=
TensorProxy
<
T
,
1
,
Vector
<
T
>>
;
using
type
=
Vector
<
T
>
;
public
:
constexpr
VectorProxy
(
T
*
data
,
UInt
n
)
:
parent
(
data
,
n
,
0
,
0
)
{}
template
<
class
Other
>
explicit
VectorProxy
(
Other
&
src
)
:
parent
(
src
)
{}
/* ---------------------------------------------------------------------- */
template
<
class
Other
>
inline
VectorProxy
<
T
>
&
operator
=
(
const
Other
&
other
)
{
parent
::
operator
=
(
other
);
return
*
this
;
}
// inline VectorProxy<T> & operator=(const VectorProxy && other) {
// parent::operator=(other);
// return *this;
// }
/* ------------------------------------------------------------------------ */
T
&
operator
()(
UInt
index
)
{
return
this
->
values
[
index
];
};
const
T
&
operator
()(
UInt
index
)
const
{
return
this
->
values
[
index
];
};
};
template
<
typename
T
>
class
MatrixProxy
:
public
TensorProxy
<
T
,
2
,
Matrix
<
T
>>
{
using
parent
=
TensorProxy
<
T
,
2
,
Matrix
<
T
>>
;
using
type
=
Matrix
<
T
>
;
public
:
MatrixProxy
(
T
*
data
,
UInt
m
,
UInt
n
)
:
parent
(
data
,
m
,
n
,
0
)
{}
template
<
class
Other
>
explicit
MatrixProxy
(
Other
&
src
)
:
parent
(
src
)
{}
/* ---------------------------------------------------------------------- */
template
<
class
Other
>
inline
MatrixProxy
<
T
>
&
operator
=
(
const
Other
&
other
)
{
parent
::
operator
=
(
other
);
return
*
this
;
}
};
template
<
typename
T
>
class
Tensor3Proxy
:
public
TensorProxy
<
T
,
3
,
Tensor3
<
T
>>
{
using
parent
=
TensorProxy
<
T
,
3
,
Tensor3
<
T
>>
;
using
type
=
Tensor3
<
T
>
;
public
:
Tensor3Proxy
(
const
T
*
data
,
UInt
m
,
UInt
n
,
UInt
k
)
:
parent
(
data
,
m
,
n
,
k
)
{}
Tensor3Proxy
(
const
Tensor3Proxy
&
src
)
:
parent
(
src
)
{}
Tensor3Proxy
(
const
Tensor3
<
T
>
&
src
)
:
parent
(
src
)
{}
/* ---------------------------------------------------------------------- */
template
<
class
Other
>
inline
Tensor3Proxy
<
T
>
&
operator
=
(
const
Other
&
other
)
{
parent
::
operator
=
(
other
);
return
*
this
;
}
};
/* -------------------------------------------------------------------------- */
/* Tensor base class */
/* -------------------------------------------------------------------------- */
template
<
typename
T
,
UInt
ndim
,
class
RetType
>
class
TensorStorage
:
public
TensorTrait
{
public
:
using
value_type
=
T
;
friend
class
Array
<
T
>
;
protected
:
template
<
class
TensorType
>
void
copySize
(
const
TensorType
&
src
)
{
for
(
UInt
d
=
0
;
d
<
ndim
;
++
d
)
this
->
n
[
d
]
=
src
.
size
(
d
);
this
->
_size
=
src
.
size
();
}
TensorStorage
()
:
values
(
nullptr
)
{
for
(
UInt
d
=
0
;
d
<
ndim
;
++
d
)
this
->
n
[
d
]
=
0
;
_size
=
0
;
}
TensorStorage
(
const
TensorProxy
<
T
,
ndim
,
RetType
>
&
proxy
)
{
this
->
copySize
(
proxy
);
this
->
values
=
proxy
.
storage
();
this
->
wrapped
=
true
;
}
public
:
TensorStorage
(
const
TensorStorage
&
src
)
=
delete
;
TensorStorage
(
const
TensorStorage
&
src
,
bool
deep_copy
)
:
values
(
nullptr
)
{
if
(
deep_copy
)
this
->
deepCopy
(
src
);
else
this
->
shallowCopy
(
src
);
}
protected
:
TensorStorage
(
UInt
m
,
UInt
n
,
UInt
p
,
const
T
&
def
)
{
static_assert
(
std
::
is_trivially_constructible
<
T
>
{},
"Cannot create a tensor on non trivial types"
);
DimHelper
<
ndim
>::
setDims
(
m
,
n
,
p
,
this
->
n
);
this
->
computeSize
();
this
->
values
=
new
T
[
this
->
_size
];
this
->
set
(
def
);
this
->
wrapped
=
false
;
}
TensorStorage
(
T
*
data
,
UInt
m
,
UInt
n
,
UInt
p
)
{
DimHelper
<
ndim
>::
setDims
(
m
,
n
,
p
,
this
->
n
);
this
->
computeSize
();
this
->
values
=
data
;
this
->
wrapped
=
true
;
}
public
:
/* ------------------------------------------------------------------------ */
template
<
class
TensorType
>
inline
void
shallowCopy
(
const
TensorType
&
src
)
{
this
->
copySize
(
src
);
if
(
!
this
->
wrapped
)
delete
[]
this
->
values
;
this
->
values
=
src
.
storage
();
this
->
wrapped
=
true
;
}
/* ------------------------------------------------------------------------ */
template
<
class
TensorType
>
inline
void
deepCopy
(
const
TensorType
&
src
)
{
this
->
copySize
(
src
);
if
(
!
this
->
wrapped
)
delete
[]
this
->
values
;
static_assert
(
std
::
is_trivially_constructible
<
T
>
{},
"Cannot create a tensor on non trivial types"
);
this
->
values
=
new
T
[
this
->
_size
];
static_assert
(
std
::
is_trivially_copyable
<
T
>
{},
"Cannot copy a tensor on non trivial types"
);
memcpy
((
void
*
)
this
->
values
,
(
void
*
)
src
.
storage
(),
this
->
_size
*
sizeof
(
T
));
this
->
wrapped
=
false
;
}
virtual
~
TensorStorage
()
{
if
(
!
this
->
wrapped
)
delete
[]
this
->
values
;
}
/* ------------------------------------------------------------------------ */
inline
TensorStorage
&
operator
=
(
const
TensorStorage
&
src
)
{
return
this
->
operator
=
(
dynamic_cast
<
RetType
&>
(
src
));
}
/* ------------------------------------------------------------------------ */
inline
TensorStorage
&
operator
=
(
const
RetType
&
src
)
{
if
(
this
!=
&
src
)
{
if
(
this
->
wrapped
)
{
static_assert
(
std
::
is_trivially_copyable
<
T
>
{},
"Cannot copy a tensor on non trivial types"
);
// this test is not sufficient for Tensor of order higher than 1
AKANTU_DEBUG_ASSERT
(
this
->
_size
==
src
.
size
(),
"Tensors of different size"
);
memcpy
((
void
*
)
this
->
values
,
(
void
*
)
src
.
storage
(),
this
->
_size
*
sizeof
(
T
));
}
else
{
deepCopy
(
src
);
}
}
return
*
this
;
}
/* ------------------------------------------------------------------------ */
template
<
class
R
>
inline
RetType
&
operator
+=
(
const
TensorStorage
<
T
,
ndim
,
R
>
&
other
)
{
T
*
a
=
this
->
storage
();
T
*
b
=
other
.
storage
();
AKANTU_DEBUG_ASSERT
(
_size
==
other
.
size
(),
"The two tensors do not have the same size, they cannot be subtracted"
);
for
(
UInt
i
=
0
;
i
<
_size
;
++
i
)
*
(
a
++
)
+=
*
(
b
++
);
return
*
(
static_cast
<
RetType
*>
(
this
));
}
/* ------------------------------------------------------------------------ */
template
<
class
R
>
inline
RetType
&
operator
-=
(
const
TensorStorage
<
T
,
ndim
,
R
>
&
other
)
{
T
*
a
=
this
->
storage
();
T
*
b
=
other
.
storage
();
AKANTU_DEBUG_ASSERT
(
_size
==
other
.
size
(),
"The two tensors do not have the same size, they cannot be subtracted"
);
for
(
UInt
i
=
0
;
i
<
_size
;
++
i
)
*
(
a
++
)
-=
*
(
b
++
);
return
*
(
static_cast
<
RetType
*>
(
this
));
}
/* ------------------------------------------------------------------------ */
inline
RetType
&
operator
+=
(
const
T
&
x
)
{
T
*
a
=
this
->
values
;
for
(
UInt
i
=
0
;
i
<
_size
;
++
i
)
*
(
a
++
)
+=
x
;
return
*
(
static_cast
<
RetType
*>
(
this
));
}
/* ------------------------------------------------------------------------ */
inline
RetType
&
operator
-=
(
const
T
&
x
)
{
T
*
a
=
this
->
values
;
for
(
UInt
i
=
0
;
i
<
_size
;
++
i
)
*
(
a
++
)
-=
x
;
return
*
(
static_cast
<
RetType
*>
(
this
));
}
/* ------------------------------------------------------------------------ */
inline
RetType
&
operator
*=
(
const
T
&
x
)
{
T
*
a
=
this
->
storage
();
for
(
UInt
i
=
0
;
i
<
_size
;
++
i
)
*
(
a
++
)
*=
x
;
return
*
(
static_cast
<
RetType
*>
(
this
));
}
/* ---------------------------------------------------------------------- */
inline
RetType
&
operator
/=
(
const
T
&
x
)
{
T
*
a
=
this
->
values
;
for
(
UInt
i
=
0
;
i
<
_size
;
++
i
)
*
(
a
++
)
/=
x
;
return
*
(
static_cast
<
RetType
*>
(
this
));
}
/// Y = \alpha X + Y
inline
RetType
&
aXplusY
(
const
TensorStorage
&
other
,
const
T
&
alpha
=
1.
)
{
AKANTU_DEBUG_ASSERT
(
_size
==
other
.
size
(),
"The two tensors do not have the same size, they cannot be subtracted"
);
Math
::
aXplusY
(
this
->
_size
,
alpha
,
other
.
storage
(),
this
->
storage
());
return
*
(
static_cast
<
RetType
*>
(
this
));
}
/* ------------------------------------------------------------------------ */
T
*
storage
()
const
{
return
values
;
}
UInt
size
()
const
{
return
_size
;
}
UInt
size
(
UInt
i
)
const
{
AKANTU_DEBUG_ASSERT
(
i
<
ndim
,
"This tensor has only "
<<
ndim
<<
" dimensions, not "
<<
(
i
+
1
));
return
n
[
i
];
};
/* ------------------------------------------------------------------------ */
inline
void
clear
()
{
memset
(
values
,
0
,
_size
*
sizeof
(
T
));
};
inline
void
set
(
const
T
&
t
)
{
std
::
fill_n
(
values
,
_size
,
t
);
};
template
<
class
TensorType
>
inline
void
copy
(
const
TensorType
&
other
)
{
AKANTU_DEBUG_ASSERT
(
_size
==
other
.
size
(),
"The two tensors do not have the same size, they cannot be copied"
);
memcpy
(
values
,
other
.
storage
(),
_size
*
sizeof
(
T
));
}
bool
isWrapped
()
const
{
return
this
->
wrapped
;
}
protected
:
inline
void
computeSize
()
{
_size
=
1
;
for
(
UInt
d
=
0
;
d
<
ndim
;
++
d
)
_size
*=
this
->
n
[
d
];
}
protected
:
template
<
typename
R
,
NormType
norm_type
>
struct
NormHelper
{
template
<
class
Ten
>
static
R
norm
(
const
Ten
&
ten
)
{
R
_norm
=
0.
;
R
*
it
=
ten
.
storage
();
R
*
end
=
ten
.
storage
()
+
ten
.
size
();
for
(;
it
<
end
;
++
it
)
_norm
+=
std
::
pow
(
std
::
abs
(
*
it
),
norm_type
);
return
std
::
pow
(
_norm
,
1.
/
norm_type
);
}
};
template
<
typename
R
>
struct
NormHelper
<
R
,
L_1
>
{
template
<
class
Ten
>
static
R
norm
(
const
Ten
&
ten
)
{
R
_norm
=
0.
;
R
*
it
=
ten
.
storage
();
R
*
end
=
ten
.
storage
()
+
ten
.
size
();
for
(;
it
<
end
;
++
it
)
_norm
+=
std
::
abs
(
*
it
);
return
_norm
;
}
};
template
<
typename
R
>
struct
NormHelper
<
R
,
L_2
>
{
template
<
class
Ten
>
static
R
norm
(
const
Ten
&
ten
)
{
R
_norm
=
0.
;
R
*
it
=
ten
.
storage
();
R
*
end
=
ten
.
storage
()
+
ten
.
size
();
for
(;
it
<
end
;
++
it
)
_norm
+=
*
it
*
*
it
;
return
sqrt
(
_norm
);
}
};
template
<
typename
R
>
struct
NormHelper
<
R
,
L_inf
>
{
template
<
class
Ten
>
static
R
norm
(
const
Ten
&
ten
)
{
R
_norm
=
0.
;
R
*
it
=
ten
.
storage
();
R
*
end
=
ten
.
storage
()
+
ten
.
size
();
for
(;
it
<
end
;
++
it
)
_norm
=
std
::
max
(
std
::
abs
(
*
it
),
_norm
);
return
_norm
;
}
};
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
<
NormType
norm_type
>
inline
T
norm
()
const
{
return
NormHelper
<
T
,
norm_type
>::
norm
(
*
this
);
}
protected
:
UInt
n
[
ndim
];
UInt
_size
;
T
*
values
;
bool
wrapped
{
false
};
};
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
namespace
types
{
namespace
details
{
template
<
typename
reference_
>
class
vector_iterator
{
public
:
using
difference_type
=
std
::
ptrdiff_t
;
using
value_type
=
std
::
decay_t
<
reference_
>
;
using
pointer
=
value_type
*
;
using
reference
=
reference_
;
using
iterator_category
=
std
::
input_iterator_tag
;
vector_iterator
(
pointer
ptr
)
:
ptr
(
ptr
)
{}
// input iterator ++it
vector_iterator
&
operator
++
()
{
++
ptr
;
return
*
this
;
}
// input iterator it++
vector_iterator
operator
++
(
int
)
{
auto
cpy
=
*
this
;
++
ptr
;
return
cpy
;
}
vector_iterator
&
operator
+=
(
int
n
)
{
ptr
+=
n
;
return
*
this
;
}
vector_iterator
operator
+
(
int
n
)
{
vector_iterator
cpy
(
*
this
);
cpy
+=
n
;
return
cpy
;
}
// input iterator it != other_it
bool
operator
!=
(
const
vector_iterator
&
other
)
const
{
return
ptr
!=
other
.
ptr
;
}
bool
operator
==
(
const
vector_iterator
&
other
)
const
{
return
ptr
==
other
.
ptr
;
}
difference_type
operator
-
(
const
vector_iterator
&
other
)
const
{
return
this
->
ptr
-
other
.
ptr
;
}
// input iterator dereference *it
reference
operator
*
()
{
return
*
ptr
;
}
private
:
pointer
ptr
;
};
}
// namespace details
}
// namespace types
/* -------------------------------------------------------------------------- */
/* Vector */
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
class
Vector
:
public
TensorStorage
<
T
,
1
,
Vector
<
T
>>
{
using
parent
=
TensorStorage
<
T
,
1
,
Vector
<
T
>>
;
public
:
using
value_type
=
typename
parent
::
value_type
;
using
proxy
=
VectorProxy
<
T
>
;
public
:
Vector
()
:
parent
()
{}
explicit
Vector
(
UInt
n
,
const
T
&
def
=
T
())
:
parent
(
n
,
0
,
0
,
def
)
{}
Vector
(
T
*
data
,
UInt
n
)
:
parent
(
data
,
n
,
0
,
0
)
{}
Vector
(
const
Vector
&
src
,
bool
deep_copy
=
true
)
:
parent
(
src
,
deep_copy
)
{}
Vector
(
const
TensorProxy
<
T
,
1
,
Vector
>
&
src
)
:
parent
(
src
)
{}
Vector
(
std
::
initializer_list
<
T
>
list
)
:
parent
(
list
.
size
(),
0
,
0
,
T
())
{
UInt
i
=
0
;
for
(
auto
val
:
list
)
{
operator
()(
i
++
)
=
val
;
}
}
public
:
using
iterator
=
types
::
details
::
vector_iterator
<
T
&>
;
using
const_iterator
=
types
::
details
::
vector_iterator
<
const
T
&>
;
iterator
begin
()
{
return
iterator
(
this
->
storage
());
}
iterator
end
()
{
return
iterator
(
this
->
storage
()
+
this
->
size
());
}
const_iterator
begin
()
const
{
return
const_iterator
(
this
->
storage
());
}
const_iterator
end
()
const
{
return
const_iterator
(
this
->
storage
()
+
this
->
size
());
}
public
:
~
Vector
()
override
=
default
;
/* ------------------------------------------------------------------------ */
inline
Vector
&
operator
=
(
const
Vector
&
src
)
{
parent
::
operator
=
(
src
);
return
*
this
;
}
/* ------------------------------------------------------------------------ */
inline
T
&
operator
()(
UInt
i
)
{
AKANTU_DEBUG_ASSERT
((
i
<
this
->
n
[
0
]),
"Access out of the vector! "
<<
"Index ("
<<
i
<<
") is out of the vector of size ("
<<
this
->
n
[
0
]
<<
")"
);
return
*
(
this
->
values
+
i
);
}
inline
const
T
&
operator
()(
UInt
i
)
const
{
AKANTU_DEBUG_ASSERT
((
i
<
this
->
n
[
0
]),
"Access out of the vector! "
<<
"Index ("
<<
i
<<
") is out of the vector of size ("
<<
this
->
n
[
0
]
<<
")"
);
return
*
(
this
->
values
+
i
);
}
inline
T
&
operator
[](
UInt
i
)
{
return
this
->
operator
()(
i
);
}
inline
const
T
&
operator
[](
UInt
i
)
const
{
return
this
->
operator
()(
i
);
}
/* ------------------------------------------------------------------------ */
inline
Vector
<
T
>
&
operator
*=
(
Real
x
)
{
return
parent
::
operator
*=
(
x
);
}
inline
Vector
<
T
>
&
operator
/=
(
Real
x
)
{
return
parent
::
operator
/=
(
x
);
}
/* ------------------------------------------------------------------------ */
inline
Vector
<
T
>
&
operator
*=
(
const
Vector
<
T
>
&
vect
)
{
AKANTU_DEBUG_ASSERT
(
this
->
_size
==
vect
.
_size
,
"The vectors have non matching sizes"
);
T
*
a
=
this
->
storage
();
T
*
b
=
vect
.
storage
();
for
(
UInt
i
=
0
;
i
<
this
->
_size
;
++
i
)
*
(
a
++
)
*=
*
(
b
++
);
return
*
this
;
}
/* ------------------------------------------------------------------------ */
inline
Real
dot
(
const
Vector
<
T
>
&
vect
)
const
{
return
Math
::
vectorDot
(
this
->
values
,
vect
.
storage
(),
this
->
_size
);
}
/* ------------------------------------------------------------------------ */
inline
Real
mean
()
const
{
Real
mean
=
0
;
T
*
a
=
this
->
storage
();
for
(
UInt
i
=
0
;
i
<
this
->
_size
;
++
i
)
mean
+=
*
(
a
++
);
return
mean
/
this
->
_size
;
}
/* ------------------------------------------------------------------------ */
inline
Vector
&
crossProduct
(
const
Vector
<
T
>
&
v1
,
const
Vector
<
T
>
&
v2
)
{
AKANTU_DEBUG_ASSERT
(
this
->
size
()
==
3
,
"crossProduct is only defined in 3D (n="
<<
this
->
size
()
<<
")"
);
AKANTU_DEBUG_ASSERT
(
this
->
size
()
==
v1
.
size
()
&&
this
->
size
()
==
v2
.
size
(),
"crossProduct is not a valid operation non matching size vectors"
);
Math
::
vectorProduct3
(
v1
.
storage
(),
v2
.
storage
(),
this
->
values
);
return
*
this
;
}
inline
Vector
crossProduct
(
const
Vector
<
T
>
&
v
)
{
Vector
<
T
>
tmp
(
this
->
size
());
tmp
.
crossProduct
(
*
this
,
v
);
return
tmp
;
}
/* ------------------------------------------------------------------------ */
inline
void
solve
(
const
Matrix
<
T
>
&
A
,
const
Vector
<
T
>
&
b
)
{
AKANTU_DEBUG_ASSERT
(
this
->
size
()
==
A
.
rows
()
&&
this
->
_size
==
A
.
cols
(),
"The size of the solution vector mismatches the size of the matrix"
);
AKANTU_DEBUG_ASSERT
(
this
->
_size
==
b
.
_size
,
"The rhs vector has a mismatch in size with the matrix"
);
Math
::
solve
(
this
->
_size
,
A
.
storage
(),
this
->
values
,
b
.
storage
());
}
/* ------------------------------------------------------------------------ */
template
<
bool
tr_A
>
inline
void
mul
(
const
Matrix
<
T
>
&
A
,
const
Vector
<
T
>
&
x
,
Real
alpha
=
1.0
);
/* ------------------------------------------------------------------------ */
inline
Real
norm
()
const
{
return
parent
::
template
norm
<
L_2
>
();
}
template
<
NormType
nt
>
inline
Real
norm
()
const
{
return
parent
::
template
norm
<
nt
>
();
}
/* ------------------------------------------------------------------------ */
inline
Vector
<
Real
>
&
normalize
()
{
Real
n
=
norm
();
operator
/=
(
n
);
return
*
this
;
}
/* ------------------------------------------------------------------------ */
/// norm of (*this - x)
inline
Real
distance
(
const
Vector
<
T
>
&
y
)
const
{
Real
*
vx
=
this
->
values
;
Real
*
vy
=
y
.
storage
();
Real
sum_2
=
0
;
for
(
UInt
i
=
0
;
i
<
this
->
_size
;
++
i
,
++
vx
,
++
vy
)
sum_2
+=
(
*
vx
-
*
vy
)
*
(
*
vx
-
*
vy
);
return
sqrt
(
sum_2
);
}
/* ------------------------------------------------------------------------ */
inline
bool
equal
(
const
Vector
<
T
>
&
v
,
Real
tolerance
=
Math
::
getTolerance
())
const
{
T
*
a
=
this
->
storage
();
T
*
b
=
v
.
storage
();
UInt
i
=
0
;
while
(
i
<
this
->
_size
&&
(
std
::
abs
(
*
(
a
++
)
-
*
(
b
++
))
<
tolerance
))
++
i
;
return
i
==
this
->
_size
;
}
/* ------------------------------------------------------------------------ */
inline
short
compare
(
const
Vector
<
T
>
&
v
,
Real
tolerance
=
Math
::
getTolerance
())
const
{
T
*
a
=
this
->
storage
();
T
*
b
=
v
.
storage
();
for
(
UInt
i
(
0
);
i
<
this
->
_size
;
++
i
,
++
a
,
++
b
)
{
if
(
std
::
abs
(
*
a
-
*
b
)
>
tolerance
)
return
(((
*
a
-
*
b
)
>
tolerance
)
?
1
:
-
1
);
}
return
0
;
}
/* ------------------------------------------------------------------------ */
inline
bool
operator
==
(
const
Vector
<
T
>
&
v
)
const
{
return
equal
(
v
);
}
inline
bool
operator
!=
(
const
Vector
<
T
>
&
v
)
const
{
return
!
operator
==
(
v
);
}
inline
bool
operator
<
(
const
Vector
<
T
>
&
v
)
const
{
return
compare
(
v
)
==
-
1
;
}
inline
bool
operator
>
(
const
Vector
<
T
>
&
v
)
const
{
return
compare
(
v
)
==
1
;
}
/* ------------------------------------------------------------------------ */
/// function to print the containt of the class
virtual
void
printself
(
std
::
ostream
&
stream
,
int
indent
=
0
)
const
{
std
::
string
space
;
for
(
Int
i
=
0
;
i
<
indent
;
i
++
,
space
+=
AKANTU_INDENT
)
;
stream
<<
"["
;
for
(
UInt
i
=
0
;
i
<
this
->
_size
;
++
i
)
{
if
(
i
!=
0
)
stream
<<
", "
;
stream
<<
this
->
values
[
i
];
}
stream
<<
"]"
;
}
// friend class ::akantu::Array<T>;
};
using
RVector
=
Vector
<
Real
>
;
/* ------------------------------------------------------------------------ */
template
<>
inline
bool
Vector
<
UInt
>::
equal
(
const
Vector
<
UInt
>
&
v
,
__attribute__
((
unused
))
Real
tolerance
)
const
{
UInt
*
a
=
this
->
storage
();
UInt
*
b
=
v
.
storage
();
UInt
i
=
0
;
while
(
i
<
this
->
_size
&&
(
*
(
a
++
)
==
*
(
b
++
)))
++
i
;
return
i
==
this
->
_size
;
}
/* ------------------------------------------------------------------------ */
/* Matrix */
/* ------------------------------------------------------------------------ */
template
<
typename
T
>
class
Matrix
:
public
TensorStorage
<
T
,
2
,
Matrix
<
T
>>
{
using
parent
=
TensorStorage
<
T
,
2
,
Matrix
<
T
>>
;
public
:
using
value_type
=
typename
parent
::
value_type
;
using
proxy
=
MatrixProxy
<
T
>
;
public
:
Matrix
()
:
parent
()
{}
Matrix
(
UInt
m
,
UInt
n
,
const
T
&
def
=
T
())
:
parent
(
m
,
n
,
0
,
def
)
{}
Matrix
(
T
*
data
,
UInt
m
,
UInt
n
)
:
parent
(
data
,
m
,
n
,
0
)
{}
Matrix
(
const
Matrix
&
src
,
bool
deep_copy
=
true
)
:
parent
(
src
,
deep_copy
)
{}
Matrix
(
const
MatrixProxy
<
T
>
&
src
)
:
parent
(
src
)
{}
Matrix
(
std
::
initializer_list
<
std
::
initializer_list
<
T
>>
list
)
{
static_assert
(
std
::
is_trivially_copyable
<
T
>
{},
"Cannot create a tensor on non trivial types"
);
std
::
size_t
n
=
0
;
std
::
size_t
m
=
list
.
size
();
for
(
auto
row
:
list
)
{
n
=
std
::
max
(
n
,
row
.
size
());
}
DimHelper
<
2
>::
setDims
(
m
,
n
,
0
,
this
->
n
);
this
->
computeSize
();
this
->
values
=
new
T
[
this
->
_size
];
this
->
set
(
0
);
UInt
i
=
0
,
j
=
0
;
for
(
auto
&
row
:
list
)
{
for
(
auto
&
val
:
row
)
{
at
(
i
,
j
++
)
=
val
;
}
++
i
;
j
=
0
;
}
}
~
Matrix
()
override
=
default
;
/* ------------------------------------------------------------------------ */
inline
Matrix
&
operator
=
(
const
Matrix
&
src
)
{
parent
::
operator
=
(
src
);
return
*
this
;
}
public
:
/* ---------------------------------------------------------------------- */
UInt
rows
()
const
{
return
this
->
n
[
0
];
}
UInt
cols
()
const
{
return
this
->
n
[
1
];
}
/* ---------------------------------------------------------------------- */
inline
T
&
at
(
UInt
i
,
UInt
j
)
{
AKANTU_DEBUG_ASSERT
(((
i
<
this
->
n
[
0
])
&&
(
j
<
this
->
n
[
1
])),
"Access out of the matrix! "
<<
"Index ("
<<
i
<<
", "
<<
j
<<
") is out of the matrix of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
")"
);
return
*
(
this
->
values
+
i
+
j
*
this
->
n
[
0
]);
}
inline
const
T
&
at
(
UInt
i
,
UInt
j
)
const
{
AKANTU_DEBUG_ASSERT
(((
i
<
this
->
n
[
0
])
&&
(
j
<
this
->
n
[
1
])),
"Access out of the matrix! "
<<
"Index ("
<<
i
<<
", "
<<
j
<<
") is out of the matrix of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
")"
);
return
*
(
this
->
values
+
i
+
j
*
this
->
n
[
0
]);
}
/* ------------------------------------------------------------------------ */
inline
T
&
operator
()(
UInt
i
,
UInt
j
)
{
return
this
->
at
(
i
,
j
);
}
inline
const
T
&
operator
()(
UInt
i
,
UInt
j
)
const
{
return
this
->
at
(
i
,
j
);
}
/// give a line vector wrapped on the column i
inline
VectorProxy
<
T
>
operator
()(
UInt
j
)
{
AKANTU_DEBUG_ASSERT
(
j
<
this
->
n
[
1
],
"Access out of the matrix! "
<<
"You are trying to access the column vector "
<<
j
<<
" in a matrix of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
")"
);
return
VectorProxy
<
T
>
(
this
->
values
+
j
*
this
->
n
[
0
],
this
->
n
[
0
]);
}
inline
const
VectorProxy
<
T
>
operator
()(
UInt
j
)
const
{
AKANTU_DEBUG_ASSERT
(
j
<
this
->
n
[
1
],
"Access out of the matrix! "
<<
"You are trying to access the column vector "
<<
j
<<
" in a matrix of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
")"
);
return
VectorProxy
<
T
>
(
this
->
values
+
j
*
this
->
n
[
0
],
this
->
n
[
0
]);
}
inline
void
block
(
const
Matrix
&
block
,
UInt
pos_i
,
UInt
pos_j
)
{
AKANTU_DEBUG_ASSERT
(
pos_i
+
block
.
rows
()
<=
rows
(),
"The block size or position are not correct"
);
AKANTU_DEBUG_ASSERT
(
pos_i
+
block
.
cols
()
<=
cols
(),
"The block size or position are not correct"
);
for
(
UInt
i
=
0
;
i
<
block
.
rows
();
++
i
)
for
(
UInt
j
=
0
;
j
<
block
.
cols
();
++
j
)
this
->
at
(
i
+
pos_i
,
j
+
pos_j
)
=
block
(
i
,
j
);
}
inline
Matrix
block
(
UInt
pos_i
,
UInt
pos_j
,
UInt
block_rows
,
UInt
block_cols
)
const
{
AKANTU_DEBUG_ASSERT
(
pos_i
+
block_rows
<=
rows
(),
"The block size or position are not correct"
);
AKANTU_DEBUG_ASSERT
(
pos_i
+
block_cols
<=
cols
(),
"The block size or position are not correct"
);
Matrix
block
(
block_rows
,
block_cols
);
for
(
UInt
i
=
0
;
i
<
block_rows
;
++
i
)
for
(
UInt
j
=
0
;
j
<
block_cols
;
++
j
)
block
(
i
,
j
)
=
this
->
at
(
i
+
pos_i
,
j
+
pos_j
);
return
block
;
}
inline
T
&
operator
[](
UInt
idx
)
{
return
*
(
this
->
values
+
idx
);
};
inline
const
T
&
operator
[](
UInt
idx
)
const
{
return
*
(
this
->
values
+
idx
);
};
/* ---------------------------------------------------------------------- */
inline
Matrix
operator
*
(
const
Matrix
&
B
)
{
Matrix
C
(
this
->
rows
(),
B
.
cols
());
C
.
mul
<
false
,
false
>
(
*
this
,
B
);
return
C
;
}
/* ----------------------------------------------------------------------- */
inline
Matrix
&
operator
*=
(
const
T
&
x
)
{
return
parent
::
operator
*=
(
x
);
}
inline
Matrix
&
operator
*=
(
const
Matrix
&
B
)
{
Matrix
C
(
*
this
);
this
->
mul
<
false
,
false
>
(
C
,
B
);
return
*
this
;
}
/* ---------------------------------------------------------------------- */
template
<
bool
tr_A
,
bool
tr_B
>
inline
void
mul
(
const
Matrix
&
A
,
const
Matrix
&
B
,
T
alpha
=
1.0
)
{
UInt
k
=
A
.
cols
();
if
(
tr_A
)
k
=
A
.
rows
();
#ifndef AKANTU_NDEBUG
if
(
tr_B
)
{
AKANTU_DEBUG_ASSERT
(
k
==
B
.
cols
(),
"matrices to multiply have no fit dimensions"
);
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
B
.
rows
(),
"matrices to multiply have no fit dimensions"
);
}
else
{
AKANTU_DEBUG_ASSERT
(
k
==
B
.
rows
(),
"matrices to multiply have no fit dimensions"
);
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
B
.
cols
(),
"matrices to multiply have no fit dimensions"
);
}
if
(
tr_A
)
{
AKANTU_DEBUG_ASSERT
(
this
->
rows
()
==
A
.
cols
(),
"matrices to multiply have no fit dimensions"
);
}
else
{
AKANTU_DEBUG_ASSERT
(
this
->
rows
()
==
A
.
rows
(),
"matrices to multiply have no fit dimensions"
);
}
#endif
// AKANTU_NDEBUG
Math
::
matMul
<
tr_A
,
tr_B
>
(
this
->
rows
(),
this
->
cols
(),
k
,
alpha
,
A
.
storage
(),
B
.
storage
(),
0.
,
this
->
storage
());
}
/* ---------------------------------------------------------------------- */
inline
void
outerProduct
(
const
Vector
<
T
>
&
A
,
const
Vector
<
T
>
&
B
)
{
AKANTU_DEBUG_ASSERT
(
A
.
size
()
==
this
->
rows
()
&&
B
.
size
()
==
this
->
cols
(),
"A and B are not compatible with the size of the matrix"
);
for
(
UInt
i
=
0
;
i
<
this
->
rows
();
++
i
)
{
for
(
UInt
j
=
0
;
j
<
this
->
cols
();
++
j
)
{
this
->
values
[
i
+
j
*
this
->
rows
()]
+=
A
[
i
]
*
B
[
j
];
}
}
}
private
:
class
EigenSorter
{
public
:
EigenSorter
(
const
Vector
<
T
>
&
eigs
)
:
eigs
(
eigs
)
{}
bool
operator
()(
const
UInt
&
a
,
const
UInt
&
b
)
const
{
return
(
eigs
(
a
)
>
eigs
(
b
));
}
private
:
const
Vector
<
T
>
&
eigs
;
};
public
:
/* ---------------------------------------------------------------------- */
inline
void
eig
(
Vector
<
T
>
&
eigenvalues
,
Matrix
<
T
>
&
eigenvectors
)
const
{
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
this
->
rows
(),
"eig is not a valid operation on a rectangular matrix"
);
AKANTU_DEBUG_ASSERT
(
eigenvalues
.
size
()
==
this
->
cols
(),
"eigenvalues should be of size "
<<
this
->
cols
()
<<
"."
);
#ifndef AKANTU_NDEBUG
if
(
eigenvectors
.
storage
()
!=
nullptr
)
AKANTU_DEBUG_ASSERT
((
eigenvectors
.
rows
()
==
eigenvectors
.
cols
())
&&
(
eigenvectors
.
rows
()
==
this
->
cols
()),
"Eigenvectors needs to be a square matrix of size "
<<
this
->
cols
()
<<
" x "
<<
this
->
cols
()
<<
"."
);
#endif
Matrix
<
T
>
tmp
=
*
this
;
Vector
<
T
>
tmp_eigs
(
eigenvalues
.
size
());
Matrix
<
T
>
tmp_eig_vects
(
eigenvectors
.
rows
(),
eigenvectors
.
cols
());
if
(
tmp_eig_vects
.
rows
()
==
0
||
tmp_eig_vects
.
cols
()
==
0
)
Math
::
matrixEig
(
tmp
.
cols
(),
tmp
.
storage
(),
tmp_eigs
.
storage
());
else
Math
::
matrixEig
(
tmp
.
cols
(),
tmp
.
storage
(),
tmp_eigs
.
storage
(),
tmp_eig_vects
.
storage
());
Vector
<
UInt
>
perm
(
eigenvalues
.
size
());
for
(
UInt
i
=
0
;
i
<
perm
.
size
();
++
i
)
perm
(
i
)
=
i
;
std
::
sort
(
perm
.
storage
(),
perm
.
storage
()
+
perm
.
size
(),
EigenSorter
(
tmp_eigs
));
for
(
UInt
i
=
0
;
i
<
perm
.
size
();
++
i
)
eigenvalues
(
i
)
=
tmp_eigs
(
perm
(
i
));
if
(
tmp_eig_vects
.
rows
()
!=
0
&&
tmp_eig_vects
.
cols
()
!=
0
)
for
(
UInt
i
=
0
;
i
<
perm
.
size
();
++
i
)
{
for
(
UInt
j
=
0
;
j
<
eigenvectors
.
rows
();
++
j
)
{
eigenvectors
(
j
,
i
)
=
tmp_eig_vects
(
j
,
perm
(
i
));
}
}
}
/* ---------------------------------------------------------------------- */
inline
void
eig
(
Vector
<
T
>
&
eigenvalues
)
const
{
Matrix
<
T
>
empty
;
eig
(
eigenvalues
,
empty
);
}
/* ---------------------------------------------------------------------- */
inline
void
eye
(
T
alpha
=
1.
)
{
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
this
->
rows
(),
"eye is not a valid operation on a rectangular matrix"
);
this
->
clear
();
for
(
UInt
i
=
0
;
i
<
this
->
cols
();
++
i
)
{
this
->
values
[
i
+
i
*
this
->
rows
()]
=
alpha
;
}
}
/* ---------------------------------------------------------------------- */
static
inline
Matrix
<
T
>
eye
(
UInt
m
,
T
alpha
=
1.
)
{
Matrix
<
T
>
tmp
(
m
,
m
);
tmp
.
eye
(
alpha
);
return
tmp
;
}
/* ---------------------------------------------------------------------- */
inline
T
trace
()
const
{
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
this
->
rows
(),
"trace is not a valid operation on a rectangular matrix"
);
T
trace
=
0.
;
for
(
UInt
i
=
0
;
i
<
this
->
rows
();
++
i
)
{
trace
+=
this
->
values
[
i
+
i
*
this
->
rows
()];
}
return
trace
;
}
/* ---------------------------------------------------------------------- */
inline
Matrix
transpose
()
const
{
Matrix
tmp
(
this
->
cols
(),
this
->
rows
());
for
(
UInt
i
=
0
;
i
<
this
->
rows
();
++
i
)
{
for
(
UInt
j
=
0
;
j
<
this
->
cols
();
++
j
)
{
tmp
(
j
,
i
)
=
operator
()(
i
,
j
);
}
}
return
tmp
;
}
/* ---------------------------------------------------------------------- */
inline
void
inverse
(
const
Matrix
&
A
)
{
AKANTU_DEBUG_ASSERT
(
A
.
cols
()
==
A
.
rows
(),
"inv is not a valid operation on a rectangular matrix"
);
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
A
.
cols
(),
"the matrix should have the same size as its inverse"
);
if
(
this
->
cols
()
==
1
)
*
this
->
values
=
1.
/
*
A
.
storage
();
else
if
(
this
->
cols
()
==
2
)
Math
::
inv2
(
A
.
storage
(),
this
->
values
);
else
if
(
this
->
cols
()
==
3
)
Math
::
inv3
(
A
.
storage
(),
this
->
values
);
else
Math
::
inv
(
this
->
cols
(),
A
.
storage
(),
this
->
values
);
}
inline
Matrix
inverse
()
{
Matrix
inv
(
this
->
rows
(),
this
->
cols
());
inv
.
inverse
(
*
this
);
return
inv
;
}
/* --------------------------------------------------------------------- */
inline
T
det
()
const
{
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
this
->
rows
(),
"inv is not a valid operation on a rectangular matrix"
);
if
(
this
->
cols
()
==
1
)
return
*
(
this
->
values
);
else
if
(
this
->
cols
()
==
2
)
return
Math
::
det2
(
this
->
values
);
else
if
(
this
->
cols
()
==
3
)
return
Math
::
det3
(
this
->
values
);
else
return
Math
::
det
(
this
->
cols
(),
this
->
values
);
}
/* --------------------------------------------------------------------- */
inline
T
doubleDot
(
const
Matrix
<
T
>
&
other
)
const
{
AKANTU_DEBUG_ASSERT
(
this
->
cols
()
==
this
->
rows
(),
"doubleDot is not a valid operation on a rectangular matrix"
);
if
(
this
->
cols
()
==
1
)
return
*
(
this
->
values
)
*
*
(
other
.
storage
());
else
if
(
this
->
cols
()
==
2
)
return
Math
::
matrixDoubleDot22
(
this
->
values
,
other
.
storage
());
else
if
(
this
->
cols
()
==
3
)
return
Math
::
matrixDoubleDot33
(
this
->
values
,
other
.
storage
());
else
AKANTU_ERROR
(
"doubleDot is not defined for other spatial dimensions"
<<
" than 1, 2 or 3."
);
return
T
();
}
/* ---------------------------------------------------------------------- */
/// function to print the containt of the class
virtual
void
printself
(
std
::
ostream
&
stream
,
int
indent
=
0
)
const
{
std
::
string
space
;
for
(
Int
i
=
0
;
i
<
indent
;
i
++
,
space
+=
AKANTU_INDENT
)
;
stream
<<
"["
;
for
(
UInt
i
=
0
;
i
<
this
->
n
[
0
];
++
i
)
{
if
(
i
!=
0
)
stream
<<
", "
;
stream
<<
"["
;
for
(
UInt
j
=
0
;
j
<
this
->
n
[
1
];
++
j
)
{
if
(
j
!=
0
)
stream
<<
", "
;
stream
<<
operator
()(
i
,
j
);
}
stream
<<
"]"
;
}
stream
<<
"]"
;
};
};
/* ------------------------------------------------------------------------ */
template
<
typename
T
>
template
<
bool
tr_A
>
inline
void
Vector
<
T
>::
mul
(
const
Matrix
<
T
>
&
A
,
const
Vector
<
T
>
&
x
,
Real
alpha
)
{
#ifndef AKANTU_NDEBUG
UInt
n
=
x
.
size
();
if
(
tr_A
)
{
AKANTU_DEBUG_ASSERT
(
n
==
A
.
rows
(),
"matrix and vector to multiply have no fit dimensions"
);
AKANTU_DEBUG_ASSERT
(
this
->
size
()
==
A
.
cols
(),
"matrix and vector to multiply have no fit dimensions"
);
}
else
{
AKANTU_DEBUG_ASSERT
(
n
==
A
.
cols
(),
"matrix and vector to multiply have no fit dimensions"
);
AKANTU_DEBUG_ASSERT
(
this
->
size
()
==
A
.
rows
(),
"matrix and vector to multiply have no fit dimensions"
);
}
#endif
Math
::
matVectMul
<
tr_A
>
(
A
.
rows
(),
A
.
cols
(),
alpha
,
A
.
storage
(),
x
.
storage
(),
0.
,
this
->
storage
());
}
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
inline
std
::
ostream
&
operator
<<
(
std
::
ostream
&
stream
,
const
Matrix
<
T
>
&
_this
)
{
_this
.
printself
(
stream
);
return
stream
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
inline
std
::
ostream
&
operator
<<
(
std
::
ostream
&
stream
,
const
Vector
<
T
>
&
_this
)
{
_this
.
printself
(
stream
);
return
stream
;
}
/* ------------------------------------------------------------------------ */
/* Tensor3 */
/* ------------------------------------------------------------------------ */
template
<
typename
T
>
class
Tensor3
:
public
TensorStorage
<
T
,
3
,
Tensor3
<
T
>>
{
using
parent
=
TensorStorage
<
T
,
3
,
Tensor3
<
T
>>
;
public
:
using
value_type
=
typename
parent
::
value_type
;
using
proxy
=
Tensor3Proxy
<
T
>
;
public
:
Tensor3
()
:
parent
(){};
Tensor3
(
UInt
m
,
UInt
n
,
UInt
p
,
const
T
&
def
=
T
())
:
parent
(
m
,
n
,
p
,
def
)
{}
Tensor3
(
T
*
data
,
UInt
m
,
UInt
n
,
UInt
p
)
:
parent
(
data
,
m
,
n
,
p
)
{}
Tensor3
(
const
Tensor3
&
src
,
bool
deep_copy
=
true
)
:
parent
(
src
,
deep_copy
)
{}
Tensor3
(
const
proxy
&
src
)
:
parent
(
src
)
{}
public
:
/* ------------------------------------------------------------------------ */
inline
Tensor3
&
operator
=
(
const
Tensor3
&
src
)
{
parent
::
operator
=
(
src
);
return
*
this
;
}
/* ---------------------------------------------------------------------- */
inline
T
&
operator
()(
UInt
i
,
UInt
j
,
UInt
k
)
{
AKANTU_DEBUG_ASSERT
(
(
i
<
this
->
n
[
0
])
&&
(
j
<
this
->
n
[
1
])
&&
(
k
<
this
->
n
[
2
]),
"Access out of the tensor3! "
<<
"You are trying to access the element "
<<
"("
<<
i
<<
", "
<<
j
<<
", "
<<
k
<<
") in a tensor of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
", "
<<
this
->
n
[
2
]
<<
")"
);
return
*
(
this
->
values
+
(
k
*
this
->
n
[
0
]
+
i
)
*
this
->
n
[
1
]
+
j
);
}
inline
const
T
&
operator
()(
UInt
i
,
UInt
j
,
UInt
k
)
const
{
AKANTU_DEBUG_ASSERT
(
(
i
<
this
->
n
[
0
])
&&
(
j
<
this
->
n
[
1
])
&&
(
k
<
this
->
n
[
2
]),
"Access out of the tensor3! "
<<
"You are trying to access the element "
<<
"("
<<
i
<<
", "
<<
j
<<
", "
<<
k
<<
") in a tensor of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
", "
<<
this
->
n
[
2
]
<<
")"
);
return
*
(
this
->
values
+
(
k
*
this
->
n
[
0
]
+
i
)
*
this
->
n
[
1
]
+
j
);
}
inline
MatrixProxy
<
T
>
operator
()(
UInt
k
)
{
AKANTU_DEBUG_ASSERT
((
k
<
this
->
n
[
2
]),
"Access out of the tensor3! "
<<
"You are trying to access the slice "
<<
k
<<
" in a tensor3 of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
", "
<<
this
->
n
[
2
]
<<
")"
);
return
MatrixProxy
<
T
>
(
this
->
values
+
k
*
this
->
n
[
0
]
*
this
->
n
[
1
],
this
->
n
[
0
],
this
->
n
[
1
]);
}
inline
const
MatrixProxy
<
T
>
operator
()(
UInt
k
)
const
{
AKANTU_DEBUG_ASSERT
((
k
<
this
->
n
[
2
]),
"Access out of the tensor3! "
<<
"You are trying to access the slice "
<<
k
<<
" in a tensor3 of size ("
<<
this
->
n
[
0
]
<<
", "
<<
this
->
n
[
1
]
<<
", "
<<
this
->
n
[
2
]
<<
")"
);
return
MatrixProxy
<
T
>
(
this
->
values
+
k
*
this
->
n
[
0
]
*
this
->
n
[
1
],
this
->
n
[
0
],
this
->
n
[
1
]);
}
inline
MatrixProxy
<
T
>
operator
[](
UInt
k
)
{
return
MatrixProxy
<
T
>
(
this
->
values
+
k
*
this
->
n
[
0
]
*
this
->
n
[
1
],
this
->
n
[
0
],
this
->
n
[
1
]);
}
inline
const
MatrixProxy
<
T
>
operator
[](
UInt
k
)
const
{
return
MatrixProxy
<
T
>
(
this
->
values
+
k
*
this
->
n
[
0
]
*
this
->
n
[
1
],
this
->
n
[
0
],
this
->
n
[
1
]);
}
};
/* -------------------------------------------------------------------------- */
// support operations for the creation of other vectors
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
Vector
<
T
>
operator
*
(
const
T
&
scalar
,
const
Vector
<
T
>
&
a
)
{
Vector
<
T
>
r
(
a
);
r
*=
scalar
;
return
r
;
}
template
<
typename
T
>
Vector
<
T
>
operator
*
(
const
Vector
<
T
>
&
a
,
const
T
&
scalar
)
{
Vector
<
T
>
r
(
a
);
r
*=
scalar
;
return
r
;
}
template
<
typename
T
>
Vector
<
T
>
operator
/
(
const
Vector
<
T
>
&
a
,
const
T
&
scalar
)
{
Vector
<
T
>
r
(
a
);
r
/=
scalar
;
return
r
;
}
template
<
typename
T
>
Vector
<
T
>
operator
*
(
const
Vector
<
T
>
&
a
,
const
Vector
<
T
>
&
b
)
{
Vector
<
T
>
r
(
a
);
r
*=
b
;
return
r
;
}
template
<
typename
T
>
Vector
<
T
>
operator
+
(
const
Vector
<
T
>
&
a
,
const
Vector
<
T
>
&
b
)
{
Vector
<
T
>
r
(
a
);
r
+=
b
;
return
r
;
}
template
<
typename
T
>
Vector
<
T
>
operator
-
(
const
Vector
<
T
>
&
a
,
const
Vector
<
T
>
&
b
)
{
Vector
<
T
>
r
(
a
);
r
-=
b
;
return
r
;
}
template
<
typename
T
>
Vector
<
T
>
operator
*
(
const
Matrix
<
T
>
&
A
,
const
Vector
<
T
>
&
b
)
{
Vector
<
T
>
r
(
b
.
size
());
r
.
template
mul
<
false
>
(
A
,
b
);
return
r
;
}
/* -------------------------------------------------------------------------- */
template
<
typename
T
>
Matrix
<
T
>
operator
*
(
const
T
&
scalar
,
const
Matrix
<
T
>
&
a
)
{
Matrix
<
T
>
r
(
a
);
r
*=
scalar
;
return
r
;
}
template
<
typename
T
>
Matrix
<
T
>
operator
*
(
const
Matrix
<
T
>
&
a
,
const
T
&
scalar
)
{
Matrix
<
T
>
r
(
a
);
r
*=
scalar
;
return
r
;
}
template
<
typename
T
>
Matrix
<
T
>
operator
/
(
const
Matrix
<
T
>
&
a
,
const
T
&
scalar
)
{
Matrix
<
T
>
r
(
a
);
r
/=
scalar
;
return
r
;
}
template
<
typename
T
>
Matrix
<
T
>
operator
+
(
const
Matrix
<
T
>
&
a
,
const
Matrix
<
T
>
&
b
)
{
Matrix
<
T
>
r
(
a
);
r
+=
b
;
return
r
;
}
template
<
typename
T
>
Matrix
<
T
>
operator
-
(
const
Matrix
<
T
>
&
a
,
const
Matrix
<
T
>
&
b
)
{
Matrix
<
T
>
r
(
a
);
r
-=
b
;
return
r
;
}
}
// namespace akantu
#include <iterator>
namespace
std
{
template
<
typename
R
>
struct
iterator_traits
<::
akantu
::
types
::
details
::
vector_iterator
<
R
>>
{
protected
:
using
iterator
=
::
akantu
::
types
::
details
::
vector_iterator
<
R
>
;
public
:
using
iterator_category
=
typename
iterator
::
iterator_category
;
using
value_type
=
typename
iterator
::
value_type
;
using
difference_type
=
typename
iterator
::
difference_type
;
using
pointer
=
typename
iterator
::
pointer
;
using
reference
=
typename
iterator
::
reference
;
};
}
// namespace std
#endif
/* __AKANTU_AKA_TYPES_HH__ */
Event Timeline
Log In to Comment