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

Rewriting scatsystem: compiles without errors now. 

Former-commit-id: cd68b0feaef7181874d94dc535fd2cc9bc89e518
This commit is contained in:
Marek Nečada 2020-01-10 11:43:35 +02:00
parent 541af5a984
commit 5a98b91b3d
4 changed files with 49 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
qpms/scatsystem.c
View File

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

12
qpms/scatsystem.h
View File

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

16
qpms/tmatrices.c
View File

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

1
qpms/tmatrices.h
View File

@ -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.
};