Using Eigen with DTSource

How to use Eigen with DTSource.

First download the Eigen distribution. Since this is a template library there is no source code to compile, just header files to include. However you need to make sure that Xcode will find the header files. Do that by going to the Build Settings, search for “Header Search Paths” and type in the full path name of the folder that contains the Eigen folder. Note that this path can not contain spaces or the compiler will not find it.

What is Eigen?

Eigen is a template library for solving linear systems, LU, QR and SVD factorizations. It handles both dense and sparse matrices.

Eigen uses its own matrix classes to handle computations. Since the default storage method for arrays matches what DTSource uses it is easy to convert back and forth and in some cases you can let Eigen share the array from DTSource so you don’t have to duplicate data.

Include Eigen in your code by adding the following line with the other include statements.

#include <Eigen/Dense>

LU factorization and solving Ax=b

As an example, consider the matrix A = [1 2;3 4] with right hand side b=[5;11] so the solution is x=[1;2]. The goal is to create the matrix A and right hand side b in DTSource, compute the solution x using Eigen and then copy it over to a DTSource object.

Create A

First create the matrix using a DTSource object.

    DTMutableDoubleArray Ad(2,2);
    Ad(0,0) = 1; Ad(0,1) = 2;
    Ad(1,0) = 3; Ad(1,1) = 4;

Now create the Eigen version of the matrix. Call this A.

Eigen::Map<const Eigen::MatrixXd,0,Eigen::OuterStride<> > A(Ad.Pointer(),Ad.m(),Ad.n(),Eigen::OuterStride<>(Ad.m()));

This needs to be broken down a little bit to explain how Eigen works. First of the Eigen:: prefix is used when you don’t specify the namespace. If you add the line “using namespace Eigen” to your source code you don’t need this and the statement becomes

Map<const MatrixXd,0,OuterStride<> > A(Ad.Pointer(),Ad.m(),Ad.n(),OuterStride<>(Ad.m()));

The Map<> template is used so that the matrix object shares the data pointer from Ad. You can create the matrix and then copy the data over, but then you have two copies of of the array and that is wasteful.

The MatrixXd means that this is a matrix of double. This is a shortcut for a much longer template definition, but look at the Eigen documentation for more information. OuterStride<> is needed because this is a two dimensional array and not a matrix. This is then used at the end to define the stride for the first dimension. This is because index (i,j) is stored at offset i + j*m, and in order for the map to know that you need to specify number of entries in each column.

Now Ad and A share the data. The Ad object owns the pointer, and you need to make sure that the Ad object lives as long or longer than the A object. Otherwise the destructor for the DTDoubleArray object will free the memory and the Eigen object might still access it.

Create b

Create the DTSource version of the right hand side.

    DTMutableDoubleArray bd(2);
    bd(0) = 5;
    bd(1) = 11;

Create the Eigen version of this array

Eigen::Map<Eigen::VectorXd> b(bd.Pointer(),bd.Length());

Note the difference betwen this and the A matrix. This is a vector and not a matrix, so you don’t need the stride information. This is also not a const object so you can change values of the b object and that will change the entries in the bd list.

Solve – #1

For the solution, first create the LU factorization and then solve the system. That is done with two lines

    auto LUFactor = A.fullPivLu();
    Eigen::MatrixXd x = LUFactor.solve(b);

The auto keyword is a handy addition to C++ that allows the compiler to automatically determine the type of the return variable based on what the fullPivLU() function returns. The first line computes the LU factorization and takes a lot longer time than the second line. Note that this does allocate memory for the L and U factors. This basically doubles the memory use because the factors have the same size as A.

The vector for x is also allocated. It is not possible to hand in an array for the solution. This is less of an issue because the memory use for x is a lot less than A. We do however need to extract this data into a DTDoubleArray if we want to return this data and use it in other functions that are based on DTSource. Do that with a simple memory copy

DTMutableDoubleArray xd(2);
memcpy(xd.Pointer(),x.data(),xd.Length()*sizeof(double));

Solve – #2

If you have 10,000 unknowns the size of A is 100 million values and the memory requirement is 800MB. You want to avoid allocating another 800MB when you compute the factorization. And the L and U factors can be saved in a single matrix because the diagonal entry of L is equal to 1 and you can store L in the strictly lower triangular part and U in the upper triangular part. To set this up in Eigen you change the definition a little bit. And in addition we show a different factorization, partial pivoting that only pivots rows and not rows and columns. Again we remove the Eigen:: prefix to make it shorter. You need to add this prefix to every entry that belongs to Eigen, i.e. the Map, OuterStride, PartialPivLU, Ref, MatrixXd

Map<MatrixXd,0,OuterStride<> > A(Ad.Pointer(),A.m(),A.n(),OuterStride<>(A.m()));
PartialPivLU<Ref<MatrixXd> > luFactor(A);
MatrixXd x = luFactor.solve(b);

Note the Ref<> wrapper around the MatrixXd type, and now the LU factorization is made by handing the A matrix into an object constructor.

This uses the data for the LU factors which means that the original value of A has been overwritten. This is often not a problem because your code might not need the matrix after you have computed the LU factor.