From 5dd93235f0c1f1bf3d87435ebf373894ad49f54d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 9 Jan 2020 10:42:49 +0200 Subject: [PATCH] New qpms_scatsys_apply_symmetry kinda done? Former-commit-id: 49a7a7984af6ad6e0e5ec1b5cc7b61ac06f81b8d --- qpms/scatsystem.c | 59 +++++++++++++++++++++++++++++++---------------- qpms/tmatrices.c | 33 ++++++++++++++++++++++++-- qpms/tmatrices.h | 10 ++++++++ 3 files changed, 80 insertions(+), 22 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 660eac7..ac9f1b8 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -176,7 +176,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, // Copy the qpms_tmatrix_fuction_t from orig ss->tmg_count = orig->tmg_count; - QPMS_CRASHING_MALLOC(ss->tmg, ss->tmg_count * sizeof(*(ss->tmg))); //TODO free + QPMS_CRASHING_MALLOC(ss->tmg, ss->tmg_count * sizeof(*(ss->tmg))); memcpy(ss->tmg, orig->tmg, ss->tmg_count * sizeof(*(ss->tmg))); // Allocate T-matrix, particle and particle orbit info arrays @@ -197,8 +197,8 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, } // Evaluate the original T-matrices at omega - qpms_tmatrix_t **tm_orig_omega; - QPMS_CRASHING_MALLOC(tm_orig_omega, orig->tmg_count * sizeof(*tm_orig_omega)); //TODO free + qpms_tmatrix_t **tm_orig_omega; + QPMS_CRASHING_MALLOC(tm_orig_omega, orig->tmg_count * sizeof(*tm_orig_omega)); for(qpms_ss_tmgi_t i = 0; i < orig->tmg_count; ++i) tm_orig_omega[i] = qpms_tmatrix_init_from_function(orig->tmg[i], omega); @@ -211,36 +211,43 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, ssw->wavenumber = qpms_wavenumber(omega, ssw->medium); // we will be using ss->tm_capacity also for ssw->tm QPMS_CRASHING_MALLOC(ssw->tm, ss->tm_capacity * sizeof(*(ssw->tm))); // returned - for (qpms_ss_tmi_t i = 0; i < ss->tm_count; ++i) - ssw->tm[i] = qpms_tmatrix_apply_operation(ss->tm[i].op, tm_orig_omega[ss->tm[i].tmgi]); - // Copy T-matrices; checking for duplicities + // Evaluate T-matrices at omega; checking for duplicities ss->max_bspecn = 0; // We'll need it later.for memory alloc estimates. 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_ss_tmi_t j; for (j = 0; j < ss->tm_count; ++j) - if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], tol->rtol, tol->atol)) { + if (qpms_tmatrix_isclose(ssw->tm[i], ssw->tm[j], tol->rtol, tol->atol)) { break; } - if (j == ss->tm_count) { // duplicity not found, copy the t-matrix - ss->tm[j] = qpms_tmatrix_copy(orig->tm[i]); - ss->max_bspecn = MAX(ss->tm[j]->spec->n, ss->max_bspecn); + if (j == ss->tm_count) { // duplicity not found, save both the "abstract" and "at omega" T-matrices + ss->tm[j].op = qpms_tmatrix_operation_copy(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); ++(ss->tm_count); } + else qpms_tmatrix_free(ti); tm_dupl_remap[i] = j; if (normalisation == QPMS_NORMALISATION_UNDEF) - normalisation = ss->tm[i]->spec->norm; + normalisation = ssw->tm[i]->spec->norm; // We expect all bspec norms to be the same. - else QPMS_ENSURE(normalisation == ss->tm[j]->spec->norm, + else QPMS_ENSURE(normalisation == ssw->tm[j]->spec->norm, "Normalisation convention must be the same for all T-matrices." - " %d != %d\n", normalisation, ss->tm[j]->spec->norm); + " %d != %d\n", normalisation, ssw->tm[j]->spec->norm); } + // Free the original T-matrices at omega + for(qpms_ss_tmgi_t i = 0; i < orig->tmg_count; ++i) + qpms_tmatrix_free(tm_orig_omega[i]); + free(tm_orig_omega); + // Copy particles, remapping the t-matrix indices for (qpms_ss_pi_t i = 0; i < orig->p_count; ++(i)) { ss->p[i] = orig->p[i]; @@ -254,18 +261,29 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, // Extend the T-matrices list by the symmetry operations for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){ - const size_t d = ss->tm[tmi]->spec->n; - complex double M[d][d]; // transformation matrix - qpms_irot3_uvswfi_dense(M[0], ss->tm[tmi]->spec, sym->rep3d[gmi]); + const size_t d = ssw->tm[tmi]->spec->n; + complex double *m; + QPMS_CRASHING_MALLOC(m, d*d*sizeof(complex double)); + 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_ss_tmi_t tmj; for (tmj = 0; tmj < ss->tm_count; ++tmj) if (qpms_tmatrix_isclose(transformed, ss->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? qpms_tmatrix_free(transformed); - } else { // MISS, save the matrix and increment the count - ss->tm[ss->tm_count] = transformed; + } 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 + 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->ops_owned[1] = true; ++(ss->tm_count); } ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix @@ -404,7 +422,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, QPMS_CRASHING_MALLOC(ss->fecv_pstarts, ss->p_count * sizeof(size_t)); for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { ss->fecv_pstarts[pi] = ss->fecv_size; - ss->fecv_size += ss->tm[ss->p[pi].tmatrix_id]->spec->n; // That's a lot of dereferencing! + ss->fecv_size += ssw->tm[ss->p[pi].tmatrix_id]->spec->n; // That's a lot of dereferencing! } QPMS_CRASHING_MALLOC(ss->saecv_sizes, sizeof(size_t) * sym->nirreps); @@ -435,13 +453,14 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, } ss->c = qpms_trans_calculator_init(lMax, normalisation); - return ss; + return ssw; } void qpms_scatsys_free(qpms_scatsys_t *ss) { if(ss) { free(ss->tm); + free(ss->tmg); free(ss->p); free(ss->fecv_pstarts); free(ss->tm_sym_map); diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 8304d06..42dc309 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1028,9 +1028,11 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { struct qpms_tmatrix_operation_compose_chain * const o = &(f->op.compose_chain); if(o->opmem) { - for(size_t i = 0; i < o->n; ++i) - qpms_tmatrix_operation_clear(&(o->opmem[i])); + for(size_t i = 0; i < o->n; ++i) + if(o->ops_owned[i]) + qpms_tmatrix_operation_clear(o->ops[i]); free(o->opmem); + free(o->ops_owned); } free(o->ops); } @@ -1081,10 +1083,12 @@ void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *dest, const qpms_tmat struct qpms_tmatrix_operation_compose_chain * const o = &(dest->op.compose_chain); QPMS_CRASHING_MALLOC(o->ops, o->n * sizeof(*(o->ops))); + QPMS_CRASHING_MALLOC(o->ops_owned, o->n * sizeof(*(o->ops_owned))); QPMS_CRASHING_MALLOC(o->opmem, o->n * sizeof(*(o->opmem))); for(size_t i = 0; i < o->n; ++i) { qpms_tmatrix_operation_copy(o->opmem + i, src->op.compose_chain.ops[i]); o->ops[i] = o->opmem + i; + o->ops_owned[i] = true; } } break; @@ -1093,6 +1097,31 @@ void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *dest, const qpms_tmat } } +void qpms_tmatrix_operation_compose_chain_init(qpms_tmatrix_operation_t *dest, + size_t nops, size_t opmem_size) { + if (nops == 0) { + QPMS_WARN("Tried to create a composed (chain) operation of zero size;" + "setting to no-op instead."); + *dest = qpms_tmatrix_operation_noop; + } + if (nops < opmem_size) + QPMS_WARN("Allocating buffer for %zu operations, in a chained operation of" + " only %zu elemens, that does not seem to make sense.", opmem_size, nops); + dest.typ = QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN; + struct qpms_tmatrix_operation_compose_chain * const o = + &(dest->op.compose_chain); + o->n = nops; + QPMS_CRASHING_MALLOC(o->ops, nops * sizeof(*(o->ops))); + if (opmem_size != 0) { + QPMS_CRASHING_CALLOC(o->ops_owned, o->n, sizeof(_Bool)); + QPMS_CRASHING_MALLOC(o->opmem, opmem_size * sizeof(*(o->opmem))); + } else { + o->ops_owned = NULL; + o->opmem = NULL; + } +} + + qpms_tmatrix_t *qpms_tmatrix_apply_operation( const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) { diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index fc54edd..98ccd50 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -579,6 +579,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.) + _Bool *ops_owned; ///< True for all sub operations owned by this and saved in opmem. NULL if opmem is NULL. }; /// Specifies an elementwise complex multiplication of type \f$ T'_{ij} = M_{ij}T_{ij} \f$ for qpms_tmatrix_operation_t. @@ -616,6 +617,8 @@ typedef struct qpms_tmatrix_operation_t { } op; ///< Operation data; actual type is determined by \a typ. } qpms_tmatrix_operation_t; +static const qpms_tmatrix_operation_t qpms_tmatrix_operation_noop = {.typ = QPMS_TMATRIX_OPERATION_NOOP}; + /// Apply an operation on a T-matrix, returning a newly allocated result. qpms_tmatrix_t *qpms_tmatrix_apply_operation(const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig); @@ -634,6 +637,13 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f); /** Makes copies of all the internal data and takes ownership over them if needed */ void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *target, const qpms_tmatrix_operation_t *src); +/// Inits a new "chain" of composed operations, some of which might be owned. +void qpms_tmatrix_operation_compose_chain_init( + qpms_tmatrix_operation_t *target, ///< The operation structure that will be set to the chain. + size_t nops, ///< Number of chained operations (length of the \a ops array) + size_t opmem_size ///< Size of the own operations buffer (length of the \a opmem array) + ); + #if 0 // Abstract types that describe T-matrix/particle/scatsystem symmetries // To be implemented later. See also the thoughts in the beginning of groups.h.