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

Cleanup: eliminate bare abort()s 

Former-commit-id: bb1e4ada19e6bcbf87d6ea1fce0897c4478fb045
This commit is contained in:
Marek Nečada 2020-03-23 23:04:26 +02:00
parent 20a13cdb2c
commit 1b4b397093
9 changed files with 61 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
qpms/bessel.c
View File

@ -160,10 +160,8 @@ static inline qpms_errno_t qpms_sbessel_calculator_ensure_lMax(qpms_sbessel_calc
if (lMax <= c->lMax) if (lMax <= c->lMax)
return QPMS_SUCCESS; return QPMS_SUCCESS;
else { else {
if ( NULL == (c->akn = realloc(c->akn, sizeof(double) * akn_index(lMax + 2, 0)))) QPMS_CRASHING_REALLOC(c->akn, sizeof(double) * akn_index(lMax + 2, 0));
abort(); // QPMS_CRASHING_REALLOC(c->bkn, sizeof(complex double) * bkn_index(lMax + 1, 0));
//if ( NULL == (c->bkn = realloc(c->bkn, sizeof(complex double) * bkn_index(lMax + 1, 0))))
// abort();
for(qpms_l_t n = c->lMax+1; n <= lMax + 1; ++n) for(qpms_l_t n = c->lMax+1; n <= lMax + 1; ++n)
for(qpms_l_t k = 0; k <= n; ++k) for(qpms_l_t k = 0; k <= n; ++k)
c->akn[akn_index(n,k)] = exp(lgamma(n + k + 1) - k*M_LN2 - lgamma(k + 1) - lgamma(n - k + 1)); c->akn[akn_index(n,k)] = exp(lgamma(n + k + 1) - k*M_LN2 - lgamma(k + 1) - lgamma(n - k + 1));
@ -174,8 +172,7 @@ static inline qpms_errno_t qpms_sbessel_calculator_ensure_lMax(qpms_sbessel_calc
} }
complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x) { complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x) {
if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, n)) QPMS_ENSURE_SUCCESS(qpms_sbessel_calculator_ensure_lMax(c, n));
abort();
complex double z = I/x; complex double z = I/x;
complex double result = 0; complex double result = 0;
for (qpms_l_t k = n; k >= 0; --k) for (qpms_l_t k = n; k >= 0; --k)
@ -188,8 +185,7 @@ complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, co
qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t * const c, qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t * const c,
const qpms_l_t lMax, const complex double x, complex double * const target) { const qpms_l_t lMax, const complex double x, complex double * const target) {
if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, lMax)) QPMS_ENSURE_SUCCESS(qpms_sbessel_calculator_ensure_lMax(c, lMax));
abort();
memset(target, 0, sizeof(complex double) * lMax); memset(target, 0, sizeof(complex double) * lMax);
complex double kahancomp[lMax]; complex double kahancomp[lMax];
memset(kahancomp, 0, sizeof(complex double) * lMax); memset(kahancomp, 0, sizeof(complex double) * lMax);

29
qpms/ewald.c
View File

@ -62,6 +62,7 @@ typedef enum {
} ewald3_constants_option; } ewald3_constants_option;
// TODO perhaps rewrite everything as agnostic.
static const ewald3_constants_option type = EWALD32_CONSTANTS_AGNOSTIC; static const ewald3_constants_option type = EWALD32_CONSTANTS_AGNOSTIC;
qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, const ewald3_constants_option type */, qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, const ewald3_constants_option type */,
@ -109,7 +110,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
/ (factorial(j) * factorial((n-m)/2-j) * factorial((n+m)/2-j)); / (factorial(j) * factorial((n-m)/2-j) * factorial((n+m)/2-j));
break; break;
default: default:
abort(); QPMS_INVALID_ENUM(type);
} }
} }
s1_constfacs_sz_cumsum += 1 + c->s1_jMaxes[y]; s1_constfacs_sz_cumsum += 1 + c->s1_jMaxes[y];
@ -139,7 +140,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
/ (factorial(j) * pow(2, 2*j) * factorial(n - 2*j)); / (factorial(j) * pow(2, 2*j) * factorial(n - 2*j));
break; break;
default: default:
abort(); // wrong type argument or not implemented QPMS_INVALID_ENUM(type); // wrong type argument or not implemented
} }
} }
s1_constfacs_1Dz_sz_cumsum += 1 + n / 2; s1_constfacs_1Dz_sz_cumsum += 1 + n / 2;
@ -154,12 +155,9 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
// N.B. here I use the GSL_SF_LEGENRE_NONE, in order to be consistent with translations.c // N.B. here I use the GSL_SF_LEGENRE_NONE, in order to be consistent with translations.c
c->legendre_normconv = GSL_SF_LEGENDRE_NONE; c->legendre_normconv = GSL_SF_LEGENDRE_NONE;
// Moreover, using this approach (i.e. gsl) takes about 64kB extra memory // Moreover, using this approach (i.e. gsl) takes about 64kB extra memory
if(GSL_SUCCESS != gsl_sf_legendre_array_e(c->legendre_normconv, lMax, 0, csphase, c->legendre0)) QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(c->legendre_normconv, lMax, 0, csphase, c->legendre0));
abort(); QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(c->legendre_normconv, lMax, +1, csphase, c->legendre_plus1));
if(GSL_SUCCESS != gsl_sf_legendre_array_e(c->legendre_normconv, lMax, +1, csphase, c->legendre_plus1)) QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(c->legendre_normconv, lMax, -1, csphase, c->legendre_minus1));
abort();
if(GSL_SUCCESS != gsl_sf_legendre_array_e(c->legendre_normconv, lMax, -1, csphase, c->legendre_minus1))
abort();
return c; return c;
} }
@ -182,11 +180,9 @@ int ewald3_sigma0(complex double *result, double *err,
{ {
qpms_csf_result gam; qpms_csf_result gam;
complex double z = -csq(k/(2*eta)); complex double z = -csq(k/(2*eta));
int retval = complex_gamma_inc_e(-0.5, z, QPMS_ENSURE_SUCCESS(complex_gamma_inc_e(-0.5, z,
// we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle] // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle]
QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, &gam); QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, &gam));
if (0 != retval)
abort();
*result = gam.val * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI; *result = gam.val * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI;
if(err) if(err)
*err = gam.err * fabs(c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI); *err = gam.err * fabs(c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI);
@ -305,7 +301,7 @@ int ewald3_21_xy_sigma_long (
int retval = complex_gamma_inc_e(0.5-j, z, int retval = complex_gamma_inc_e(0.5-j, z,
// we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle] // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle]
QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j); QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j);
if(!(retval==0 || retval==GSL_EUNDRFLW)) abort(); QPMS_ENSURE_SUCCESS_OR(retval, GSL_EUNDRFLW);
} }
if (latdim & LAT1D) if (latdim & LAT1D)
factor1d = M_SQRT1_2 * .5 * k * gamma_pq; factor1d = M_SQRT1_2 * .5 * k * gamma_pq;
@ -446,7 +442,7 @@ int ewald3_1_z_sigma_long (
int retval = complex_gamma_inc_e(0.5-j, z, int retval = complex_gamma_inc_e(0.5-j, z,
// we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle] // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle]
QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j); QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j);
if(!(retval==0 || retval==GSL_EUNDRFLW)) abort(); QPMS_ENSURE_SUCCESS_OR(retval, GSL_EUNDRFLW);
} }
// R-DEPENDENT END // R-DEPENDENT END
// TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
@ -713,8 +709,7 @@ int ewald3_sigma_short(
legendre_array = c->legendre0; legendre_array = c->legendre0;
break; break;
default: default:
if(GSL_SUCCESS != gsl_sf_legendre_array_e(c->legendre_normconv, lMax, cos(Rpq_shifted_theta), c->legendre_csphase, legendre_buf)) QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(c->legendre_normconv, lMax, cos(Rpq_shifted_theta), c->legendre_csphase, legendre_buf));
abort();
legendre_array = legendre_buf; legendre_array = legendre_buf;
break; break;
} }
@ -732,7 +727,7 @@ int ewald3_sigma_short(
} else } else
retval = ewald32_sr_integral_ck(r_pq_shifted, k, n, eta, retval = ewald32_sr_integral_ck(r_pq_shifted, k, n, eta,
cintres+n, interr + n, workspace); cintres+n, interr + n, workspace);
if (retval) abort(); QPMS_ENSURE_SUCCESS(retval);
} // otherwise recycle the integrals } // otherwise recycle the integrals
for (qpms_m_t m = -n; m <= n; ++m){ for (qpms_m_t m = -n; m <= n; ++m){
complex double e_imf; complex double e_imf;

10
qpms/latticegens.c
View File

@ -320,7 +320,7 @@ PGen PGen_1D_new_minMaxR(double period, double offset, double minR, bool inc_min
s->ptindex = ptindex_dec(s->ptindex); s->ptindex = ptindex_dec(s->ptindex);
break; break;
default: default:
abort(); // invalid argument / not implemented QPMS_WTF;
} }
s->a = period; s->a = period;
@ -357,7 +357,7 @@ PGenZReturnData PGen_1D_next_z(PGen *g) {
} else theEnd = true; } else theEnd = true;
break; break;
default: default:
abort(); // invalid value QPMS_INVALID_ENUM(s->incdir);
} }
if (!theEnd) { if (!theEnd) {
const PGenZReturnData retval = {PGEN_NOTDONE | PGEN_AT_Z, zval}; const PGenZReturnData retval = {PGEN_NOTDONE | PGEN_AT_Z, zval};
@ -391,7 +391,7 @@ PGenSphReturnData PGen_1D_next_sph(PGen *g) {
} else theEnd = true; } else theEnd = true;
break; break;
default: default:
abort(); // invalid value QPMS_INVALID_ENUM(s->incdir);
} }
if (!theEnd) { if (!theEnd) {
const PGenSphReturnData retval = {PGEN_NOTDONE | PGEN_AT_Z | PGEN_COORDS_SPH, const PGenSphReturnData retval = {PGEN_NOTDONE | PGEN_AT_Z | PGEN_COORDS_SPH,
@ -512,7 +512,7 @@ PGenCart2ReturnData PGen_xyWeb_next_cart2(PGen *g) {
if(s->j >= 0) ++s->phase; if(s->j >= 0) ++s->phase;
break; break;
default: default:
abort(); QPMS_WTF;
} }
if(s->phase == 4) { // phase overflow, start new layer if(s->phase == 4) { // phase overflow, start new layer
++s->layer; ++s->layer;
@ -549,7 +549,7 @@ PGenCart2ReturnData PGen_xyWeb_next_cart2(PGen *g) {
if(s->j >= 0) ++s->phase; if(s->j >= 0) ++s->phase;
break; break;
default: default:
abort(); QPMS_WTF;
} }
if(s->phase == 6) { // phase overflow, start next layer if(s->phase == 6) { // phase overflow, start next layer
++s->layer; ++s->layer;

49
qpms/lattices2d.c
View File

@ -362,14 +362,9 @@ static int trilatgen_pointlist_extend_capacity(triangular_lattice_gen_t *g, size
trilatgen_pointlist_linearise(g); trilatgen_pointlist_linearise(g);
intcoord2_t *newmem = realloc(p->pointlist_base, newcapacity * sizeof(intcoord2_t)); QPMS_CRASHING_REALLOC(p->pointlist_base, newcapacity * sizeof(intcoord2_t));
if (newmem != NULL) { p->pointlist_capacity = newcapacity;
p->pointlist_base = newmem; return 0;
p->pointlist_capacity = newcapacity;
return 0;
} else
abort();
} }
// lower estimate for the number of lattice points inside the circumscribed hexagon, but outside the circle // lower estimate for the number of lattice points inside the circumscribed hexagon, but outside the circle
@ -404,16 +399,8 @@ static int trilatgen_ensure_ps_rs_capacity(triangular_lattice_gen_t *g, int maxs
if (needed_capacity <= g->priv->ps_rs_capacity) if (needed_capacity <= g->priv->ps_rs_capacity)
return 0; // probably does not happen, but fuck it return 0; // probably does not happen, but fuck it
double *newmem = realloc(g->ps.rs, needed_capacity * sizeof(double)); QPMS_CRASHING_REALLOC(g->ps.rs, needed_capacity * sizeof(double));
if (newmem != NULL) QPMS_CRASHING_REALLOC(g->ps.r_offsets, (needed_capacity + 1) * sizeof(ptrdiff_t));
g->ps.rs = newmem;
else
abort();
ptrdiff_t *newmem2 = realloc(g->ps.r_offsets, (needed_capacity + 1) * sizeof(ptrdiff_t));
if (newmem2 != NULL)
g->ps.r_offsets = newmem2;
else
abort();
g->priv->ps_rs_capacity = needed_capacity; g->priv->ps_rs_capacity = needed_capacity;
return 0; return 0;
@ -426,12 +413,7 @@ static int trilatgen_ensure_ps_points_capacity(triangular_lattice_gen_t *g, int
if(needed_capacity <= g->priv->ps_points_capacity) if(needed_capacity <= g->priv->ps_points_capacity)
return 0; return 0;
point2d *newmem = realloc(g->ps.base, needed_capacity * sizeof(point2d)); QPMS_CRASHING_REALLOC(g->ps.base, needed_capacity * sizeof(point2d));
if (newmem != NULL)
g->ps.base = newmem;
else
abort();
g->priv->ps_points_capacity = needed_capacity; g->priv->ps_points_capacity = needed_capacity;
return 0; return 0;
} }
@ -468,7 +450,7 @@ static void trilatgen_sort_pointlist(triangular_lattice_gen_t *g) {
compar = trilat_cmp_intcoord2_by_3r2_plus1s; compar = trilat_cmp_intcoord2_by_3r2_plus1s;
break; break;
default: default:
abort(); QPMS_WTF;
} }
qsort(p->pointlist_base + p->pointlist_beg, p->pointlist_n, sizeof(intcoord2_t), compar); qsort(p->pointlist_base + p->pointlist_beg, p->pointlist_n, sizeof(intcoord2_t), compar);
} }
@ -580,7 +562,7 @@ int triangular_lattice_gen_extend_to_steps(triangular_lattice_gen_t * g, int max
thepoint = point2d_fromxy(-(M_SQRT3_2*coord.j + g->hexshift*M_1_SQRT3)*g->a, (coord.i+.5*coord.j)*g->a); thepoint = point2d_fromxy(-(M_SQRT3_2*coord.j + g->hexshift*M_1_SQRT3)*g->a, (coord.i+.5*coord.j)*g->a);
break; break;
default: default:
abort(); QPMS_INVALID_ENUM(g->orientation);
} }
g->ps.base[g->ps.r_offsets[g->ps.nrs+1]] = thepoint; g->ps.base[g->ps.r_offsets[g->ps.nrs+1]] = thepoint;
++(g->ps.r_offsets[g->ps.nrs+1]); ++(g->ps.r_offsets[g->ps.nrs+1]);
@ -624,18 +606,9 @@ int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, const int
return 0; return 0;
triangular_lattice_gen_extend_to_steps(g->tg, maxsteps); triangular_lattice_gen_extend_to_steps(g->tg, maxsteps);
double *newmem = realloc(g->ps.rs, g->tg->ps.nrs * sizeof(double)); QPMS_CRASHING_REALLOC(g->ps.rs, g->tg->ps.nrs * sizeof(double));
if (NULL != newmem) QPMS_CRASHING_REALLOC(g->ps.r_offsets, (g->tg->ps.nrs+1) * sizeof(ptrdiff_t));
g->ps.rs = newmem; QPMS_CRASHING_REALLOC(g->ps.base, 2 * (g->tg->ps.r_offsets[g->tg->ps.nrs]) * sizeof(point2d));
else abort();
ptrdiff_t *newmem2 = realloc(g->ps.r_offsets, (g->tg->ps.nrs+1) * sizeof(ptrdiff_t));
if (NULL != newmem2)
g->ps.r_offsets = newmem2;
else abort();
point2d *newmem3 = realloc(g->ps.base, 2 * (g->tg->ps.r_offsets[g->tg->ps.nrs]) * sizeof(point2d));
if (NULL != newmem3)
g->ps.base = newmem3;
else abort();
// Now copy (new) contents of g->tg->ps into g->ps, but with inverse copy of each point // Now copy (new) contents of g->tg->ps into g->ps, but with inverse copy of each point
for (size_t ri = g->ps.nrs; ri <= g->tg->ps.nrs; ++ri) for (size_t ri = g->ps.nrs; ri <= g->tg->ps.nrs; ++ri)

2
qpms/scatsystem.c
View File

@ -749,7 +749,7 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
// Get the projector (also workspace for right sg. vect.) // Get the projector (also workspace for right sg. vect.)
complex double *projector = qpms_orbit_irrep_projector_matrix(NULL, complex double *projector = qpms_orbit_irrep_projector_matrix(NULL,
ot, bspec, sym, iri); ot, bspec, sym, iri);
if(!projector) abort(); QPMS_ENSURE(projector != NULL, "got NULL from qpms_orbit_irrep_projector_matrix()");
// Workspace for the right singular vectors. // Workspace for the right singular vectors.
complex double *V_H; QPMS_CRASHING_MALLOC(V_H, n*n*N*N*sizeof(complex double)); complex double *V_H; QPMS_CRASHING_MALLOC(V_H, n*n*N*N*sizeof(complex double));
// THIS SHOULD NOT BE NECESSARY // THIS SHOULD NOT BE NECESSARY

16
qpms/symmetries.c
View File

@ -54,7 +54,7 @@ complex double *qpms_zflip_uvswi_dense(
qpms_vswf_type_t ct; qpms_vswf_type_t ct;
qpms_l_t cl; qpms_l_t cl;
qpms_m_t cm; qpms_m_t cm;
if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort(); QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl));
if (rl == cl && rm == cm && rt == ct) if (rl == cl && rm == cm && rt == ct)
switch(rt) { switch(rt) {
case QPMS_VSWF_ELECTRIC: case QPMS_VSWF_ELECTRIC:
@ -65,7 +65,7 @@ complex double *qpms_zflip_uvswi_dense(
target[n*row + col] = -min1pow(cm + cl); target[n*row + col] = -min1pow(cm + cl);
break; break;
default: default:
abort(); QPMS_INVALID_ENUM(rt);
} }
else target[n*row + col] = 0; else target[n*row + col] = 0;
} }
@ -91,7 +91,7 @@ complex double *qpms_yflip_uvswi_dense(
qpms_vswf_type_t ct; qpms_vswf_type_t ct;
qpms_l_t cl; qpms_l_t cl;
qpms_m_t cm; qpms_m_t cm;
if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort(); QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl));
if (rl == cl && rm == -cm && rt == ct) if (rl == cl && rm == -cm && rt == ct)
switch(rt) { switch(rt) {
case QPMS_VSWF_ELECTRIC: case QPMS_VSWF_ELECTRIC:
@ -102,7 +102,7 @@ complex double *qpms_yflip_uvswi_dense(
target[n*row + col] = -min1pow(rm); target[n*row + col] = -min1pow(rm);
break; break;
default: default:
abort(); QPMS_INVALID_ENUM(rt);
} }
else target[n*row + col] = 0; else target[n*row + col] = 0;
} }
@ -128,7 +128,7 @@ complex double *qpms_xflip_uvswi_dense(
qpms_vswf_type_t ct; qpms_vswf_type_t ct;
qpms_l_t cl; qpms_l_t cl;
qpms_m_t cm; qpms_m_t cm;
if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort(); QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl));
if (rl == cl && rm == -cm && rt == ct) if (rl == cl && rm == -cm && rt == ct)
switch(rt) { switch(rt) {
case QPMS_VSWF_ELECTRIC: case QPMS_VSWF_ELECTRIC:
@ -139,7 +139,7 @@ complex double *qpms_xflip_uvswi_dense(
target[n*row + col] = -1; target[n*row + col] = -1;
break; break;
default: default:
abort(); QPMS_INVALID_ENUM(rt);
} }
else target[n*row + col] = 0; else target[n*row + col] = 0;
} }
@ -169,7 +169,7 @@ complex double *qpms_zrot_uvswi_dense(
qpms_vswf_type_t ct; qpms_vswf_type_t ct;
qpms_l_t cl; qpms_l_t cl;
qpms_m_t cm; qpms_m_t cm;
if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort(); QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl));
if (rl == cl && rm == cm && rt == ct) // TODO COMPARE WITH PYTHON if (rl == cl && rm == cm && rt == ct) // TODO COMPARE WITH PYTHON
target[n*row + col] = cexp(/* - ?*/I * rm * phi); target[n*row + col] = cexp(/* - ?*/I * rm * phi);
else target[n*row + col] = 0; else target[n*row + col] = 0;
@ -212,7 +212,7 @@ complex double *qpms_irot3_uvswfi_dense(
qpms_vswf_type_t ct; qpms_vswf_type_t ct;
qpms_l_t cl; qpms_l_t cl;
qpms_m_t cm; qpms_m_t cm;
if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort(); QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl));
if (rl == cl && rt == ct) if (rl == cl && rt == ct)
// TODO qpms_vswf_irot_elem_from_irot3 might be slow and not too accurate for large l // TODO qpms_vswf_irot_elem_from_irot3 might be slow and not too accurate for large l
target[n*row + col] = // Checkme rm and cm order target[n*row + col] = // Checkme rm and cm order

42
qpms/tmatrices.c
View File

@ -41,15 +41,12 @@
#define MIN(x,y) ((x) > (y) ? (y) : (x)) #define MIN(x,y) ((x) > (y) ? (y) : (x))
qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) { qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t)); qpms_tmatrix_t *t;
if (!t) abort(); QPMS_CRASHING_MALLOC(t, sizeof(qpms_tmatrix_t));
else { t->spec = bspec;
t->spec = bspec; size_t n = bspec->n;
size_t n = bspec->n; QPMS_CRASHING_CALLOC(t->m, n*n, sizeof(complex double));
t->m = calloc(n*n, sizeof(complex double)); t->owns_m = true;
if (!t->m) abort();
t->owns_m = true;
}
return t; return t;
} }
@ -296,11 +293,11 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incou
// check if all matrices have the same bspec // check if all matrices have the same bspec
for (size_t i = 0; i < incount; ++i) for (size_t i = 0; i < incount; ++i)
if (!qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec)) QPMS_ENSURE(qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec),
abort(); "All T-matrices for interpolation must have the same basis!");
if (!(ip->splines_real = calloc(n*n,sizeof(gsl_spline *)))) abort(); QPMS_CRASHING_CALLOC(ip->splines_real, n*n, sizeof(gsl_spline *));
if (!(ip->splines_imag = calloc(n*n,sizeof(gsl_spline *)))) abort(); QPMS_CRASHING_CALLOC(ip->splines_imag, n*n,sizeof(gsl_spline *));
for (size_t row = 0; row < n; ++row) for (size_t row = 0; row < n; ++row)
for (size_t col = 0; col < n; ++col) { for (size_t col = 0; col < n; ++col) {
double y_real[incount], y_imag[incount]; double y_real[incount], y_imag[incount];
@ -313,13 +310,13 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incou
if (n0_real) { if (n0_real) {
gsl_spline *s = gsl_spline *s =
ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount); ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount);
if (gsl_spline_init(s, freqs, y_real, incount) != 0 /*GSL_SUCCESS*/) abort(); QPMS_ENSURE_SUCCESS(gsl_spline_init(s, freqs, y_real, incount));
} }
else ip->splines_real[n * row + col] = NULL; else ip->splines_real[n * row + col] = NULL;
if (n0_imag) { if (n0_imag) {
gsl_spline *s = gsl_spline *s =
ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount); ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount);
if (gsl_spline_init(s, freqs, y_imag, incount) != 0 /*GSL_SUCCESS*/) abort(); QPMS_ENSURE_SUCCESS(gsl_spline_init(s, freqs, y_imag, incount));
} }
else ip->splines_imag[n * row + col] = NULL; else ip->splines_imag[n * row + col] = NULL;
} }
@ -403,21 +400,18 @@ qpms_errno_t qpms_read_scuff_tmatrix(
if (linebuf[0] == '#') continue; if (linebuf[0] == '#') continue;
int Alpha, LAlpha, MAlpha, PAlpha, Beta, LBeta, MBeta, PBeta; int Alpha, LAlpha, MAlpha, PAlpha, Beta, LBeta, MBeta, PBeta;
double currentfreq_su, tr, ti; double currentfreq_su, tr, ti;
if (11 != sscanf(linebuf, "%lf %d %d %d %d %d %d %d %d %lf %lf", QPMS_ENSURE(11 == sscanf(linebuf, "%lf %d %d %d %d %d %d %d %d %lf %lf",
&currentfreq_su, &Alpha, &LAlpha, &MAlpha, &PAlpha, &currentfreq_su, &Alpha, &LAlpha, &MAlpha, &PAlpha,
&Beta, &LBeta, &MBeta, &PBeta, &tr, &ti)) &Beta, &LBeta, &MBeta, &PBeta, &tr, &ti),
abort(); // Malformed T-matrix file "Malformed T-matrix file");
if (currentfreq_su != lastfreq_su) { // New frequency -> new T-matrix if (currentfreq_su != lastfreq_su) { // New frequency -> new T-matrix
++*n; ++*n;
lastfreq_su = currentfreq_su; lastfreq_su = currentfreq_su;
if(*n > n_alloc) { if(*n > n_alloc) {
n_alloc *= 2; n_alloc *= 2;
*freqs = realloc(*freqs, n_alloc * sizeof(double)); QPMS_CRASHING_REALLOC(*freqs, n_alloc * sizeof(double));
if (freqs_su) *freqs_su = realloc(*freqs_su, n_alloc * sizeof(double)); if (freqs_su) {QPMS_CRASHING_REALLOC(*freqs_su, n_alloc * sizeof(double));}
*tmdata = realloc(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc); QPMS_CRASHING_REALLOC(*tmdata, sizeof(complex double) * bs->n * bs->n * n_alloc);
if (!*freqs || (!freqs_su != !*freqs_su) || !*tmdata)
qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
"realloc() failed.");
} }
if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su; if (freqs_su) (*freqs_su)[*n-1] = currentfreq_su;
(*freqs)[*n-1] = qpms_SU2SI(currentfreq_su); (*freqs)[*n-1] = qpms_SU2SI(currentfreq_su);

3
qpms/translations_python.c
View File

@ -10,7 +10,6 @@
#include "tiny_inlines.h" #include "tiny_inlines.h"
#include "assert_cython_workaround.h" #include "assert_cython_workaround.h"
#include "kahansum.h" #include "kahansum.h"
#include <stdlib.h> //abort()
#include <gsl/gsl_sf_coupling.h> #include <gsl/gsl_sf_coupling.h>
#include "qpms_error.h" #include "qpms_error.h"
#include "normalisation.h" #include "normalisation.h"
@ -110,7 +109,7 @@ int qpms_cython_trans_calculator_get_AB_arrays_loop(
// these casts are not needed // these casts are not needed
*((double *) r_p), *((double *) theta_p), *((double *)phi_p), *((double *) r_p), *((double *) theta_p), *((double *)phi_p),
(int)(*((npy_bool *) r_ge_d_p)), J); (int)(*((npy_bool *) r_ge_d_p)), J);
if (errval_local) abort(); QPMS_ENSURE_SUCCESS(errval_local);
// increment the last index 'digit' (ax is now resnd-1; we don't have do-while loop in python) // increment the last index 'digit' (ax is now resnd-1; we don't have do-while loop in python)
++local_indices[ax]; ++local_indices[ax];

2
qpms/wigner.c
View File

@ -9,7 +9,7 @@
complex double qpms_wignerD_elem(const qpms_quat_t R, complex double qpms_wignerD_elem(const qpms_quat_t R,
const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m) { const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m) {
// TODO do some optimisations... The combinatoric coeffs could be precomputed. // TODO do some optimisations... The combinatoric coeffs could be precomputed.
if (abs(m) > l || abs(mp) > l) abort();//return 0; // It's safer to crash. :D QPMS_ENSURE(abs(m) <= l || abs(mp) <= l, "Got invalid values l = %d, m = %d", (int)l, (int)m);
const double mRa = cabs(R.a), mRb = cabs(R.b); const double mRa = cabs(R.a), mRb = cabs(R.b);
if (mRa < WIGNER_ATOL) { if (mRa < WIGNER_ATOL) {
if (m != -mp) return 0; if (m != -mp) return 0;