新手问题-使用Eigen

如何解决新手问题-使用Eigen

我是C ++的新手，正在尝试将以下MATLAB代码转换为C ++：

x = [1,2,3];
y = [4,5,7];
px = 0.5;
py = 0.5;

m = -1:1:1;

v = [1+i,2+i,3+i;1+2i,2+2i,3+2i;1+3i,2+3i,3+3i];

Tr5 = real(exp(2i*pi/px*x(1)*m)*v*exp(2i*pi/py*m'*y(1)));

fprintf('Tr5 = %f\n',Tr5);

返回Tr5 = 18。

这是我在C ++中所做的尝试，在这里我在每一行上都做了尽可能少的尝试来查找错误：

#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>
#include <math.h>
#include <iterator>
#include <vector>
#include <algorithm>
#include <complex>
#include <iostream>

int main()
{

double x [3] = {1,3};
double y [3] = {4,7};
const std::complex<double> im(0.0,1.0);
double px = 0.5;
double py = 0.5;

Eigen::MatrixXcd Tr1;
Eigen::MatrixXcd Tr2;
Eigen::MatrixXcd Tr3;
Eigen::MatrixXcd Tr4;

Eigen::MatrixXd m = Eigen::VectorXd::LinSpaced(2*1+1,-1,1);
Eigen::MatrixXcd v4(3,3);
v4(0,0) = 1.0d + 1.0 * im;
v4(0,1) = 2.0d + 1.0 * im;
v4(0,2) = 3.0d + 1.0 * im;
v4(1,0) = 1.0d + 2.0 * im;
v4(1,1) = 2.0d + 2.0 * im;
v4(1,2) = 3.0d + 2.0 * im;
v4(2,0) = 1.0d + 3.0 * im;
v4(2,1) = 2.0d + 3.0 * im;
v4(2,2) = 3.0d + 3.0 * im;

Tr1 = 2.0*im*M_PI/px*x[1]*m;
Tr1 = Tr1.exp();

Tr2 = 2.0*im*M_PI/py*y[1]*m;
Tr2 = Tr2.exp();

Tr3 = Tr2.transpose();
Tr4 = Tr1 * v4 * Tr3;

Eigen::MatrixXcd Tr5 = Tr4.real(); // This should Now be "18"

std::cout << "Tr5 = " << Tr5 << std::endl;

}

但是，当我尝试编译并运行时，出现错误：

/usr/local/include/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:436: const Eigen::MatrixExponentialReturnValue<Derived> Eigen::MatrixBase<Derived>::exp() const [with Derived = Eigen::Matrix<std::complex<double>,-1>]: Assertion `rows() == cols()' Failed.

有人可以帮我提供代码吗？

非常感谢，汤姆·怀特

解决方法

似乎您没有使用按系数计算的指数，而是矩阵指数（包含在中）。矩阵指数（请参阅https://en.wikipedia.org/wiki/Matrix_exponential）是平方矩阵的矩阵函数。您会看到该断言错误，因为矩阵'm'不是正方形。

您应该使用逐系数系数。在Eigen中，为数组定义了逐系数函数，但是也可以通过将其首先转换为数组来用于矩阵：m.array()（请参阅https://eigen.tuxfamily.org/dox/group__CoeffwiseMathFunctions.html）但是，您需要将数组重新投影执行矩阵乘法之前先移至矩阵：leftSide.array().exp().matrix() * v * rightSide.array().exp().matrix()

最后，您可以使用逗号初始值设定项语法填充矩阵v（我相信它更容易阅读）。我将以一个示例为例，当然，您的初始化也可以。

示例（语法类似于MATLAB）

#define _USE_MATH_DEFINES

#include "Eigen/Core"
#include <iostream>
#include <math.h>

std::complex<double> operator"" _i(long double x ) 
{ 
    return std::complex<double>(0.0,x);
} 

int main()
{
    double x [3] = {1,2,3};
    double y [3] = {4,5,7};

    double px = 0.5;
    double py = 0.5;

    Eigen::RowVectorXd m = Eigen::RowVectorXd::LinSpaced(3,-1,1);

    Eigen::MatrixXcd v(3,3);
    v << (1.0 + 1.0_i),(2.0 + 1.0_i),(3.0 + 1.0_i),(1.0 + 2.0_i),(2.0 + 2.0_i),(3.0 + 2.0_i),(1.0 + 3.0_i),(2.0 + 3.0_i),(3.0 + 3.0_i);

    double Tr5 = ((2.0_i * M_PI / px * x[0] * m).array().exp().matrix() * v * (2.0_i * M_PI / py * y[0] * m).array().exp().matrix().transpose()).real()(0);

    std::cout << "Tr5 = " << Tr5 << std::endl;

    return 0;
}

这将产生所需的输出：

Tr5 = 18

您忘记了#include <iostream>
请使用std::cout和std::endl而不是cout，endl
您错误地认为.real()成员函数返回实数（double）标量。它返回一个矩阵，请参见https://eigen.tuxfamily.org/dox/namespaceEigen.html#title49。类似地，在Matlab中（实际上，我在Octave中进行了检查，https://octave-online.net/）real(matrix)返回与Eigen完全一致的元素“即”矩阵的“实部”。所以更换

double Tr5;

使用

Eigen::MatrixXcd Tr5;

实际上，在C ++中，将对象定义推迟到真正需要时是一个好习惯，因此您可以考虑将声明与初始化结合起来

Eigen::MatrixXcd Tr5 = Tr4.real(); // This should now be a real double