Are there any memory issues with this Eigen QR wrapper?

I am writing a wrapper to Eigen QR for my personal use and I am wondering if there are any memory leaks or undocumented behavior in my implementation, especially in the function get_QR().

The answer is as expected. This is related to my previous question here.

using std::cout;
using std::endl;
using namespace Eigen;

/*!
Obtains the QR decomposition as A=QR, where all the matrices are in Eigen MatrixXd format.
*/
void get_QR(MatrixXd A, MatrixXd& Q, MatrixXd& R) {

int m           =       A.rows();
int n           =       A.cols();
int minmn       =       min(m,n);

//      A_E     =       Q_E*R_E.
HouseholderQR<MatrixXd> qr(A);
Q  =   qr.householderQ()*(MatrixXd::Identity(m, minmn));
R  =   qr.matrixQR().block(0, 0, minmn, n).triangularView<Upper>();
}

/*!
Obtains the QR decomposition as A=QR, where all the matrices are in double format.
*/
void get_QR(double* A, int m, int n, double*& Q, double*& R) {
MatrixXd Q_E, R_E;
int minmn       =       min(m,n);

//      Maps the double to MatrixXd.
Map<MatrixXd> A_E(A, m, n);

get_QR(A_E, Q_E, R_E);

Q       =       (double*)realloc(Q_E.data(), m*minmn*sizeof(double));
R       =       (double*)realloc(R_E.data(), minmn*n*sizeof(double));
}

int main(int argc, char* argv[]) {
srand(time(NULL));

int m   =       atoi(argv[1]);
int n   =       atoi(argv[2]);

//      Check the double version.
int minmn       =       min(m,n);
double* A       =       (double*)malloc(m*n*sizeof(double));
double* Q       =       (double*)malloc(m*minmn*sizeof(double));
double* R       =       (double*)malloc(minmn*n*sizeof(double));

double RANDMAX  =       double(RAND_MAX);
//      Initialize A as a random matrix.
for (int index=0; index<m*n; ++index) {
A[index]        =       rand()/RANDMAX;
}

get_QR(A, m, n, Q, R);
std::cout << Q[0] << std::endl;

//      Check the MatrixXd version.
Map<MatrixXd> A_E(A, m, n);
MatrixXd Q_E, R_E;
get_QR(A_E, Q_E, R_E);

cout << Q[0] << endl;
cout << Q_E(0,0) << endl;

free(A);
free(Q);
free(R);
}


For instance, I get the output as

-0.360995
-0.360995
-0.360995

• You mean apart from using malloc in C++ code? – Martin York Mar 1 '14 at 18:19

First check that there are the enough command line parameters

// If there are not enough parameters
//
// If there are no parameters passed then
// Then argv[1] is NULL
//      argv[2] is random
int m   =       atoi(argv[1]);
int n   =       atoi(argv[2]);


If not enough parameters are passed then this can get you funny number in m and n which will affect how much memory you allocate.

Don't use malloc in C++ code.
You can not mix malloc/free with new/delete. So by only using one system of memory management you will not mix up the two different types of memory.

    double* A       =       (double*)malloc(m*n*sizeof(double));

// Also the C++ interface is simpler.

double* A       =       new double[m*n];


Secondly don't use pointer.
In modern C++ it is very rare to see RAW pointers (as they have no concept of ownership). And they are prone to leaking. you should be using either smart pointers or containers to hold your dynamically allocated memory.

In this case a container would have been better.

   std::vector<double> A(m*n);


All memory management now handled correctly.
This can then be used in your C code with

   double* Ap = &A[0];  // get the address of the first element.


This is not going to go well:

    Q       =       (double*)realloc(Q_E.data(), m*minmn*sizeof(double));

// Q_E.data() better not only return you a pointer allocated with malloc.
//            But it better set the internal pointers to NULL or some
//            other safe value. As soon as the realloc() finishes the
//            old pointer is no longer valid.


Also this is not the correct way to use realloc.

  size_t    size = /* Calc Size */;
double*   data = (double*)malloc(size);
...
size_t    newSiz = /* New Space needed */;
double*   resize = (double*)realloc(data, newSize);
if (resize != NULL)
{
data = resize;  // if realloc() returns NULL the original array
size = newSiz;  // is unaffected. You should check this after
}                   // the call and only reassign the original
// value if the resize worked correctly.


Major Points:

I don't know what is happening inside MatrixXd but it looks like it is exposing its memory management outside the class. That is not a good idea. An object should know how to handle its memory in all situation (and do its own resizing as appropriate).

Two don't use RAW pointers (use C++ ownership symantics).

Learn the Rule of Three/Five so that objects manage their own memory.

Don't write C code and definately don't write C memory management code.

• A reason why I do not want to use vector is it actually slows down the computation by a few factors. That is why I want to do all the memory management myself and also give the user some responsibility on the memory management front. – John Smith Mar 2 '14 at 17:29
• @LeslieFaulkner: I don't believe that. We have shown several times on SO that vector is as fast (and faster when you use realloc) than manual memory management. If there is a slow down then it is because you doing something non idiomatic and causing unnecessary copies. If you want to publish your version with std::vector we can show you where you are going wrong. – Martin York Mar 2 '14 at 17:55
• Anyway fast and broken is not really as useful as slow and correct. – Martin York Mar 2 '14 at 17:58

@Loki Astari has hit on all the important points, but I have just a few things to add:

• Make sure you're including all of the libraries here, especially if someone wants to test this code. You're missing <iostream>, <ctime>, and <cstdlib>.

• As Loki has mentioned, you're using a lot of C code in this C++ program. Another instance of this is some of your casting.

Instead of casting the C way (for this instance):

double RANDMAX  =       double(RAND_MAX);


cast the C++ way:

double RANDMAX  =       static_cast<double>(RAND_MAX);

• Your use of whitespace for assignments is a little unusual. Although it is okay to align a list of initialized variables (at least when you must have a list), I don't see the need in other instance such as this:

for (int index=0; index<m*n; ++index) {
A[index]        =       rand()/RANDMAX;
}


That makes the loop harder to read, and there's just no point to it. In addition, you should also add a space between the operator and operands in the loop statement.

You'll then have this:

for (int index = 0; index < m*n; ++index) {
A[index] = rand() / RANDMAX;
}

• In get_QR():

void get_QR(double* A, int m, int n, double*& Q, double*& R) {}


I'd group the arguments and parameters by types, making it less-likely to mismatch the arguments when making the function call. Perhaps something like this:

void get_QR(double* A, double*& Q, double*& R, int m, int n) {}


I am wondering if there is any memory leak or undocumented behavior in my implementation

I don't see a leak, but it's difficult to confirm there isn't one.

This statement is dubious:

Q       =       (double*)realloc(Q_E.data(), m*minmn*sizeof(double));


Q       =       (double*)realloc(Q, m*minmn*sizeof(double));


... then it would have been less dubious; however if realloc fails it will return NULL, therefore overwriting a value with realloc can leak the previously-allocated value, so even that should be more like the following using a temp variable to test the result of realloc before using it to overwrite the previous value ...

double* temp       =       (double*)realloc(Q, m*minmn*sizeof(double));
if (!temp) {
... handle error here and don't forget to free Q ...
return;
}
// else realloc successful
Q = temp;


However in your code I don't see how your get_QR function can be making use of the pre-allocated, passed-in double* Q and double* R arrays: and therefore I think that they're being leaked; maybe there wouldn't be a leak if you didn't allocate/initialize Q and R in main.

Depending on which compiler you're using, there may be built-in tools to detect memory leaks: e.g. the Debug Heap for Microsoft's compiler); or for Linux I recommend using valgrind.