QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

csphvec_reldiff version with tolerance 

Former-commit-id: f2ddb1fd9c55aff96b1e9f4145cb93bf255b6318
This commit is contained in:
Marek Nečada 2018-05-12 06:09:13 +00:00
parent 9bf0736dd4
commit 75e0d135d3
2 changed files with 17 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
qpms/tests/vswf_errs.c
View File

@ -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)
);

8
qpms/vectors.h
View File

@ -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