From 75e0d135d3972c2f24110878300cf19269165a12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sat, 12 May 2018 06:09:13 +0000 Subject: [PATCH] csphvec_reldiff version with tolerance Former-commit-id: f2ddb1fd9c55aff96b1e9f4145cb93bf255b6318 --- qpms/tests/vswf_errs.c | 20 +++++++++++--------- qpms/vectors.h | 8 ++++++-- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/qpms/tests/vswf_errs.c b/qpms/tests/vswf_errs.c index b390ab0..0058c94 100644 --- a/qpms/tests/vswf_errs.c +++ b/qpms/tests/vswf_errs.c @@ -29,6 +29,8 @@ char *normstr(qpms_normalisation_t norm) { } } +#define DIFFTOL (1e-13) + int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavetype, cart3_t o2minuso1, int npoints, cart3_t *o1points); //int test_planewave_decomposition(cart3_t k, ccart3_t E, qpms_l_t lMax, qpms_normalisation_t norm, int npoints, cart3_t *points); @@ -72,7 +74,7 @@ int main() { ++i){ qpms_normalisation_t norm = i | (use_csbit ? QPMS_NORMALISATION_T_CSBIT : 0); qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, norm); - for(int J = 2; J <= 2; ++J) + for(int J = 3; J <= 3; ++J) test_sphwave_translation(c, J, o2minuso1, npoints, points); qpms_trans_calculator_free(c); } @@ -140,10 +142,10 @@ int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavet csphvec_t M2s2_sgAB_regw = qpms_eval_vswf(w2s, NULL, A_sg, B_sg, lMax,QPMS_BESSEL_REGULAR, c->normalisation); csphvec_t M2s2_sgAB_sgw = qpms_eval_vswf(w2s, NULL, A_sg, B_sg, lMax,wavetype, c->normalisation); printf("Merr:\tRC_RW %.2e\tRC_SW %.2e\tSC_RW %.2e\tSC_SW %.2e\n", - csphvec_reldiff(M1s2, M2s2_regAB_regw), - csphvec_reldiff(M1s2, M2s2_regAB_sgw), - csphvec_reldiff(M1s2, M2s2_sgAB_regw), - csphvec_reldiff(M1s2, M2s2_sgAB_sgw) + csphvec_reldiff_abstol(M1s2, M2s2_regAB_regw, DIFFTOL), + csphvec_reldiff_abstol(M1s2, M2s2_regAB_sgw, DIFFTOL), + csphvec_reldiff_abstol(M1s2, M2s2_sgAB_regw, DIFFTOL), + csphvec_reldiff_abstol(M1s2, M2s2_sgAB_sgw, DIFFTOL) ); @@ -168,10 +170,10 @@ int test_sphwave_translation(const qpms_trans_calculator *c, qpms_bessel_t wavet csphvec_t N2s2_sgAB_regw = qpms_eval_vswf(w2s, NULL, B_sg, A_sg, lMax,QPMS_BESSEL_REGULAR, c->normalisation); csphvec_t N2s2_sgAB_sgw = qpms_eval_vswf(w2s, NULL, B_sg, A_sg, lMax,wavetype, c->normalisation); printf("Nerr:\tRC_RW %.2e\tRC_SW %.2e\tSC_RW %.2e\tSC_SW %.2e\n", - csphvec_reldiff(N1s2, N2s2_regAB_regw), - csphvec_reldiff(N1s2, N2s2_regAB_sgw), - csphvec_reldiff(N1s2, N2s2_sgAB_regw), - csphvec_reldiff(N1s2, N2s2_sgAB_sgw) + csphvec_reldiff_abstol(N1s2, N2s2_regAB_regw, DIFFTOL), + csphvec_reldiff_abstol(N1s2, N2s2_regAB_sgw, DIFFTOL), + csphvec_reldiff_abstol(N1s2, N2s2_sgAB_regw, DIFFTOL), + csphvec_reldiff_abstol(N1s2, N2s2_sgAB_sgw, DIFFTOL) ); diff --git a/qpms/vectors.h b/qpms/vectors.h index 137a703..dda691a 100644 --- a/qpms/vectors.h +++ b/qpms/vectors.h @@ -159,11 +159,15 @@ static inline double csphvec_norm(const csphvec_t a) { return sqrt(creal(a.rc * conj(a.rc) + a.thetac * conj(a.thetac) + a.phic * conj(a.phic))); } -static inline double csphvec_reldiff(const csphvec_t a, const csphvec_t b) { +static inline double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance) { double anorm = csphvec_norm(a); double bnorm = csphvec_norm(b); - if (anorm == 0 && bnorm == 0) return 0; + if (anorm <= tolerance && bnorm <= tolerance) return 0; return csphvec_norm(csphvec_substract(a,b)) / (anorm + bnorm); } +static inline double csphvec_reldiff(const csphvec_t a, const csphvec_t b) { + return csphvec_reldiff_abstol(a, b, 0); +} + #endif //VECTORS_H