update dcel to avoid merge conflict
This commit is contained in:
James E McClure
commit
a49915ccbc
@ -35,10 +35,6 @@ IF ( NOT CMAKE_CXX_STANDARD )
|
||||
ENDIF()
|
||||
|
||||
|
||||
# Disable Link Time Optimization until we have a chance to test
|
||||
SET( DISABLE_LTO )
|
||||
|
||||
|
||||
# Set source/install paths
|
||||
SET( ${PROJ}_SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}" )
|
||||
SET( ${PROJ}_BUILD_DIR "${CMAKE_CURRENT_BINARY_DIR}" )
|
||||
|
||||
@ -11,6 +11,9 @@
|
||||
#include "IO/Reader.h"
|
||||
#include "IO/Writer.h"
|
||||
|
||||
#include "ProfilerApp.h"
|
||||
|
||||
|
||||
#define PI 3.14159265359
|
||||
|
||||
// Constructor
|
||||
@ -55,11 +58,11 @@ double Minkowski::X(){
|
||||
return Xi_global;
|
||||
}
|
||||
|
||||
void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
|
||||
void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue)
|
||||
{
|
||||
PROFILE_START("ComputeScalar");
|
||||
Xi = Ji = Ai = 0.0;
|
||||
DECL object;
|
||||
Point P1,P2,P3;
|
||||
int e1,e2,e3;
|
||||
double s,s1,s2,s3;
|
||||
double a1,a2,a3;
|
||||
@ -75,13 +78,13 @@ void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
|
||||
e1 = object.Face(idx);
|
||||
e2 = object.halfedge.next(e1);
|
||||
e3 = object.halfedge.next(e2);
|
||||
P1 = object.vertex.coords(object.halfedge.v1(e1));
|
||||
P2 = object.vertex.coords(object.halfedge.v1(e2));
|
||||
P3 = object.vertex.coords(object.halfedge.v1(e3));
|
||||
auto P1 = object.vertex.coords(object.halfedge.v1(e1));
|
||||
auto P2 = object.vertex.coords(object.halfedge.v1(e2));
|
||||
auto P3 = object.vertex.coords(object.halfedge.v1(e3));
|
||||
// Surface area
|
||||
s1 = sqrt((P1.x-P2.x)*(P1.x-P2.x)+(P1.y-P2.y)*(P1.y-P2.y)+(P1.z-P2.z)*(P1.z-P2.z));
|
||||
s2 = sqrt((P2.x-P3.x)*(P2.x-P3.x)+(P2.y-P3.y)*(P2.y-P3.y)+(P2.z-P3.z)*(P2.z-P3.z));
|
||||
s3 = sqrt((P1.x-P3.x)*(P1.x-P3.x)+(P1.y-P3.y)*(P1.y-P3.y)+(P1.z-P3.z)*(P1.z-P3.z));
|
||||
s1 = Distance( P1, P2 );
|
||||
s2 = Distance( P2, P3 );
|
||||
s3 = Distance( P1, P3 );
|
||||
s = 0.5*(s1+s2+s3);
|
||||
Ai += sqrt(s*(s-s1)*(s-s2)*(s-s3));
|
||||
// Mean curvature based on half edge angle
|
||||
@ -122,7 +125,7 @@ void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
|
||||
MPI_Allreduce(&Ai,&Ai_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
|
||||
MPI_Allreduce(&Ji,&Ji_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
|
||||
MPI_Barrier(Dm->Comm);
|
||||
|
||||
PROFILE_STOP("ComputeScalar");
|
||||
}
|
||||
|
||||
void Minkowski::PrintAll()
|
||||
|
||||
@ -48,7 +48,7 @@ public:
|
||||
//...........................................................................
|
||||
Minkowski(std::shared_ptr <Domain> Dm);
|
||||
~Minkowski();
|
||||
void ComputeScalar(const DoubleArray Field, const double isovalue);
|
||||
void ComputeScalar(const DoubleArray& Field, const double isovalue);
|
||||
void PrintAll();
|
||||
|
||||
};
|
||||
|
||||
@ -1,5 +1,7 @@
|
||||
#include "analysis/analysis.h"
|
||||
#include "ProfilerApp.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
|
||||
@ -76,10 +76,6 @@ int Halfedge::size(){
|
||||
|
||||
|
||||
DECL::DECL(){
|
||||
cellvertices=DTMutableList<Point>(16);
|
||||
FaceData.resize(16);
|
||||
Triangles.resize(3,16);
|
||||
halfedge.data.resize(6,48);
|
||||
}
|
||||
|
||||
DECL::~DECL(){
|
||||
@ -99,7 +95,17 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const
|
||||
|
||||
int CubeIndex;
|
||||
int nTris = 0;
|
||||
int nVert = 0;
|
||||
int nVert =0;
|
||||
|
||||
Point VertexList[12];
|
||||
Point NewVertexList[12];
|
||||
int LocalRemap[12];
|
||||
|
||||
DTMutableList<Point> cellvertices = DTMutableList<Point>(20);
|
||||
IntArray Triangles = IntArray(3,20);
|
||||
|
||||
// Values from array 'A' at the cube corners
|
||||
double CubeValues[8];
|
||||
|
||||
// Points corresponding to cube corners
|
||||
C0.x = 0.0; C0.y = 0.0; C0.z = 0.0;
|
||||
@ -231,8 +237,10 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const
|
||||
|
||||
// Now add the local values to the DECL data structure
|
||||
if (nTris>0){
|
||||
//FaceData.resize(TriangleCount);
|
||||
//halfedge.data.resize(6,nTris*3);
|
||||
FaceData.resize(TriangleCount);
|
||||
//printf("Construct halfedge structure... \n");
|
||||
//printf(" Construct %i triangles \n",nTris);
|
||||
halfedge.data.resize(6,nTris*3);
|
||||
int idx_edge=0;
|
||||
for (int idx=0; idx<TriangleCount; idx++){
|
||||
int V1 = Triangles(0,idx);
|
||||
|
||||
@ -46,11 +46,12 @@ public:
|
||||
DECL();
|
||||
~DECL();
|
||||
|
||||
int Face();
|
||||
int face();
|
||||
Vertex vertex;
|
||||
Halfedge halfedge;
|
||||
void LocalIsosurface(const DoubleArray A, double value, int i, int j, int k);
|
||||
void LocalIsosurface(const DoubleArray A, double value, int i, int j, int k);
|
||||
int Face(int index);
|
||||
|
||||
double origin(int edge);
|
||||
double EdgeAngle(int edge);
|
||||
Point TriNormal(int edge);
|
||||
@ -59,13 +60,4 @@ public:
|
||||
|
||||
private:
|
||||
Array <int> FaceData;
|
||||
|
||||
Point VertexList[12];
|
||||
Point NewVertexList[12];
|
||||
int LocalRemap[12];
|
||||
double CubeValues[8];
|
||||
|
||||
DTMutableList<Point> cellvertices;// = DTMutableList<Point>(20);
|
||||
IntArray Triangles;// = IntArray(3,20);
|
||||
|
||||
};
|
||||
|
||||
@ -1,6 +1,7 @@
|
||||
#include "analysis/distance.h"
|
||||
|
||||
|
||||
|
||||
/******************************************************************
|
||||
* A fast distance calculation *
|
||||
******************************************************************/
|
||||
|
||||
@ -2,6 +2,7 @@
|
||||
#define Distance_H_INC
|
||||
|
||||
#include "common/Domain.h"
|
||||
#include "common/Array.hpp"
|
||||
|
||||
|
||||
struct Vec {
|
||||
|
||||
117
common/Array.cpp
Normal file
117
common/Array.cpp
Normal file
@ -0,0 +1,117 @@
|
||||
#include "common/Array.h"
|
||||
#include "common/Array.hpp"
|
||||
|
||||
#include <complex>
|
||||
|
||||
|
||||
/********************************************************
|
||||
* ArraySize *
|
||||
********************************************************/
|
||||
ArraySize::ArraySize( const std::vector<size_t>& N )
|
||||
{
|
||||
d_ndim = N.size();
|
||||
d_N[0] = 0;
|
||||
d_N[1] = 1;
|
||||
d_N[2] = 1;
|
||||
d_N[3] = 1;
|
||||
d_N[4] = 1;
|
||||
for ( size_t i = 0; i < d_ndim; i++ )
|
||||
d_N[i] = N[i];
|
||||
d_length = 1;
|
||||
for ( unsigned long i : d_N )
|
||||
d_length *= i;
|
||||
if ( d_ndim == 0 )
|
||||
d_length = 0;
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Explicit instantiations of Array *
|
||||
********************************************************/
|
||||
template class Array<char, FunctionTable>;
|
||||
template class Array<uint8_t, FunctionTable>;
|
||||
template class Array<uint16_t, FunctionTable>;
|
||||
template class Array<uint32_t, FunctionTable>;
|
||||
template class Array<uint64_t, FunctionTable>;
|
||||
template class Array<int8_t, FunctionTable>;
|
||||
template class Array<int16_t, FunctionTable>;
|
||||
template class Array<int32_t, FunctionTable>;
|
||||
template class Array<int64_t, FunctionTable>;
|
||||
template class Array<float, FunctionTable>;
|
||||
template class Array<double, FunctionTable>;
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Explicit instantiations of Array<bool> *
|
||||
********************************************************/
|
||||
// clang-format off
|
||||
template Array<bool, FunctionTable>::Array();
|
||||
template Array<bool, FunctionTable>::~Array();
|
||||
template Array<bool, FunctionTable>::Array( size_t );
|
||||
template Array<bool, FunctionTable>::Array( size_t, size_t );
|
||||
template Array<bool, FunctionTable>::Array( size_t, size_t, size_t );
|
||||
template Array<bool, FunctionTable>::Array( size_t, size_t, size_t, size_t );
|
||||
template Array<bool, FunctionTable>::Array( size_t, size_t, size_t, size_t, size_t );
|
||||
template Array<bool, FunctionTable>::Array( const std::vector<size_t>&, const bool* );
|
||||
template Array<bool, FunctionTable>::Array( std::string );
|
||||
template Array<bool, FunctionTable>::Array( std::initializer_list<bool> );
|
||||
template Array<bool, FunctionTable>::Array( const Array<bool, FunctionTable>& );
|
||||
template Array<bool, FunctionTable>::Array( Array<bool, FunctionTable>&& );
|
||||
template Array<bool, FunctionTable>& Array<bool, FunctionTable>::operator=( const Array<bool, FunctionTable>& );
|
||||
template Array<bool, FunctionTable>& Array<bool, FunctionTable>::operator=( Array<bool, FunctionTable>&& );
|
||||
template Array<bool, FunctionTable>& Array<bool, FunctionTable>::operator=( const std::vector<bool>& );
|
||||
template void Array<bool, FunctionTable>::fill(bool const&);
|
||||
template void Array<bool, FunctionTable>::clear();
|
||||
template bool Array<bool, FunctionTable>::operator==(Array<bool, FunctionTable> const&) const;
|
||||
template void Array<bool, FunctionTable>::resize( ArraySize const& );
|
||||
// clang-format on
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Explicit instantiations of Array<std::complex> *
|
||||
********************************************************/
|
||||
// clang-format off
|
||||
template Array<std::complex<double>, FunctionTable>::Array();
|
||||
template Array<std::complex<double>, FunctionTable>::~Array();
|
||||
template Array<std::complex<double>, FunctionTable>::Array( size_t );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( size_t, size_t );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( size_t, size_t, size_t );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( size_t, size_t, size_t, size_t );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( size_t, size_t, size_t, size_t, size_t );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( const std::vector<size_t>&, const std::complex<double>* );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( std::initializer_list<std::complex<double>> );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( const Range<std::complex<double>>& range );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( const Array<std::complex<double>, FunctionTable>& );
|
||||
template Array<std::complex<double>, FunctionTable>::Array( Array<std::complex<double>, FunctionTable>&& );
|
||||
template Array<std::complex<double>, FunctionTable>& Array<std::complex<double>, FunctionTable>::operator=( const Array<std::complex<double>, FunctionTable>& );
|
||||
template Array<std::complex<double>, FunctionTable>& Array<std::complex<double>, FunctionTable>::operator=( Array<std::complex<double>, FunctionTable>&& );
|
||||
template Array<std::complex<double>, FunctionTable>& Array<std::complex<double>, FunctionTable>::operator=( const std::vector<std::complex<double>>& );
|
||||
template void Array<std::complex<double>, FunctionTable>::resize( ArraySize const& );
|
||||
template void Array<std::complex<double>, FunctionTable>::viewRaw( ArraySize const&, std::complex<double>*, bool, bool );
|
||||
template void Array<std::complex<double>, FunctionTable>::fill(std::complex<double> const&);
|
||||
template void Array<std::complex<double>, FunctionTable>::clear();
|
||||
template bool Array<std::complex<double>, FunctionTable>::operator==(Array<std::complex<double>, FunctionTable> const&) const;
|
||||
template Array<std::complex<double>, FunctionTable> Array<std::complex<double>, FunctionTable>::repmat(std::vector<unsigned long, std::allocator<unsigned long> > const&) const;
|
||||
// clang-format on
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Explicit instantiations of Array<std::string> *
|
||||
********************************************************/
|
||||
// clang-format off
|
||||
template Array<std::string, FunctionTable>::Array();
|
||||
template Array<std::string, FunctionTable>::~Array();
|
||||
template Array<std::string, FunctionTable>::Array( size_t );
|
||||
template Array<std::string, FunctionTable>::Array( size_t, size_t );
|
||||
template Array<std::string, FunctionTable>::Array( size_t, size_t, size_t );
|
||||
template Array<std::string, FunctionTable>::Array( size_t, size_t, size_t, size_t );
|
||||
template Array<std::string, FunctionTable>::Array( size_t, size_t, size_t, size_t, size_t );
|
||||
template Array<std::string, FunctionTable>::Array( const std::vector<size_t>&, const std::string* );
|
||||
template Array<std::string, FunctionTable>::Array( std::initializer_list<std::string> );
|
||||
template Array<std::string, FunctionTable>::Array( const Array<std::string, FunctionTable>& );
|
||||
template Array<std::string, FunctionTable>::Array( Array<std::string, FunctionTable>&& );
|
||||
template Array<std::string, FunctionTable>& Array<std::string, FunctionTable>::operator=( const Array<std::string, FunctionTable>& );
|
||||
template Array<std::string, FunctionTable>& Array<std::string, FunctionTable>::operator=( Array<std::string, FunctionTable>&& );
|
||||
template Array<std::string, FunctionTable>& Array<std::string, FunctionTable>::operator=( const std::vector<std::string>& );
|
||||
template void Array<std::string, FunctionTable>::resize( ArraySize const& );
|
||||
// clang-format on
|
||||
610
common/Array.h
610
common/Array.h
@ -1,247 +1,20 @@
|
||||
#ifndef included_ArrayClass
|
||||
#define included_ArrayClass
|
||||
|
||||
#include "common/ArraySize.h"
|
||||
|
||||
#include <array>
|
||||
#include <cstring>
|
||||
#include <functional>
|
||||
#include <initializer_list>
|
||||
#include <iostream>
|
||||
#include <memory>
|
||||
#include <vector>
|
||||
|
||||
#include "Utilities.h"
|
||||
|
||||
|
||||
#if defined( __CUDA_ARCH__ )
|
||||
#include <cuda.h>
|
||||
#define HOST_DEVICE __host__ __device__
|
||||
#else
|
||||
#define HOST_DEVICE
|
||||
#endif
|
||||
#if defined( USING_GCC ) || defined( USING_CLANG )
|
||||
#define ATTRIBUTE_INLINE __attribute__( ( always_inline ) )
|
||||
#else
|
||||
#define ATTRIBUTE_INLINE
|
||||
#endif
|
||||
|
||||
|
||||
#if ( defined( DEBUG ) || defined( _DEBUG ) ) && !defined( NDEBUG )
|
||||
#define CHECK_ARRAY_LENGTH( i ) \
|
||||
do { \
|
||||
if ( i >= d_length ) \
|
||||
throw std::length_error( "Index exceeds array bounds" ); \
|
||||
} while ( 0 )
|
||||
#else
|
||||
#define CHECK_ARRAY_LENGTH( i ) \
|
||||
do { \
|
||||
} while ( 0 )
|
||||
#endif
|
||||
|
||||
|
||||
// Forward decleration
|
||||
class FunctionTable;
|
||||
|
||||
|
||||
//! Simple range class
|
||||
template<class TYPE = size_t>
|
||||
class Range final
|
||||
{
|
||||
public:
|
||||
//! Empty constructor
|
||||
Range() : i( 0 ), j( -1 ), k( 1 ) {}
|
||||
|
||||
/*!
|
||||
* Create a range i:k:j (or i:j)
|
||||
* @param i_ Starting value
|
||||
* @param j_ Ending value
|
||||
* @param k_ Increment value
|
||||
*/
|
||||
Range( TYPE i_, TYPE j_, TYPE k_ = 1 ) : i( i_ ), j( j_ ), k( k_ ) {}
|
||||
|
||||
TYPE i, j, k;
|
||||
};
|
||||
|
||||
|
||||
//! Simple class to store the array dimensions
|
||||
class ArraySize final
|
||||
{
|
||||
public:
|
||||
//! Empty constructor
|
||||
inline ArraySize();
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
*/
|
||||
inline ArraySize( size_t N1 );
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
*/
|
||||
inline ArraySize( size_t N1, size_t N2 );
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
* @param N3 Number of elements in the third dimension
|
||||
*/
|
||||
inline ArraySize( size_t N1, size_t N2, size_t N3 );
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
* @param N3 Number of elements in the third dimension
|
||||
* @param N4 Number of elements in the fourth dimension
|
||||
*/
|
||||
inline ArraySize( size_t N1, size_t N2, size_t N3, size_t N4 );
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
* @param N3 Number of elements in the third dimension
|
||||
* @param N4 Number of elements in the fourth dimension
|
||||
* @param N5 Number of elements in the fifth dimension
|
||||
*/
|
||||
inline ArraySize( size_t N1, size_t N2, size_t N3, size_t N4, size_t N5 );
|
||||
|
||||
/*!
|
||||
* Create from initializer list
|
||||
* @param N Size of the array
|
||||
*/
|
||||
inline ArraySize( std::initializer_list<size_t> N );
|
||||
|
||||
/*!
|
||||
* Create from raw pointer
|
||||
* @param ndim Number of dimensions
|
||||
* @param ndim Dimensions
|
||||
*/
|
||||
inline ArraySize( size_t ndim, const size_t *dims );
|
||||
|
||||
/*!
|
||||
* Create from std::vector
|
||||
* @param N Size of the array
|
||||
*/
|
||||
inline ArraySize( const std::vector<size_t> &N );
|
||||
|
||||
/*!
|
||||
* Copy constructor
|
||||
* @param rhs Array to copy
|
||||
*/
|
||||
inline ArraySize( const ArraySize &rhs );
|
||||
|
||||
/*!
|
||||
* Move constructor
|
||||
* @param rhs Array to copy
|
||||
*/
|
||||
inline ArraySize( ArraySize &&rhs );
|
||||
|
||||
/*!
|
||||
* Assignment operator
|
||||
* @param rhs Array to copy
|
||||
*/
|
||||
inline ArraySize &operator=( const ArraySize &rhs );
|
||||
|
||||
/*!
|
||||
* Move assignment operator
|
||||
* @param rhs Array to copy
|
||||
*/
|
||||
inline ArraySize &operator=( ArraySize &&rhs );
|
||||
|
||||
/*!
|
||||
* Access the ith dimension
|
||||
* @param i Index to access
|
||||
*/
|
||||
inline size_t operator[]( size_t i ) const { return d_N[i]; }
|
||||
|
||||
//! Sum the elements
|
||||
inline uint8_t ndim() const ATTRIBUTE_INLINE { return d_ndim; }
|
||||
|
||||
//! Sum the elements
|
||||
inline size_t size() const ATTRIBUTE_INLINE { return d_ndim; }
|
||||
|
||||
//! Sum the elements
|
||||
inline size_t length() const ATTRIBUTE_INLINE { return d_length; }
|
||||
|
||||
//! Sum the elements
|
||||
inline void resize( uint8_t dim, size_t N );
|
||||
|
||||
//! Returns an iterator to the beginning
|
||||
inline const size_t *begin() const ATTRIBUTE_INLINE { return d_N; }
|
||||
|
||||
//! Returns an iterator to the end
|
||||
inline const size_t *end() const ATTRIBUTE_INLINE { return d_N + d_ndim; }
|
||||
|
||||
// Check if two matrices are equal
|
||||
inline bool operator==( const ArraySize &rhs ) const ATTRIBUTE_INLINE
|
||||
{
|
||||
return d_ndim == rhs.d_ndim && memcmp( d_N, rhs.d_N, sizeof( d_N ) ) == 0;
|
||||
}
|
||||
|
||||
//! Check if two matrices are not equal
|
||||
inline bool operator!=( const ArraySize &rhs ) const ATTRIBUTE_INLINE
|
||||
{
|
||||
return d_ndim != rhs.d_ndim || memcmp( d_N, rhs.d_N, sizeof( d_N ) ) != 0;
|
||||
}
|
||||
|
||||
//! Maximum supported dimension
|
||||
constexpr static uint8_t maxDim() ATTRIBUTE_INLINE { return 5u; }
|
||||
|
||||
//! Get the index
|
||||
inline size_t index( size_t i ) const ATTRIBUTE_INLINE
|
||||
{
|
||||
CHECK_ARRAY_LENGTH( i );
|
||||
return i;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
inline size_t index( size_t i1, size_t i2 ) const ATTRIBUTE_INLINE
|
||||
{
|
||||
size_t index = i1 + i2 * d_N[0];
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
inline size_t index( size_t i1, size_t i2, size_t i3 ) const ATTRIBUTE_INLINE
|
||||
{
|
||||
size_t index = i1 + d_N[0] * ( i2 + d_N[1] * i3 );
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
inline size_t index( size_t i1, size_t i2, size_t i3, size_t i4 ) const ATTRIBUTE_INLINE
|
||||
{
|
||||
size_t index = i1 + d_N[0] * ( i2 + d_N[1] * ( i3 + d_N[2] * i4 ) );
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
inline size_t index(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const ATTRIBUTE_INLINE
|
||||
{
|
||||
size_t index = i1 + d_N[0] * ( i2 + d_N[1] * ( i3 + d_N[2] * ( i4 + d_N[3] * i5 ) ) );
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
private:
|
||||
uint8_t d_ndim;
|
||||
size_t d_length;
|
||||
size_t d_N[5];
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
* Class Array is a multi-dimensional array class written by Mark Berrill
|
||||
*/
|
||||
template<class TYPE, class FUN = FunctionTable>
|
||||
template<class TYPE, class FUN, class Allocator>
|
||||
class Array final
|
||||
{
|
||||
public: // Constructors / assignment operators
|
||||
@ -309,6 +82,12 @@ public: // Constructors / assignment operators
|
||||
*/
|
||||
explicit Array( const Range<TYPE> &range );
|
||||
|
||||
/*!
|
||||
* Create a 1D Array using a string that mimic's MATLAB
|
||||
* @param range Range of the data
|
||||
*/
|
||||
explicit Array( std::string range );
|
||||
|
||||
/*!
|
||||
* Create a 1D Array with the given initializer list
|
||||
* @param data Input data
|
||||
@ -346,74 +125,34 @@ public: // Constructors / assignment operators
|
||||
*/
|
||||
Array &operator=( const std::vector<TYPE> &rhs );
|
||||
|
||||
//! Is copyable?
|
||||
inline bool isCopyable() const { return d_isCopyable; }
|
||||
|
||||
//! Set is copyable
|
||||
inline void setCopyable( bool flag ) { d_isCopyable = flag; }
|
||||
|
||||
//! Is fixed size?
|
||||
inline bool isFixedSize() const { return d_isFixedSize; }
|
||||
|
||||
//! Set is copyable
|
||||
inline void setFixedSize( bool flag ) { d_isFixedSize = flag; }
|
||||
|
||||
|
||||
public: // Views/copies/subset
|
||||
/*!
|
||||
* Create a 1D Array view to a raw block of data
|
||||
* @param N Number of elements in the array
|
||||
* Create a multi-dimensional Array view to a raw block of data
|
||||
* @param N Number of elements in each dimension
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<Array> view( size_t N, std::shared_ptr<TYPE> const &data );
|
||||
static std::unique_ptr<Array> view( const ArraySize &N, std::shared_ptr<TYPE> &data );
|
||||
|
||||
/*!
|
||||
* Create a new 2D Array with the given number of rows and columns
|
||||
* @param N_rows Number of rows
|
||||
* @param N_columns Number of columns
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<Array> view(
|
||||
size_t N_rows, size_t N_columns, std::shared_ptr<TYPE> const &data );
|
||||
|
||||
/*!
|
||||
* Create a new 3D Array view to a raw block of data
|
||||
* @param N1 Number of rows
|
||||
* @param N2 Number of columns
|
||||
* @param N3 Number of elements in the third dimension
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<Array> view(
|
||||
size_t N1, size_t N2, size_t N3, std::shared_ptr<TYPE> const &data );
|
||||
|
||||
/*!
|
||||
* Create a multi-dimensional Array view to a raw block of data
|
||||
* @param N Number of elements in each dimension
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<Array> view( const ArraySize &N, std::shared_ptr<TYPE> const &data );
|
||||
|
||||
|
||||
/*!
|
||||
* Create a 1D Array view to a raw block of data
|
||||
* @param N Number of elements in the array
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<const Array> constView(
|
||||
size_t N, std::shared_ptr<const TYPE> const &data );
|
||||
|
||||
/*!
|
||||
* Create a new 2D Array with the given number of rows and columns
|
||||
* @param N_rows Number of rows
|
||||
* @param N_columns Number of columns
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<const Array> constView(
|
||||
size_t N_rows, size_t N_columns, std::shared_ptr<const TYPE> const &data );
|
||||
|
||||
/*!
|
||||
* Create a new 3D Array view to a raw block of data
|
||||
* @param N1 Number of rows
|
||||
* @param N2 Number of columns
|
||||
* @param N3 Number of elements in the third dimension
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<const Array> constView(
|
||||
size_t N1, size_t N2, size_t N3, std::shared_ptr<const TYPE> const &data );
|
||||
|
||||
/*!
|
||||
* Create a multi-dimensional Array view to a raw block of data
|
||||
* @param N Number of elements in each dimension
|
||||
* @param data Pointer to the data
|
||||
*/
|
||||
static std::shared_ptr<const Array> constView(
|
||||
static std::unique_ptr<const Array> constView(
|
||||
const ArraySize &N, std::shared_ptr<const TYPE> const &data );
|
||||
|
||||
|
||||
@ -432,7 +171,7 @@ public: // Views/copies/subset
|
||||
|
||||
/*!
|
||||
* Make this object a view of the raw data (expert use only).
|
||||
* Use view2( N, std::shared_ptr(data,[](TYPE*){}) ) instead.
|
||||
* Use view2( N, shared_ptr(data,[](TYPE*){}) ) instead.
|
||||
* Note: this interface is not recommended as it does not protect from
|
||||
* the src data being deleted while still being used by the Array.
|
||||
* Additionally for maximum performance it does not set the internal shared_ptr
|
||||
@ -440,27 +179,28 @@ public: // Views/copies/subset
|
||||
* @param ndim Number of dimensions
|
||||
* @param dims Number of elements in each dimension
|
||||
* @param data Pointer to the data
|
||||
* @param isCopyable Once the view is created, can the array be copied
|
||||
* @param isFixedSize Once the view is created, is the array size fixed
|
||||
*/
|
||||
void viewRaw( int ndim, const size_t *dims, TYPE *data );
|
||||
inline void viewRaw(
|
||||
int ndim, const size_t *dims, TYPE *data, bool isCopyable = true, bool isFixedSize = true )
|
||||
{
|
||||
viewRaw( ArraySize( ndim, dims ), data, isCopyable, isFixedSize );
|
||||
}
|
||||
|
||||
/*!
|
||||
* Make this object a view of the raw data (expert use only).
|
||||
* Use view2( N, std::shared_ptr(data,[](TYPE*){}) ) instead.
|
||||
* Use view2( N, shared_ptr(data,[](TYPE*){}) ) instead.
|
||||
* Note: this interface is not recommended as it does not protect from
|
||||
* the src data being deleted while still being used by the Array.
|
||||
* Additionally for maximum performance it does not set the internal shared_ptr
|
||||
* so functions like getPtr and resize will not work correctly.
|
||||
* @param N Number of elements in each dimension
|
||||
* @param data Pointer to the data
|
||||
* @param isCopyable Once the view is created, can the array be copied
|
||||
* @param isFixedSize Once the view is created, is the array size fixed
|
||||
*/
|
||||
void viewRaw( const ArraySize &N, TYPE *data );
|
||||
|
||||
/*!
|
||||
* Convert an array of one type to another. This may or may not allocate new memory.
|
||||
* @param array Input array
|
||||
*/
|
||||
template<class TYPE2>
|
||||
static std::shared_ptr<Array<TYPE2>> convert( std::shared_ptr<Array<TYPE, FUN>> array );
|
||||
void viewRaw( const ArraySize &N, TYPE *data, bool isCopyable = true, bool isFixedSize = true );
|
||||
|
||||
|
||||
/*!
|
||||
@ -468,8 +208,27 @@ public: // Views/copies/subset
|
||||
* @param array Input array
|
||||
*/
|
||||
template<class TYPE2>
|
||||
static std::shared_ptr<const Array<TYPE2>> convert(
|
||||
std::shared_ptr<const Array<TYPE, FUN>> array );
|
||||
static inline std::unique_ptr<Array<TYPE2, FUN, Allocator>> convert(
|
||||
std::shared_ptr<Array<TYPE, FUN, Allocator>> array )
|
||||
{
|
||||
auto array2 = std::make_unique<Array<TYPE2>>( array->size() );
|
||||
array2.copy( *array );
|
||||
return array2;
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* Convert an array of one type to another. This may or may not allocate new memory.
|
||||
* @param array Input array
|
||||
*/
|
||||
template<class TYPE2>
|
||||
static inline std::unique_ptr<const Array<TYPE2, FUN, Allocator>> convert(
|
||||
std::shared_ptr<const Array<TYPE, FUN, Allocator>> array )
|
||||
{
|
||||
auto array2 = std::make_unique<Array<TYPE2>>( array->size() );
|
||||
array2.copy( *array );
|
||||
return array2;
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
@ -477,7 +236,11 @@ public: // Views/copies/subset
|
||||
* @param array Source array
|
||||
*/
|
||||
template<class TYPE2>
|
||||
void copy( const Array<TYPE2> &array );
|
||||
void inline copy( const Array<TYPE2, FUN, Allocator> &array )
|
||||
{
|
||||
resize( array.size() );
|
||||
copy( array.data() );
|
||||
}
|
||||
|
||||
/*!
|
||||
* Copy and convert data from a raw vector to this array.
|
||||
@ -485,22 +248,36 @@ public: // Views/copies/subset
|
||||
* @param array Source array
|
||||
*/
|
||||
template<class TYPE2>
|
||||
void copy( const TYPE2 *array );
|
||||
void inline copy( const TYPE2 *data )
|
||||
{
|
||||
for ( size_t i = 0; i < d_size.length(); i++ )
|
||||
d_data[i] = static_cast<TYPE>( data[i] );
|
||||
}
|
||||
|
||||
/*!
|
||||
* Copy and convert data from this array to a raw vector.
|
||||
* @param array Source array
|
||||
*/
|
||||
template<class TYPE2>
|
||||
void copyTo( TYPE2 *array ) const;
|
||||
void inline copyTo( TYPE2 *data ) const
|
||||
{
|
||||
for ( size_t i = 0; i < d_size.length(); i++ )
|
||||
data[i] = static_cast<TYPE2>( d_data[i] );
|
||||
}
|
||||
|
||||
/*!
|
||||
* Copy and convert data from this array to a raw vector.
|
||||
* @param array Source array
|
||||
* Copy and convert data from this array to a new array
|
||||
*/
|
||||
template<class TYPE2>
|
||||
Array<TYPE2, FUN> cloneTo() const;
|
||||
Array<TYPE2, FUN, Allocator> inline cloneTo() const
|
||||
{
|
||||
Array<TYPE2, FUN> dst( this->size() );
|
||||
copyTo( dst.data() );
|
||||
return dst;
|
||||
}
|
||||
|
||||
/*! swap the raw data pointers for the Arrays after checking for compatibility */
|
||||
void swap( Array &other );
|
||||
|
||||
/*!
|
||||
* Fill the array with the given value
|
||||
@ -519,7 +296,7 @@ public: // Views/copies/subset
|
||||
* @param base Base array
|
||||
* @param exp Exponent value
|
||||
*/
|
||||
void pow( const Array<TYPE, FUN> &base, const TYPE &exp );
|
||||
void pow( const Array &base, const TYPE &exp );
|
||||
|
||||
//! Destructor
|
||||
~Array();
|
||||
@ -534,11 +311,7 @@ public: // Views/copies/subset
|
||||
|
||||
|
||||
//! Return the size of the Array
|
||||
inline ArraySize &size() { return d_size; }
|
||||
|
||||
|
||||
//! Return the size of the Array
|
||||
inline ArraySize size() const { return d_size; }
|
||||
inline const ArraySize &size() const { return d_size; }
|
||||
|
||||
|
||||
//! Return the size of the Array
|
||||
@ -557,14 +330,15 @@ public: // Views/copies/subset
|
||||
* Resize the Array
|
||||
* @param N NUmber of elements
|
||||
*/
|
||||
void resize( size_t N );
|
||||
inline void resize( size_t N ) { resize( ArraySize( N ) ); }
|
||||
|
||||
|
||||
/*!
|
||||
* Resize the Array
|
||||
* @param N_rows Number of rows
|
||||
* @param N_columns Number of columns
|
||||
* @param N_row Number of rows
|
||||
* @param N_col Number of columns
|
||||
*/
|
||||
void resize( size_t N_rows, size_t N_columns );
|
||||
inline void resize( size_t N_row, size_t N_col ) { resize( ArraySize( N_row, N_col ) ); }
|
||||
|
||||
/*!
|
||||
* Resize the Array
|
||||
@ -572,7 +346,7 @@ public: // Views/copies/subset
|
||||
* @param N2 Number of columns
|
||||
* @param N3 Number of elements in the third dimension
|
||||
*/
|
||||
void resize( size_t N1, size_t N2, size_t N3 );
|
||||
inline void resize( size_t N1, size_t N2, size_t N3 ) { resize( ArraySize( N1, N2, N3 ) ); }
|
||||
|
||||
/*!
|
||||
* Resize the Array
|
||||
@ -598,19 +372,25 @@ public: // Views/copies/subset
|
||||
|
||||
|
||||
/*!
|
||||
* Subset the Array (total size of array will not change)
|
||||
* Reshape the Array so that the number of dimensions is the
|
||||
* max of ndim and the largest dim>1.
|
||||
* @param ndim Desired number of dimensions
|
||||
*/
|
||||
inline void setNdim( int ndim ) { d_size.setNdim( ndim ); }
|
||||
|
||||
|
||||
/*!
|
||||
* Subset the Array
|
||||
* @param index Index to subset (imin,imax,jmin,jmax,kmin,kmax,...)
|
||||
*/
|
||||
template<class TYPE2 = TYPE>
|
||||
Array<TYPE2, FUN> subset( const std::vector<size_t> &index ) const;
|
||||
Array subset( const std::vector<size_t> &index ) const;
|
||||
|
||||
|
||||
/*!
|
||||
* Subset the Array (total size of array will not change)
|
||||
* Subset the Array
|
||||
* @param index Index to subset (ix:kx:jx,iy:ky:jy,...)
|
||||
*/
|
||||
template<class TYPE2 = TYPE>
|
||||
Array<TYPE2, FUN> subset( const std::vector<Range<size_t>> &index ) const;
|
||||
Array subset( const std::vector<Range<size_t>> &index ) const;
|
||||
|
||||
|
||||
/*!
|
||||
@ -618,30 +398,28 @@ public: // Views/copies/subset
|
||||
* @param index Index of the subset (imin,imax,jmin,jmax,kmin,kmax,...)
|
||||
* @param subset The subset array to copy from
|
||||
*/
|
||||
template<class TYPE2>
|
||||
void copySubset( const std::vector<size_t> &index, const Array<TYPE2, FUN> &subset );
|
||||
void copySubset( const std::vector<size_t> &index, const Array &subset );
|
||||
|
||||
/*!
|
||||
* Copy data from an array into a subset of this array
|
||||
* @param index Index of the subset
|
||||
* @param subset The subset array to copy from
|
||||
*/
|
||||
template<class TYPE2>
|
||||
void copySubset( const std::vector<Range<size_t>> &index, const Array<TYPE2, FUN> &subset );
|
||||
void copySubset( const std::vector<Range<size_t>> &index, const Array &subset );
|
||||
|
||||
/*!
|
||||
* Add data from an array into a subset of this array
|
||||
* @param index Index of the subset (imin,imax,jmin,jmax,kmin,kmax,...)
|
||||
* @param subset The subset array to add from
|
||||
*/
|
||||
void addSubset( const std::vector<size_t> &index, const Array<TYPE, FUN> &subset );
|
||||
void addSubset( const std::vector<size_t> &index, const Array &subset );
|
||||
|
||||
/*!
|
||||
* Add data from an array into a subset of this array
|
||||
* @param index Index of the subset
|
||||
* @param subset The subset array to add from
|
||||
*/
|
||||
void addSubset( const std::vector<Range<size_t>> &index, const Array<TYPE, FUN> &subset );
|
||||
void addSubset( const std::vector<Range<size_t>> &index, const Array &subset );
|
||||
|
||||
|
||||
public: // Accessors
|
||||
@ -649,16 +427,13 @@ public: // Accessors
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
*/
|
||||
HOST_DEVICE inline TYPE &operator()( size_t i ) ATTRIBUTE_INLINE
|
||||
{
|
||||
return d_data[d_size.index( i )];
|
||||
}
|
||||
ARRAY_ATTRIBUTE inline TYPE &operator()( size_t i ) { return d_data[d_size.index( i )]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
*/
|
||||
HOST_DEVICE inline const TYPE &operator()( size_t i ) const ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline const TYPE &operator()( size_t i ) const
|
||||
{
|
||||
return d_data[d_size.index( i )];
|
||||
}
|
||||
@ -668,7 +443,7 @@ public: // Accessors
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
HOST_DEVICE inline TYPE &operator()( size_t i, size_t j ) ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline TYPE &operator()( size_t i, size_t j )
|
||||
{
|
||||
return d_data[d_size.index( i, j )];
|
||||
}
|
||||
@ -678,7 +453,7 @@ public: // Accessors
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
HOST_DEVICE inline const TYPE &operator()( size_t i, size_t j ) const ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline const TYPE &operator()( size_t i, size_t j ) const
|
||||
{
|
||||
return d_data[d_size.index( i, j )];
|
||||
}
|
||||
@ -689,7 +464,7 @@ public: // Accessors
|
||||
* @param j The column index
|
||||
* @param k The third index
|
||||
*/
|
||||
HOST_DEVICE inline TYPE &operator()( size_t i, size_t j, size_t k ) ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline TYPE &operator()( size_t i, size_t j, size_t k )
|
||||
{
|
||||
return d_data[d_size.index( i, j, k )];
|
||||
}
|
||||
@ -700,7 +475,7 @@ public: // Accessors
|
||||
* @param j The column index
|
||||
* @param k The third index
|
||||
*/
|
||||
HOST_DEVICE inline const TYPE &operator()( size_t i, size_t j, size_t k ) const ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline const TYPE &operator()( size_t i, size_t j, size_t k ) const
|
||||
{
|
||||
return d_data[d_size.index( i, j, k )];
|
||||
}
|
||||
@ -712,8 +487,7 @@ public: // Accessors
|
||||
* @param i3 The third index
|
||||
* @param i4 The fourth index
|
||||
*/
|
||||
HOST_DEVICE inline TYPE &operator()(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4 ) ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline TYPE &operator()( size_t i1, size_t i2, size_t i3, size_t i4 )
|
||||
{
|
||||
return d_data[d_size.index( i1, i2, i3, i4 )];
|
||||
}
|
||||
@ -725,8 +499,8 @@ public: // Accessors
|
||||
* @param i3 The third index
|
||||
* @param i4 The fourth index
|
||||
*/
|
||||
HOST_DEVICE inline const TYPE &operator()(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4 ) const ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline const TYPE &operator()(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4 ) const
|
||||
{
|
||||
return d_data[d_size.index( i1, i2, i3, i4 )];
|
||||
}
|
||||
@ -739,8 +513,7 @@ public: // Accessors
|
||||
* @param i4 The fourth index
|
||||
* @param i5 The fifth index
|
||||
*/
|
||||
HOST_DEVICE inline TYPE &operator()(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline TYPE &operator()( size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 )
|
||||
{
|
||||
return d_data[d_size.index( i1, i2, i3, i4, i5 )];
|
||||
}
|
||||
@ -753,8 +526,8 @@ public: // Accessors
|
||||
* @param i4 The fourth index
|
||||
* @param i5 The fifth index
|
||||
*/
|
||||
HOST_DEVICE inline const TYPE &operator()(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline const TYPE &operator()(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const
|
||||
{
|
||||
return d_data[d_size.index( i1, i2, i3, i4, i5 )];
|
||||
}
|
||||
@ -763,7 +536,7 @@ public: // Accessors
|
||||
* Access the desired element as a raw pointer
|
||||
* @param i The global index
|
||||
*/
|
||||
HOST_DEVICE inline TYPE *ptr( size_t i ) ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline TYPE *ptr( size_t i )
|
||||
{
|
||||
return i >= d_size.length() ? nullptr : &d_data[i];
|
||||
}
|
||||
@ -772,34 +545,34 @@ public: // Accessors
|
||||
* Access the desired element as a raw pointer
|
||||
* @param i The global index
|
||||
*/
|
||||
HOST_DEVICE inline const TYPE *ptr( size_t i ) const ATTRIBUTE_INLINE
|
||||
ARRAY_ATTRIBUTE inline const TYPE *ptr( size_t i ) const
|
||||
{
|
||||
return i >= d_size.length() ? nullptr : &d_data[i];
|
||||
}
|
||||
|
||||
//! Get iterator to beginning of data
|
||||
inline TYPE *begin() ATTRIBUTE_INLINE { return d_data; }
|
||||
inline TYPE *begin() { return d_data; }
|
||||
|
||||
//! Get iterator to beginning of data
|
||||
inline const TYPE *begin() const ATTRIBUTE_INLINE { return d_data; }
|
||||
inline const TYPE *begin() const { return d_data; }
|
||||
|
||||
//! Get iterator to beginning of data
|
||||
inline TYPE *end() ATTRIBUTE_INLINE { return d_data + d_size.length(); }
|
||||
inline TYPE *end() { return d_data + d_size.length(); }
|
||||
|
||||
//! Get iterator to beginning of data
|
||||
inline const TYPE *end() const ATTRIBUTE_INLINE { return d_data + d_size.length(); }
|
||||
inline const TYPE *end() const { return d_data + d_size.length(); }
|
||||
|
||||
//! Return the pointer to the raw data
|
||||
inline std::shared_ptr<TYPE> getPtr() ATTRIBUTE_INLINE { return d_ptr; }
|
||||
inline std::shared_ptr<TYPE> getPtr() { return d_ptr; }
|
||||
|
||||
//! Return the pointer to the raw data
|
||||
inline std::shared_ptr<const TYPE> getPtr() const ATTRIBUTE_INLINE { return d_ptr; }
|
||||
inline std::shared_ptr<const TYPE> getPtr() const { return d_ptr; }
|
||||
|
||||
//! Return the pointer to the raw data
|
||||
HOST_DEVICE inline TYPE *data() ATTRIBUTE_INLINE { return d_data; }
|
||||
ARRAY_ATTRIBUTE inline TYPE *data() { return d_data; }
|
||||
|
||||
//! Return the pointer to the raw data
|
||||
HOST_DEVICE inline const TYPE *data() const ATTRIBUTE_INLINE { return d_data; }
|
||||
ARRAY_ATTRIBUTE inline const TYPE *data() const { return d_data; }
|
||||
|
||||
|
||||
public: // Operator overloading
|
||||
@ -834,52 +607,52 @@ public: // Math operations
|
||||
void rand();
|
||||
|
||||
//! Return true if NaNs are present
|
||||
inline bool NaNs() const;
|
||||
bool NaNs() const;
|
||||
|
||||
//! Return the smallest value
|
||||
inline TYPE min() const;
|
||||
TYPE min() const;
|
||||
|
||||
//! Return the largest value
|
||||
inline TYPE max() const;
|
||||
TYPE max() const;
|
||||
|
||||
//! Return the sum of all elements
|
||||
inline TYPE sum() const;
|
||||
TYPE sum() const;
|
||||
|
||||
//! Return the mean of all elements
|
||||
inline TYPE mean() const;
|
||||
TYPE mean() const;
|
||||
|
||||
//! Return the min of all elements in a given direction
|
||||
Array<TYPE, FUN> min( int dir ) const;
|
||||
Array min( int dir ) const;
|
||||
|
||||
//! Return the max of all elements in a given direction
|
||||
Array<TYPE, FUN> max( int dir ) const;
|
||||
Array max( int dir ) const;
|
||||
|
||||
//! Return the sum of all elements in a given direction
|
||||
Array<TYPE, FUN> sum( int dir ) const;
|
||||
Array sum( int dir ) const;
|
||||
|
||||
//! Return the smallest value
|
||||
inline TYPE min( const std::vector<size_t> &index ) const;
|
||||
TYPE min( const std::vector<size_t> &index ) const;
|
||||
|
||||
//! Return the largest value
|
||||
inline TYPE max( const std::vector<size_t> &index ) const;
|
||||
TYPE max( const std::vector<size_t> &index ) const;
|
||||
|
||||
//! Return the sum of all elements
|
||||
inline TYPE sum( const std::vector<size_t> &index ) const;
|
||||
TYPE sum( const std::vector<size_t> &index ) const;
|
||||
|
||||
//! Return the mean of all elements
|
||||
inline TYPE mean( const std::vector<size_t> &index ) const;
|
||||
TYPE mean( const std::vector<size_t> &index ) const;
|
||||
|
||||
//! Return the smallest value
|
||||
inline TYPE min( const std::vector<Range<size_t>> &index ) const;
|
||||
TYPE min( const std::vector<Range<size_t>> &index ) const;
|
||||
|
||||
//! Return the largest value
|
||||
inline TYPE max( const std::vector<Range<size_t>> &index ) const;
|
||||
TYPE max( const std::vector<Range<size_t>> &index ) const;
|
||||
|
||||
//! Return the sum of all elements
|
||||
inline TYPE sum( const std::vector<Range<size_t>> &index ) const;
|
||||
TYPE sum( const std::vector<Range<size_t>> &index ) const;
|
||||
|
||||
//! Return the mean of all elements
|
||||
inline TYPE mean( const std::vector<Range<size_t>> &index ) const;
|
||||
TYPE mean( const std::vector<Range<size_t>> &index ) const;
|
||||
|
||||
//! Find all elements that match the operator
|
||||
std::vector<size_t> find(
|
||||
@ -894,17 +667,17 @@ public: // Math operations
|
||||
static Array multiply( const Array &a, const Array &b );
|
||||
|
||||
//! Transpose an array
|
||||
Array<TYPE, FUN> reverseDim() const;
|
||||
Array reverseDim() const;
|
||||
|
||||
//! Replicate an array a given number of times in each direction
|
||||
Array<TYPE, FUN> repmat( const std::vector<size_t> &N ) const;
|
||||
Array repmat( const std::vector<size_t> &N ) const;
|
||||
|
||||
//! Coarsen an array using the given filter
|
||||
Array<TYPE, FUN> coarsen( const Array<TYPE, FUN> &filter ) const;
|
||||
Array coarsen( const Array &filter ) const;
|
||||
|
||||
//! Coarsen an array using the given filter
|
||||
Array<TYPE, FUN> coarsen( const std::vector<size_t> &ratio,
|
||||
std::function<TYPE( const Array<TYPE, FUN> & )> filter ) const;
|
||||
Array coarsen(
|
||||
const std::vector<size_t> &ratio, std::function<TYPE( const Array & )> filter ) const;
|
||||
|
||||
/*!
|
||||
* Perform a element-wise operation y = f(x)
|
||||
@ -928,21 +701,31 @@ public: // Math operations
|
||||
* @param[in] x x
|
||||
* @param[in] beta beta
|
||||
*/
|
||||
void axpby( const TYPE &alpha, const Array<TYPE, FUN> &x, const TYPE &beta );
|
||||
void axpby( const TYPE &alpha, const Array &x, const TYPE &beta );
|
||||
|
||||
/*!
|
||||
* Linear interpolation
|
||||
* @param[in] x Position as a decimal index
|
||||
*/
|
||||
TYPE interp( const std::vector<double> &x ) const;
|
||||
|
||||
/**
|
||||
* \fn equals (Array & const rhs, TYPE tol )
|
||||
* \brief Determine if two Arrays are equal using an absolute tolerance
|
||||
* \param[in] rhs Vector to compare to
|
||||
* \param[in] tol Tolerance of comparison
|
||||
* \return True iff \f$||\mathit{rhs} - x||_\infty < \mathit{tol}\f$
|
||||
*/
|
||||
bool equals( const Array &rhs, TYPE tol = 0.000001 ) const;
|
||||
|
||||
private:
|
||||
bool d_isCopyable; // Can the array be copied
|
||||
bool d_isFixedSize; // Can the array be resized
|
||||
ArraySize d_size; // Size of each dimension
|
||||
TYPE *d_data; // Raw pointer to data in array
|
||||
std::shared_ptr<TYPE> d_ptr; // Shared pointer to data in array
|
||||
void allocate( const ArraySize &N );
|
||||
|
||||
public:
|
||||
template<class TYPE2, class FUN2>
|
||||
inline bool sizeMatch( const Array<TYPE2, FUN2> &rhs ) const
|
||||
{
|
||||
return d_size == rhs.d_size;
|
||||
}
|
||||
|
||||
private:
|
||||
inline void checkSubsetIndex( const std::vector<Range<size_t>> &range ) const;
|
||||
inline std::vector<Range<size_t>> convert( const std::vector<size_t> &index ) const;
|
||||
@ -952,11 +735,60 @@ private:
|
||||
};
|
||||
|
||||
|
||||
/********************************************************
|
||||
* ostream operator *
|
||||
********************************************************/
|
||||
inline std::ostream &operator<<( std::ostream &out, const ArraySize &s )
|
||||
{
|
||||
out << "[" << s[0];
|
||||
for ( size_t i = 1; i < s.ndim(); i++ )
|
||||
out << "," << s[i];
|
||||
out << "]";
|
||||
return out;
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Math operations *
|
||||
********************************************************/
|
||||
template<class TYPE, class FUN, class Allocator>
|
||||
inline Array<TYPE, FUN, Allocator> operator+(
|
||||
const Array<TYPE, FUN, Allocator> &a, const Array<TYPE, FUN, Allocator> &b )
|
||||
{
|
||||
Array<TYPE, FUN, Allocator> c;
|
||||
const auto &fun = []( const TYPE &a, const TYPE &b ) { return a + b; };
|
||||
FUN::transform( fun, a, b, c );
|
||||
return c;
|
||||
}
|
||||
template<class TYPE, class FUN, class Allocator>
|
||||
inline Array<TYPE, FUN, Allocator> operator-(
|
||||
const Array<TYPE, FUN, Allocator> &a, const Array<TYPE, FUN, Allocator> &b )
|
||||
{
|
||||
Array<TYPE, FUN, Allocator> c;
|
||||
const auto &fun = []( const TYPE &a, const TYPE &b ) { return a - b; };
|
||||
FUN::transform( fun, a, b, c );
|
||||
return c;
|
||||
}
|
||||
template<class TYPE, class FUN, class Allocator>
|
||||
inline Array<TYPE, FUN, Allocator> operator*(
|
||||
const Array<TYPE, FUN, Allocator> &a, const Array<TYPE, FUN, Allocator> &b )
|
||||
{
|
||||
return Array<TYPE, FUN, Allocator>::multiply( a, b );
|
||||
}
|
||||
template<class TYPE, class FUN, class Allocator>
|
||||
inline Array<TYPE, FUN, Allocator> operator*(
|
||||
const Array<TYPE, FUN, Allocator> &a, const std::vector<TYPE> &b )
|
||||
{
|
||||
Array<TYPE, FUN, Allocator> b2;
|
||||
b2.viewRaw( { b.size() }, const_cast<TYPE *>( b.data() ) );
|
||||
return Array<TYPE, FUN, Allocator>::multiply( a, b2 );
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Convience typedefs *
|
||||
********************************************************/
|
||||
typedef Array<double> DoubleArray;
|
||||
typedef Array<float> FloatArray;
|
||||
typedef Array<int> IntArray;
|
||||
|
||||
|
||||
#include "common/Array.hpp"
|
||||
|
||||
#endif
|
||||
|
||||
967
common/Array.hpp
967
common/Array.hpp
File diff suppressed because it is too large
Load Diff
294
common/ArraySize.h
Normal file
294
common/ArraySize.h
Normal file
@ -0,0 +1,294 @@
|
||||
#ifndef included_ArraySizeClass
|
||||
#define included_ArraySizeClass
|
||||
|
||||
|
||||
#include <array>
|
||||
#include <cstring>
|
||||
#include <initializer_list>
|
||||
#include <vector>
|
||||
|
||||
|
||||
#if defined( __CUDA_ARCH__ )
|
||||
#include <cuda.h>
|
||||
#define HOST_DEVICE __host__ __device__
|
||||
#else
|
||||
#define HOST_DEVICE
|
||||
#endif
|
||||
#if defined( USING_GCC ) || defined( USING_CLANG )
|
||||
#define ARRAY_ATTRIBUTE HOST_DEVICE __attribute__( ( always_inline ) )
|
||||
#else
|
||||
#define ARRAY_ATTRIBUTE HOST_DEVICE
|
||||
#endif
|
||||
|
||||
|
||||
#if ( defined( DEBUG ) || defined( _DEBUG ) ) && !defined( NDEBUG )
|
||||
#define CHECK_ARRAY_LENGTH( i ) \
|
||||
do { \
|
||||
if ( i >= d_length ) \
|
||||
throw std::out_of_range( "Index exceeds array bounds" ); \
|
||||
} while ( 0 )
|
||||
#else
|
||||
#define CHECK_ARRAY_LENGTH( i ) \
|
||||
do { \
|
||||
} while ( 0 )
|
||||
#endif
|
||||
|
||||
|
||||
// Forward declerations
|
||||
class FunctionTable;
|
||||
template<class TYPE, class FUN = FunctionTable, class Allocator = std::nullptr_t>
|
||||
class Array;
|
||||
|
||||
|
||||
//! Simple range class
|
||||
template<class TYPE = size_t>
|
||||
class Range final
|
||||
{
|
||||
public:
|
||||
//! Empty constructor
|
||||
constexpr Range() : i( 0 ), j( -1 ), k( 1 ) {}
|
||||
|
||||
/*!
|
||||
* Create a range i:k:j (or i:j)
|
||||
* @param i_ Starting value
|
||||
* @param j_ Ending value
|
||||
* @param k_ Increment value
|
||||
*/
|
||||
constexpr Range( TYPE i_, TYPE j_, TYPE k_ = 1 ) : i( i_ ), j( j_ ), k( k_ ) {}
|
||||
|
||||
TYPE i, j, k;
|
||||
};
|
||||
|
||||
|
||||
//! Simple class to store the array dimensions
|
||||
class ArraySize final
|
||||
{
|
||||
public:
|
||||
//! Empty constructor
|
||||
constexpr ArraySize() : d_ndim( 1 ), d_length( 0 ), d_N{ 0, 1, 1, 1, 1 } {}
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
*/
|
||||
constexpr ArraySize( size_t N1 ) : d_ndim( 1 ), d_length( N1 ), d_N{ N1, 1, 1, 1, 1 } {}
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
*/
|
||||
constexpr ArraySize( size_t N1, size_t N2 )
|
||||
: d_ndim( 2 ), d_length( N1 * N2 ), d_N{ N1, N2, 1, 1, 1 }
|
||||
{
|
||||
}
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
* @param N3 Number of elements in the third dimension
|
||||
*/
|
||||
constexpr ArraySize( size_t N1, size_t N2, size_t N3 )
|
||||
: d_ndim( 3 ), d_length( N1 * N2 * N3 ), d_N{ N1, N2, N3, 1, 1 }
|
||||
{
|
||||
}
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
* @param N3 Number of elements in the third dimension
|
||||
* @param N4 Number of elements in the fourth dimension
|
||||
*/
|
||||
constexpr ArraySize( size_t N1, size_t N2, size_t N3, size_t N4 )
|
||||
: d_ndim( 4 ), d_length( N1 * N2 * N3 * N4 ), d_N{ N1, N2, N3, N4, 1 }
|
||||
{
|
||||
}
|
||||
|
||||
/*!
|
||||
* Create the vector size
|
||||
* @param N1 Number of elements in the first dimension
|
||||
* @param N2 Number of elements in the second dimension
|
||||
* @param N3 Number of elements in the third dimension
|
||||
* @param N4 Number of elements in the fourth dimension
|
||||
* @param N5 Number of elements in the fifth dimension
|
||||
*/
|
||||
constexpr ArraySize( size_t N1, size_t N2, size_t N3, size_t N4, size_t N5 )
|
||||
: d_ndim( 5 ), d_length( N1 * N2 * N3 * N4 * N5 ), d_N{ N1, N2, N3, N4, N5 }
|
||||
{
|
||||
}
|
||||
|
||||
/*!
|
||||
* Create from initializer list
|
||||
* @param N Size of the array
|
||||
*/
|
||||
constexpr ArraySize( std::initializer_list<size_t> N )
|
||||
: d_ndim( N.size() ), d_length( 0 ), d_N{ 0, 1, 1, 1, 1 }
|
||||
{
|
||||
if ( d_ndim > maxDim() )
|
||||
throw std::out_of_range( "Maximum number of dimensions exceeded" );
|
||||
auto it = N.begin();
|
||||
for ( size_t i = 0; i < d_ndim; i++, ++it )
|
||||
d_N[i] = *it;
|
||||
d_length = 1;
|
||||
for ( unsigned long i : d_N )
|
||||
d_length *= i;
|
||||
if ( d_ndim == 0 )
|
||||
d_length = 0;
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* Create from raw pointer
|
||||
* @param ndim Number of dimensions
|
||||
* @param dims Dimensions
|
||||
*/
|
||||
constexpr ArraySize( size_t ndim, const size_t *dims )
|
||||
: d_ndim( ndim ), d_length( 0 ), d_N{ 0, 1, 1, 1, 1 }
|
||||
{
|
||||
if ( d_ndim > maxDim() )
|
||||
throw std::out_of_range( "Maximum number of dimensions exceeded" );
|
||||
for ( size_t i = 0; i < ndim; i++ )
|
||||
d_N[i] = dims[i];
|
||||
d_length = 1;
|
||||
for ( unsigned long i : d_N )
|
||||
d_length *= i;
|
||||
if ( d_ndim == 0 )
|
||||
d_length = 0;
|
||||
}
|
||||
|
||||
/*!
|
||||
* Create from std::vector
|
||||
* @param N Size of the array
|
||||
*/
|
||||
ArraySize( const std::vector<size_t> &N );
|
||||
|
||||
// Copy/assignment constructors
|
||||
constexpr ArraySize( ArraySize &&rhs ) = default;
|
||||
constexpr ArraySize( const ArraySize &rhs ) = default;
|
||||
constexpr ArraySize &operator=( ArraySize &&rhs ) = default;
|
||||
constexpr ArraySize &operator=( const ArraySize &rhs ) = default;
|
||||
|
||||
/*!
|
||||
* Access the ith dimension
|
||||
* @param i Index to access
|
||||
*/
|
||||
constexpr ARRAY_ATTRIBUTE size_t operator[]( size_t i ) const { return d_N[i]; }
|
||||
|
||||
//! Return the number of dimensions
|
||||
constexpr ARRAY_ATTRIBUTE uint8_t ndim() const { return d_ndim; }
|
||||
|
||||
//! Return the number of dimensions
|
||||
constexpr ARRAY_ATTRIBUTE size_t size() const { return d_ndim; }
|
||||
|
||||
//! Return the total number of elements in the array
|
||||
constexpr ARRAY_ATTRIBUTE size_t length() const { return d_length; }
|
||||
|
||||
//! Resize the dimension
|
||||
constexpr void resize( uint8_t dim, size_t N )
|
||||
{
|
||||
if ( dim >= d_ndim )
|
||||
throw std::out_of_range( "Invalid dimension" );
|
||||
d_N[dim] = N;
|
||||
d_length = 1;
|
||||
for ( unsigned long i : d_N )
|
||||
d_length *= i;
|
||||
}
|
||||
|
||||
/*!
|
||||
* Reshape the Array so that the number of dimensions is the
|
||||
* max of ndim and the largest dim>1.
|
||||
* @param ndim Desired number of dimensions
|
||||
*/
|
||||
constexpr void setNdim( uint8_t ndim ) { d_ndim = std::max( ndim, d_ndim ); }
|
||||
|
||||
//! Returns an iterator to the beginning
|
||||
constexpr const size_t *begin() const { return d_N; }
|
||||
|
||||
//! Returns an iterator to the end
|
||||
constexpr const size_t *end() const { return d_N + d_ndim; }
|
||||
|
||||
// Check if two array sizes are equal
|
||||
constexpr ARRAY_ATTRIBUTE bool operator==( const ArraySize &rhs ) const
|
||||
{
|
||||
return d_ndim == rhs.d_ndim && memcmp( d_N, rhs.d_N, sizeof( d_N ) ) == 0;
|
||||
}
|
||||
|
||||
// Check if two array sizes are equal (ignoring the dimension)
|
||||
constexpr ARRAY_ATTRIBUTE bool approxEqual( const ArraySize &rhs ) const
|
||||
{
|
||||
return ( length() == 0 && rhs.length() == 0 ) || memcmp( d_N, rhs.d_N, sizeof( d_N ) ) == 0;
|
||||
}
|
||||
|
||||
//! Check if two matrices are not equal
|
||||
constexpr ARRAY_ATTRIBUTE bool operator!=( const ArraySize &rhs ) const
|
||||
{
|
||||
return d_ndim != rhs.d_ndim || memcmp( d_N, rhs.d_N, sizeof( d_N ) ) != 0;
|
||||
}
|
||||
|
||||
//! Maximum supported dimension
|
||||
constexpr ARRAY_ATTRIBUTE static uint8_t maxDim() { return 5u; }
|
||||
|
||||
//! Get the index
|
||||
constexpr ARRAY_ATTRIBUTE size_t index( size_t i ) const
|
||||
{
|
||||
CHECK_ARRAY_LENGTH( i );
|
||||
return i;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
constexpr ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2 ) const
|
||||
{
|
||||
size_t index = i1 + i2 * d_N[0];
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
constexpr ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2, size_t i3 ) const
|
||||
{
|
||||
size_t index = i1 + d_N[0] * ( i2 + d_N[1] * i3 );
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
constexpr ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2, size_t i3, size_t i4 ) const
|
||||
{
|
||||
size_t index = i1 + d_N[0] * ( i2 + d_N[1] * ( i3 + d_N[2] * i4 ) );
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
//! Get the index
|
||||
constexpr ARRAY_ATTRIBUTE size_t index(
|
||||
size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const
|
||||
{
|
||||
size_t index = i1 + d_N[0] * ( i2 + d_N[1] * ( i3 + d_N[2] * ( i4 + d_N[3] * i5 ) ) );
|
||||
CHECK_ARRAY_LENGTH( index );
|
||||
return index;
|
||||
}
|
||||
|
||||
private:
|
||||
uint8_t d_ndim;
|
||||
size_t d_length;
|
||||
size_t d_N[5];
|
||||
};
|
||||
|
||||
|
||||
// Function to concatenate dimensions of two array sizes
|
||||
constexpr ArraySize cat( const ArraySize &x, const ArraySize &y )
|
||||
{
|
||||
if ( x.ndim() + y.ndim() > ArraySize::maxDim() )
|
||||
throw std::out_of_range( "Maximum number of dimensions exceeded" );
|
||||
size_t N[ArraySize::maxDim()] = { 0 };
|
||||
for ( int i = 0; i < x.ndim(); i++ )
|
||||
N[i] = x[i];
|
||||
for ( int i = 0; i < y.ndim(); i++ )
|
||||
N[i + x.ndim()] = y[i];
|
||||
return ArraySize( x.ndim() + y.ndim(), N );
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
@ -2,7 +2,7 @@
|
||||
#define included_FunctionTable
|
||||
|
||||
|
||||
#include "common/Array.h"
|
||||
#include "common/ArraySize.h"
|
||||
|
||||
#include <functional>
|
||||
|
||||
@ -65,17 +65,43 @@ public:
|
||||
* @param[out] c The output array
|
||||
*/
|
||||
template<class TYPE, class FUN>
|
||||
static void multiply(
|
||||
static inline void multiply(
|
||||
const Array<TYPE, FUN> &a, const Array<TYPE, FUN> &b, Array<TYPE, FUN> &c );
|
||||
|
||||
/*!
|
||||
* Perform dgemv/dgemm equavalent operation ( C = alpha*A*B + beta*C )
|
||||
* @param[in] alpha The scalar value alpha
|
||||
* @param[in] A The first array
|
||||
* @param[in] B The second array
|
||||
* @param[in] beta The scalar value alpha
|
||||
* @param[in,out] c The output array C
|
||||
*/
|
||||
template<class TYPE, class FUN>
|
||||
static void gemm( const TYPE alpha, const Array<TYPE, FUN> &A, const Array<TYPE, FUN> &B,
|
||||
const TYPE beta, Array<TYPE, FUN> &C );
|
||||
|
||||
/*!
|
||||
* Perform axpy equavalent operation ( y = alpha*x + y )
|
||||
* @param[in] alpha The scalar value alpha
|
||||
* @param[in] x The input array x
|
||||
* @param[in,out] y The output array y
|
||||
*/
|
||||
template<class TYPE, class FUN>
|
||||
static void axpy( const TYPE alpha, const Array<TYPE, FUN> &x, Array<TYPE, FUN> &y );
|
||||
|
||||
/*!
|
||||
* Check if two arrays are approximately equal
|
||||
* @param[in] A The first array
|
||||
* @param[in] B The second array
|
||||
* @param[in] tol The tolerance
|
||||
*/
|
||||
template<class TYPE, class FUN>
|
||||
static bool equals( const Array<TYPE, FUN> &A, const Array<TYPE, FUN> &B, TYPE tol );
|
||||
|
||||
|
||||
private:
|
||||
FunctionTable();
|
||||
|
||||
template<class T>
|
||||
static inline void rand( size_t N, T *x );
|
||||
};
|
||||
|
||||
#include "common/FunctionTable.hpp"
|
||||
|
||||
#endif
|
||||
|
||||
@ -13,37 +13,30 @@
|
||||
/********************************************************
|
||||
* Random number initialization *
|
||||
********************************************************/
|
||||
template<class TYPE>
|
||||
static inline typename std::enable_if<std::is_integral<TYPE>::value>::type genRand(
|
||||
size_t N, TYPE* x )
|
||||
{
|
||||
std::random_device rd;
|
||||
std::mt19937 gen( rd() );
|
||||
std::uniform_int_distribution<TYPE> dis;
|
||||
for ( size_t i = 0; i < N; i++ )
|
||||
x[i] = dis( gen );
|
||||
}
|
||||
template<class TYPE>
|
||||
static inline typename std::enable_if<std::is_floating_point<TYPE>::value>::type genRand(
|
||||
size_t N, TYPE* x )
|
||||
{
|
||||
std::random_device rd;
|
||||
std::mt19937 gen( rd() );
|
||||
std::uniform_real_distribution<TYPE> dis( 0, 1 );
|
||||
for ( size_t i = 0; i < N; i++ )
|
||||
x[i] = dis( gen );
|
||||
}
|
||||
template<class TYPE, class FUN>
|
||||
void FunctionTable::rand( Array<TYPE, FUN> &x )
|
||||
inline void FunctionTable::rand( Array<TYPE, FUN>& x )
|
||||
{
|
||||
FunctionTable::rand<TYPE>( x.length(), x.data() );
|
||||
}
|
||||
template<>
|
||||
inline void FunctionTable::rand<double>( size_t N, double *x )
|
||||
{
|
||||
std::random_device rd;
|
||||
std::mt19937 gen( rd() );
|
||||
std::uniform_real_distribution<> dis( 0, 1 );
|
||||
for ( size_t i = 0; i < N; i++ )
|
||||
x[i] = dis( gen );
|
||||
}
|
||||
template<>
|
||||
inline void FunctionTable::rand<float>( size_t N, float *x )
|
||||
{
|
||||
std::random_device rd;
|
||||
std::mt19937 gen( rd() );
|
||||
std::uniform_real_distribution<> dis( 0, 1 );
|
||||
for ( size_t i = 0; i < N; i++ )
|
||||
x[i] = dis( gen );
|
||||
}
|
||||
template<>
|
||||
inline void FunctionTable::rand<int>( size_t N, int *x )
|
||||
{
|
||||
std::random_device rd;
|
||||
std::mt19937 gen( rd() );
|
||||
std::uniform_int_distribution<> dis;
|
||||
for ( size_t i = 0; i < N; i++ )
|
||||
x[i] = dis( gen );
|
||||
genRand<TYPE>( x.length(), x.data() );
|
||||
}
|
||||
|
||||
|
||||
@ -51,11 +44,11 @@ inline void FunctionTable::rand<int>( size_t N, int *x )
|
||||
* Reduction *
|
||||
********************************************************/
|
||||
template<class TYPE, class FUN, typename LAMBDA>
|
||||
inline TYPE FunctionTable::reduce( LAMBDA &op, const Array<TYPE, FUN> &A )
|
||||
inline TYPE FunctionTable::reduce( LAMBDA& op, const Array<TYPE, FUN>& A )
|
||||
{
|
||||
if ( A.length() == 0 )
|
||||
return TYPE();
|
||||
const TYPE *x = A.data();
|
||||
const TYPE* x = A.data();
|
||||
TYPE y = x[0];
|
||||
const size_t N = A.length();
|
||||
for ( size_t i = 1; i < N; i++ )
|
||||
@ -68,7 +61,7 @@ inline TYPE FunctionTable::reduce( LAMBDA &op, const Array<TYPE, FUN> &A )
|
||||
* Unary transformation *
|
||||
********************************************************/
|
||||
template<class TYPE, class FUN, typename LAMBDA>
|
||||
inline void FunctionTable::transform( LAMBDA &fun, const Array<TYPE, FUN> &x, Array<TYPE, FUN> &y )
|
||||
inline void FunctionTable::transform( LAMBDA& fun, const Array<TYPE, FUN>& x, Array<TYPE, FUN>& y )
|
||||
{
|
||||
y.resize( x.size() );
|
||||
const size_t N = x.length();
|
||||
@ -77,9 +70,9 @@ inline void FunctionTable::transform( LAMBDA &fun, const Array<TYPE, FUN> &x, Ar
|
||||
}
|
||||
template<class TYPE, class FUN, typename LAMBDA>
|
||||
inline void FunctionTable::transform(
|
||||
LAMBDA &fun, const Array<TYPE, FUN> &x, const Array<TYPE, FUN> &y, Array<TYPE, FUN> &z )
|
||||
LAMBDA& fun, const Array<TYPE, FUN>& x, const Array<TYPE, FUN>& y, Array<TYPE, FUN>& z )
|
||||
{
|
||||
if ( !x.sizeMatch( y ) )
|
||||
if ( x.size() != y.size() )
|
||||
throw std::logic_error( "Sizes of x and y do not match" );
|
||||
z.resize( x.size() );
|
||||
const size_t N = x.length();
|
||||
@ -88,29 +81,157 @@ inline void FunctionTable::transform(
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* axpy *
|
||||
********************************************************/
|
||||
template<class TYPE>
|
||||
inline void call_axpy( size_t N, const TYPE alpha, const TYPE* x, TYPE* y );
|
||||
template<>
|
||||
inline void call_axpy<float>( size_t, const float, const float*, float* )
|
||||
{
|
||||
throw std::logic_error( "LapackWrappers not configured" );
|
||||
}
|
||||
template<>
|
||||
inline void call_axpy<double>( size_t, const double, const double*, double* )
|
||||
{
|
||||
throw std::logic_error( "LapackWrappers not configured" );
|
||||
}
|
||||
template<class TYPE>
|
||||
inline void call_axpy( size_t N, const TYPE alpha, const TYPE* x, TYPE* y )
|
||||
{
|
||||
for ( size_t i = 0; i < N; i++ )
|
||||
y[i] += alpha * x[i];
|
||||
}
|
||||
template<class TYPE, class FUN>
|
||||
void FunctionTable::axpy( const TYPE alpha, const Array<TYPE, FUN>& x, Array<TYPE, FUN>& y )
|
||||
{
|
||||
if ( x.size() != y.size() )
|
||||
throw std::logic_error( "Array sizes do not match" );
|
||||
call_axpy( x.length(), alpha, x.data(), y.data() );
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Multiply two arrays *
|
||||
********************************************************/
|
||||
template<class TYPE>
|
||||
inline void call_gemv( size_t M, size_t N, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* x, TYPE* y );
|
||||
template<>
|
||||
inline void call_gemv<double>(
|
||||
size_t, size_t, double, double, const double*, const double*, double* )
|
||||
{
|
||||
throw std::logic_error( "LapackWrappers not configured" );
|
||||
}
|
||||
template<>
|
||||
inline void call_gemv<float>( size_t, size_t, float, float, const float*, const float*, float* )
|
||||
{
|
||||
throw std::logic_error( "LapackWrappers not configured" );
|
||||
}
|
||||
template<class TYPE>
|
||||
inline void call_gemv(
|
||||
size_t M, size_t N, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* x, TYPE* y )
|
||||
{
|
||||
for ( size_t i = 0; i < M; i++ )
|
||||
y[i] = beta * y[i];
|
||||
for ( size_t j = 0; j < N; j++ ) {
|
||||
for ( size_t i = 0; i < M; i++ )
|
||||
y[i] += alpha * A[i + j * M] * x[j];
|
||||
}
|
||||
}
|
||||
template<class TYPE>
|
||||
inline void call_gemm(
|
||||
size_t M, size_t N, size_t K, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* B, TYPE* C );
|
||||
template<>
|
||||
inline void call_gemm<double>( size_t, size_t, size_t, double, double, const double*, const double*, double* )
|
||||
{
|
||||
throw std::logic_error( "LapackWrappers not configured" );
|
||||
}
|
||||
template<>
|
||||
inline void call_gemm<float>( size_t, size_t, size_t, float, float, const float*, const float*, float* )
|
||||
{
|
||||
throw std::logic_error( "LapackWrappers not configured" );
|
||||
}
|
||||
template<class TYPE>
|
||||
inline void call_gemm(
|
||||
size_t M, size_t N, size_t K, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* B, TYPE* C )
|
||||
{
|
||||
for ( size_t i = 0; i < K * M; i++ )
|
||||
C[i] = beta * C[i];
|
||||
for ( size_t k = 0; k < K; k++ ) {
|
||||
for ( size_t j = 0; j < N; j++ ) {
|
||||
for ( size_t i = 0; i < M; i++ )
|
||||
C[i + k * M] += alpha * A[i + j * M] * B[j + k * N];
|
||||
}
|
||||
}
|
||||
}
|
||||
template<class TYPE, class FUN>
|
||||
void FunctionTable::gemm( const TYPE alpha, const Array<TYPE, FUN>& a, const Array<TYPE, FUN>& b,
|
||||
const TYPE beta, Array<TYPE, FUN>& c )
|
||||
{
|
||||
if ( a.ndim() == 2 && b.ndim() == 1 ) {
|
||||
if ( a.size( 1 ) != b.size( 0 ) )
|
||||
throw std::logic_error( "Inner dimensions must match" );
|
||||
call_gemv<TYPE>( a.size( 0 ), a.size( 1 ), alpha, beta, a.data(), b.data(), c.data() );
|
||||
} else if ( a.ndim() <= 2 && b.ndim() <= 2 ) {
|
||||
if ( a.size( 1 ) != b.size( 0 ) )
|
||||
throw std::logic_error( "Inner dimensions must match" );
|
||||
call_gemm<TYPE>(
|
||||
a.size( 0 ), a.size( 1 ), b.size( 1 ), alpha, beta, a.data(), b.data(), c.data() );
|
||||
} else {
|
||||
throw std::logic_error( "Not finished yet" );
|
||||
}
|
||||
}
|
||||
template<class TYPE, class FUN>
|
||||
void FunctionTable::multiply(
|
||||
const Array<TYPE, FUN> &a, const Array<TYPE, FUN> &b, Array<TYPE, FUN> &c )
|
||||
const Array<TYPE, FUN>& a, const Array<TYPE, FUN>& b, Array<TYPE, FUN>& c )
|
||||
{
|
||||
if ( a.ndim() <= 2 && b.ndim() <= 2 ) {
|
||||
if ( a.ndim() == 2 && b.ndim() == 1 ) {
|
||||
if ( a.size( 1 ) != b.size( 0 ) )
|
||||
throw std::logic_error( "Inner dimensions must match" );
|
||||
c.resize( a.size( 0 ) );
|
||||
call_gemv<TYPE>( a.size( 0 ), a.size( 1 ), 1, 0, a.data(), b.data(), c.data() );
|
||||
} else if ( a.ndim() <= 2 && b.ndim() <= 2 ) {
|
||||
if ( a.size( 1 ) != b.size( 0 ) )
|
||||
throw std::logic_error( "Inner dimensions must match" );
|
||||
c.resize( a.size( 0 ), b.size( 1 ) );
|
||||
c.fill( 0 );
|
||||
for ( size_t k = 0; k < b.size( 1 ); k++ ) {
|
||||
for ( size_t j = 0; j < a.size( 1 ); j++ ) {
|
||||
for ( size_t i = 0; i < a.size( 0 ); i++ ) {
|
||||
c( i, k ) += a( i, j ) * b( j, k );
|
||||
}
|
||||
}
|
||||
}
|
||||
call_gemm<TYPE>(
|
||||
a.size( 0 ), a.size( 1 ), b.size( 1 ), 1, 0, a.data(), b.data(), c.data() );
|
||||
} else {
|
||||
throw std::logic_error( "Not finished yet" );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Check if two arrays are equal *
|
||||
********************************************************/
|
||||
template<class TYPE, class FUN>
|
||||
inline typename std::enable_if<!std::is_floating_point<TYPE>::value, bool>::type
|
||||
FunctionTableCompare( const Array<TYPE, FUN>& a, const Array<TYPE, FUN>& b, TYPE )
|
||||
{
|
||||
bool pass = true;
|
||||
if ( a.size() != b.size() )
|
||||
throw std::logic_error( "Sizes of x and y do not match" );
|
||||
for ( size_t i = 0; i < a.length(); i++ )
|
||||
pass = pass && a( i ) == b( i );
|
||||
return pass;
|
||||
}
|
||||
template<class TYPE, class FUN>
|
||||
inline typename std::enable_if<std::is_floating_point<TYPE>::value, bool>::type
|
||||
FunctionTableCompare( const Array<TYPE, FUN>& a, const Array<TYPE, FUN>& b, TYPE tol )
|
||||
{
|
||||
bool pass = true;
|
||||
if ( a.size() != b.size() )
|
||||
throw std::logic_error( "Sizes of x and y do not match" );
|
||||
for ( size_t i = 0; i < a.length(); i++ )
|
||||
pass = pass && ( std::abs( a( i ) - b( i ) ) < tol );
|
||||
return pass;
|
||||
}
|
||||
template<class TYPE, class FUN>
|
||||
bool FunctionTable::equals( const Array<TYPE, FUN>& a, const Array<TYPE, FUN>& b, TYPE tol )
|
||||
{
|
||||
return FunctionTableCompare( a, b, tol );
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
@ -4,13 +4,13 @@
|
||||
#include "common/Domain.h"
|
||||
#include "common/SpherePack.h"
|
||||
|
||||
using namespace std;
|
||||
#include "ProfilerApp.h"
|
||||
|
||||
|
||||
/*
|
||||
* Compare the measured and analytical curvature for a sphere
|
||||
*
|
||||
*/
|
||||
|
||||
std::shared_ptr<Database> loadInputs( )
|
||||
{
|
||||
//auto db = std::make_shared<Database>( "Domain.in" );
|
||||
@ -38,7 +38,7 @@ int main(int argc, char **argv)
|
||||
int Nx = db->getVector<int>( "n" )[0];
|
||||
int Ny = db->getVector<int>( "n" )[1];
|
||||
int Nz = db->getVector<int>( "n" )[2];
|
||||
std::shared_ptr<Domain> Dm = std::shared_ptr<Domain>(new Domain(db,comm));
|
||||
auto Dm = std::make_shared<Domain>( db, comm );
|
||||
|
||||
Nx+=2; Ny+=2; Nz+=2;
|
||||
DoubleArray Phase(Nx,Ny,Nz);
|
||||
@ -98,6 +98,7 @@ int main(int argc, char **argv)
|
||||
printf(" Euler characteristic = %f (analytical = 2.0) \n",sphere.Xi);
|
||||
|
||||
}
|
||||
PROFILE_SAVE("test_dcel_minkowski");
|
||||
MPI_Barrier(comm);
|
||||
MPI_Finalize();
|
||||
return toReturn;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user