[eigen] Memory fault when using a function returning a SelfAdjointView |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: [eigen] Memory fault when using a function returning a SelfAdjointView
- From: Douglas Bates <bates@xxxxxxxxxxxxx>
- Date: Wed, 26 Oct 2011 13:20:34 -0500
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:sender:date:x-google-sender-auth:message-id:subject :from:to:content-type; bh=C4daexdm5CagGWkm/UV7IeEZnLIIkhenRNWuMW2X52k=; b=ZomHgI8wmNzjP3FZnMEI30LgX++ChxjJQp41DA4HnXEPnWK7NiTW4UU+k00xJ/oBQQ npdNktb49ns62q73iEkF9S8+1AwG8NcQODeERA8wSJjNNWwdDtkt3VGFe6SpeFmzRL37 i0osqib+0DCqxsXb7Mw6qTOiHOCpNhleN1cC4=
In statistical computing I will frequently want to calculate X'X from
an m by n model matrix X where m >= n. I use the rankUpdate method
for the SelfAdjointView<MatrixXd, Lower> class to calculate the
symmetric result.
I can do this inline but if I try to calculate the value in a function
returning the SelfAdjointView, rather than a MatrixXd, I encounter a
bad_alloc error under Ubuntu 11.10 (64-bit). I enclose an example
program.
The reason that I want to return the SelfAdjointView is because I want
to create different types of decompositions, such as LLT, LDLT and
SelfAdjointEigenSolver from the result. If I return the value as a
matrix I think I need to apply the selfadjointView method before
creating the decomposition, which seems redundant.
#include <iostream>
#include <eigen3/Eigen/Dense>
using namespace std;
using namespace Eigen;
SelfAdjointView<MatrixXd, Lower> XtX(const MatrixXd& X) {
return MatrixXd(X.cols(), X.cols()).setZero().selfadjointView<Lower>().rankUpdate(X.adjoint());
}
MatrixXd XtXmat(const MatrixXd& X) {
return MatrixXd(MatrixXd(X.cols(), X.cols()).setZero().
selfadjointView<Lower>().rankUpdate(X.adjoint()));
}
int main() {
MatrixXd X(4, 2);
X <<
1, 1,
1, 2,
1, 3,
1, 4;
// inline calculation of XtX is okay
MatrixXd xtx(MatrixXd(X.cols(), X.cols()).setZero().
selfadjointView<Lower>().rankUpdate(X.adjoint()));
cout << "XtX(inline):\n" << xtx << endl;
// returning XtX as a MatrixXd is okay
cout << "XtX(mat):\n" << XtXmat(X) << endl;
// but not returning a SelfAdjointView
MatrixXd XtXres(XtX(X));
cout << "XtX(function):\n" << XtXres << endl;
}