Roundoff error cleaning for zeros.

Former-commit-id: 1efdd13b64022a93ae396bbc92345917de014ed8

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;
+    }
+  }
+}
 

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