Re: [eigen] Convert Eigen::Tensor<double, 2> to Eigen::Matrix<double, -1, -1>

[Benoit Steiner <bsteiner@xxxxxxxxxx>]

You can use a matrix map to wrap the coefficients of the output tensor and use them in the solver:

f = g.contract(D, ctr);

Eigen::Map<Eigen::MatrixXd> F(f.data(), 5, 5);

Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, -1, -1> > eigensolver(F);

...

On Tue, May 17, 2016 at 7:11 PM, Andreas Copan <avcopan@xxxxxxxxx> wrote:

I am writing a code in which I need to contract a four-dimensional tensor (Eigen::Tensor<double, 4>) with a two-dimensional tensor (Eigen::Tensor<double, 2>) and then diagonalize it.  It looks like there currently isn't a .asMatrix() method for converting an Eigen::Tensor to a matrix for diagonalization, is that right?  If not, what would people recommend as the best way to do this using the functionality currently available?  Even just a simple method of filling an Eigen::Matrix with values from an Eigen::Tensor would help me out.  Here is a minimal example of what I would like to do.

#include <iostream>
#include <array>
#include <Eigen/Dense>                    // Eigen::Matrix
#include <unsupported/Eigen/CXX11/Tensor> // Eigen::Tensor

using IndexPair = Eigen::IndexPair<int>;
template<int naxes>
using Contraction = std::array<IndexPair, naxes>;

int main() {
  Eigen::Tensor<double, 4> g(5, 5, 5, 5);
  Eigen::Tensor<double, 2> D(5, 5);
  Eigen::Tensor<double, 2> f(5, 5);
  g.setRandom();
  D.setRandom();
  Contraction<2> ctr({IndexPair(1, 1), IndexPair(3, 0)});
  f = g.contract(D, ctr); // f_ij = sum_kl g_ikjl * D_lk

  /* this is what I would like to avoid: */
  Eigen::Matrix<double, -1, -1> F(5, 5);
  for(int i=0; i<5; ++i)
    for(int j=0; j<5; ++j)
      F(i,j) = f(i,j);

  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, -1, -1> > eigensolver(F);
  Eigen::Matrix<double, -1, -1> eigenvectors = eigensolver.eigenvectors();
  Eigen::Matrix<double, -1,  1> eigenvalues  = eigensolver.eigenvalues();
  std::cout << eigenvectors << std::endl;
  std::cout << eigenvalues  << std::endl;
}