0

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.

3
  • They seem to be some floating points going towards zero, which often happens in some numerical algorithms. You could e.g. threshold to have your own interpretation of what's zero.
    – keyser
    Jan 12, 2015 at 15:36
  • I have not done any algorithms yet, just converted from complex to MatrixXcd in the Map function. Jan 12, 2015 at 15:39
  • Good point. Look at how that constructor makes us of the matrix data then.
    – keyser
    Jan 12, 2015 at 15:42

1 Answer 1

2

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.

1
  • Thanks alot. This helps. I cannot use Eigen directly that would force everyone to use Eigen and they don't want to. Jan 12, 2015 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.