From 43e53b0bbc44f987de4504bea130112b13100b07 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Thu, 11 Aug 2022 16:02:15 +0200
Subject: [PATCH] optional threshold for sparsemat::simplify

---
 sparsemat.cc | 4 ++--
 sparsemat.h  | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/sparsemat.cc b/sparsemat.cc
index 6507ff6..8966488 100644
--- a/sparsemat.cc
+++ b/sparsemat.cc
@@ -315,7 +315,7 @@ return nonzero; //number of (in principle) nonzero elements
 
 
 template <class T>
-void SparseMat<T>::simplify()
+void SparseMat<T>::simplify(const double sparseepsilon)
 {
 unsigned int n;
 if(!list) return;
@@ -352,7 +352,7 @@ for(i=1; i<n;i++)
 
 //check if  summed to zero
 for(i=0; i<n;i++) if(rowsorted[i] &&
-	std::abs(rowsorted[i]->elem)<=SPARSEEPSILON
+	std::abs(rowsorted[i]->elem)<=sparseepsilon
 	)	 {delete rowsorted[i]; rowsorted[i]=NULL;}
 
 //restore connectivity
diff --git a/sparsemat.h b/sparsemat.h
index 26c89cd..6b83567 100644
--- a/sparsemat.h
+++ b/sparsemat.h
@@ -23,7 +23,7 @@ namespace LA {
 
 //threshold for neglecting elements, if not defined, no tests are done except exact zero test in simplify - might be even faster
 //seems to perform better with a threshold, in spite of abs() tests
-const double SPARSEEPSILON=1e-35;
+const double SPARSEEPSILON=1e-15;
 
 typedef unsigned int SPMatindex;
 typedef int SPMatindexdiff; //more clear would be to use traits
@@ -140,7 +140,7 @@ public:
 	unsigned int length() const;
 	void copyonwrite(bool detachonly=false);
 	void clear() {copyonwrite(LA_traits<T>::is_plaindata()); if(count) {delete count; count=NULL;}}
-	void simplify();
+	void simplify(const double sparseepsilon=SPARSEEPSILON);
 	const T trace() const;
 	const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first
 	unsigned int sort(int type) const;