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

New qpms_scatsys_apply_symmetry kinda done? 

Former-commit-id: 49a7a7984af6ad6e0e5ec1b5cc7b61ac06f81b8d
This commit is contained in:
Marek Nečada 2020-01-09 10:42:49 +02:00
parent 7e57c3cc81
commit 5dd93235f0
3 changed files with 80 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
qpms/scatsystem.c
View File

@ -176,7 +176,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig,
// Copy the qpms_tmatrix_fuction_t from orig // Copy the qpms_tmatrix_fuction_t from orig
ss->tmg_count = orig->tmg_count; 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))); memcpy(ss->tmg, orig->tmg, ss->tmg_count * sizeof(*(ss->tmg)));
// Allocate T-matrix, particle and particle orbit info arrays // Allocate T-matrix, particle and particle orbit info arrays
@ -198,7 +198,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig,
// Evaluate the original T-matrices at omega // Evaluate the original T-matrices at omega
qpms_tmatrix_t **tm_orig_omega; qpms_tmatrix_t **tm_orig_omega;
QPMS_CRASHING_MALLOC(tm_orig_omega, orig->tmg_count * sizeof(*tm_orig_omega)); //TODO free 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) 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); 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); ssw->wavenumber = qpms_wavenumber(omega, ssw->medium);
// we will be using ss->tm_capacity also for ssw->tm // we will be using ss->tm_capacity also for ssw->tm
QPMS_CRASHING_MALLOC(ssw->tm, ss->tm_capacity * sizeof(*(ssw->tm))); // returned 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. 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! 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; ss->tm_count = 0;
for (qpms_ss_tmi_t i = 0; i < orig->tm_count; ++i) { 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; qpms_ss_tmi_t j;
for (j = 0; j < ss->tm_count; ++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; break;
} }
if (j == ss->tm_count) { // duplicity not found, copy the t-matrix if (j == ss->tm_count) { // duplicity not found, save both the "abstract" and "at omega" T-matrices
ss->tm[j] = qpms_tmatrix_copy(orig->tm[i]); ss->tm[j].op = qpms_tmatrix_operation_copy(orig->tm[j].op);
ss->max_bspecn = MAX(ss->tm[j]->spec->n, ss->max_bspecn); 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, ss->tm[j]->spec->lMax);
++(ss->tm_count); ++(ss->tm_count);
} }
else qpms_tmatrix_free(ti);
tm_dupl_remap[i] = j; tm_dupl_remap[i] = j;
if (normalisation == QPMS_NORMALISATION_UNDEF) 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. // 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." "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 // Copy particles, remapping the t-matrix indices
for (qpms_ss_pi_t i = 0; i < orig->p_count; ++(i)) { for (qpms_ss_pi_t i = 0; i < orig->p_count; ++(i)) {
ss->p[i] = orig->p[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 // Extend the T-matrices list by the symmetry operations
for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi)
for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){ for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){
const size_t d = ss->tm[tmi]->spec->n; const size_t d = ssw->tm[tmi]->spec->n;
complex double M[d][d]; // transformation matrix complex double *m;
qpms_irot3_uvswfi_dense(M[0], ss->tm[tmi]->spec, sym->rep3d[gmi]); 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_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ss->tm[tmi], M[0]);
qpms_ss_tmi_t tmj; qpms_ss_tmi_t tmj;
for (tmj = 0; tmj < ss->tm_count; ++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, ss->tm[tmj], tol->rtol, tol->atol))
break; break;
if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists 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); qpms_tmatrix_free(transformed);
} else { // MISS, save the matrix and increment the count } else { // MISS, save the matrix (also the "abstract" one)
ss->tm[ss->tm_count] = transformed; 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_count);
} }
ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix 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)); QPMS_CRASHING_MALLOC(ss->fecv_pstarts, ss->p_count * sizeof(size_t));
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
ss->fecv_pstarts[pi] = ss->fecv_size; 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); 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); ss->c = qpms_trans_calculator_init(lMax, normalisation);
return ss; return ssw;
} }
void qpms_scatsys_free(qpms_scatsys_t *ss) { void qpms_scatsys_free(qpms_scatsys_t *ss) {
if(ss) { if(ss) {
free(ss->tm); free(ss->tm);
free(ss->tmg);
free(ss->p); free(ss->p);
free(ss->fecv_pstarts); free(ss->fecv_pstarts);
free(ss->tm_sym_map); free(ss->tm_sym_map);

31
qpms/tmatrices.c
View File

@ -1029,8 +1029,10 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
&(f->op.compose_chain); &(f->op.compose_chain);
if(o->opmem) { if(o->opmem) {
for(size_t i = 0; i < o->n; ++i) for(size_t i = 0; i < o->n; ++i)
qpms_tmatrix_operation_clear(&(o->opmem[i])); if(o->ops_owned[i])
qpms_tmatrix_operation_clear(o->ops[i]);
free(o->opmem); free(o->opmem);
free(o->ops_owned);
} }
free(o->ops); 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 = struct qpms_tmatrix_operation_compose_chain * const o =
&(dest->op.compose_chain); &(dest->op.compose_chain);
QPMS_CRASHING_MALLOC(o->ops, o->n * sizeof(*(o->ops))); 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))); QPMS_CRASHING_MALLOC(o->opmem, o->n * sizeof(*(o->opmem)));
for(size_t i = 0; i < o->n; ++i) { for(size_t i = 0; i < o->n; ++i) {
qpms_tmatrix_operation_copy(o->opmem + i, src->op.compose_chain.ops[i]); qpms_tmatrix_operation_copy(o->opmem + i, src->op.compose_chain.ops[i]);
o->ops[i] = o->opmem + i; o->ops[i] = o->opmem + i;
o->ops_owned[i] = true;
} }
} }
break; 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( qpms_tmatrix_t *qpms_tmatrix_apply_operation(
const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) { const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) {

10
qpms/tmatrices.h
View File

@ -579,6 +579,7 @@ struct qpms_tmatrix_operation_compose_chain {
size_t n; ///< Number of operations in ops; size_t n; ///< Number of operations in ops;
const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers owned by this.) 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.) 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. /// 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. } op; ///< Operation data; actual type is determined by \a typ.
} qpms_tmatrix_operation_t; } 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. /// 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); 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 */ /** 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); 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 #if 0
// Abstract types that describe T-matrix/particle/scatsystem symmetries // Abstract types that describe T-matrix/particle/scatsystem symmetries
// To be implemented later. See also the thoughts in the beginning of groups.h. // To be implemented later. See also the thoughts in the beginning of groups.h.