[eigen] Dynamic allocations in SelfAdjointEigenSolver and HouseholderQR

[Ricard Marxer <ricardmp@xxxxxxxxx>]

Is it possible to avoid dynamic allocations in
SelfAdjointEigenSolver.compute() and HouseholderQR.solve(), when using
dynamic size matrices?

In the examples the matrices are small and their size is known at
compile time.  But in my use case they are dynamic and their size is
only known at runtime.

To try the attachment tests simply run:
g++ -I../eigen dynalloc_selfadjoint.cpp -o dynalloc_selfadjoint
&&./dynalloc_selfadjointg++ -I../eigen dynalloc_qr.cpp -o dynalloc_qr
&&  ./dynalloc_qr
The backtraces show the following:
(gdb) bt#0  0x00007ffff72ef405 in *__GI_raise (sig=<optimized out>)
at../nptl/sysdeps/unix/sysv/linux/raise.c:64#1  0x00007ffff72f2680 in
*__GI_abort () at abort.c:92#2  0x00007ffff72e85b1 in
*__GI___assert_fail (   assertion=0x4247c8 "is_malloc_allowed() &&
\"heap allocation isforbidden (EIGEN_RUNTIME_NO_MALLOC is defined and
g_is_malloc_allowedis false)\"", file=<optimized out>, line=197,
function=0x425500
"voidEigen::internal::check_that_malloc_is_allowed()") at
assert.c:81#3  0x0000000000402b7c in
Eigen::internal::check_that_malloc_is_allowed() ()#4
0x0000000000402b8f in Eigen::internal::aligned_malloc(unsigned long)
()#5  0x0000000000409092 in
void*Eigen::internal::conditional_aligned_malloc<true>(unsigned long)
()#6  0x0000000000406b0b in
float*Eigen::internal::conditional_aligned_new_auto<float,
true>(unsignedlong) ()#7  0x00000000004069e8 in
Eigen::DenseStorage<float, -1, -1, 1,0>::DenseStorage(long, long,
long) ()#8  0x0000000000405796 in
Eigen::PlainObjectBase<Eigen::Matrix<float,-1, 1, 0, -1, 1>
>::PlainObjectBase(long, long, long) ()#9  0x0000000000404415 in
Eigen::Matrix<float, -1, 1, 0, -1, 1>::Matrix(long) ()#10
0x0000000000405e3d in
voidEigen::internal::tridiagonalization_inplace_selector<Eigen::Matrix<float,-1,
-1, 0, -1, -1>, -1, false>::run<Eigen::Matrix<float, -1, 1, 0, -1,1>,
Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::Matrix<float, -1,-1, 0,
-1, -1>&, Eigen::Matrix<float, -1, 1, 0, -1, 1>&,Eigen::Matrix<float,
-1, 1, 0, -1, 1>&, bool) ()#11 0x0000000000404b1e in
voidEigen::internal::tridiagonalization_inplace<Eigen::Matrix<float,
-1,-1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1,
1>,Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::Matrix<float, -1,
-1,0, -1, -1>&, Eigen::Matrix<float, -1, 1, 0, -1,
1>&,Eigen::Matrix<float, -1, 1, 0, -1, 1>&, bool) ()#12
0x0000000000403926
inEigen::SelfAdjointEigenSolver<Eigen::Matrix<float, -1, -1, 0, -1,
-1>>::compute(Eigen::Matrix<float, -1, -1, 0, -1, -1> const&, int)
()#13 0x000000000040142a in main ()
and
(gdb) bt#0  0x00007ffff72ef405 in *__GI_raise (sig=<optimized out>)
at../nptl/sysdeps/unix/sysv/linux/raise.c:64#1  0x00007ffff72f2680 in
*__GI_abort () at abort.c:92#2  0x00007ffff72e85b1 in
*__GI___assert_fail (   assertion=0x42dcc8 "is_malloc_allowed() &&
\"heap allocation isforbidden (EIGEN_RUNTIME_NO_MALLOC is defined and
g_is_malloc_allowedis false)\"", file=<optimized out>, line=197,
function=0x42ecc0
"voidEigen::internal::check_that_malloc_is_allowed()") at
assert.c:81#3  0x000000000040180a in
Eigen::internal::check_that_malloc_is_allowed() ()#4
0x000000000040181d in Eigen::internal::aligned_malloc(unsigned long)
()#5  0x0000000000406576 in
void*Eigen::internal::conditional_aligned_malloc<true>(unsigned long)
()#6  0x00000000004050ed in
float*Eigen::internal::conditional_aligned_new_auto<float,
true>(unsignedlong) ()#7  0x00000000004043b0 in
Eigen::DenseStorage<float, -1, -1, 1,0>::DenseStorage(long, long,
long) ()#8  0x0000000000403564 in
Eigen::PlainObjectBase<Eigen::Matrix<float,-1, 1, 0, -1, 1>
>::PlainObjectBase(long, long, long) ()#9  0x0000000000402c51 in
Eigen::Matrix<float, -1, 1, 0, -1,1>::Matrix(Eigen::Matrix<float, -1,
1, 0, -1, 1> const&) ()#10 0x000000000040a775 in
voidEigen::internal::solve_retval<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1,
1>>::evalTo<Eigen::Matrix<float, -1, 1, 0, -1, 1>
>(Eigen::Matrix<float,-1, 1, 0, -1, 1>&) const ()#11
0x0000000000408505 in
voidEigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1,
1>>::evalTo<Eigen::Matrix<float, -1, 1, 0, -1, 1>
>(Eigen::Matrix<float,-1, 1, 0, -1, 1>&) const ()#12
0x00000000004063e1 in
voidEigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1, 1>
>>::evalTo<Eigen::Matrix<float, -1, 1, 0, -1, 1>
>(Eigen::Matrix<float,-1, 1, 0, -1, 1>&) const ()#13
0x0000000000404d9e in Eigen::Matrix<float, -1, 1, 0, -1,
1>&Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1,
1>>::operator=<Eigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1, 1>
>>(Eigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1, 1> > >const&)
()#14 0x0000000000403a8c in Eigen::Matrix<float, -1, 1, 0, -1,
1>&Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1,
1>>::operator=<Eigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1, 1>
>>(Eigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1, 1> > >const&)
()#15 0x0000000000402e3d in Eigen::Matrix<float, -1, 1, 0, -1,
1>&Eigen::Matrix<float, -1, 1, 0,
-1,1>::operator=<Eigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1, 1>
>>(Eigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::HouseholderQR<Eigen::Matrix<float,-1,
-1, 0, -1, -1> >, Eigen::Matrix<float, -1, 1, 0, -1, 1> > >const&)
()#16 0x000000000040148d in main ()

Thanks
-- 
ricard
http://twitter.com/ricardmp
http://www.ricardmarxer.com
http://www.caligraft.com

#define EIGEN_RUNTIME_NO_MALLOC

#include <Eigen/Eigen>
#include <Eigen/QR>
#include <iostream>

using namespace Eigen;
using namespace std;

int main(int, char *[])
{ 
    HouseholderQR<MatrixXf> es(4, 4);
    MatrixXf A = MatrixXf::Random(4,4);
    VectorXf y = VectorXf::Random(4);
    VectorXf p(4);
    
    Eigen::internal::set_is_malloc_allowed(false);
    es.compute(A);
    Eigen::internal::set_is_malloc_allowed(true);
    
    cout << "Performed QR decomposition." << endl;

    Eigen::internal::set_is_malloc_allowed(false);
    p = es.solve(y);
    Eigen::internal::set_is_malloc_allowed(true);

    cout << "The solution to Ap=y is p: " << p << endl;
}

#define EIGEN_RUNTIME_NO_MALLOC

#include <Eigen/Eigen>
#include <Eigen/Eigenvalues>
#include <iostream>

using namespace Eigen;
using namespace std;

int main(int, char *[])
{ 
    SelfAdjointEigenSolver<MatrixXf> es(4);
    MatrixXf X = MatrixXf::Random(4,4);
    MatrixXf A = X + X.transpose();
    
    Eigen::internal::set_is_malloc_allowed(false);
    es.compute(A);
    Eigen::internal::set_is_malloc_allowed(true);
    
    cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl;

    Eigen::internal::set_is_malloc_allowed(false);
    es.compute(A + MatrixXf::Identity(4,4)); // re-use es to compute eigenvalues of A+I
    Eigen::internal::set_is_malloc_allowed(true);
    
    cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl;
}