From 5a98b91b3dc8823c4ba3a91b45d57db1f70ae726 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 10 Jan 2020 11:43:35 +0200 Subject: [PATCH] Rewriting scatsystem: compiles without errors now. Former-commit-id: cd68b0feaef7181874d94dc535fd2cc9bc89e518 --- qpms/scatsystem.c | 60 +++++++++++++++++++++++------------------------ qpms/scatsystem.h | 12 +++++----- qpms/tmatrices.c | 16 ++++++++++--- qpms/tmatrices.h | 1 + 4 files changed, 49 insertions(+), 40 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 2d41036..8aa43f6 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -220,18 +220,18 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities; VLA! ss->tm_count = 0; for (qpms_ss_tmi_t i = 0; i < orig->tm_count; ++i) { - qpms_tmatrix_t *ti = qpms_tmatrix_apply_operation(orig->tm[i].op, tm_orig_omega[orig->tm[i].tmgi]); + qpms_tmatrix_t *ti = qpms_tmatrix_apply_operation(&(orig->tm[i].op), tm_orig_omega[orig->tm[i].tmgi]); qpms_ss_tmi_t j; for (j = 0; j < ss->tm_count; ++j) if (qpms_tmatrix_isclose(ssw->tm[i], ssw->tm[j], tol->rtol, tol->atol)) { break; } if (j == ss->tm_count) { // duplicity not found, save both the "abstract" and "at omega" T-matrices - qpms_tmatrix_operation_copy(&(ss->tm[j].op), orig->tm[j].op); + qpms_tmatrix_operation_copy(&ss->tm[j].op, &orig->tm[j].op); ss->tm[j].tmgi = orig->tm[j].tmgi; // T-matrix functions are preserved. ssw->tm[j] = ti; ss->max_bspecn = MAX(ssw->tm[j]->spec->n, ss->max_bspecn); - lMax = MAX(lMax, ss->tm[j].spec->lMax); + lMax = MAX(lMax, ssw->tm[j]->spec->lMax); ++(ss->tm_count); } else qpms_tmatrix_free(ti); @@ -264,12 +264,12 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){ const size_t d = ssw->tm[tmi]->spec->n; complex double *m; - QPMS_CRASHING_MALLOC(m, d*d*sizeof(complex double)); + QPMS_CRASHING_MALLOC(m, d*d*sizeof(complex double)); // ownership passed to ss->tm[ss->tm_count].op qpms_irot3_uvswfi_dense(m, ssw->tm[tmi]->spec, sym->rep3d[gmi]); - qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ss->tm[tmi], M[0]); + qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ssw->tm[tmi], m); qpms_ss_tmi_t tmj; for (tmj = 0; tmj < ss->tm_count; ++tmj) - if (qpms_tmatrix_isclose(transformed, ss->tm[tmj], tol->rtol, tol->atol)) + if (qpms_tmatrix_isclose(transformed, ssw->tm[tmj], tol->rtol, tol->atol)) break; if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists //TODO some "rounding error cleanup" symmetrisation could be performed here? @@ -277,13 +277,13 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, } else { // MISS, save the matrix (also the "abstract" one) ssw->tm[ss->tm_count] = transformed; qpms_tmatrix_operation_compose_chain_init(&(ss->tm[ss->tm_count].op), 2, 1); - struct qpms_tmatrix_operation_compose_chain * const o = &(ss->tm[ss->tm_count].op.compose_chain); - o->ops[0] = &(ss->tm[tmj].op); // Let's just borrow this + struct qpms_tmatrix_operation_compose_chain * const o = &(ss->tm[ss->tm_count].op.op.compose_chain); + o->ops[0] = & ss->tm[tmj].op; // Let's just borrow this o->ops_owned[0] = false; - o->ops[1] = o->opmem[0]; - o->ops[1]->typ = QPMS_TMATRIX_OPERATION_LRMATRIX; - o->ops[1]->op.lrmatrix.m = m; - o->ops[1]->op.lrmatrix.m_size = d * d; + o->opmem[0].typ = QPMS_TMATRIX_OPERATION_LRMATRIX; + o->opmem[0].op.lrmatrix.m = m; + o->opmem[0].op.lrmatrix.m_size = d * d; + o->ops[1] = o->opmem; o->ops_owned[1] = true; ++(ss->tm_count); } @@ -1025,11 +1025,11 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( { // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC) size_t fullvec_offsetR = 0; for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) { - const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecR = qpms_ss_bspec_pi(ss, piR); const cart3_t posR = ss->p[piR].pos; size_t fullvec_offsetC = 0; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { - const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecC = qpms_ss_bspec_pi(ss, piC); if(piC != piR) { // The diagonal will be dealt with later. const cart3_t posC = ss->p[piC].pos; QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, @@ -1053,7 +1053,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( complex double *qpms_scatsys_at_omega_build_modeproblem_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, - const qpms_scatsys_at_omega_t *ssw, + const qpms_scatsys_at_omega_t *ssw ) { const complex double k = ssw->wavenumber; @@ -1105,7 +1105,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. complex double *target_packed, const qpms_scatsys_at_omega_t *ssw, - qpms_iri_t iri, + qpms_iri_t iri ) { const qpms_scatsys_t *ss = ssw->ss; @@ -1426,7 +1426,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, Sblock, // Sblock is S(piR->piC) bspecR, bspecC->n, bspecC, 1, - a->k, posR, posC, QPMS_HANKEL_PLUS)); + k, posR, posC, QPMS_HANKEL_PLUS)); SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, @@ -1471,11 +1471,12 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread // this differs from the ...build_modeproblem_matrix... only by the `J` // maybe I should use this one there as well to save lines... TODO struct qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thread_arg{ - const qpms_scatsys_at_omega_t *ssw; + const qpms_scatsys_t *ss; qpms_ss_pi_t *opistartR_ptr; pthread_mutex_t *opistartR_mutex; qpms_iri_t iri; complex double *target_packed; + complex double k; qpms_bessel_t J; }; @@ -1483,9 +1484,7 @@ static void *qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thre { const struct qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thread_arg *a = arg; - const qpms_scatsys_at_omega_t *ssw = a->ssw; - const qpms_scatsys_t *ss = ssw->ss; - const complex double k = ssw->wavenumber; + const qpms_scatsys_t *ss = a->ss; const qpms_iri_t iri = a->iri; const size_t packedlen = ss->saecv_sizes[iri]; const qpms_bessel_t J = a->J; @@ -1524,7 +1523,7 @@ static void *qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thre if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n - const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecR = qpms_ss_bspec_pi(ss, orbitstartpiR); // This is the orbit-level matrix projecting the whole orbit onto the irrep. const complex double *omR = otR->irbases + otR->irbase_offsets[iri]; // Orbit coeff vector's full size: @@ -1548,7 +1547,7 @@ static void *qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thre const size_t packed_orbit_offsetC = ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC] + osnC * otC->irbase_sizes[iri]; - const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecC = qpms_ss_bspec_pi(ss, piC); // Orbit coeff vector's full size: const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n @@ -1607,7 +1606,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( complex double *target_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, - const complex double k; + const complex double k, qpms_bessel_t J ) { @@ -1623,7 +1622,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( pthread_mutex_t opistartR_mutex; QPMS_ENSURE_SUCCESS(pthread_mutex_init(&opistartR_mutex, NULL)); const struct qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thread_arg - arg = {ssw, &opistartR, &opistartR_mutex, iri, target_packed, J}; + arg = {ss, &opistartR, &opistartR_mutex, iri, target_packed, J}; // FIXME THIS IS NOT PORTABLE: long nthreads; @@ -1743,7 +1742,6 @@ complex double *qpms_scatsys_apply_Tmatrices_full( for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { complex double *ptarget = target_full + ss->fecv_pstarts[pi]; const complex double *psrc = inc_full + ss->fecv_pstarts[pi]; - const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi); // TODO check whether T-matrix is non-virtual after virtual t-matrices are implemented. const qpms_tmatrix_t *T = ssw->tm[ss->p[pi].tmatrix_id]; qpms_apply_tmatrix(ptarget, psrc, T); @@ -1789,7 +1787,7 @@ void qpms_ss_LU_free(qpms_ss_LU lu) { free(lu.ipiv); } -qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(complex double *mpmatrix_full, +qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmatrix_full, int *target_piv, const qpms_scatsys_at_omega_t *ssw) { const qpms_scatsys_t *ss = ssw->ss; QPMS_ENSURE(mpmatrix_full, "A non-NULL pointer to the pre-calculated mode matrix is required"); @@ -1805,7 +1803,7 @@ qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(complex double *mpmatr return lu; } -qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed, +qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed, int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) { QPMS_ENSURE(mpmatrix_packed, "A non-NULL pointer to the pre-calculated mode matrix is required"); size_t n = ssw->ss->saecv_sizes[iri]; @@ -1821,11 +1819,11 @@ qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(complex double return lu; } -qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU( +qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( complex double *target, int *target_piv, const qpms_scatsys_at_omega_t *ssw){ - target = qpms_scatsys_build_modeproblem_matrix_full(target, ssw); - return qpms_scatsys_modeproblem_matrix_full_factorise(target, target_piv, ssw); + target = qpms_scatsysw_build_modeproblem_matrix_full(target, ssw); + return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, ssw); } qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index cdf1883..7490b02 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -125,7 +125,7 @@ typedef struct qpms_ss_particle_orbitinfo { /// Auxillary type used in qpms_scatsys_t: A recepy to create another T-matrices by symmetry operations. typedef struct qpms_ss_derived_tmatrix { qpms_ss_tmgi_t tmgi; ///< Index of the corresponding qpms_scatsys_t::tm element. - struct qpms_tmatrix_operation_t *op; ///< Operation to derive this particular T-matrix. + struct qpms_tmatrix_operation_t op; ///< Operation to derive this particular T-matrix. } qpms_ss_derived_tmatrix_t; @@ -188,12 +188,12 @@ typedef struct qpms_scatsys_t { } qpms_scatsys_t; /// Retrieve the bspec of \a tmi'th element of \a ss->tm. -static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) { +static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) { return ss->tmg[ss->tm[tmi].tmgi].spec; } /// Retrieve the bspec of \a pi'th particle in \a ss->p. -static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(qpms_scatsys_t *ss, qpms_ss_pi_t pi) { +static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) { return ss->tmg[ss->tm[ss->p[pi].tmatrix_id].tmgi].spec; } @@ -330,7 +330,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( ); /// Creates the full \f$ (I - TS) \f$ matrix of the scattering system. -complex double *qpms_scatsys_at_omega_build_modeproblem_matrix_full( +complex double *qpms_scatsysw_build_modeproblem_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_at_omega_t *ssw @@ -370,7 +370,7 @@ typedef struct qpms_ss_LU { void qpms_ss_LU_free(qpms_ss_LU); /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. -qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU( +qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_at_omega_t *ssw @@ -385,7 +385,7 @@ qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( ); /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents. -qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise( +qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise( complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory. int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). const qpms_scatsys_at_omega_t *ssw diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 11ad54c..96516ba 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1038,8 +1038,17 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { &(f->op.compose_chain); if(o->opmem) { for(size_t i = 0; i < o->n; ++i) - if(o->ops_owned[i]) - qpms_tmatrix_operation_clear(o->ops[i]); + if(o->ops_owned[i]) { + // A bit complicated yet imperfect sanity/bound check, + // but this way we at least get rid of the const discard warning. + ptrdiff_t d = o->ops[i] - o->opmem; + QPMS_ENSURE(d >= 0 && d < o->opmem_size, + "Bound check failed; is there a bug in the" + " construction of compose_chain operation?\n" + "opmem_size = %zu, o->ops[%zu] - o->opmem = %td.", + o->opmem_size, i, d); + qpms_tmatrix_operation_clear(o->opmem + d); + } free(o->opmem); free(o->ops_owned); } @@ -1128,6 +1137,7 @@ void qpms_tmatrix_operation_compose_chain_init(qpms_tmatrix_operation_t *dest, o->ops_owned = NULL; o->opmem = NULL; } + o->opmem_size = opmem_size; } @@ -1175,7 +1185,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace( const qpms_tmatrix_operation_t *f, qpms_tmatrix_t *T) { switch(f->typ) { case QPMS_TMATRIX_OPERATION_NOOP: - return f; + return T; case QPMS_TMATRIX_OPERATION_LRMATRIX: return qpms_tmatrix_apply_symop_inplace(T, f->op.lrmatrix.m); case QPMS_TMATRIX_OPERATION_IROT3: diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index 82c537d..8b4d2e6 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -576,6 +576,7 @@ struct qpms_tmatrix_operation_compose_chain { size_t n; ///< Number of operations in ops; const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers owned by this.) struct qpms_tmatrix_operation_t *opmem; ///< (Optional) operations buffer into which elements of \a ops point. (Owned by this or NULL.) + size_t opmem_size; ///< Length of the opmem array. _Bool *ops_owned; ///< True for all sub operations owned by this and saved in opmem. NULL if opmem is NULL. };