mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
171 lines
4.3 KiB
C++
171 lines
4.3 KiB
C++
#include <gsl/gsl_math.h>
|
|
#include <gsl/gsl_eigen.h>
|
|
|
|
#include <vec3.hh>
|
|
|
|
template<typename T>
|
|
class mat3{
|
|
protected:
|
|
std::array<T,9> data_;
|
|
gsl_matrix_view m_;
|
|
gsl_vector *eval_;
|
|
gsl_matrix *evec_;
|
|
gsl_eigen_symmv_workspace * wsp_;
|
|
bool bdid_alloc_gsl_;
|
|
|
|
void init_gsl(){
|
|
// allocate memory for GSL operations if we haven't done so yet
|
|
if( !bdid_alloc_gsl_ )
|
|
{
|
|
m_ = gsl_matrix_view_array (&data_[0], 3, 3);
|
|
eval_ = gsl_vector_alloc (3);
|
|
evec_ = gsl_matrix_alloc (3, 3);
|
|
wsp_ = gsl_eigen_symmv_alloc (3);
|
|
bdid_alloc_gsl_ = true;
|
|
}
|
|
}
|
|
|
|
void free_gsl(){
|
|
// free memory for GSL operations if it was allocated
|
|
if( bdid_alloc_gsl_ )
|
|
{
|
|
gsl_eigen_symmv_free (wsp_);
|
|
gsl_vector_free (eval_);
|
|
gsl_matrix_free (evec_);
|
|
}
|
|
}
|
|
|
|
public:
|
|
|
|
mat3()
|
|
: bdid_alloc_gsl_(false)
|
|
{}
|
|
|
|
//! copy constructor
|
|
mat3( const mat3<T> &m)
|
|
: data_(m.data_), bdid_alloc_gsl_(false)
|
|
{}
|
|
|
|
//! move constructor
|
|
mat3( mat3<T> &&m)
|
|
: data_(std::move(m.data_)), bdid_alloc_gsl_(false)
|
|
{}
|
|
|
|
//! construct mat3 from initializer list
|
|
template<typename ...E>
|
|
mat3(E&&...e)
|
|
: data_{{std::forward<E>(e)...}}, bdid_alloc_gsl_(false)
|
|
{}
|
|
|
|
mat3<T>& operator=(const mat3<T>& m) noexcept{
|
|
data_ = m.data_;
|
|
return *this;
|
|
}
|
|
|
|
mat3<T>& operator=(const mat3<T>&& m) noexcept{
|
|
data_ = std::move(m.data_);
|
|
return *this;
|
|
}
|
|
|
|
//! destructor
|
|
~mat3(){
|
|
this->free_gsl();
|
|
}
|
|
|
|
//! bracket index access to vector components
|
|
T &operator[](size_t i) noexcept { return data_[i];}
|
|
|
|
//! const bracket index access to vector components
|
|
const T &operator[](size_t i) const noexcept { return data_[i]; }
|
|
|
|
//! matrix 2d index access
|
|
T &operator()(size_t i, size_t j) noexcept { return data_[3*i+j]; }
|
|
|
|
//! const matrix 2d index access
|
|
const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; }
|
|
|
|
//! in-place addition
|
|
mat3<T>& operator+=( const mat3<T>& rhs ) noexcept{
|
|
for (size_t i = 0; i < 9; ++i) {
|
|
(*this)[i] += rhs[i];
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
//! in-place subtraction
|
|
mat3<T>& operator-=( const mat3<T>& rhs ) noexcept{
|
|
for (size_t i = 0; i < 9; ++i) {
|
|
(*this)[i] -= rhs[i];
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
void zero() noexcept{
|
|
for (size_t i = 0; i < 9; ++i) data_[i]=0;
|
|
}
|
|
|
|
void eigen( vec3<T>& evals, vec3<T>& evec1, vec3<T>& evec2, vec3<T>& evec3 )
|
|
{
|
|
// for( auto x : data_ ){
|
|
// std::cerr << x << " " ;
|
|
// }
|
|
// std::cerr << std::endl;
|
|
// resort into symmetrix matrix
|
|
// data_[8] = data_[5];
|
|
// data_[7] = data_[4];
|
|
// data_[6] = data_[2];
|
|
// data_[5] = data_[4];
|
|
// data_[4] = data_[3];
|
|
// data_[3] = data_[1];
|
|
|
|
this->init_gsl();
|
|
|
|
gsl_eigen_symmv (&m_.matrix, eval_, evec_, wsp_);
|
|
gsl_eigen_symmv_sort (eval_, evec_, GSL_EIGEN_SORT_VAL_ASC);
|
|
|
|
for( int i=0; i<3; ++i ){
|
|
evals[i] = gsl_vector_get( eval_, i );
|
|
evec1[i] = gsl_matrix_get( evec_, i, 0 );
|
|
evec2[i] = gsl_matrix_get( evec_, i, 1 );
|
|
evec3[i] = gsl_matrix_get( evec_, i, 2 );
|
|
}
|
|
|
|
// std::cerr << "(" << evals[0] << " " << evals[1] << " " << evals[2] << ")" << std::endl;
|
|
}
|
|
};
|
|
|
|
template<typename T>
|
|
constexpr const mat3<T> operator+(const mat3<T> &lhs, const mat3<T> &rhs) noexcept
|
|
{
|
|
mat3<T> result;
|
|
for (size_t i = 0; i < 9; ++i) {
|
|
result[i] = lhs[i] + rhs[i];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
// matrix - vector multiplication
|
|
template<typename T>
|
|
vec3<T> operator*( const mat3<T> &A, const vec3<T> &v ) noexcept
|
|
{
|
|
vec3<T> result;
|
|
for( int mu=0; mu<3; ++mu ){
|
|
result[mu] = 0.0;
|
|
for( int nu=0; nu<3; ++nu ){
|
|
result[mu] += A(mu,nu)*v[nu];
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
// template<typename T>
|
|
// vec3<T> operator*( const vec3<T> &v, const mat3<T> &A ) noexcept
|
|
// {
|
|
// vec3<T> result = 0.0;
|
|
// for( int mu=0; mu<3; ++mu ){
|
|
// for( int nu=0; nu<3; ++nu ){
|
|
// result[nu] += v[mu]*A(mu,nu);
|
|
// }
|
|
// }
|
|
// return result;
|
|
// }
|