Hi

I am trying to port some code from MKL BLAS to Eigen. It works with fairly large dense matrices, i.e. it's feasible to have only 1 copy of some matrix A in memory (to put it differently - keeping 2 copies of A will considerably limit volume of data that can be handled by the program).

To make it work I need a way to do LU and QR decompositions in place or at least with minimum amount of extra memory allocation. With LAPACK I can do it easily, because its primitives for LU and QR are designed as in-place transformations. With Eigen - not so much. The best I could do is build my own QR using some existing Eigen primitives (from namespace Eigen::internal::). It works, but it's extremely ugly (and I suspect can be broken with future Eigen updates).

Is there any trick to it that I am missing? (It doesn't seem so, just want to be sure.)

Also, there are things that look like they *could* be doable, but interface doesn't allow it, e.g.:

Code: Select all

FullPivLU<MatrixXf> lu(100, 100);
auto& m = lu.matrixLU();
prepareMatrix(m);
lu.compute();

Here Eigen considers the call to matrixLU() as error, because lu object was not initialized. (Also there is no compute() method without args.) It kind of makes sense, but then I would expect some code like this would work:

Code: Select all

FullPivLU<MatrixXf> lu;
MatrixXf m(100, 100);
prepareMatrix(m);
lu.compute(std::move(m));

But currently it doesn't. Is there any reason *not* to support it?

yes, that's a recurrent question for which there is no clean solution yet.

A workaround, that I suggested several times, is to create a PartialPivLU object:

PartialPivLU<MatrixXd> lu(rows,cols);

and then grab its own matrix:

MatrixXd& mat = const_cast<MatrixXd&>(lu.matrixLU());

Then you can fill mat and work inplace by calling lu.compute(mat);

For a cleaner solver, I'm still not sure what would be the best API. Perhaps C++11 move semantic could be used for dynamic objects:

lu.compute(std::move(A));

This approach is straightforward to implement as it only requires to add a compute() method taking a rvalue-ref.... but it's C++11 only (could be emulated with a Eigen::move()).

Another approach could be to use Ref<>: PartialPivLU<Ref<MatrixXd> >. Actually, it could be that this approach already works in the devel branch! (or requires only minor fixes).

interesting, the following works:

Code: Select all

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

int main()
{
  MatrixXd A0 = MatrixXd::Random(10,10);
  MatrixXd A1 = A0;
  VectorXd b = VectorXd::Random(10), x(10);
  PartialPivLU<Ref<MatrixXd> > lu(A0);
  cout << "Check that A0 has been modified: " << (A0-A1).norm()/A0.norm() << endl;
  x = lu.solve(b);
  cout << "Error: " << (A1*x-b).norm()/b.norm() << endl;
}


after changing only one line:

Code: Select all

 PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix)
-  : m_lu(matrix.rows(), matrix.rows()),
+  : m_lu(matrix.const_cast_derived()),
     m_p(matrix.rows()),


We could thus add a non-const ctor for that precise purpose....

2 ggael

Thanks.

PartialPivLU<MatrixXd> lu(rows,cols);

and then grab its own matrix:

MatrixXd& mat = const_cast<MatrixXd&>(lu.matrixLU());

This was my plan so far. But 1) I was somewhat discouraged by eigen_assert() inside matrixLU(), 2) I can't see how to make it work if I need to do e.g. LU followed by QR for the same matrix (I could of course dump the matrix to HDD and then read it back it into another transform object, but it seems rather suboptimal).

This approach is straightforward to implement as it only requires to add a compute() method taking a rvalue-ref.... but it's C++11 only (could be emulated with a Eigen::move()).

C++11 is not a problem for me. Indeed, it seems easy to modify Eigen source to handle my case using std::move() and rvalue argument. However, I am reluctant to introduce incompatible changes to our copy of Eigen. If there is a native solution coming, I'd rather wait for it.

right, I forgot about this assert. We have a related bug entry: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=707. The Ref<> approach seems to be quite reasonable and more general than std::move as Ref<> would allow working inplace within sub-matrices, whereas we can only move a MatrixXd within a MatrixXd.

It's in: http://eigen.tuxfamily.org/dox-devel/gr ... ition.html