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