From a339af5c09039c4ae7dca71b59b4f39139ea146c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Wed, 6 Mar 2019 17:50:11 +0000
Subject: [PATCH] Roundoff error cleaning for zeros.

Former-commit-id: 1efdd13b64022a93ae396bbc92345917de014ed8
---
 qpms/symmetries.c | 25 ++++++++++++++++++++++++-
 qpms/symmetries.h | 17 +++++++++++++++++
 2 files changed, 41 insertions(+), 1 deletion(-)

diff --git a/qpms/symmetries.c b/qpms/symmetries.c
index b418a83..3ee2826 100644
--- a/qpms/symmetries.c
+++ b/qpms/symmetries.c
@@ -207,4 +207,27 @@ complex double *qpms_irot3_uvswfi_dense(
   return target;
 }
 
- 
+size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol) {
+  size_t changed = 0;
+  for(size_t i = 0; i < nmemb; ++i)
+    if(fabs(arr[i]) <= atol) {
+      arr[i] = 0;
+      ++changed;
+    }
+  return changed;
+}
+
+size_t qpms_czero_roundoff_clean(complex double *arr, size_t nmemb, double atol) {
+  size_t changed = 0;
+  for(size_t i = 0; i < nmemb; ++i) {
+    if(fabs(creal(arr[i])) <= atol) {
+      arr[i] = I*cimag(arr[i]);
+      ++changed;
+    }
+    if(fabs(cimag(arr[i])) <= atol) {
+      arr[i] = creal(arr[i]);
+      ++changed;
+    }
+  }
+}
+
diff --git a/qpms/symmetries.h b/qpms/symmetries.h
index 9e4266a..97e15cd 100644
--- a/qpms/symmetries.h
+++ b/qpms/symmetries.h
@@ -58,4 +58,21 @@ complex double *qpms_irot3_uvswfi_dense(
 		const qpms_vswf_set_spec_t *bspec,
 		const qpms_irot3_t transf);
 
+// Some auxilliary functions for "numerical cleaning" of the roundoff error
+
+/// Cleans the roundoff error of an array.
+/**
+ * Replaces all the array members \a x for which |\a x| <= \a atol
+ * with exact zeros.
+ *
+ * Returns the count of modified array members.
+ */
+size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol);
+
+/// Cleans the roundoff error of an array.
+/**
+ * Works on real and imaginary parts separately.
+ * TODO doc.
+ */
+size_t qpms_czero_roundoff_clean(complex double *arr, size_t nmemb, double atol);
 #endif // SYMMETRIES_H