Loading...
Searching...
No Matches
tensor-exp.cc
Go to the documentation of this file.
// exp(tensor) : 2D is explicit while 3D uses pade' approx
//
#include "rheolef/tensor.h"
#include "rheolef/compiler_eigen.h"
#undef _RHEOLEF_EIGEN_HAVE_EXPM // not yet supported officially, but faster than via eig(a)
#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
#include <unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h>
#endif // _RHEOLEF_EIGEN_HAVE_EXPM
namespace rheolef {
// =========================================================================
// 2D case: explicit
// =========================================================================
//
// see tensor_exp_tst.mac
// see also maxima: matrixexp.usage allows explicit expm() computation in 2D
// http://www.ma.utexas.edu/pipermail/maxima/2006/003031.html
// /usr/share/maxima/5.27.0/share/linearalgebra/matrixexp.usage
// refs. KanGueFor-2009, p. 50 (is buggy: missing 2 factor in A00 & A11...)
template <class T>
static
tensor_basic<T>
exp2d (const tensor_basic<T>& chi) {
static T eps = 1e3*std::numeric_limits<T>::epsilon();
T a = chi(0,0), b = chi(1,0), c = chi(1,1);
T b2=sqr(b);
T d2 = sqr(a-c)+4*b2;
T d = sqrt(d2);
tensor_basic<T> A;
if (fabs(d) < eps) { // chi = a*I
A(0,0) = A(1,1) = exp(a);
A(0,1) = A(1,0) = 0.;
return A;
}
T a2=sqr(a), c2=sqr(c);
T ed = exp(d);
T k = exp((a+c-d)/2)/(2*d);
T x1 = (ed+1)*d;
T x2 = (ed-1)*(a-c);
T x3 = 2*(ed-1)*b;
A(0,0) = k*(x1 + x2);
A(1,1) = k*(x1 - x2);
A(1,0) = A(0,1) = k*x3;
return A;
}
template<typename T, int N>
Eigen::Matrix<T,N,N>
expm_eig (const Eigen::Matrix<T,N,N>& a) {
using namespace Eigen;
// eig of sym matrix, then exp
SelfAdjointEigenSolver<Matrix<T,N,N> > eig (a);
Matrix<T,N,1> e;
for (size_t i = 0; i < N; ++i) {
e(i) = exp(eig.eigenvalues()(i));
}
return eig.eigenvectors()*e.asDiagonal()*eig.eigenvectors().transpose();
}
#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
template<typename T, int N>
Eigen::Matrix<T,N,N>
expm_pade (const Eigen::Matrix<T,N,N>& a) {
return a.exp(); // in "eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h"
}
#endif // _RHEOLEF_EIGEN_HAVE_EXPM
// =========================================================================
// tensor interface
// =========================================================================
template <class T>
tensor_basic<T>
exp (const tensor_basic<T>& a, size_t d) {
// check that a is symmetric ; otherwise, the result would be complex
check_macro (a.is_symmetric(d), "exp: tensor should be symmetric");
using namespace std;
using namespace Eigen;
if (d == 1) return tensor(exp(a(0,0)));
if (d == 2) return exp2d (a);
// 3D case:
Matrix<T,3,3> a1 (d,d);
for (size_t i = 0; i < d; ++i)
for (size_t j = 0; j < d; ++j) a1(i,j) = a(i,j);
#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
Matrix<T,3,3> b1 = expm_pade(a1);
#else // _RHEOLEF_EIGEN_HAVE_EXPM
Matrix<T,3,3> b1 = expm_eig (a1);
#endif // _RHEOLEF_EIGEN_HAVE_EXPM
tensor b;
for (size_t i = 0; i < d; ++i)
for (size_t j = 0; j < d; ++j) b(i,j) = b1(i,j);
return b;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template tensor_basic<T> exp (const tensor_basic<T>& a, size_t d); \
_RHEOLEF_instanciation(Float)
}// namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
Float
see the Float page for the full documentation
rheolef::tensor_basic
Definition tensor.h:90
tensor
see the tensor page for the full documentation
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Eigen
Definition field_eigen.h:35
mkgeo_ball.b
int b
Definition mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition mkgeo_ball.sh:153
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
mkgeo_contraction.c2
int c2
Definition mkgeo_contraction.sh:199
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::exp
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition tensor-exp.cc:92
rheolef::expm_eig
Eigen::Matrix< T, N, N > expm_eig(const Eigen::Matrix< T, N, N > &a)
Definition tensor-exp.cc:70
std
STL namespace.
N
size_t N
Definition neumann-laplace-check.cc:24
d
size_t d
Definition neumann-laplace-lambda.cc:28