Move some Cholmodwrapper into cpp file and changed approx equal (of floats) for more consitency

include/xerus/cholmod_wrapper.h

@@ -39,24 +39,6 @@ namespace xerus {
 	namespace internal {
 		
 		
-		///@brief wrapper object for the cholmod_common struct to automatically call the constructor and destructor
-		struct CholmodCommon final {
-			struct RestrictedAccess final {
-				cholmod_common* const c;
-				std::mutex &lock;
-				RestrictedAccess(cholmod_common* const _c, std::mutex &_lock);
-				operator cholmod_common*() const;
-				~RestrictedAccess();
-			};
-			
-			std::unique_ptr<cholmod_common> c;
-			std::mutex lock;
-			CholmodCommon();
-			~CholmodCommon();
-			RestrictedAccess get();
-			std::function<void(cholmod_sparse*)> get_deleter();
-		};
-		
 		///@brief wrapper class for the cholmod sparse matrix objects
 		class CholmodSparse final {
 		public:
@@ -137,8 +119,5 @@ namespace xerus {
 				bool _fullrank
 			);
 		};
-		
-		
-		extern thread_local CholmodCommon cholmodObject;
 	}
 }

include/xerus/misc/math.h

@@ -46,6 +46,7 @@ namespace xerus {
 		
 		
 		///@brief: Calculates _base^_exp by binary exponentiation
+		///@note @a _exp is assumed to be a built-in integral type and thus taken by value, while @a _base may be a matrix, rational number etc. and is thus taken as reference
 		template<class T, class I, typename std::enable_if<std::is_integral<I>::value && std::is_unsigned<I>::value, bool>::type = true> 
 		constexpr T pow(const T &_base, const I _exp) noexcept {
 			return _exp==0 ? T(1) : (_exp%2==0 ? pow(T(_base*_base), _exp/2) : T(_base*pow(_base, _exp-1)));
@@ -55,11 +56,11 @@ namespace xerus {
 		///@brief: Calculates _base^_exp by binary exponentiation
 		template<class T, class I, typename std::enable_if<std::is_integral<I>::value && !std::is_unsigned<I>::value, bool>::type = true> 
 		constexpr T pow(const T &_base, const I _exp) noexcept {
-			return _exp==0 ? 
+			return _exp==0 ?
 						T(1) :
 						(
 							_exp<0 ? 
-								T( 1/pow(_base, -_exp)) 
+								T( 1/pow(_base, -_exp))
 							: 
 								( _exp%2==0 ? pow(T(_base*_base), _exp/2) : T(_base*pow(_base, _exp-1)) )
 						);
@@ -69,7 +70,7 @@ namespace xerus {
 		///@brief: Checks whether the relative difference between @a _a and @a _b (i.e. |a-b|/(|a|/2+|b|/2)) is smaller than @a _eps.
 		template<class T, typename std::enable_if<std::is_floating_point<T>::value, bool>::type = true>
 		bool approx_equal(T _a, T _b, T _eps = 4*std::numeric_limits<T>::epsilon()) noexcept {
-			return std::abs(_a-_b) <= _eps*0.5*(std::abs(_a)+std::abs(_b));
+			return std::abs(_a-_b) <= _eps*0.5*(std::abs(_a+_b)) || std::abs(_a-_b) < std::numeric_limits<T>::min();
 		}
 		
 		

src/xerus/cholmod_wrapper.cpp

@@ -33,6 +33,27 @@
 
 namespace xerus { namespace internal {
 	
+	
+		
+	///@brief Wrapper object for the cholmod_common struct to automatically call the constructor and destructor (cholmod MAY use its own allocator)
+	/// @note any CholmodSparse object stored outside (eg in Tensor) needs the corresponding deleter function defined here to be able to free it (regardles of the thread freeing it)
+	struct CholmodCommon final {
+		struct RestrictedAccess final {
+			cholmod_common* const c;
+			std::mutex &lock;
+			RestrictedAccess(cholmod_common* const _c, std::mutex &_lock);
+			operator cholmod_common*() const;
+			~RestrictedAccess();
+		};
+		
+		std::unique_ptr<cholmod_common> c;
+		std::mutex lock;
+		CholmodCommon();
+		~CholmodCommon();
+		RestrictedAccess get();
+		std::function<void(cholmod_sparse*)> get_deleter();
+	};
+	
 	CholmodCommon::RestrictedAccess::RestrictedAccess(cholmod_common* const _c, std::mutex& _lock) 
 		: c(_c), lock(_lock)
 	{
@@ -73,7 +94,7 @@ namespace xerus { namespace internal {
 	}
 
 	std::function<void(cholmod_sparse*)> CholmodCommon::get_deleter() {
-		return [&](cholmod_sparse* _toDelete) {
+		return [this](cholmod_sparse* _toDelete) {
 			cholmod_l_free_sparse(&_toDelete, this->get());
 		};
 	}