CODEWITHSUNDEEP
c++matrixeigen
user1876942
user1876942

Reputation: 1501

Mapping a complex matrix with Eigen

I have the following matrix:

int N = 3;
complex<double>* complexvector = new complex<double>[N];
for(int i=0; i<N; i++)
        {
            complexvector[i]=complex<double>(i+1,i+1);
        }



complex<double>** complexMatrix = new complex<double>*[N];
    for(int i=0; i<N; i++)
    {
        complexMatrix[i] = new complex<double>[N];// allocate an array to each pointer
    }



complexMatrix[0][0] = complex<double>(1,1);
complexMatrix[0][1] = complex<double>(2,2);
complexMatrix[0][2] = complex<double>(3,3);

complexMatrix[1][0] = complex<double>(4,4);
complexMatrix[1][1] = complex<double>(5,5);
complexMatrix[1][2] = complex<double>(6,6);

complexMatrix[2][0] = complex<double>(7,7);
complexMatrix[2][1] = complex<double>(8,8);
complexMatrix[2][2] = complex<double>(9,9);


complex<double>* returnvector = new complex<double>[N];
returnvector = preformcholesky(complexMatrix, complexvector, N);

And I do the Cholesky here:

complex<double>* preformcholesky(complex<double>** matrix, complex<double>* vector, 
int  size)
{
std::cout << "start cholesky" << std::endl;
Map<MatrixXcd, RowMajor> mat(*matrix,size,size);
Map<VectorXcd> vec(vector,size);
printMatrix(matrix,size);
std::cout << "mat" << mat << std::endl;
std::cout << "vec" << vec << std::endl;
//mat.llt();
VectorXcd x = mat.llt().solve(vec);
std::cout << "coeff" << x << std::endl;
return x.data();
}

The trouble is when printing out "mat", it goes like this:

mat           (1,1) (0,3.21143e-322)            (6,6)
       (2,2)            (4,4) (0,3.21143e-322)
       (3,3)            (5,5)            (7,7)

Where do those (0,3.21143e-322) come from? If I print "matrix" then that is OK, so I think that the map has gone wrong. If you need some more code then just let me know. I am new to Eigen so maybe some basic mistake. I am using C++ and Linux.

Upvotes: 0

Views: 1781

Answers (1)

Wintermute
Wintermute

Reputation: 44023

There are several things wrong with that code. The return value of performcholesky, for example, is a pointer to a local variable that does not exist anymore after performcholesky has returned, and you overwrite a pointer to allocated memory with this dangling pointer, leaking that memory. Why don't you work with Eigen::MatrixXcd directly? You'd save yourself a world of pain debugging all that later.

Anyway, the problem you ran into this time is that the Eigen::Map constructor expects a flat array, not a pointer to pointer. It is to be used like this:

// It would also work with naked new, but naked new leads to pain, fear, anger,
// and suffering. Also to long hours of debugging.
std::vector<std::complex<double> > data(N * N);

for(int i = 0; i < N * N; ++i) {
  data[i] = std::complex<double>(i, i);
}

Eigen::Map<Eigen::MatrixXcd, Eigen::RowMajor> mat(&data[0], N, N);

An interesting question is why the two arrays that you did not pass to the constructor nevertheless show up in the map, if somewhat shifted. The answer is: By pure happenstance, new gave you three slabs of memory that were right behind each other on the heap. This also explains the strange artifacts you see in the map; those are bookkeeping data of your heap implementation for the corresponding allocated blocks, reinterpreted as doubles. In essence, the memory around complexmatrix[0] looks like this:

  the area that Eigen::Map uses
+-----------------------------------------------------------------+
|(1,1)|(2,2)|(3,3)|bookkeeping|(4,4)|(5,5)|(6,6)|bookkeeping|(7,7)|(8,8)|(9,9)|
+-----------------+           +-----------------+           +-----------------+
 ^-complexmatrix[0]            ^-complexmatrix[1]            ^-complexmatrix[2]

Naturally, this is not behavior you can depend upon. The code could just as well show any random data or crash with a segfault if the blocks new gives you are arranged in a different manner.

Upvotes: 2

Related Questions