* bug #206: correctly forward computationOptions and work towards avoiding mallocs after preallocation, with unit test.

* added EIGEN_RUNTIME_NO_MALLOC and new set_is_malloc_allowed() function to implement that test
This commit is contained in:
Benoit Jacob 2011-03-06 20:59:25 -05:00
parent b464fc19bc
commit bfcad536e8
3 changed files with 114 additions and 11 deletions

View File

@ -167,14 +167,36 @@ inline void* generic_aligned_realloc(void* ptr, size_t size, size_t old_size)
*** Implementation of portable aligned versions of malloc/free/realloc ***
*****************************************************************************/
#ifdef EIGEN_NO_MALLOC
inline void check_that_malloc_is_allowed()
{
eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
}
#elif defined EIGEN_RUNTIME_NO_MALLOC
inline bool is_malloc_allowed_impl(bool update, bool new_value = false)
{
static bool value = true;
if (update == 1)
value = new_value;
return value;
}
inline bool is_malloc_allowed() { return is_malloc_allowed_impl(false); }
inline bool set_is_malloc_allowed(bool new_value) { return is_malloc_allowed_impl(true, new_value); }
inline void check_that_malloc_is_allowed()
{
eigen_assert(is_malloc_allowed() && "heap allocation is forbidden (EIGEN_RUNTIME_NO_MALLOC is defined and g_is_malloc_allowed is false)");
}
#else
inline void check_that_malloc_is_allowed()
{}
#endif
/** \internal Allocates \a size bytes. The returned pointer is guaranteed to have 16 bytes alignment.
* On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown.
*/
inline void* aligned_malloc(size_t size)
{
#ifdef EIGEN_NO_MALLOC
eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
#endif
check_that_malloc_is_allowed();
void *result;
#if !EIGEN_ALIGN
@ -268,9 +290,7 @@ template<bool Align> inline void* conditional_aligned_malloc(size_t size)
template<> inline void* conditional_aligned_malloc<false>(size_t size)
{
#ifdef EIGEN_NO_MALLOC
eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
#endif
check_that_malloc_is_allowed();
void *result = std::malloc(size);
#ifdef EIGEN_EXCEPTIONS

View File

@ -376,7 +376,12 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
* The default constructor is useful in cases in which the user intends to
* perform decompositions via JacobiSVD::compute(const MatrixType&).
*/
JacobiSVD() : m_isInitialized(false) {}
JacobiSVD()
: m_isInitialized(false),
m_isAllocated(false),
m_computationOptions(0),
m_rows(-1), m_cols(-1)
{}
/** \brief Default Constructor with memory preallocation
@ -386,6 +391,10 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
* \sa JacobiSVD()
*/
JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
: m_isInitialized(false),
m_isAllocated(false),
m_computationOptions(0),
m_rows(-1), m_cols(-1)
{
allocate(rows, cols, computationOptions);
}
@ -401,11 +410,15 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
* available with the (non-default) FullPivHouseholderQR preconditioner.
*/
JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
: m_isInitialized(false),
m_isAllocated(false),
m_computationOptions(0),
m_rows(-1), m_cols(-1)
{
compute(matrix, computationOptions);
}
/** \brief Method performing the decomposition of given matrix.
/** \brief Method performing the decomposition of given matrix using custom options.
*
* \param matrix the matrix to decompose
* \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
@ -415,7 +428,18 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
* Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
* available with the (non-default) FullPivHouseholderQR preconditioner.
*/
JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0);
JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
/** \brief Method performing the decomposition of given matrix using current options.
*
* \param matrix the matrix to decompose
*
* This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int).
*/
JacobiSVD& compute(const MatrixType& matrix)
{
return compute(matrix, m_computationOptions);
}
/** \returns the \a U matrix.
*
@ -494,16 +518,17 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
inline Index cols() const { return m_cols; }
private:
void allocate(Index rows, Index cols, unsigned int computationOptions = 0);
void allocate(Index rows, Index cols, unsigned int computationOptions);
protected:
MatrixUType m_matrixU;
MatrixVType m_matrixV;
SingularValuesType m_singularValues;
WorkMatrixType m_workMatrix;
bool m_isInitialized;
bool m_isInitialized, m_isAllocated;
bool m_computeFullU, m_computeThinU;
bool m_computeFullV, m_computeThinV;
unsigned int m_computationOptions;
Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
@ -515,9 +540,21 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
template<typename MatrixType, int QRPreconditioner>
void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
{
eigen_assert(rows >= 0 && cols >= 0);
if (m_isAllocated &&
rows == m_rows &&
cols == m_cols &&
computationOptions == m_computationOptions)
{
return;
}
m_rows = rows;
m_cols = cols;
m_isInitialized = false;
m_isAllocated = true;
m_computationOptions = computationOptions;
m_computeFullU = (computationOptions & ComputeFullU) != 0;
m_computeThinU = (computationOptions & ComputeThinU) != 0;
m_computeFullV = (computationOptions & ComputeFullV) != 0;

View File

@ -23,6 +23,9 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
// discard stack allocation as that too bypasses malloc
#define EIGEN_STACK_ALLOCATION_LIMIT 0
#define EIGEN_RUNTIME_NO_MALLOC
#include "main.h"
#include <Eigen/SVD>
@ -241,6 +244,46 @@ void jacobisvd_inf_nan()
svd.compute(m, ComputeFullU | ComputeFullV);
}
void jacobisvd_preallocate()
{
Vector3f v(3.f, 2.f, 1.f);
MatrixXf m = v.asDiagonal();
internal::set_is_malloc_allowed(false);
VERIFY_RAISES_ASSERT(VectorXf v(10);)
JacobiSVD<MatrixXf> svd;
internal::set_is_malloc_allowed(true);
svd.compute(m);
VERIFY_IS_APPROX(svd.singularValues(), v);
JacobiSVD<MatrixXf> svd2(3,3);
internal::set_is_malloc_allowed(false);
svd2.compute(m);
internal::set_is_malloc_allowed(true);
VERIFY_IS_APPROX(svd2.singularValues(), v);
VERIFY_RAISES_ASSERT(svd2.matrixU());
VERIFY_RAISES_ASSERT(svd2.matrixV());
svd2.compute(m, ComputeFullU | ComputeFullV);
VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
internal::set_is_malloc_allowed(false);
svd2.compute(m);
internal::set_is_malloc_allowed(true);
JacobiSVD<MatrixXf> svd3(3,3,ComputeFullU|ComputeFullV);
internal::set_is_malloc_allowed(false);
svd2.compute(m);
internal::set_is_malloc_allowed(true);
VERIFY_IS_APPROX(svd2.singularValues(), v);
VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
internal::set_is_malloc_allowed(false);
svd2.compute(m, ComputeFullU|ComputeFullV);
internal::set_is_malloc_allowed(true);
}
void test_jacobisvd()
{
CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) ));
@ -290,4 +333,7 @@ void test_jacobisvd()
// Test problem size constructors
CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) );
// Check that preallocation avoids subsequent mallocs
CALL_SUBTEST_9( jacobisvd_preallocate() );
}