This Eigen sequence

Map<MatrixXd> W(W_ptr, count, count);
MatrixXd T(count, count);
...
T.computeInverse(&W);

gives an error for the last line stating

no matching function for call to 'Eigen::Matrix<double, 10000, 10000, 3, 10000, 10000>::computeInverse(Eigen::Map<Eigen::Matrix<double, 10000, 10000, 3, 10000, 10000>, 1>*)

So, how does one properly pass

1) a reference or
2) a pointer

to an Eigen method that expects the *non-mapped* equivalent? Without the Map<>, everything is fine but I want the Map so as to avoid excessive malloc() calls [by pre-allocating a big buffer].

TIA.

Hm this was a limitation of Eigen 2.0. computeInverse() expected that the destination type is the plain-matrix version of the current type.

I've fixed that for you in the 2.0 branch, enjoy! Meanwhile the default branch has been fixed long ago.

Perhaps I misunderstood your reply ("2.0 branch").

I went to the Eigen site and re-downloaded the 2.0.11 "latest stable release" (.zip) but nothing has changed. It still gives the error cited above.

I did not, and never have, downloaded source code or anything but the stable-release link.

The 2.0 branch is the branch that will eventually be released as 2.0.12. It's the stable branch. In order to get it:
http://bitbucket.org/eigen/eigen/get/2.0.tar.bz2
(other formats possible, see http://eigen.tuxfamily.org/index.php?ti ... e#Download)

After using the modified version of Eigen, the error with computeInverse() went away but I got a similar one elsewhere, viz.,

-(void)aFunc:(MatrixXd&)M
{
...
}

Map<MatrixXd> M(M_ptr, count, count);
[self afunc:M];

gives

warning: cannot pass objects of non-POD type 'class Eigen::Map<Eigen::Matrix<double, 10000, 10000, 3, 10000, 10000>, 1>' through '...'; call will abort at runtime

The code that you pasted isn't c++

But since that seems to be the question, if you want to write a function that can take seamlessly either a matrix or a map, you have to let it take a MatrixBase<Derived>, template in Derived.

See the code example there:
http://eigen.tuxfamily.org/dox/classEig ... xBase.html

"The code that you pasted isn't c++ "

but it *was* a header for a .mm file so it should have been OK.

On the other hand, I cannot begin to fathom what

"...you have to let it take a MatrixBase<Derived>, template in Derived"

is supposed to mean and the supplied link has so many examples of so many kinds that I cannot tell to which you refer. Could you give a line number on that page? Thanks.

Well, there are indeed many examples for the specific methods in this class, but there's only one example in the class' own description,

Code: Select all

    template<typename Derived>
    void printFirstRow(const Eigen::MatrixBase<Derived>& x)
    {
      cout << x.row(0) << endl;
    }


Sorry, I don't know what .mm files are?

Hi,

I have question that is similar to the OP. I would like to do a LDLt factorization on a Mapped matrix, but the LDLt function seems to only like plain Eigen matrices, even though

const LDLT< typename MatrixBase< Derived >::PlainMatrixType > ldlt ( ) const [inline]

is specified in

http://eigen.tuxfamily.org/dox/classEig ... xBase.html

I think I am misunderstanding Benoit's recommendation.

Here is my function (I am using the devel branch)

Code: Select all

#include "/tmp/eigen/Eigen/Cholesky"

using namespace Eigen;

extern "C" {
    void foo(const double *A, const int *n, double *d) {
//      Map<MatrixXd> AA(A, n[0], n[0]);
//      I would like to do the above but only the below compiles
        int count = 0;
        MatrixXd AA(n[0], n[0]);
        for(int i=0; i<n[0]; i++) {
            for(int j=0; j<n[0]; j++) {
                AA(i,j) = A[count];
                count++;
            }
        }

//      LDLt method does not like Map, only Matrix
        VectorXd D = AA.ldlt().vectorD();
        for(int i=0; i<n[0]; i++) {
            d[i] = D(i);
        }
    }
}


and here is the compiler error I get when I try to do it with Map

Code: Select all

g++ -I/usr/share/R/include      -fpic  -g -O2 -c LDLt.cc -o LDLt.o
In file included from /tmp/eigen/Eigen/Cholesky:36,
                 from LDLt.cc:1:
/tmp/eigen/Eigen/src/Cholesky/LDLT.h: In member function ‘const Eigen::LDLT<typename Eigen::MatrixBase<Derived>::PlainObject> Eigen::MatrixBase<Derived>::ldlt() const [with Derived = Eigen::Map<Eigen::Matrix<double, 33331, 33331, 0, 33331, 33331>, 0>]’:
LDLt.cc:20:   instantiated from here
/tmp/eigen/Eigen/src/Cholesky/LDLT.h:325: error: conversion from ‘const Eigen::Map<Eigen::Matrix<double, 33331, 33331, 0, 33331, 33331>, 0>’ to non-scalar type ‘const Eigen::LDLT<Eigen::Matrix<double, 33331, 33331, 0, 33331, 33331> >’ requested
make: *** [LDLt.o] Error 1


indeed, this is a bug in the ldlt() function. Problem solved en both branches.

Also note that Map objects are writable, so you can do:

Map<VectorXd>(d,n[0]) = AA.ldlt().vectorD();

Hi Gael,

Thank you for fixing this compilation error. It is working now, apart from an additional bug. When the matrix whose LDLt decomposition is sought is positive semi-definite but rank-deficient, the vectorD() method (sometimes) puts 1.0 where 0.0 should be. Thus, when you multiply the all factors together, it does not equal the original matrix.

Here is an example of a 5x5 matrix whose rank is 2 where vectorD() produces the wrong answer, namely

Code: Select all

Here is D:
1
0.814626
1.3018e-15
7.04623e-16
1<- 1 is wrong

Here is the difference:
           0            0            0            0            0
           0            0            0            0            0
           0            0            0            0            0
           0            0            0            0            0
           0            0            0            0           -1<- -1 is wrong


Code: Select all

#include "/tmp/eigen/Eigen/Cholesky"
#include<iostream>

using namespace Eigen;
using namespace std;

extern "C" {

    int main() {
        MatrixXd AA(5,5);
        AA << 1.000000000000000, -0.430551015593539, -0.231799415339481, -0.608640054687203, -0.399959174714208,
             -0.430551015593539,  1.000000000000000, -0.778182118510235,  0.978188583485658,  0.999434535106194,
             -0.231799415339481, -0.778182118510235,  1.000000000000000, -0.630753429260956, -0.798859598981201,
             -0.608640054687203,  0.978188583485658, -0.630753429260956,  1.000000000000000,  0.970651002394712,
             -0.399959174714208,  0.999434535106194, -0.798859598981201,  0.970651002394712,  1.000000000000000;

        LDLT<MatrixXd> AA_factored = AA.ldlt();
        MatrixXd L = AA_factored.matrixL();

        cout << "Here is D:\n" << AA_factored.vectorD() << "<- 1 is wrong\n" << endl;


        cout << "Here is the difference:\n" << AA - (L * AA_factored.vectorD().asDiagonal() * L.transpose()) <<
                "<- -1 is wrong\n" << endl;
        return 0;
    }
}


Thanks for the great test case. Bug fixed, see the mailing list.