Integration of a matrix function over time with Boost / calculate the Controllability Gramian in C++

Question

I have Matrices A and B. I want to calculate the Controllability Gramian like defined on Wikipedia:

Luckily, Boost already provides integration methods, so I wrote following function:

#include <boost/math/quadrature/gauss.hpp>
#include <unsupported/Eigen/MatrixFunctions>
#include <Eigen/Geometry>

typedef Eigen::Matrix<std::complex<double>, 12, 12> GramianMatrix;

GramianMatrix calcControllabilityGramian(Eigen::Matrix<double, 12, 12> A, Eigen::Matrix<double, 12, 4> B) {
  // the integral function
  auto controllabilityGramianIntegral = [&A, &B](const double & t) -> GramianMatrix {
      return (A * t).exp() * B * B.transpose() * (A.transpose() * t).exp();
  };
  // now do the integration
  boost::math::quadrature::gauss<double, 10> integrator;
  GramianMatrix W_t = integrator.integrate(controllabilityGramianIntegral, 0.0, 1.0);

  return W_t;
}

It seems that Boost expects the output of this integral to be a double value and not a matrix. The compiler gives following error: /usr/local/include/boost/math/quadrature/gauss.hpp:1183:9: No viable conversion from 'double' to 'K' (aka 'Matrix<std::complex<double>, 12, 12>'). So is there a way to do this with a matrix type output? Should I use Boost or are there better libraries for this issue? I would like to hear your suggestions!