0

I am writing a wrapper to Eigen QR for my personal use and I am wondering if there is any memory leak or undocumented behavior in my implementation especially in the function void get_QR(double* A, int m, int n, double*& Q, double*& R). 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
Community
  • 1
  • 1
John Smith
  • 161
  • 1
  • 6

2 Answers2

1

realloc is not guaranteed to be compatible with operator new[] and operator delete[], only with malloc and free.

I would just replace all your double* pointers with std::vector<double> objects.

aschepler
  • 70,891
  • 9
  • 107
  • 161
  • However there is placement new – Sebastian Hoffmann Mar 01 '14 at 16:28
  • Yes, right. I have got ride of the `new` now and replaced with `malloc`. I do not want to use `vector`. If so, is my implementation, leak free and undocumented behavior free. – John Smith Mar 01 '14 at 16:29
  • `I do not want to use vector` Why not? – PaulMcKenzie Mar 01 '14 at 16:31
  • @PaulMcKenzie Because this will be used by other codes and those rely purely on double* to denote array or matrices. – John Smith Mar 01 '14 at 16:33
  • 1
    @LeslieFaulkner - a vector is compatible with an array by specifying the address of the first element. So I see no reason why you can't use vector. – PaulMcKenzie Mar 01 '14 at 16:34
  • @PaulMcKenzie Thanks, that would be useful, though could you elaborate on it. For instance, if a user passes an array to a function, how would I use a `vector` without copying the array contents (i.e., just by manipulating the memory of the first element)? Thanks. – John Smith Mar 01 '14 at 16:37
  • @LeslieFaulkner - Yes. A vector basically takes what you did and puts it into a class. That's basically all it does. All you need to learn is how to get into the internals of this class, and taking the address of the first element is the way to do it. – PaulMcKenzie Mar 01 '14 at 16:39
  • @PaulMcKenzie Thanks. For instance, if I have a pointer `a` how should I put it in a vector? Can you add that one line of code? Thanks. – John Smith Mar 01 '14 at 16:41
  • @PaulMcKenzie Could you also answer my original question, as to if my code is right or has some memory leaks or UB? Thanks again and sorry to bug you repeatedly. – John Smith Mar 01 '14 at 16:42
  • @LeslieFaulkner - See my answer. Your original code does not check if malloc() or realloc() returns NULL. – PaulMcKenzie Mar 01 '14 at 16:45
1

You can replace your calls to malloc and free with a vector. Here is sample code of how you would use vector with existing code that takes a pointer to double.

#include <vector>
//...
void foo(double* d, int size)
{
   d[0] = 10;
}

int main()
{
   std::vector<double> d(10);  // 10 doubles
   foo(&d[0], d.size());
}

The foo() function is an "old C" function that expects a pointer to double to denote a double array. Note that to do this magic, you pass the address of the first element.

So your original code could look like this, given the above:

#include <vector>
#include <algorithm>

typedef std::vector<double> VectDouble;

void get_QR(VectDouble& A, int m, int n, VectDouble& Q, VectDouble& R) 
{
    MatrixXd Q_E, R_E;
    int minmn       =       min(m,n);

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

    get_QR(A_E, Q_E, R_E);

    Q.resize(m * minmn);
    R.resize(minmn *n, 0);
}

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

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

    int minmn       =       min(m,n);
    VectDouble A(m*n);
    VectDouble Q(m*minmn);
    VectDouble R(minmn*n);

    double RANDMAX  =       double(RAND_MAX);
    //      Initialize A as a random matrix.
    std::fill(A.begin(), A.end(), rand()/RANDMAX);

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

    Map<MatrixXd> A_E(&A[0], 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;
}
PaulMcKenzie
  • 34,698
  • 4
  • 24
  • 45