mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
371 lines
12 KiB
C++
371 lines
12 KiB
C++
#ifndef GRID_MATH_TENSORS_H
|
|
#define GRID_MATH_TENSORS_H
|
|
|
|
namespace Grid {
|
|
|
|
///////////////////////////////////////////////////
|
|
// Scalar, Vector, Matrix objects.
|
|
// These can be composed to form tensor products of internal indices.
|
|
///////////////////////////////////////////////////
|
|
|
|
// It is useful to NOT have any constructors
|
|
// so that these classes assert "is_pod<class> == true"
|
|
// because then the standard C++ valarray container eliminates fill overhead on new allocation and
|
|
// non-move copying.
|
|
//
|
|
// However note that doing this eliminates some syntactical sugar such as
|
|
// calling the constructor explicitly or implicitly
|
|
//
|
|
class GridTensorBase {};
|
|
|
|
template<class vtype> class iScalar
|
|
{
|
|
public:
|
|
vtype _internal;
|
|
|
|
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
|
|
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
|
|
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
|
|
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
|
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
|
|
typedef iScalar<recurse_scalar_object> scalar_object;
|
|
|
|
// substitutes a real or complex version with same tensor structure
|
|
typedef iScalar<typename GridTypeMapper<vtype>::Complexified > Complexified;
|
|
typedef iScalar<typename GridTypeMapper<vtype>::Realified > Realified;
|
|
|
|
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
|
|
|
|
// Scalar no action
|
|
// template<int Level> using tensor_reduce_level = typename iScalar<GridTypeMapper<vtype>::tensor_reduce_level<Level> >;
|
|
iScalar() = default;
|
|
/*
|
|
iScalar(const iScalar<vtype> ©me)=default;
|
|
iScalar(iScalar<vtype> &©me)=default;
|
|
iScalar<vtype> & operator= (const iScalar<vtype> ©me) = default;
|
|
iScalar<vtype> & operator= (iScalar<vtype> &©me) = default;
|
|
*/
|
|
iScalar(scalar_type s) : _internal(s) {};// recurse down and hit the constructor for vector_type
|
|
iScalar(const Zero &z){ *this = zero; };
|
|
|
|
iScalar<vtype> & operator= (const Zero &hero){
|
|
zeroit(*this);
|
|
return *this;
|
|
}
|
|
friend strong_inline void vstream(iScalar<vtype> &out,const iScalar<vtype> &in){
|
|
vstream(out._internal,in._internal);
|
|
}
|
|
friend strong_inline void zeroit(iScalar<vtype> &that){
|
|
zeroit(that._internal);
|
|
}
|
|
friend strong_inline void prefetch(iScalar<vtype> &that){
|
|
prefetch(that._internal);
|
|
}
|
|
friend strong_inline void permute(iScalar<vtype> &out,const iScalar<vtype> &in,int permutetype){
|
|
permute(out._internal,in._internal,permutetype);
|
|
}
|
|
|
|
// Unary negation
|
|
friend strong_inline iScalar<vtype> operator -(const iScalar<vtype> &r) {
|
|
iScalar<vtype> ret;
|
|
ret._internal= -r._internal;
|
|
return ret;
|
|
}
|
|
// *=,+=,-= operators inherit from corresponding "*,-,+" behaviour
|
|
strong_inline iScalar<vtype> &operator *=(const iScalar<vtype> &r) {
|
|
*this = (*this)*r;
|
|
return *this;
|
|
}
|
|
strong_inline iScalar<vtype> &operator -=(const iScalar<vtype> &r) {
|
|
*this = (*this)-r;
|
|
return *this;
|
|
}
|
|
strong_inline iScalar<vtype> &operator +=(const iScalar<vtype> &r) {
|
|
*this = (*this)+r;
|
|
return *this;
|
|
}
|
|
strong_inline vtype & operator ()(void) {
|
|
return _internal;
|
|
}
|
|
strong_inline const vtype & operator ()(void) const {
|
|
return _internal;
|
|
}
|
|
|
|
// Type casts meta programmed, must be pure scalar to match TensorRemove
|
|
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator ComplexF () const { return(TensorRemove(_internal)); };
|
|
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator ComplexD () const { return(TensorRemove(_internal)); };
|
|
// template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator RealD () const { return(real(TensorRemove(_internal))); }
|
|
template<class U=vtype,class V=scalar_type,IfReal<V> = 0,IfNotSimd<U> = 0> operator RealD () const { return TensorRemove(_internal); }
|
|
template<class U=vtype,class V=scalar_type,IfInteger<V> = 0,IfNotSimd<U> = 0> operator Integer () const { return Integer(TensorRemove(_internal)); }
|
|
|
|
// convert from a something to a scalar via constructor of something arg
|
|
template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline iScalar<vtype> operator = (T arg)
|
|
{
|
|
_internal = arg;
|
|
return *this;
|
|
}
|
|
|
|
friend std::ostream& operator<< (std::ostream& stream, const iScalar<vtype> &o){
|
|
stream<< "S {"<<o._internal<<"}";
|
|
return stream;
|
|
};
|
|
};
|
|
///////////////////////////////////////////////////////////
|
|
// Allows to turn scalar<scalar<scalar<double>>>> back to double.
|
|
///////////////////////////////////////////////////////////
|
|
template<class T> strong_inline typename std::enable_if<!isGridTensor<T>::value, T>::type TensorRemove(T arg) { return arg;}
|
|
template<class vtype> strong_inline auto TensorRemove(iScalar<vtype> arg) -> decltype(TensorRemove(arg._internal))
|
|
{
|
|
return TensorRemove(arg._internal);
|
|
}
|
|
|
|
template<class vtype,int N> class iVector
|
|
{
|
|
public:
|
|
vtype _internal[N];
|
|
|
|
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
|
|
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
|
|
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
|
|
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
|
|
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
|
typedef iVector<recurse_scalar_object,N> scalar_object;
|
|
|
|
// substitutes a real or complex version with same tensor structure
|
|
typedef iVector<typename GridTypeMapper<vtype>::Complexified,N > Complexified;
|
|
typedef iVector<typename GridTypeMapper<vtype>::Realified,N > Realified;
|
|
|
|
template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline auto operator = (T arg) -> iVector<vtype,N>
|
|
{
|
|
zeroit(*this);
|
|
for(int i=0;i<N;i++)
|
|
_internal[i] = arg;
|
|
return *this;
|
|
}
|
|
|
|
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
|
|
iVector(const Zero &z){ *this = zero; };
|
|
iVector() =default;
|
|
/*
|
|
iVector(const iVector<vtype,N> ©me)=default;
|
|
iVector(iVector<vtype,N> &©me)=default;
|
|
iVector<vtype,N> & operator= (const iVector<vtype,N> ©me) = default;
|
|
iVector<vtype,N> & operator= (iVector<vtype,N> &©me) = default;
|
|
*/
|
|
|
|
iVector<vtype,N> & operator= (const Zero &hero){
|
|
zeroit(*this);
|
|
return *this;
|
|
}
|
|
friend strong_inline void zeroit(iVector<vtype,N> &that){
|
|
for(int i=0;i<N;i++){
|
|
zeroit(that._internal[i]);
|
|
}
|
|
}
|
|
friend strong_inline void prefetch(iVector<vtype,N> &that){
|
|
for(int i=0;i<N;i++) prefetch(that._internal[i]);
|
|
}
|
|
friend strong_inline void vstream(iVector<vtype,N> &out,const iVector<vtype,N> &in){
|
|
for(int i=0;i<N;i++){
|
|
vstream(out._internal[i],in._internal[i]);
|
|
}
|
|
}
|
|
friend strong_inline void permute(iVector<vtype,N> &out,const iVector<vtype,N> &in,int permutetype){
|
|
for(int i=0;i<N;i++){
|
|
permute(out._internal[i],in._internal[i],permutetype);
|
|
}
|
|
}
|
|
// Unary negation
|
|
friend strong_inline iVector<vtype,N> operator -(const iVector<vtype,N> &r) {
|
|
iVector<vtype,N> ret;
|
|
for(int i=0;i<N;i++) ret._internal[i]= -r._internal[i];
|
|
return ret;
|
|
}
|
|
// *=,+=,-= operators inherit from corresponding "*,-,+" behaviour
|
|
strong_inline iVector<vtype,N> &operator *=(const iScalar<vtype> &r) {
|
|
*this = (*this)*r;
|
|
return *this;
|
|
}
|
|
strong_inline iVector<vtype,N> &operator -=(const iVector<vtype,N> &r) {
|
|
*this = (*this)-r;
|
|
return *this;
|
|
}
|
|
strong_inline iVector<vtype,N> &operator +=(const iVector<vtype,N> &r) {
|
|
*this = (*this)+r;
|
|
return *this;
|
|
}
|
|
strong_inline vtype & operator ()(int i) {
|
|
return _internal[i];
|
|
}
|
|
strong_inline const vtype & operator ()(int i) const {
|
|
return _internal[i];
|
|
}
|
|
friend std::ostream& operator<< (std::ostream& stream, const iVector<vtype,N> &o){
|
|
stream<< "V<"<<N<<">{";
|
|
for(int i=0;i<N;i++) {
|
|
stream<<o._internal[i];
|
|
if (i<N-1) stream<<",";
|
|
}
|
|
stream<<"}";
|
|
return stream;
|
|
};
|
|
// strong_inline vtype && operator ()(int i) {
|
|
// return _internal[i];
|
|
// }
|
|
};
|
|
|
|
template<class vtype,int N> class iMatrix
|
|
{
|
|
public:
|
|
vtype _internal[N][N];
|
|
|
|
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
|
|
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
|
|
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
|
|
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
|
|
|
|
// substitutes a real or complex version with same tensor structure
|
|
typedef iMatrix<typename GridTypeMapper<vtype>::Complexified,N > Complexified;
|
|
typedef iMatrix<typename GridTypeMapper<vtype>::Realified,N > Realified;
|
|
|
|
// Tensure removal
|
|
typedef iScalar<tensor_reduced_v> tensor_reduced;
|
|
typedef iMatrix<recurse_scalar_object,N> scalar_object;
|
|
|
|
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
|
|
|
|
|
|
iMatrix(const Zero &z){ *this = zero; };
|
|
iMatrix() =default;
|
|
|
|
iMatrix& operator=(const iMatrix& rhs){
|
|
for(int i=0;i<N;i++)
|
|
for(int j=0;j<N;j++)
|
|
vstream(_internal[i][j],rhs._internal[i][j]);
|
|
return *this;
|
|
};
|
|
|
|
|
|
|
|
iMatrix(scalar_type s) { (*this) = s ;};// recurse down and hit the constructor for vector_type
|
|
|
|
/*
|
|
iMatrix(const iMatrix<vtype,N> ©me)=default;
|
|
iMatrix(iMatrix<vtype,N> &©me)=default;
|
|
iMatrix<vtype,N> & operator= (const iMatrix<vtype,N> ©me) = default;
|
|
iMatrix<vtype,N> & operator= (iMatrix<vtype,N> &©me) = default;
|
|
*/
|
|
|
|
|
|
|
|
iMatrix<vtype,N> & operator= (const Zero &hero){
|
|
zeroit(*this);
|
|
return *this;
|
|
}
|
|
template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline auto operator = (T arg) -> iMatrix<vtype,N>
|
|
{
|
|
zeroit(*this);
|
|
for(int i=0;i<N;i++)
|
|
_internal[i][i] = arg;
|
|
return *this;
|
|
}
|
|
|
|
friend strong_inline void zeroit(iMatrix<vtype,N> &that){
|
|
for(int i=0;i<N;i++){
|
|
for(int j=0;j<N;j++){
|
|
zeroit(that._internal[i][j]);
|
|
}}
|
|
}
|
|
friend strong_inline void prefetch(iMatrix<vtype,N> &that){
|
|
for(int i=0;i<N;i++)
|
|
for(int j=0;j<N;j++)
|
|
prefetch(that._internal[i][j]);
|
|
}
|
|
friend strong_inline void vstream(iMatrix<vtype,N> &out,const iMatrix<vtype,N> &in){
|
|
for(int i=0;i<N;i++){
|
|
for(int j=0;j<N;j++){
|
|
vstream(out._internal[i][j],in._internal[i][j]);
|
|
}}
|
|
}
|
|
friend strong_inline void permute(iMatrix<vtype,N> &out,const iMatrix<vtype,N> &in,int permutetype){
|
|
for(int i=0;i<N;i++){
|
|
for(int j=0;j<N;j++){
|
|
permute(out._internal[i][j],in._internal[i][j],permutetype);
|
|
}}
|
|
}
|
|
// Unary negation
|
|
friend strong_inline iMatrix<vtype,N> operator -(const iMatrix<vtype,N> &r) {
|
|
iMatrix<vtype,N> ret;
|
|
for(int i=0;i<N;i++){
|
|
for(int j=0;j<N;j++){
|
|
ret._internal[i][j]= -r._internal[i][j];
|
|
}}
|
|
return ret;
|
|
}
|
|
// *=,+=,-= operators inherit from corresponding "*,-,+" behaviour
|
|
template<class T>
|
|
strong_inline iMatrix<vtype,N> &operator *=(const T &r) {
|
|
*this = (*this)*r;
|
|
return *this;
|
|
}
|
|
template<class T>
|
|
strong_inline iMatrix<vtype,N> &operator -=(const T &r) {
|
|
*this = (*this)-r;
|
|
return *this;
|
|
}
|
|
template<class T>
|
|
strong_inline iMatrix<vtype,N> &operator +=(const T &r) {
|
|
*this = (*this)+r;
|
|
return *this;
|
|
}
|
|
|
|
// returns an lvalue reference
|
|
strong_inline vtype & operator ()(int i,int j) {
|
|
return _internal[i][j];
|
|
}
|
|
strong_inline const vtype & operator ()(int i,int j) const {
|
|
return _internal[i][j];
|
|
}
|
|
friend std::ostream& operator<< (std::ostream& stream, const iMatrix<vtype,N> &o){
|
|
stream<< "M<"<<N<<">{";
|
|
for(int i=0;i<N;i++) {
|
|
stream<< "{";
|
|
for(int j=0;j<N;j++) {
|
|
stream<<o._internal[i][j];
|
|
if (i<N-1) stream<<",";
|
|
}
|
|
stream<<"}";
|
|
if(i!=N-1) stream<<"\n\t\t";
|
|
}
|
|
stream<<"}";
|
|
return stream;
|
|
};
|
|
|
|
// strong_inline vtype && operator ()(int i,int j) {
|
|
// return _internal[i][j];
|
|
// }
|
|
|
|
};
|
|
|
|
template<class v> void vprefetch(const iScalar<v> &vv)
|
|
{
|
|
vprefetch(vv._internal);
|
|
}
|
|
template<class v,int N> void vprefetch(const iVector<v,N> &vv)
|
|
{
|
|
for(int i=0;i<N;i++){
|
|
vprefetch(vv._internal[i]);
|
|
}
|
|
}
|
|
template<class v,int N> void vprefetch(const iMatrix<v,N> &vv)
|
|
{
|
|
for(int i=0;i<N;i++){
|
|
for(int j=0;j<N;j++){
|
|
vprefetch(vv._internal[i][j]);
|
|
}}
|
|
}
|
|
|
|
|
|
}
|
|
#endif
|