使用本征时的怪异行为 [英] Weird Behavior when using Eigen

问题描述

我正在给Eigen写一个包装供个人使用，遇到以下奇怪的行为：

I am writing a wrapper to Eigen for my personal use and I encountered the following weird behavior:

void get_QR(MatrixXd A, MatrixXd& Q, MatrixXd& R) {
        HouseholderQR<MatrixXd> qr(A);
        Q  =   qr.householderQ()*(MatrixXd::Identity(A.rows(),A.cols()));
        R  =   qr.matrixQR().block(0,0,A.cols(),A.cols()).triangularView<Upper>();
}

void get_QR(double* A, int m, int n, double*& Q, double*& R) {
        //      Maps the double to MatrixXd.
        Map<MatrixXd> A_E(A, m, n);

        //      Obtains the QR of A_E.
        MatrixXd Q_E, R_E;
        get_QR(A_E, Q_E, R_E);

        //      Maps the MatrixXd to double.
        Q       =       Q_E.data();
        R       =       R_E.data();
}

以下是测试：

int main(int argc, char* argv[]) {
    srand(time(NULL));

    int m   =       atoi(argv[1]);
    int n   =       atoi(argv[2]);

    //      Check the double version.

    double* A       =       new double[m*n];
    double* Q;
    double* R;

    double RANDMAX  =       double(RAND_MAX);

    //      Initialize A as a random matrix.

    for (int index=0; index<m*n; ++index) {
            A[index]        =       rand()/RANDMAX;
    }

    get_QR(A, m, n, Q, R);
    std::cout << Q[0] << std::endl;

    //      Check the MatrixXd version.

    Map<MatrixXd> A_E(A, m, n);
    MatrixXd Q_E(m,n), R_E(n,n);
    get_QR(A_E, Q_E, R_E);
    std::cout << Q[0] << std::endl;
}

我得到了不同的Q [0]值。例如，我得到 -0.421857和 -1.49563。

I get different values of Q[0]. For instance, I get "-0.421857" and "-1.49563".

谢谢

推荐答案

乔治的答案是正确的，但存在不必要的副本。更好的解决方案包括映射Q和R：

The answer of George is correct but suffers from unnecessary copies. A better solution consists in mapping Q and R:

void get_QR(const double* A, int m, int n, double*& Q, double*& R) {
    Map<const MatrixXd> A_E(A, m, n);
    Map<MatrixXd> Q_E(Q, m, n);
    Map<MatrixXd> R_E(Q, n, n);

    HouseholderQR<MatrixXd> qr(A_E);
    Q_E = qr.householderQ()*(MatrixXd::Identity(m,n));
    R_E = qr.matrixQR().block(0,0,n,n).triangularView<Upper>();
}

为了能够重用 get_QR 函数获取Eigen的对象，然后使用 Ref< MatrixXd> 代替 MatrixXd ：

In order to be able to reuse the get_QR function taking Eigen's object, then use Ref<MatrixXd> instead of MatrixXd:

void get_QR(Ref<const MatrixXd> A, Ref<MatrixXd> Q, Ref<MatrixXd> R) {
    HouseholderQR<MatrixXd> qr(A);
    Q = qr.householderQ()*(MatrixXd::Identity(A.rows(),A.cols()));
    R = qr.matrixQR().block(0,0,A.cols(),A.cols()).triangularView<Upper>();
}

void get_QR(const double* A, int m, int n, double* Q, double* R) {
    Map<const MatrixXd> A_E(A, m, n);
    Map<MatrixXd> Q_E(Q, m, n);
    Map<MatrixXd> R_E(R, n, n);
    get_QR(A_E, Q_E, R_E);
}

Ref< MatrixXd> 可以包装类似于 MatrixXd 的任何本征对象，而无需任何副本。其中包括 MatrixXd 本身，以及 Map 和 Block 表达式。这需要Eigen 3.2。

The Ref<MatrixXd> can wrap any Eigen's object that is similar to a MatrixXd without any copy. This include MatrixXd itself, as well as Map and Block expressions. This requires Eigen 3.2.

这篇关于使用本征时的怪异行为的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！