From d53f2964f00ac03951aa5096450622afea2fa63a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 19 Dec 2019 13:50:12 +0200 Subject: [PATCH 01/29] WIP abstract t-matrices Former-commit-id: 8c573ac3a62bf92195246d6eb95f95df240c48a1 --- qpms/scatsystem.h | 6 ++++-- qpms/tmatrices.h | 21 +++++++++++++++++++++ 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 16578a6..4b0fb9f 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -124,8 +124,8 @@ typedef struct qpms_ss_particle_orbitinfo { struct qpms_trans_calculator; typedef struct qpms_scatsys_t { - // TODO does bspec belong here? - qpms_tmatrix_t **tm; ///< T-matrices in the system + qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. + qpms_abstract_tmatrix_t **tm; ///< T-matrices in the system qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[]. qpms_particle_tid_t *p; ///< Particles. @@ -167,6 +167,8 @@ typedef struct qpms_scatsys_t { struct qpms_trans_calculator *c; } qpms_scatsys_t; +typedef struct qpms_scatsys_at_omega_t { + /// Convenience function to access pi'th particle's bspec. static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) { return ss->tm[ss->p[pi].tmatrix_id]->spec; diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index c208aa0..ba29e48 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -8,6 +8,8 @@ #include "materials.h" #include + + struct qpms_finite_group_t; typedef struct qpms_finite_group_t qpms_finite_group_t; @@ -503,6 +505,25 @@ qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target, ); +/// An "abstract" t-matrix, contains a T-matrix generator instead of actual data. +typedef struct qpms_tmatrix_function_t { + /** \brief VSWF basis specification, NOT owned by qpms_tmatrix_t by default. + * + * Usually not checked for meaningfulness by the functions (methods), + * so the caller should take care that \a spec->ilist does not + * contain any duplicities and that for each wave with order \a m + * there is also one with order \a −m. + */ + const qpms_vswf_set_spec_t *spec; + const qpms_tmatrix_generator_t *gen; ///< A T-matrix generator function. +} qpms_tmatrix_function_t; + +/// A recepy to create another T-matrices from qpms_tmatrix_fuction_t by symmetry operations. +typedef struct qpms_derived_tmatrix_function_t { + const qpms_tmatrix_function_t *t; +} + + #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. From d17a5e5eea4cb7e91856185a31449008533ce1e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sat, 21 Dec 2019 11:36:01 +0200 Subject: [PATCH 02/29] WIP data structure redefinitions. Former-commit-id: 5e2baffb4a47657233e792407630507ba611b129 --- qpms/scatsystem.h | 15 +++++++++++++++ qpms/tmatrices.h | 2 +- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 4b0fb9f..602ac4a 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -167,7 +167,22 @@ typedef struct qpms_scatsys_t { struct qpms_trans_calculator *c; } qpms_scatsys_t; + typedef struct qpms_scatsys_at_omega_t { + const qpms_scatsys_t *ss; ///< Parent scattering system. + /// T-matrices from \a ss, evaluated at \a omega. + /** The T-matrices are in the same order as in \a ss, + * i.e in the order corresponding to TODO WHAT E7ACTLY?? + */ + qpms_tmatrix_t **tm; + complex double omega; ///< Angular frequency + complex double eps; ///< Background medium relative electric permittivity + complex double mu; ///< Background medium relative magnetic permeability + complex double wavenumber; ///< Background medium wavenumber + complex double refractive_index; ///< Background medium refractive index +} qpms_scatsys_at_omega_t; + +void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *); /// Convenience function to access pi'th particle's bspec. static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) { diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index ba29e48..f583aa5 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -505,7 +505,7 @@ qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target, ); -/// An "abstract" t-matrix, contains a T-matrix generator instead of actual data. +/// An "abstract" T-matrix, contains a T-matrix generator instead of actual data. typedef struct qpms_tmatrix_function_t { /** \brief VSWF basis specification, NOT owned by qpms_tmatrix_t by default. * From dff8293e6d3f10704a8488bec82db430bb144138 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 6 Jan 2020 02:15:45 +0200 Subject: [PATCH 03/29] T-matrix general operation type definitions and destructor. Former-commit-id: 491a4d8ad602a7252aa9f4446b55c7c905102de9 --- qpms/tmatrices.c | 46 +++++++++++++++++++++++ qpms/tmatrices.h | 98 +++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 142 insertions(+), 2 deletions(-) diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 2811b87..1e2c529 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -991,3 +991,49 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, return QPMS_SUCCESS; } +void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { + switch(f->typ) { + case QPMS_TMATRIX_OPERATION_LRMATRIX: + if(f->op.lrmatrix.owns_m) + free(f->op.lrmatrix.m); + break; + case QPMS_TMATRIX_OPERATION_SCMULZ: + if(f->op.scmulz.owns_m) + free(f->op.scmulz.m); + break; + case QPMS_TMATRIX_OPERATION_IROT3: + break; + case QPMS_TMATRIX_OPERATION_IROT3ARR: + if(f->op.irot3arr.owns_ops) + free(f->op.irot3arr.ops); + break; + case QPMS_TMATRIX_OPERATION_FINITE_GROUP: // Group never owned + break; + case QPMS_TMATRIX_OPERATION_COMPOSE_SUM: + { + struct qpms_tmatrix_operation_compose_sum * const o = + &(f->op.compose_sum); + if(o->opmem) { + for(size_t i = 0; i < o->n; ++i) + qpms_tmatrix_operation_clear(o->ops[i]); + free(o->opmem); + } + free(o->ops); + } + break; + case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN: + { + 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->ops[i]); + free(o->opmem); + } + free(o->ops); + } + break; + default: + QPMS_WTF; + } +} diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index f583aa5..83aa71b 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -518,10 +518,104 @@ typedef struct qpms_tmatrix_function_t { const qpms_tmatrix_generator_t *gen; ///< A T-matrix generator function. } qpms_tmatrix_function_t; -/// A recepy to create another T-matrices from qpms_tmatrix_fuction_t by symmetry operations. +/// Specifies different kinds of operations done on T-matrices +typedef enum { + QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f[ T' = MTM^\dagger \f]; see @ref qpms_tmatrix_operation_lrmatrix. + QPMS_TMATRIX_OPERATION_IROT3, ///< Single rotoreflection specified by a qpms_irot3_t. + QPMS_TMATRIX_OPERATION_IROT3ARR, ///< Symmetrise using an array of rotoreflection; see @ref qpms_tmatrix_operation_irot3arr. + QPMS_TMATRIX_OPERATION_COMPOSE_SUM, ///< Apply several transformations and sum the results, see @ref qpms_tmatrix_operation_compose_sum. + QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN, ///< Chain several different transformations; see @ref qpms_tmatrix_operation_compose_chain. + QPMS_TMATRIX_OPERATION_SCMULZ, ///< Elementwise scalar multiplication with a complex matrix; see @ref qpms_tmatrix_operation_scmulz. + //QPMS_TMATRIX_OPERATION_POINTGROUP, ///< TODO the new point group in pointgroup.h + QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE ///< Legacy qpms_finite_group_t with filled rep3d; see @ref qpms_tmatrix_operation_finite_group. +} qpms_tmatrix_operation_kind_t; + +/// General matrix transformation \f[ T' = MTM^\dagger \f] spec for qpms_tmatrix_operation_t. +struct qpms_tmatrix_operation_lrmatrix { + /// Raw matrix data of \a M in row-major order. + /** The matrix must be taylored for the given bspec! */ + const complex double *m; + bool owns_m; ///< Whether \a m is owned by this; +}; + +struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations. + +/// Specifies a composed operation of type \f$ T' = c\Sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t. +struct qpms_tmatrix_operation_compose_sum { + size_t n; ///< Number of operations in ops; + const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers array \a ops[] always owned by this.) + double factor; ///< Overall factor \a c. + /// (Optional) operations buffer into which elements of \a ops point. + /** Owned by this or NULL. If not NULL, all \a ops members are assumed to point into + * the memory held by \a opmem and to be properly initialised. + * Each \a ops member has to point to _different_ elements of \a opmem. + */ + qpms_tmatrix_operation_t *opmem; +}; + +/// Specifies a composed operation of type \f$ T' = f_{n-1}(f_{n-2}(\dots f_0(T)\dots))) \f$ for qpms_tmatrix_operation_t. +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.) + qpms_tmatrix_operation_t *opmem; ///< (Optional) operations buffer into which elements of \a ops point. (Owned by this or NULL.) +}; + +/// Specifies an elementwise complex multiplication of type \f$ T'_{ij} = M_{ij}T_{ij} \f$ for qpms_tmatrix_operation_t. +struct qpms_tmatrix_operation_scmulz { + /// Raw matrix data of \a M in row-major order. + /** The matrix must be taylored for the given bspec! */ + const complex double *m; + bool owns_m; ///< Whether \a m is owned by this. +}; + +/// Specifies a symmetrisation using a set of rotoreflections (with equal weights) for qpms_tmatrix_operation_t. +/** Internally, this is evaluated using a call to qpms_symmetrise_tmdata_irot3arr() + * or qpms_symmetrise_tmdata_irot3arr_inplace(). */ +struct qpms_tmatrix_operation_irot3arr { + size_t n; ///< Number of rotoreflections; + const qpms_irot3_t *ops; ///< Rotoreflection array of size \a n. + bool owns_ops; ///< Whether \a ops array is owned by this. +}; + +/// A generic T-matrix transformation operator. +typedef struct qpms_tmatrix_operation_t { + /// Specifies the operation kind to be performed and which type \op actually contains. + qpms_tmatrix_operation_kind_t typ; + union { + struct qpms_tmatrix_operation_lrmatrix lrmatrix; + struct qpms_tmatrix_operation_scmulz scmulz; + qpms_irot3 irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3 + struct qpms_tmatrix_operation_finite_group finite_group; + struct qpms_tmatrix_operation_irot3arr irot3arr; + struct qpms_tmatrix_operation_compose_sum compose_sum; + struct qpms_tmatrix_operation_compose_chain compose_chain; + /// Finite group for QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE. + /** Not owned by this; \a rep3d must be filled. */ + const qpms_finite_group_t *finitegroup; + } op; ///< Operation data; actual type is determined by \a typ. +} qpms_tmatrix_operation_t; + +/// Apply an operation on a T-matrix, returning a newly allocated result. +// TODO IMPLEMENT +qpms_tmatrix_t *qpms_tmatrix_apply_operation(const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig); + +/// Apply an operation on a T-matrix and replace it with the result. +// TODO IMPLEMENT +qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(const qpms_tmatrix_operation_t *op, qpms_tmatrix_t *orig); + +/// (Recursively) deallocates internal data of qpms_tmatrix_operation_t. +/** It does NOT clear the memory pointed to by op it self, but only + * heap-allocated data of certain operations, if labeled as owned by it. + * In case of compose operations, the recursion stops if the children are + * not owned by them (the opmem pointer is NULL). + */ +void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f); + +/// A recepy to create another T-matrices from qpms_tmatrix_fuction_t by symmetry (or other) operations. typedef struct qpms_derived_tmatrix_function_t { const qpms_tmatrix_function_t *t; -} + const qpms_tmatrix_operation_t *op; +} qpms_derived_tmatrix_function_t; #if 0 From c2b4787cd5c5181d6c891228e493f6d9ce2b242e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 6 Jan 2020 17:01:40 +0200 Subject: [PATCH 04/29] Implementation of qpms_tmatrix_apply_operation(). Former-commit-id: 6773f0e1d02d5f929c2039f99338f08c25d0ccab --- qpms/CMakeLists.txt | 3 +++ qpms/scatsystem.h | 3 ++- qpms/tmatrices.c | 65 ++++++++++++++++++++++++++++++++++++++++++--- qpms/tmatrices.h | 23 ++++++++-------- 4 files changed, 78 insertions(+), 16 deletions(-) diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 2a63a3f..d90239c 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -3,6 +3,9 @@ find_package(GSL 2.0 REQUIRED) find_package(BLAS REQUIRED) find_package(LAPACK REQUIRED) +# disable an annoying warning that gives false positives probably due to a bug in gcc +# and other not very relevant warnings +set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wno-int-in-bool-context -Wno-comment") #includes set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS}) diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 602ac4a..5c1f23f 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -122,9 +122,10 @@ typedef struct qpms_ss_particle_orbitinfo { } qpms_ss_particle_orbitinfo_t; struct qpms_trans_calculator; +struct qpms_epsmu_generator_t; typedef struct qpms_scatsys_t { - qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. + struct qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. qpms_abstract_tmatrix_t **tm; ///< T-matrices in the system qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[]. diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 1e2c529..eed05e9 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1007,7 +1007,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { if(f->op.irot3arr.owns_ops) free(f->op.irot3arr.ops); break; - case QPMS_TMATRIX_OPERATION_FINITE_GROUP: // Group never owned + case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: // Group never owned break; case QPMS_TMATRIX_OPERATION_COMPOSE_SUM: { @@ -1015,7 +1015,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { &(f->op.compose_sum); if(o->opmem) { for(size_t i = 0; i < o->n; ++i) - qpms_tmatrix_operation_clear(o->ops[i]); + qpms_tmatrix_operation_clear(&(o->opmem[i])); free(o->opmem); } free(o->ops); @@ -1027,7 +1027,7 @@ 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) - qpms_tmatrix_operation_clear(o->ops[i]); + qpms_tmatrix_operation_clear(&(o->opmem[i])); free(o->opmem); } free(o->ops); @@ -1037,3 +1037,62 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { QPMS_WTF; } } + +qpms_tmatrix_t *qpms_tmatrix_apply_operation( + const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) { + // Certain operations could be optimized, but the effect would be marginal. + qpms_tmatrix_t *res = qpms_tmatrix_copy(orig); + return qpms_tmatrix_apply_operation_inplace(f, res); +} + +static qpms_tmatrix_t *qtao_compose_sum_inplace(qpms_tmatrix_t *T, + const struct qpms_tmatrix_operation_compose_sum *cs) { + qpms_tmatrix_t *tmp_target = qpms_tmatrix_init(T->spec); + qpms_tmatrix_t *sum = qpms_tmatrix_init(T->spec); + for (size_t i = 0; i < cs->n; ++i) { + memcpy(tmp_target->m, T->m, SQ(T->spec->n) * sizeof(complex double)); + QPMS_ENSURE(qpms_tmatrix_apply_operation_inplace(cs->ops[i] , tmp_target), + "Got NULL pointer from qpms_tmatrix_apply_operation_inplace, hupsis!"); + for (size_t j = 0; j < SQ(T->spec->n); ++j) + sum->m[j] += tmp_target->m[j]; + } + for(size_t j = 0; j < SQ(T->spec->n); ++j) + T->m[j] = sum->m[j] * cs->factor; + qpms_tmatrix_free(sum); + qpms_tmatrix_free(tmp_target); + return T; +} + +static qpms_tmatrix_t *qtao_compose_chain_inplace(qpms_tmatrix_t *T, + const struct qpms_tmatrix_operation_compose_chain *cc) { + for(size_t i = 0; i < cc->n; ++i) + qpms_tmatrix_apply_operation_inplace(cc->ops[i], T); + return T; +} + +static qpms_tmatrix_t *qtao_scmulz_inplace(qpms_tmatrix_t *T, + const struct qpms_tmatrix_operation_scmulz *s) { + for(size_t i = 0; i < SQ(T->spec->n); ++i) + T->m[i] *= s->m[i]; + return T; +} + +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_LRMATRIX: + return qpms_tmatrix_apply_symop_inplace(T, f->op.lrmatrix.m); + case QPMS_TMATRIX_OPERATION_IROT3: + return qpms_tmatrix_symmetrise_irot3arr_inplace(T, 1, &(f->op.irot3)); + case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: + return qpms_tmatrix_symmetrise_finite_group_inplace(T, f->op.finite_group); + case QPMS_TMATRIX_OPERATION_COMPOSE_SUM: + return qtao_compose_sum_inplace(T, &(f->op.compose_sum)); + case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN: + return qtao_compose_chain_inplace(T, &(f->op.compose_chain)); + case QPMS_TMATRIX_OPERATION_SCMULZ: + return qtao_scmulz_inplace(T, &(f->op.scmulz)); + default: + QPMS_WTF; + } +} diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index 83aa71b..163c6d3 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -235,7 +235,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( complex double *qpms_apply_tmatrix( complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated. const complex double *a, ///< Incident field coefficient array of size T->spec->n. - const qpms_tmatrix_t *T + const qpms_tmatrix_t *T ///< T-matrix \a T to apply. ); /// Generic T-matrix generator function that fills a pre-initialised qpms_tmatrix_t given a frequency. @@ -429,7 +429,7 @@ qpms_arc_function_retval_t qpms_arc_cylinder(double theta, ); /// Arc parametrisation of spherical particle; for qpms_arc_function_t. -/** Useful mostly only for benchmarks, as one can use the Mie-Lorentz solution. */ +/** Useful mostly only for benchmarks or debugging, as one can use the Mie-Lorentz solution. */ qpms_arc_function_retval_t qpms_arc_sphere(double theta, const void *R; ///< Points to double containing particle's radius. ); @@ -520,7 +520,7 @@ typedef struct qpms_tmatrix_function_t { /// Specifies different kinds of operations done on T-matrices typedef enum { - QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f[ T' = MTM^\dagger \f]; see @ref qpms_tmatrix_operation_lrmatrix. + QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f$ T' = MTM^\dagger \f$; see @ref qpms_tmatrix_operation_lrmatrix. QPMS_TMATRIX_OPERATION_IROT3, ///< Single rotoreflection specified by a qpms_irot3_t. QPMS_TMATRIX_OPERATION_IROT3ARR, ///< Symmetrise using an array of rotoreflection; see @ref qpms_tmatrix_operation_irot3arr. QPMS_TMATRIX_OPERATION_COMPOSE_SUM, ///< Apply several transformations and sum the results, see @ref qpms_tmatrix_operation_compose_sum. @@ -534,13 +534,13 @@ typedef enum { struct qpms_tmatrix_operation_lrmatrix { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ - const complex double *m; + complex double *m; bool owns_m; ///< Whether \a m is owned by this; }; struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations. -/// Specifies a composed operation of type \f$ T' = c\Sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t. +/// Specifies a composed operation of type \f$ T' = c\sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t. struct qpms_tmatrix_operation_compose_sum { size_t n; ///< Number of operations in ops; const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers array \a ops[] always owned by this.) @@ -550,21 +550,21 @@ struct qpms_tmatrix_operation_compose_sum { * the memory held by \a opmem and to be properly initialised. * Each \a ops member has to point to _different_ elements of \a opmem. */ - qpms_tmatrix_operation_t *opmem; + struct qpms_tmatrix_operation_t *opmem; }; /// Specifies a composed operation of type \f$ T' = f_{n-1}(f_{n-2}(\dots f_0(T)\dots))) \f$ for qpms_tmatrix_operation_t. 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.) - 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.) }; /// Specifies an elementwise complex multiplication of type \f$ T'_{ij} = M_{ij}T_{ij} \f$ for qpms_tmatrix_operation_t. struct qpms_tmatrix_operation_scmulz { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ - const complex double *m; + complex double *m; bool owns_m; ///< Whether \a m is owned by this. }; @@ -573,7 +573,7 @@ struct qpms_tmatrix_operation_scmulz { * or qpms_symmetrise_tmdata_irot3arr_inplace(). */ struct qpms_tmatrix_operation_irot3arr { size_t n; ///< Number of rotoreflections; - const qpms_irot3_t *ops; ///< Rotoreflection array of size \a n. + qpms_irot3_t *ops; ///< Rotoreflection array of size \a n. bool owns_ops; ///< Whether \a ops array is owned by this. }; @@ -584,14 +584,13 @@ typedef struct qpms_tmatrix_operation_t { union { struct qpms_tmatrix_operation_lrmatrix lrmatrix; struct qpms_tmatrix_operation_scmulz scmulz; - qpms_irot3 irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3 - struct qpms_tmatrix_operation_finite_group finite_group; + qpms_irot3_t irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3 struct qpms_tmatrix_operation_irot3arr irot3arr; struct qpms_tmatrix_operation_compose_sum compose_sum; struct qpms_tmatrix_operation_compose_chain compose_chain; /// Finite group for QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE. /** Not owned by this; \a rep3d must be filled. */ - const qpms_finite_group_t *finitegroup; + const qpms_finite_group_t *finite_group; } op; ///< Operation data; actual type is determined by \a typ. } qpms_tmatrix_operation_t; From d31d8737b85fccfa65a87c827596f700a845b957 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 7 Jan 2020 16:57:59 +0200 Subject: [PATCH 05/29] WIP scatsystem update for "abstract" T-matrices. Former-commit-id: 7f723a0f459f263e12282edfb1e8deb440650880 --- qpms/qpms_types.h | 3 +++ qpms/scatsystem.h | 54 +++++++++++++++++++++++++++++++++++++---------- qpms/tmatrices.c | 4 ++++ qpms/tmatrices.h | 7 +----- 4 files changed, 51 insertions(+), 17 deletions(-) diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index fd73ce0..0f599b6 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -297,6 +297,9 @@ typedef struct qpms_vswf_set_spec_t { /// T-matrix index used in qpms_scatsys_t and related structures. typedef int32_t qpms_ss_tmi_t; +/// T-matrix generator index used in qpms_scatsys_t and related structures. +typedef int32_t qpms_ss_tmgi_t; + /// Particle index used in qpms_scatsys_t and related structures. typedef int32_t qpms_ss_pi_t; diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 5c1f23f..9ab7177 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -121,14 +121,32 @@ typedef struct qpms_ss_particle_orbitinfo { qpms_ss_orbit_pi_t p; ///< Order (sija, ei rankki) of the particle inside that orbit type. } qpms_ss_particle_orbitinfo_t; +/// 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. + qpms_tmatrix_operation_t *op; ///< Operation to derive this particular T-matrix. +} qpms_ss_derived_tmatrix_t; + + struct qpms_trans_calculator; struct qpms_epsmu_generator_t; typedef struct qpms_scatsys_t { struct qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. - qpms_abstract_tmatrix_t **tm; ///< T-matrices in the system - qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices + + /// (Template) T-matrix functions in the system. + /** The qpms_abstract_tmatrix_t objects (onto which this array member point) + * are NOT owned by this and must be kept alive for the whole lifetime + * of all qpms_scatsys_t objects that are built upon them. + */ + qpms_abstract_tmatrix_t **tmg; + qpms_ss_tmgi_t tmg_count; ///< Number of all different original T-matrix generators in the system. + + /// All the different T-matrix functions in the system, including those derived from \a tmg elements by symmetries. + qpms_ss_derived_tmatrix_t *tm; + qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices in the system (length of tm[]). qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[]. + qpms_particle_tid_t *p; ///< Particles. qpms_ss_pi_t p_count; ///< Size of particles array. qpms_ss_pi_t p_capacity; ///< Capacity of p[]. @@ -173,7 +191,7 @@ typedef struct qpms_scatsys_at_omega_t { const qpms_scatsys_t *ss; ///< Parent scattering system. /// T-matrices from \a ss, evaluated at \a omega. /** The T-matrices are in the same order as in \a ss, - * i.e in the order corresponding to TODO WHAT E7ACTLY?? + * i.e in the order corresponding to \a ss->tm. */ qpms_tmatrix_t **tm; complex double omega; ///< Angular frequency @@ -190,18 +208,32 @@ static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t return ss->tm[ss->p[pi].tmatrix_id]->spec; } -/// Creates a new scatsys by applying a symmetry group, copying particles if needed. -/** In fact, it copies everything except the vswf set specs, so keep them alive until scatsys is destroyed. - * The following fields must be filled: - * orig->tm - * orig->tm_count - * orig->p - * orig->p_count +#if 0 //TODO!!! +/// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed. +/** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances, + * so keep them alive until scatsys is destroyed. + * + * The following fields must be filled in the "proto- scattering system" \a orig: + * + * * orig->tmg – The pointers are copied to the new qpms_scatsys_t instance; + * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting + * qpms_scatsys_t instances are properly destroyed. The pointers from orig->tmg, however, are copied. + * * orig->tmg_count + * * orig->tm – Must be filled, although the operations will typically be identities + * (QPMS_TMATRIX_OPERATION_NOOP). N.B. these NOOPs might be replaced with some symmetrisation operation + * in the resulting "full" qpms_scatsys_t instance. + * * orig->tm_count + * * orig->p + * * orig->p_count */ -qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym); +qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym); +#endif /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load. void qpms_scatsys_free(qpms_scatsys_t *s); +/// Destroys the result of qpms_scatsys_eval_at_omega() +void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega *so); + /// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis. /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient. * diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index eed05e9..e13ec1c 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -993,6 +993,8 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { switch(f->typ) { + case QPMS_TMATRIX_OPERATION_NOOP: + break; case QPMS_TMATRIX_OPERATION_LRMATRIX: if(f->op.lrmatrix.owns_m) free(f->op.lrmatrix.m); @@ -1080,6 +1082,8 @@ static qpms_tmatrix_t *qtao_scmulz_inplace(qpms_tmatrix_t *T, 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; 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 163c6d3..dfa282f 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -520,6 +520,7 @@ typedef struct qpms_tmatrix_function_t { /// Specifies different kinds of operations done on T-matrices typedef enum { + QPMS_TMATRIX_OPERATION_NOOP, ///< Identity operation. QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f$ T' = MTM^\dagger \f$; see @ref qpms_tmatrix_operation_lrmatrix. QPMS_TMATRIX_OPERATION_IROT3, ///< Single rotoreflection specified by a qpms_irot3_t. QPMS_TMATRIX_OPERATION_IROT3ARR, ///< Symmetrise using an array of rotoreflection; see @ref qpms_tmatrix_operation_irot3arr. @@ -610,12 +611,6 @@ qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(const qpms_tmatrix_operatio */ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f); -/// A recepy to create another T-matrices from qpms_tmatrix_fuction_t by symmetry (or other) operations. -typedef struct qpms_derived_tmatrix_function_t { - const qpms_tmatrix_function_t *t; - const qpms_tmatrix_operation_t *op; -} qpms_derived_tmatrix_function_t; - #if 0 // Abstract types that describe T-matrix/particle/scatsystem symmetries From e1a63892321c8ba6d71b1fed481522e9cfda8437 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 8 Jan 2020 14:39:29 +0200 Subject: [PATCH 06/29] Some new convenience functions and types. Former-commit-id: 7701cd8ee779e06ba18d6e19bfe650bd9465487f --- qpms/materials.h | 6 ++++++ qpms/qpms_types.h | 1 + qpms/tmatrices.h | 21 +++++++++++++++++++-- qpms/tolerances.h | 10 ++++++++++ 4 files changed, 36 insertions(+), 2 deletions(-) create mode 100644 qpms/tolerances.h diff --git a/qpms/materials.h b/qpms/materials.h index 0db7372..ce27d76 100644 --- a/qpms/materials.h +++ b/qpms/materials.h @@ -24,6 +24,12 @@ typedef struct qpms_epsmu_generator_t { const void *params; } qpms_epsmu_generator_t; +/// Convenience function for generating material properties at given frequency. +static inline qpms_epsmu_t qpms_epsmu_generator_eval( + qpms_epsmu_generator_t gen, complex double omega) { + return gen.function(omega, gen.params); +} + /// Constant optical property "generator" for qpms_epsmu_generator_t. qpms_epsmu_t qpms_epsmu_const_g(complex double omega, ///< Frequency ignored. const void *epsmu ///< Points to the qpms_epsmu_t to be returned. diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 0f599b6..c4c2bdf 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -410,6 +410,7 @@ typedef struct qpms_epsmu_t { complex double mu; ///< Relative permeability. } qpms_epsmu_t; +struct qpms_tolerance_spec_t; // See tolerances.c #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l)) #endif // QPMS_TYPES diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index dfa282f..c9f819f 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -19,8 +19,11 @@ static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){ } /// Initialises a zero T-matrix. +/** \sa qpms_tmatrix_init_from_generator() + * \sa qpms_tmatrix_init_from_function() */ qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec); + /// Copies a T-matrix, allocating new array for the T-matrix data. qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T); @@ -253,6 +256,16 @@ typedef struct qpms_tmatrix_generator_t { const void *params; ///< Parameter pointer passed to the function. } qpms_tmatrix_generator_t; +/// Initialises and evaluates a new T-matrix using a generator. +static inline qpms_tmatrix_t *qpms_tmatrix_init_from_generator( + const qpms_vswf_set_spec_t *bspec, + qpms_tmatrix_generator_t gen, + complex double omega) { + qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); + QPMS_ENSURE_SUCCESS(gen.function(t, omega, gen.params)); + return t; +} + /// Implementation of qpms_matrix_generator_t that just copies a constant matrix. /** N.B. this does almost no checks at all, so use it only for t-matrices with * the same base spec. @@ -518,6 +531,12 @@ typedef struct qpms_tmatrix_function_t { const qpms_tmatrix_generator_t *gen; ///< A T-matrix generator function. } qpms_tmatrix_function_t; +/// Convenience function to create a new T-matrix from qpms_tmatrix_function_t. +// FIXME the name is not very intuitive. +static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, complex double omega) { + return qpms_tmatrix_init_from_generator(f.spec, omega, f.gen); +} + /// Specifies different kinds of operations done on T-matrices typedef enum { QPMS_TMATRIX_OPERATION_NOOP, ///< Identity operation. @@ -596,11 +615,9 @@ typedef struct qpms_tmatrix_operation_t { } qpms_tmatrix_operation_t; /// Apply an operation on a T-matrix, returning a newly allocated result. -// TODO IMPLEMENT qpms_tmatrix_t *qpms_tmatrix_apply_operation(const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig); /// Apply an operation on a T-matrix and replace it with the result. -// TODO IMPLEMENT qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(const qpms_tmatrix_operation_t *op, qpms_tmatrix_t *orig); /// (Recursively) deallocates internal data of qpms_tmatrix_operation_t. diff --git a/qpms/tolerances.h b/qpms/tolerances.h new file mode 100644 index 0000000..e3259a9 --- /dev/null +++ b/qpms/tolerances.h @@ -0,0 +1,10 @@ +/*! \file tolerances.h */ +#ifndef QPMS_TOLERANCES_H +#define QPMS_TOLERANCES_H + +typedef struct qpms_tolerance_spec_t { + double atol; + double rtol; +} qpms_tolerance_spec_t; + +#endif // QPMS_TOLERANCES_H From 3bf263c4f34621294ebd97c7d8aaedd156fe4d69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 8 Jan 2020 15:18:25 +0200 Subject: [PATCH 07/29] Copying of T-matrix operations Former-commit-id: dea91f97e5e72146039868ab5f0c8ac5e7ea7a57 --- qpms/qpms_error.h | 9 ++++++++ qpms/tmatrices.c | 54 +++++++++++++++++++++++++++++++++++++++++++++++ qpms/tmatrices.h | 7 +++++- 3 files changed, 69 insertions(+), 1 deletion(-) diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 10fb7b8..393157c 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -56,6 +56,15 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types); (size_t) (size));\ } +#define QPMS_CRASHING_MALLOCPY(dest, src, size) {\ + (dest) = malloc(size);\ + if(QPMS_UNLIKELY(!(dest) && (size)))\ + qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__,\ + "Allocation of %zd bytes for " #dest " failed.",\ + (size_t) (size));\ + memcpy((dest), (src), (size));\ +} + #define QPMS_CRASHING_CALLOC(pointer, nmemb, size) {\ (pointer) = calloc((nmemb), (size));\ if(QPMS_UNLIKELY(!(pointer) && (nmemb) && (size)))\ diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index e13ec1c..8304d06 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1040,6 +1040,60 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) { } } +void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *dest, const qpms_tmatrix_operation_t *src) { + memcpy(dest, src, sizeof(qpms_tmatrix_operation_t)); + switch(dest->typ) { + case QPMS_TMATRIX_OPERATION_NOOP: + break; + case QPMS_TMATRIX_OPERATION_LRMATRIX: + QPMS_CRASHING_MALLOCPY(dest->op.lrmatrix.m, src->op.lrmatrix.m, + sizeof(complex double) * src->op.lrmatrix.m_size); + dest->op.lrmatrix.owns_m = true; + break; + case QPMS_TMATRIX_OPERATION_SCMULZ: + QPMS_CRASHING_MALLOCPY(dest->op.scmulz.m, src->op.scmulz.m, + sizeof(complex double) * src->op.scmulz.m_size); + dest->op.scmulz.owns_m = true; + break; + case QPMS_TMATRIX_OPERATION_IROT3: + break; + case QPMS_TMATRIX_OPERATION_IROT3ARR: + QPMS_CRASHING_MALLOCPY(dest->op.irot3arr.ops, src->op.irot3arr.ops, + src->op.irot3arr.n * sizeof(qpms_irot3_t)); + dest->op.irot3arr.owns_ops = true; + break; + case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: // Group never owned + break; + case QPMS_TMATRIX_OPERATION_COMPOSE_SUM: + { + struct qpms_tmatrix_operation_compose_sum * const o = + &(dest->op.compose_sum); + QPMS_CRASHING_MALLOC(o->ops, o->n * sizeof(*(o->ops))); + 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_sum.ops[i]); + o->ops[i] = o->opmem + i; + } + } + break; + case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN: + { + 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->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; + } + } + break; + default: + QPMS_WTF; + } +} + + qpms_tmatrix_t *qpms_tmatrix_apply_operation( const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) { // Certain operations could be optimized, but the effect would be marginal. diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index c9f819f..fc54edd 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -555,6 +555,7 @@ struct qpms_tmatrix_operation_lrmatrix { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ complex double *m; + size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double). bool owns_m; ///< Whether \a m is owned by this; }; @@ -584,7 +585,8 @@ struct qpms_tmatrix_operation_compose_chain { struct qpms_tmatrix_operation_scmulz { /// Raw matrix data of \a M in row-major order. /** The matrix must be taylored for the given bspec! */ - complex double *m; + complex double *m; + size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double). bool owns_m; ///< Whether \a m is owned by this. }; @@ -628,6 +630,9 @@ qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(const qpms_tmatrix_operatio */ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f); +/// (Recursively) copies an qpms_tmatrix_operation_t. +/** 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); #if 0 // Abstract types that describe T-matrix/particle/scatsystem symmetries From 7e57c3cc81ab2ae020d38f7ad6395689ac572e5f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 8 Jan 2020 16:00:09 +0200 Subject: [PATCH 08/29] WIP new scatsystem (keskeytetty n. scatsystem.c:214) Former-commit-id: 5d205f46f38f7b9e988bd03b8a9bff70b3986808 --- qpms/scatsystem.c | 133 +++++++++++++++++++++++++--------------------- qpms/scatsystem.h | 31 +++++++---- 2 files changed, 95 insertions(+), 69 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index a697d0f..660eac7 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -80,7 +80,7 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu qpms_ss_orbit_type_t * const ot_new = & (ss->orbit_types[ss->orbit_type_count]); ot_new->size = ot_current->size; - const qpms_vswf_set_spec_t *bspec = ss->tm[ot_current->tmatrices[0]]->spec; + const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_tmi(ss, ot_current->tmatrices[0]); const size_t bspecn = bspec->n; ot_new->bspecn = bspecn; @@ -143,7 +143,9 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu // Almost 200 lines. This whole thing deserves a rewrite! -qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) { +qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, + const qpms_finite_group_t *sym, complex double omega, + const qpms_tolerance_spec_t *tol) { // TODO check data sanity qpms_l_t lMax = 0; // the overall lMax of all base specs. @@ -171,32 +173,57 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp for (qpms_ss_pi_t j = 0; j < i; ++j) assert(!cart3_isclose(orig->p[i].pos, orig->p[j].pos, 0, QPMS_SCATSYS_LEN_RTOL * lenscale)); + + // 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 + memcpy(ss->tmg, orig->tmg, ss->tmg_count * sizeof(*(ss->tmg))); + // Allocate T-matrix, particle and particle orbit info arrays - qpms_scatsys_t *ss = malloc(sizeof(qpms_scatsys_t)); + qpms_scatsys_t *ss; + ss = QPMS_CRASHING_MALLOC(ss, sizeof(qpms_scatsys_t)); ss->lenscale = lenscale; ss->sym = sym; ss->tm_capacity = sym->order * orig->tm_count; - ss->tm = malloc(ss->tm_capacity * sizeof(qpms_tmatrix_t *)); + QPMS_CRASHING_MALLOC(ss->tm, ss->tm_capacity * sizeof(*(ss->tm))); ss->p_capacity = sym->order * orig->p_count; - ss->p = malloc(ss->p_capacity * sizeof(qpms_particle_tid_t)); - ss->p_orbitinfo = malloc(ss->p_capacity * sizeof(qpms_ss_particle_orbitinfo_t)); + QPMS_CRASHING_MALLOC(ss->p, ss->p_capacity * sizeof(qpms_particle_tid_t)); + QPMS_CRASHING_MALLOC(ss->p_orbitinfo, ss->p_capacity * sizeof(qpms_ss_particle_orbitinfo_t)); for (qpms_ss_pi_t pi = 0; pi < ss->p_capacity; ++pi) { ss->p_orbitinfo[pi].t = QPMS_SS_P_ORBITINFO_UNDEF; ss->p_orbitinfo[pi].p = QPMS_SS_P_ORBITINFO_UNDEF; } + // 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 + 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); + + // Evaluate the medium and derived T-matrices at omega. + qpms_scatsys_at_omega_t *ssw; + QPMS_CRASHING_MALLOC(ssw, sizeof(*ssw)); // returned + ssw->ss = ss; + ssw->omega = omega; + ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega); + 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 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 + 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_ss_tmi_t j; for (j = 0; j < ss->tm_count; ++j) - if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) { + if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], tol->rtol, tol->atol)) { break; } if (j == ss->tm_count) { // duplicity not found, copy the t-matrix @@ -222,7 +249,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp ss->p_count = orig->p_count; // allocate t-matrix symmetry map - ss->tm_sym_map = malloc(sizeof(qpms_ss_tmi_t) * sym->order * sym->order * ss->tm_count); + QPMS_CRASHING_MALLOC(ss->tm_sym_map, sizeof(qpms_ss_tmi_t) * sym->order * sym->order * ss->tm_count); // Extend the T-matrices list by the symmetry operations for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) @@ -233,7 +260,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp 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], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) + if (qpms_tmatrix_isclose(transformed, ss->tm[tmj], tol->rtol, tol->atol)) break; if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists qpms_tmatrix_free(transformed); @@ -244,22 +271,23 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix } // Possibly free some space using the new ss->tm_count instead of (old) ss->tm_count*sym->order - ss->tm_sym_map = realloc(ss->tm_sym_map, sizeof(qpms_ss_tmi_t) * sym->order * ss->tm_count); + QPMS_CRASHING_REALLOC(ss->tm_sym_map, sizeof(qpms_ss_tmi_t) * sym->order * ss->tm_count); // tm could be realloc'd as well, but those are just pointers, not likely many. // allocate particle symmetry map - ss->p_sym_map = malloc(sizeof(qpms_ss_pi_t) * sym->order * sym->order * ss->p_count); + QPMS_CRASHING_MALLOC(ss->p_sym_map, sizeof(qpms_ss_pi_t) * sym->order * sym->order * ss->p_count); // allocate orbit type array (TODO realloc later if too long) ss->orbit_type_count = 0; - ss->orbit_types = calloc(ss->p_count, sizeof(qpms_ss_orbit_type_t)); + QPMS_CRASHING_CALLOC(ss->orbit_types, ss->p_count, sizeof(qpms_ss_orbit_type_t)); - ss->otspace_end = ss->otspace = malloc( // reallocate later + QPMS_CRASHING_MALLOC(ss->otspace, // reallocate later (sizeof(qpms_ss_orbit_pi_t) * sym->order * sym->order + sizeof(qpms_ss_tmi_t) * sym->order + 3*sizeof(size_t) * sym->nirreps + sizeof(complex double) * SQ(sym->order * ss->max_bspecn)) * ss->p_count ); + ss->otspace_end = ss->otspace; // Workspace for the orbit type determination qpms_ss_orbit_type_t ot_current; @@ -348,15 +376,15 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp } } // Possibly free some space using the new ss->p_count instead of (old) ss->p_count*sym->order - ss->p_sym_map = realloc(ss->p_sym_map, sizeof(qpms_ss_pi_t) * sym->order * ss->p_count); - ss->p = realloc(ss->p, sizeof(qpms_particle_tid_t) * ss->p_count); - ss->p_orbitinfo = realloc(ss->p_orbitinfo, sizeof(qpms_ss_particle_orbitinfo_t)*ss->p_count); + QPMS_CRASHING_REALLOC(ss->p_sym_map, sizeof(qpms_ss_pi_t) * sym->order * ss->p_count); + QPMS_CRASHING_REALLOC(ss->p, sizeof(qpms_particle_tid_t) * ss->p_count); + QPMS_CRASHING_REALLOC(ss->p_orbitinfo, sizeof(qpms_ss_particle_orbitinfo_t)*ss->p_count); ss->p_capacity = ss->p_count; { // Reallocate the orbit type data space and update the pointers if needed. size_t otspace_sz = ss->otspace_end - ss->otspace; char *old_otspace = ss->otspace; - ss->otspace = realloc(ss->otspace, otspace_sz); + QPMS_CRASHING_REALLOC(ss->otspace, otspace_sz); ptrdiff_t shift = ss->otspace - old_otspace; if(shift) { for (size_t oi = 0; oi < ss->orbit_type_count; ++oi) { @@ -373,15 +401,14 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp // Set ss->fecv_size and ss->fecv_pstarts ss->fecv_size = 0; - ss->fecv_pstarts = malloc(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) { 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->saecv_sizes = malloc(sizeof(size_t) * sym->nirreps); if (!ss->saecv_sizes) abort(); - ss->saecv_ot_offsets = malloc(sizeof(size_t) * sym->nirreps * ss->orbit_type_count); - if (!ss->saecv_ot_offsets) abort(); + QPMS_CRASHING_MALLOC(ss->saecv_sizes, sizeof(size_t) * sym->nirreps); + QPMS_CRASHING_MALLOC(ss->saecv_ot_offsets, sizeof(size_t) * sym->nirreps * ss->orbit_type_count); for(qpms_iri_t iri = 0; iri < sym->nirreps; ++iri) { ss->saecv_sizes[iri] = 0; for(qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) { @@ -442,7 +469,7 @@ static inline void check_norm_compat(const qpms_vswf_set_spec_t *s) case QPMS_NORMALISATION_NORM_SPHARM: break; default: - abort(); // Only SPHARM and POWER norms are supported right now. + QPMS_WTF; // Only SPHARM and POWER norms are supported right now. } } @@ -455,9 +482,8 @@ complex double *qpms_orbit_action_matrix(complex double *target, // check_norm_compat(bspec); not needed here, the qpms_irot3_uvswfi_dense should complain if necessary const size_t n = bspec->n; const qpms_gmi_t N = ot->size; - if (target == NULL) - target = malloc(n*n*N*N*sizeof(complex double)); - if (target == NULL) abort(); + if (target == NULL) + QPMS_CRASHING_MALLOC(target, n*n*N*N*sizeof(complex double)); memset(target, 0, n*n*N*N*sizeof(complex double)); complex double tmp[n][n]; // this is the 'per-particle' action qpms_irot3_uvswfi_dense(tmp[0], bspec, sym->rep3d[g]); @@ -500,8 +526,7 @@ complex double *qpms_orbit_irrep_projector_matrix(complex double *target, const size_t n = bspec->n; const qpms_gmi_t N = ot->size; if (target == NULL) - target = malloc(n*n*N*N*sizeof(complex double)); - if (target == NULL) abort(); + QPMS_CRASHING_MALLOC(target, n*n*N*N*sizeof(complex double)); memset(target, 0, n*n*N*N*sizeof(complex double)); // Workspace for the orbit group action matrices complex double *tmp = malloc(n*n*N*N*sizeof(complex double)); @@ -556,7 +581,6 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size, const bool newtarget = (target == NULL); if (newtarget) QPMS_CRASHING_MALLOC(target,n*n*N*N*sizeof(complex double)); - if (target == NULL) abort(); memset(target, 0, n*n*N*N*sizeof(complex double)); // Get the projector (also workspace for right sg. vect.) @@ -564,12 +588,10 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size, ot, bspec, sym, iri); if(!projector) abort(); // Workspace for the right singular vectors. - complex double *V_H = malloc(n*n*N*N*sizeof(complex double)); - if(!V_H) abort(); + complex double *V_H; QPMS_CRASHING_MALLOC(V_H, n*n*N*N*sizeof(complex double)); // THIS SHOULD NOT BE NECESSARY - complex double *U = malloc(n*n*N*N*sizeof(complex double)); - if(!U) abort(); - double *s = malloc(n*N*sizeof(double)); if(!s) abort(); + complex double *U; QPMS_CRASHING_MALLOC(U, n*n*N*N*sizeof(complex double)); + double *s; QPMS_CRASHING_MALLOC(s, n*N*sizeof(double)); int info = LAPACKE_zgesdd(LAPACK_ROW_MAJOR, 'A', // jobz; overwrite projector with left sg.vec. and write right into V_H @@ -646,8 +668,7 @@ complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_pac return target_packed; const size_t full_len = ss->fecv_size; if (target_packed == NULL) - target_packed = malloc(SQ(packedlen)*sizeof(complex double)); - if (target_packed == NULL) abort(); + QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); // Workspace for the intermediate matrix @@ -684,8 +705,7 @@ complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_f const size_t packedlen = ss->saecv_sizes[iri]; const size_t full_len = ss->fecv_size; if (target_full == NULL) - target_full = malloc(SQ(full_len)*sizeof(complex double)); - if (target_full == NULL) abort(); + QPMS_CRASHING_MALLOC(target_full, SQ(full_len)*sizeof(complex double)); if(!add) memset(target_full, 0, SQ(full_len)*sizeof(complex double)); if(!packedlen) return target_full; // Empty irrep, do nothing. @@ -724,13 +744,12 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed, return target_packed; const size_t full_len = ss->fecv_size; if (target_packed == NULL) - target_packed = malloc(SQ(packedlen)*sizeof(complex double)); - if (target_packed == NULL) abort(); + QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); // Workspace for the intermediate particle-orbit matrix result - complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn) - * ss->sym->order); if (!tmp) abort(); + complex double *tmp; + QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); const complex double one = 1, zero = 0; @@ -806,15 +825,14 @@ complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full, const size_t packedlen = ss->saecv_sizes[iri]; const size_t full_len = ss->fecv_size; if (target_full == NULL) - target_full = malloc(SQ(full_len)*sizeof(complex double)); - if (target_full == NULL) abort(); + QPMS_CRASHING_MALLOC(target_full, SQ(full_len)*sizeof(complex double)); if(!add) memset(target_full, 0, SQ(full_len)*sizeof(complex double)); if(!packedlen) return target_full; // Empty irrep, do nothing. // Workspace for the intermediate particle-orbit matrix result - complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn) - * ss->sym->order); if (!tmp) abort(); + complex double *tmp; + QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); const complex double one = 1, zero = 0; @@ -889,8 +907,7 @@ complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed, const size_t packedlen = ss->saecv_sizes[iri]; if (!packedlen) return target_packed; // Empty irrep if (target_packed == NULL) - target_packed = malloc(packedlen*sizeof(complex double)); - if (target_packed == NULL) abort(); + QPMS_CRASHING_MALLOC(target_packed, packedlen*sizeof(complex double)); memset(target_packed, 0, packedlen*sizeof(complex double)); const complex double one = 1; @@ -928,8 +945,7 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full, const qpms_iri_t iri, bool add) { const size_t full_len = ss->fecv_size; if (target_full == NULL) - target_full = malloc(full_len*sizeof(complex double)); - if (target_full == NULL) abort(); + QPMS_CRASHING_MALLOC(target_full, full_len*sizeof(complex double)); if (!add) memset(target_full, 0, full_len*sizeof(complex double)); const complex double one = 1; @@ -1075,8 +1091,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; if (target_packed == NULL) - target_packed = malloc(SQ(packedlen)*sizeof(complex double)); - if (target_packed == NULL) abort(); + QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); // some of the following workspaces are probably redundant; TODO optimize later. @@ -1185,8 +1200,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; if (target_packed == NULL) - target_packed = malloc(SQ(packedlen)*sizeof(complex double)); - if (target_packed == NULL) abort(); + QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); // some of the following workspaces are probably redundant; TODO optimize later. @@ -1576,8 +1590,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; if (target_packed == NULL) - target_packed = malloc(SQ(packedlen)*sizeof(complex double)); - if (target_packed == NULL) abort(); + QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); qpms_ss_pi_t opistartR = 0; @@ -1629,8 +1642,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; if (target_packed == NULL) - target_packed = malloc(SQ(packedlen)*sizeof(complex double)); - if (target_packed == NULL) abort(); + QPMS_CRASHING_MALLOC(target_packed,SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); qpms_ss_pi_t opistartR = 0; @@ -1698,8 +1710,9 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed( complex double *qpms_scatsys_apply_Tmatrices_full( complex double *target_full, const complex double *inc_full, - const qpms_scatsys_t *ss) { + const qpms_scatsys_at_omega_t *ssw) { QPMS_UNTESTED; + const qpms_scatsys_t *ss = ssw->ss; if (!target_full) QPMS_CRASHING_CALLOC(target_full, ss->fecv_size, sizeof(complex double)); for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { @@ -1707,7 +1720,7 @@ complex double *qpms_scatsys_apply_Tmatrices_full( 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 = ss->tm[ss->p[pi].tmatrix_id]; + const qpms_tmatrix_t *T = ssw->tm[ss->p[pi].tmatrix_id]; qpms_apply_tmatrix(ptarget, psrc, T); } return target_full; diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 9ab7177..5690bdd 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -132,14 +132,14 @@ struct qpms_trans_calculator; struct qpms_epsmu_generator_t; typedef struct qpms_scatsys_t { - struct qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. + struct qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. /// (Template) T-matrix functions in the system. /** The qpms_abstract_tmatrix_t objects (onto which this array member point) * are NOT owned by this and must be kept alive for the whole lifetime * of all qpms_scatsys_t objects that are built upon them. */ - qpms_abstract_tmatrix_t **tmg; + qpms_tmatrix_function_t *tmg; qpms_ss_tmgi_t tmg_count; ///< Number of all different original T-matrix generators in the system. /// All the different T-matrix functions in the system, including those derived from \a tmg elements by symmetries. @@ -186,6 +186,10 @@ typedef struct qpms_scatsys_t { struct qpms_trans_calculator *c; } 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) { + return ss->tmg[ss->tm[tmi].tmgi]->spec; +} typedef struct qpms_scatsys_at_omega_t { const qpms_scatsys_t *ss; ///< Parent scattering system. @@ -195,10 +199,8 @@ typedef struct qpms_scatsys_at_omega_t { */ qpms_tmatrix_t **tm; complex double omega; ///< Angular frequency - complex double eps; ///< Background medium relative electric permittivity - complex double mu; ///< Background medium relative magnetic permeability + qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency complex double wavenumber; ///< Background medium wavenumber - complex double refractive_index; ///< Background medium refractive index } qpms_scatsys_at_omega_t; void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *); @@ -208,13 +210,13 @@ static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t return ss->tm[ss->p[pi].tmatrix_id]->spec; } -#if 0 //TODO!!! /// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed. /** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances, * so keep them alive until scatsys is destroyed. * * The following fields must be filled in the "proto- scattering system" \a orig: - * + * * orig->medium – The pointer is copied to the new qpms_scatsys_t instance; + * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting * * orig->tmg – The pointers are copied to the new qpms_scatsys_t instance; * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting * qpms_scatsys_t instances are properly destroyed. The pointers from orig->tmg, however, are copied. @@ -225,9 +227,20 @@ static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t * * orig->tm_count * * orig->p * * orig->p_count + * + * The resulting qpms_scatsys_t is obtained by actually evaluating the T-matrices + * at the given frequency \a omega and where applicable, these are compared + * by their values with given tolerances. The T-matrix generators are expected + * to preserve the point group symmetries for all frequencies. + * + * \returns An instance \a sso of qpms_scatsys_omega_t. Note that \a sso->ss + * must be saved by the caller before destroying \a sso + * (with qpms_scatsys_at_omega_free(), and destroyed only afterwards with + * qpms_scatsys_free() when not needed anymore. */ -qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym); -#endif +qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym, + complex double omega, const qpms_tolerance_spec_t *tol); + /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load. void qpms_scatsys_free(qpms_scatsys_t *s); 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 09/29] 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. From 541af5a984be888b671cdaa3a7a83834ad58d7c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 9 Jan 2020 16:57:30 +0200 Subject: [PATCH 10/29] WIP Rewriting scatsystem. Former-commit-id: 17f0f48ab54b84c4701b17846f941dd0142eb668 --- qpms/qpms_types.h | 2 +- qpms/scatsystem.c | 115 ++++++++++++++++++++++++---------------------- qpms/scatsystem.h | 56 +++++++++++----------- qpms/tmatrices.c | 11 ++++- qpms/tmatrices.h | 11 ++--- 5 files changed, 104 insertions(+), 91 deletions(-) diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index c4c2bdf..1f708de 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -410,7 +410,7 @@ typedef struct qpms_epsmu_t { complex double mu; ///< Relative permeability. } qpms_epsmu_t; -struct qpms_tolerance_spec_t; // See tolerances.c +struct qpms_tolerance_spec_t; // See tolerances.h #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l)) #endif // QPMS_TYPES diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index ac9f1b8..2d41036 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -21,6 +21,7 @@ #include "tmatrices.h" #include #include "kahansum.h" +#include "tolerances.h" #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS #include "qpmsblas.h" @@ -143,7 +144,7 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu // Almost 200 lines. This whole thing deserves a rewrite! -qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, +qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym, complex double omega, const qpms_tolerance_spec_t *tol) { // TODO check data sanity @@ -174,17 +175,17 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, assert(!cart3_isclose(orig->p[i].pos, orig->p[j].pos, 0, QPMS_SCATSYS_LEN_RTOL * lenscale)); + // Allocate T-matrix, particle and particle orbit info arrays + qpms_scatsys_t *ss; + QPMS_CRASHING_MALLOC(ss, sizeof(qpms_scatsys_t)); + ss->lenscale = lenscale; + ss->sym = sym; + // 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))); memcpy(ss->tmg, orig->tmg, ss->tmg_count * sizeof(*(ss->tmg))); - // Allocate T-matrix, particle and particle orbit info arrays - qpms_scatsys_t *ss; - ss = QPMS_CRASHING_MALLOC(ss, sizeof(qpms_scatsys_t)); - ss->lenscale = lenscale; - ss->sym = sym; - ss->tm_capacity = sym->order * orig->tm_count; QPMS_CRASHING_MALLOC(ss->tm, ss->tm_capacity * sizeof(*(ss->tm))); @@ -226,11 +227,11 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, break; } 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); + 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, ss->tm[j].spec->lMax); ++(ss->tm_count); } else qpms_tmatrix_free(ti); @@ -1049,13 +1050,14 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( -complex double *qpms_scatsys_build_modeproblem_matrix_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_t *ss, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, ) { + const complex double k = ssw->wavenumber; + const qpms_scatsys_t *ss = ssw->ss; const size_t full_len = ss->fecv_size; if(!target) QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); @@ -1066,13 +1068,13 @@ complex double *qpms_scatsys_build_modeproblem_matrix_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 = ssw->tm[ss->p[piR].tmatrix_id]->spec; const cart3_t posR = ss->p[piR].pos; size_t fullvec_offsetC = 0; // dest particle T-matrix - const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m; + const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m; 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 = ssw->tm[ss->p[piC].tmatrix_id]->spec; 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, @@ -1102,10 +1104,12 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full( 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_t *ss, qpms_iri_t iri, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, + qpms_iri_t iri, ) { + const qpms_scatsys_t *ss = ssw->ss; + const complex double k = ssw->wavenumber; const size_t packedlen = ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; @@ -1136,7 +1140,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( const size_t packed_orbit_offsetR = ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiR] + osnR * otR->irbase_sizes[iri]; - const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[piR].tmatrix_id]->spec; // Orbit coeff vector's full size: const size_t orbit_fullsizeR = otR->size * otR->bspecn; const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n @@ -1146,7 +1150,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( const cart3_t posR = ss->p[piR].pos; if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep // dest particle T-matrix - const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m; + const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t; const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC; @@ -1156,7 +1160,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( 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 = ss->tm[ss->p[piC].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec; // Orbit coeff vector's full size: const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n @@ -1211,10 +1215,11 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. complex double *target_packed, - const qpms_scatsys_t *ss, qpms_iri_t iri, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ) { + const qpms_scatsys_t *ss = ssw->ss; + const complex double k = ssw->wavenumber; const size_t packedlen = ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; @@ -1248,7 +1253,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( 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 = ss->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; // 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: @@ -1264,7 +1269,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( assert(ss->p_orbitinfo[piR].osn == osnR); const cart3_t posR = ss->p[piR].pos; // dest particle T-matrix - const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m; + const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t; const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC; @@ -1274,7 +1279,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( 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 = ss->tm[ss->p[piC].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec; // Orbit coeff vector's full size: const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n @@ -1328,19 +1333,20 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( } struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{ - const qpms_scatsys_t *ss; + const qpms_scatsys_at_omega_t *ssw; qpms_ss_pi_t *opistartR_ptr; pthread_mutex_t *opistartR_mutex; qpms_iri_t iri; complex double *target_packed; - complex double k; }; static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg) { const struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg *a = arg; - const qpms_scatsys_t *ss = a->ss; + const qpms_scatsys_at_omega_t *ssw = a->ssw; + const complex double k = ssw->wavenumber; + const qpms_scatsys_t *ss = ssw->ss; const qpms_iri_t iri = a->iri; const size_t packedlen = ss->saecv_sizes[iri]; @@ -1380,7 +1386,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread 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 = ss->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; // 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: @@ -1396,7 +1402,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread assert(ss->p_orbitinfo[piR].osn == osnR); const cart3_t posR = ss->p[piR].pos; // dest particle T-matrix - const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m; + const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t; const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC; @@ -1406,7 +1412,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread 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 = ss->tm[ss->p[piC].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec; // Orbit coeff vector's full size: const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n @@ -1465,12 +1471,11 @@ 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_t *ss; + const qpms_scatsys_at_omega_t *ssw; 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; }; @@ -1478,7 +1483,9 @@ 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_t *ss = a->ss; + const qpms_scatsys_at_omega_t *ssw = a->ssw; + const qpms_scatsys_t *ss = ssw->ss; + const complex double k = ssw->wavenumber; const qpms_iri_t iri = a->iri; const size_t packedlen = ss->saecv_sizes[iri]; const qpms_bessel_t J = a->J; @@ -1517,7 +1524,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 = ss->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; // 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: @@ -1541,7 +1548,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 = ss->tm[ss->p[piC].tmatrix_id]->spec; + const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec; // Orbit coeff vector's full size: const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n @@ -1600,7 +1607,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( complex double *target_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, - complex double k, ///< Wave number to use in the translation matrix. + const complex double k; qpms_bessel_t J ) { @@ -1616,7 +1623,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 = {ss, &opistartR, &opistartR_mutex, iri, target_packed, k, J}; + arg = {ssw, &opistartR, &opistartR_mutex, iri, target_packed, J}; // FIXME THIS IS NOT PORTABLE: long nthreads; @@ -1653,11 +1660,10 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. complex double *target_packed, - const qpms_scatsys_t *ss, qpms_iri_t iri, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ) { - const size_t packedlen = ss->saecv_sizes[iri]; + const size_t packedlen = ssw->ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; if (target_packed == NULL) @@ -1668,7 +1674,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( pthread_mutex_t opistartR_mutex; QPMS_ENSURE_SUCCESS(pthread_mutex_init(&opistartR_mutex, NULL)); const struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg - arg = {ss, &opistartR, &opistartR_mutex, iri, target_packed, k}; + arg = {ssw, &opistartR, &opistartR_mutex, iri, target_packed}; // FIXME THIS IS NOT PORTABLE: long nthreads; @@ -1784,7 +1790,8 @@ void qpms_ss_LU_free(qpms_ss_LU lu) { } qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(complex double *mpmatrix_full, - int *target_piv, const qpms_scatsys_t *ss) { + 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"); if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, ss->fecv_size * sizeof(int)); QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, ss->fecv_size, ss->fecv_size, @@ -1792,23 +1799,23 @@ qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(complex double *mpmatr qpms_ss_LU lu; lu.a = mpmatrix_full; lu.ipiv = target_piv; - lu.ss = ss; + lu.ssw = ssw; lu.full = true; lu.iri = -1; return lu; } qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed, - int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri) { + 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 = ss->saecv_sizes[iri]; + size_t n = ssw->ss->saecv_sizes[iri]; if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, n * sizeof(int)); QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, mpmatrix_packed, n, target_piv)); qpms_ss_LU lu; lu.a = mpmatrix_packed; lu.ipiv = target_piv; - lu.ss = ss; + lu.ssw = ssw; lu.full = false; lu.iri = iri; return lu; @@ -1816,21 +1823,21 @@ qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(complex double qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU( complex double *target, int *target_piv, - const qpms_scatsys_t *ss, complex double k){ - target = qpms_scatsys_build_modeproblem_matrix_full(target, ss, k); - return qpms_scatsys_modeproblem_matrix_full_factorise(target, target_piv, ss); + 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); } qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( complex double *target, int *target_piv, - const qpms_scatsys_t *ss, qpms_iri_t iri, complex double k){ - target = qpms_scatsys_build_modeproblem_matrix_irrep_packed(target, ss, iri, k); - return qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ss, iri); + const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri){ + target = qpms_scatsys_build_modeproblem_matrix_irrep_packed(target, ssw, iri); + return qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ssw, iri); } complex double *qpms_scatsys_scatter_solve( complex double *f, const complex double *a_inc, qpms_ss_LU lu) { - const size_t n = lu.full ? lu.ss->fecv_size : lu.ss->saecv_sizes[lu.iri]; + const size_t n = lu.full ? lu.ssw->ss->fecv_size : lu.ssw->ss->saecv_sizes[lu.iri]; if (!f) QPMS_CRASHING_MALLOC(f, n * sizeof(complex double)); memcpy(f, a_inc, n*sizeof(complex double)); // It will be rewritten by zgetrs QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/, n /*n*/, 1 /*nrhs number of right hand sides*/, diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 5690bdd..cdf1883 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -12,6 +12,7 @@ #define QPMS_SCATSYSTEM_H #include "qpms_types.h" #include "vswf.h" +#include "tmatrices.h" #include /// Overrides the number of threads spawned by the paralellized functions. @@ -124,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. - 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; @@ -132,7 +133,7 @@ struct qpms_trans_calculator; struct qpms_epsmu_generator_t; typedef struct qpms_scatsys_t { - struct qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium. + struct qpms_epsmu_generator_t medium; ///< Optical properties of the background medium. /// (Template) T-matrix functions in the system. /** The qpms_abstract_tmatrix_t objects (onto which this array member point) @@ -187,8 +188,13 @@ 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) { - return ss->tmg[ss->tm[tmi].tmgi]->spec; +static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(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) { + return ss->tmg[ss->tm[ss->p[pi].tmatrix_id].tmgi].spec; } typedef struct qpms_scatsys_at_omega_t { @@ -205,11 +211,6 @@ typedef struct qpms_scatsys_at_omega_t { void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *); -/// Convenience function to access pi'th particle's bspec. -static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) { - return ss->tm[ss->p[pi].tmatrix_id]->spec; -} - /// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed. /** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances, * so keep them alive until scatsys is destroyed. @@ -239,13 +240,13 @@ static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t * qpms_scatsys_free() when not needed anymore. */ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym, - complex double omega, const qpms_tolerance_spec_t *tol); + complex double omega, const struct qpms_tolerance_spec_t *tol); /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load. void qpms_scatsys_free(qpms_scatsys_t *s); /// Destroys the result of qpms_scatsys_eval_at_omega() -void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega *so); +void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw); /// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis. /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient. @@ -329,37 +330,36 @@ 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_build_modeproblem_matrix_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_t *ss, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw ); /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form. complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, - const qpms_scatsys_t *ss, qpms_iri_t iri, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, + qpms_iri_t iri ); /// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, - const qpms_scatsys_t *ss, qpms_iri_t iri, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, + qpms_iri_t iri ); /// Alternative (serial reference) implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorder_serial( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, - const qpms_scatsys_t *ss, qpms_iri_t iri, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, + qpms_iri_t iri ); /// LU factorisation (LAPACKE_zgetrf) result holder. typedef struct qpms_ss_LU { - const qpms_scatsys_t *ss; + const qpms_scatsys_at_omega_t *ssw; bool full; ///< true if full matrix; false if irrep-packed. qpms_iri_t iri; ///< Irrep index if `full == false`. /// LU decomposition array. @@ -373,30 +373,30 @@ void qpms_ss_LU_free(qpms_ss_LU); qpms_ss_LU qpms_scatsys_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_t *ss, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw ); /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch. qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_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_t *ss, qpms_iri_t iri, - complex double k ///< Wave number to use in the translation matrix. + const qpms_scatsys_at_omega_t *ssw, + qpms_iri_t iri ); /// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents. qpms_ss_LU qpms_scatsys_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_t *ss + const qpms_scatsys_at_omega_t *ssw ); /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents. qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise( complex double *modeproblem_matrix_irrep_packed, ///< 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_t *ss, qpms_iri_t iri + const qpms_scatsys_at_omega_t *ssw, + qpms_iri_t iri ); /// Solves a (possibly partial, irrep-packed) scattering problem \f$ (I-TS)f = Ta_\mathrm{inc} \f$ using a pre-factorised \f$ (I-TS) \f$. @@ -488,7 +488,7 @@ complex double *qpms_scatsys_incident_field_vector_full( complex double *qpms_scatsys_apply_Tmatrices_full( complex double *target_full, /// Target vector array. If NULL, a new one is allocated. const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL. - const qpms_scatsys_t *ss + const qpms_scatsys_at_omega_t *ssw ); diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 42dc309..11ad54c 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -53,6 +53,15 @@ qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) { return t; } +qpms_tmatrix_t *qpms_tmatrix_init_from_generator( + const qpms_vswf_set_spec_t *bspec, + qpms_tmatrix_generator_t gen, + complex double omega) { + qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); + QPMS_ENSURE_SUCCESS(gen.function(t, omega, gen.params)); + return t; +} + qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) { qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec); size_t n = T->spec->n; @@ -1107,7 +1116,7 @@ void qpms_tmatrix_operation_compose_chain_init(qpms_tmatrix_operation_t *dest, 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; + dest->typ = QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN; struct qpms_tmatrix_operation_compose_chain * const o = &(dest->op.compose_chain); o->n = nops; diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index 98ccd50..82c537d 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -257,14 +257,11 @@ typedef struct qpms_tmatrix_generator_t { } qpms_tmatrix_generator_t; /// Initialises and evaluates a new T-matrix using a generator. -static inline qpms_tmatrix_t *qpms_tmatrix_init_from_generator( +qpms_tmatrix_t *qpms_tmatrix_init_from_generator( const qpms_vswf_set_spec_t *bspec, qpms_tmatrix_generator_t gen, - complex double omega) { - qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); - QPMS_ENSURE_SUCCESS(gen.function(t, omega, gen.params)); - return t; -} + complex double omega); + /// Implementation of qpms_matrix_generator_t that just copies a constant matrix. /** N.B. this does almost no checks at all, so use it only for t-matrices with @@ -534,7 +531,7 @@ typedef struct qpms_tmatrix_function_t { /// Convenience function to create a new T-matrix from qpms_tmatrix_function_t. // FIXME the name is not very intuitive. static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, complex double omega) { - return qpms_tmatrix_init_from_generator(f.spec, omega, f.gen); + return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega); } /// Specifies different kinds of operations done on T-matrices 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 11/29] 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. }; From c86b88108878670b082aa2b8f308ffbb75038fdb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 10 Jan 2020 12:10:08 +0200 Subject: [PATCH 12/29] Forgotten function renames Former-commit-id: 89fa50cd8cfbdcf055e57f54093464f1e028c5bc --- qpms/scatsystem.c | 24 ++++++++++++------------ qpms/scatsystem.h | 18 ++++++++---------- 2 files changed, 20 insertions(+), 22 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 8aa43f6..ef29643 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1101,7 +1101,7 @@ complex double *qpms_scatsys_at_omega_build_modeproblem_matrix_full( } // Serial reference implementation. -complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( +complex double *qpms_scatsysw_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, @@ -1212,7 +1212,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( return target_packed; } -complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( +complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( /// 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 @@ -1332,7 +1332,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( return target_packed; } -struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{ +struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{ const qpms_scatsys_at_omega_t *ssw; qpms_ss_pi_t *opistartR_ptr; pthread_mutex_t *opistartR_mutex; @@ -1340,9 +1340,9 @@ struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{ complex double *target_packed; }; -static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg) +static void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg) { - const struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg + const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg *a = arg; const qpms_scatsys_at_omega_t *ssw = a->ssw; const complex double k = ssw->wavenumber; @@ -1656,7 +1656,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( // Parallel implementation, now default -complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( +complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( /// 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 @@ -1672,7 +1672,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( qpms_ss_pi_t opistartR = 0; pthread_mutex_t opistartR_mutex; QPMS_ENSURE_SUCCESS(pthread_mutex_init(&opistartR_mutex, NULL)); - const struct qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg + const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg arg = {ssw, &opistartR, &opistartR_mutex, iri, target_packed}; // FIXME THIS IS NOT PORTABLE: @@ -1694,7 +1694,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( pthread_t thread_ids[nthreads]; for(long thi = 0; thi < nthreads; ++thi) QPMS_ENSURE_SUCCESS(pthread_create(thread_ids + thi, NULL, - qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread, + qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread, (void *) &arg)); for(long thi = 0; thi < nthreads; ++thi) { void *retval; @@ -1732,7 +1732,7 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed( #endif -complex double *qpms_scatsys_apply_Tmatrices_full( +complex double *qpms_scatsysw_apply_Tmatrices_full( complex double *target_full, const complex double *inc_full, const qpms_scatsys_at_omega_t *ssw) { QPMS_UNTESTED; @@ -1826,11 +1826,11 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, ssw); } -qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( +qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( complex double *target, int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri){ - target = qpms_scatsys_build_modeproblem_matrix_irrep_packed(target, ssw, iri); - return qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ssw, iri); + target = qpms_scatsysw_build_modeproblem_matrix_irrep_packed(target, ssw, iri); + return qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ssw, iri); } complex double *qpms_scatsys_scatter_solve( diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 7490b02..b701dc8 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -209,8 +209,6 @@ typedef struct qpms_scatsys_at_omega_t { complex double wavenumber; ///< Background medium wavenumber } qpms_scatsys_at_omega_t; -void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *); - /// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed. /** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances, * so keep them alive until scatsys is destroyed. @@ -336,21 +334,21 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full( const qpms_scatsys_at_omega_t *ssw ); /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form. -complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( +complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( /// 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, qpms_iri_t iri ); -/// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). -complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( +/// Alternative implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed(). +complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( /// 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, qpms_iri_t iri ); -/// Alternative (serial reference) implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). -complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorder_serial( +/// Alternative (serial reference) implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed(). +complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorder_serial( /// 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, @@ -377,7 +375,7 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( ); /// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch. -qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( +qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_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, @@ -392,7 +390,7 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise( ); /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents. -qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise( +qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise( complex double *modeproblem_matrix_irrep_packed, ///< 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, @@ -485,7 +483,7 @@ complex double *qpms_scatsys_incident_field_vector_full( ); /// Applies T-matrices onto an incident field vector in the full basis. -complex double *qpms_scatsys_apply_Tmatrices_full( +complex double *qpms_scatsysw_apply_Tmatrices_full( complex double *target_full, /// Target vector array. If NULL, a new one is allocated. const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL. const qpms_scatsys_at_omega_t *ssw From 6cf1f667deedde59136a0ad46d1a6ff4fb7337c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 10 Jan 2020 15:20:45 +0200 Subject: [PATCH 13/29] Default tolerance constant. Former-commit-id: 19951825d21f94455da0228430a674eec37e7653 --- qpms/tolerances.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/qpms/tolerances.h b/qpms/tolerances.h index e3259a9..48e0a13 100644 --- a/qpms/tolerances.h +++ b/qpms/tolerances.h @@ -2,9 +2,14 @@ #ifndef QPMS_TOLERANCES_H #define QPMS_TOLERANCES_H +// TODO DOC + typedef struct qpms_tolerance_spec_t { double atol; double rtol; } qpms_tolerance_spec_t; +/// A rather arbitrary default tolerance. +static const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT = {.atol = 1e-15, .rtol = 1e-8}; + #endif // QPMS_TOLERANCES_H From be8f55eb1f17c7af8f083691d2efcbcf8393f4b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 10 Jan 2020 15:45:42 +0200 Subject: [PATCH 14/29] Update qpms_cdefs.pyx Former-commit-id: 7e1126b1bb594ffb1a8e5e9ed7a91839c0833b8f --- qpms/qpms_cdefs.pxd | 81 +++++++++++++++++++++++++++++++++------------ 1 file changed, 60 insertions(+), 21 deletions(-) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 56ac441..bf083ad 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -102,6 +102,7 @@ cdef extern from "qpms_types.h": QPMS_VSWF_ELECTRIC QPMS_VSWF_MAGNETIC QPMS_VSWF_LONGITUDINAL + ctypedef int32_t qpms_ss_tmgi_t ctypedef int32_t qpms_ss_tmi_t ctypedef int32_t qpms_ss_pi_t ctypedef int qpms_gmi_t @@ -150,6 +151,11 @@ cdef extern from "qpms_error.h": qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types) qpms_dbgmsg_flags qpms_dbgmsg_disable(qpms_dbgmsg_flags types) +cdef extern from "qpms/tolerances.h": + struct qpms_tolerance_spec_t: + pass # TODO + qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT + # This is due to the fact that cython apparently cannot nest the unnamed struct/unions in an obvious way ctypedef union qpms_incfield_planewave_params_k: @@ -475,6 +481,23 @@ cdef extern from "tmatrices.h": double omega, cdouble epsilon_fg, cdouble epsilon_bg) qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(cdouble *target, cdouble omega, const qpms_tmatrix_generator_axialsym_param_t *param, qpms_normalisation_t norm, qpms_bessel_t J) + struct qpms_tmatrix_function_t: + const qpms_vswf_set_spec_t *spec + const qpms_tmatrix_generator_t *gen + ctypedef enum qpms_tmatrix_operation_kind_t: + QPMS_TMATRIX_OPERATION_NOOP + QPMS_TMATRIX_OPERATION_LRMATRIX + QPMS_TMATRIX_OPERATION_IROT3 + QPMS_TMATRIX_OPERATION_IROT3ARR + QPMS_TMATRIX_OPERATION_COMPOSE_SUM + QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN + QPMS_TMATRIX_OPERATION_SCMULZ + QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE + + struct qpms_tmatrix_operation_t: + qpms_tmatrix_operation_kind_t typ + pass # TODO add the op union later if needed + const qpms_tmatrix_operation_t qpms_tmatrix_operation_noop cdef extern from "pointgroups.h": bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls) @@ -499,8 +522,14 @@ cdef extern from "scatsystem.h": struct qpms_particle_tid_t: cart3_t pos qpms_ss_tmi_t tmatrix_id + struct qpms_ss_derived_tmatrix_t: + qpms_ss_tmgi_t tmgi + qpms_tmatrix_operation_t op struct qpms_scatsys_t: - qpms_tmatrix_t **tm + qpms_epsmu_generator_t medium + qpms_tmatrix_function_t *tmg + qpms_ss_tmgi_t tmg_count + qpms_ss_derived_tmatrix_t **tm qpms_ss_tmi_t tm_count qpms_particle_tid_t *p qpms_ss_pi_t p_count @@ -508,10 +537,18 @@ cdef extern from "scatsystem.h": size_t fecv_size size_t *saecv_sizes const qpms_finite_group_t *sym - qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) void qpms_scatsys_free(qpms_scatsys_t *s) qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI qpms_scatsys_t *qpms_scatsys_load(char *path) #NI + struct qpms_scatsys_at_omega_t: + const qpms_scatsys_t *ss + qpms_tmatrix_t **tm, + cdouble omega + qpms_epsmu_t medium + cdouble wavenumber + qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym, + cdouble omega, const qpms_tolerance_spec_t *tol) + void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw) cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed, const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full, @@ -520,41 +557,43 @@ cdef extern from "scatsystem.h": const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full, const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add) - cdouble *qpms_scatsys_build_modeproblem_matrix_full(cdouble *target, - const qpms_scatsys_t *ss, cdouble k) + cdouble *qpms_scatsysw_build_modeproblem_matrix_full(cdouble *target, + const qpms_scatsys_at_omega_t *ssw) cdouble *qpms_scatsys_build_translation_matrix_full(cdouble *target, const qpms_scatsys_t *ss, cdouble k) cdouble *qpms_scatsys_build_translation_matrix_e_full(cdouble *target, const qpms_scatsys_t *ss, cdouble k, qpms_bessel_t J) - cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed(cdouble *target, - const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil + cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(cdouble *target, + const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil cdouble *qpms_scatsys_build_translation_matrix_e_irrep_packed(cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k, qpms_bessel_t J) nogil - cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( - cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil - cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( - cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil + cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( + cdouble *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil + cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( + cdouble *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil cdouble *qpms_scatsys_incident_field_vector_full(cdouble *target_full, const qpms_scatsys_t *ss, qpms_incfield_t field_at_point, const void *args, bint add) - cdouble *qpms_scatsys_apply_Tmatrices_full(cdouble *target_full, const cdouble *inc_full, - const qpms_scatsys_t *ss) + cdouble *qpms_scatsysw_apply_Tmatrices_full(cdouble *target_full, const cdouble *inc_full, + const qpms_scatsys_at_omega_t *ssw) struct qpms_ss_LU: - const qpms_scatsys_t *ss + const qpms_scatsys_at_omega_t *ssw bint full qpms_iri_t iri cdouble *a int *ipiv void qpms_ss_LU_free(qpms_ss_LU lu) - qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(cdouble *target, - int *target_piv, const qpms_scatsys_t *ss, cdouble k) - qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(cdouble *target, - int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) - qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full, - int *target_piv, const qpms_scatsys_t *ss) - qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full, - int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri) + qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(cdouble *target, + int *target_piv, const qpms_scatsys_at_omega_t *ssw) + qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(cdouble *target, + int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) + qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full, + int *target_piv, const qpms_scatsys_at_omega_t *ssw) + qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full, + int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) cdouble *qpms_scatsys_scatter_solve(cdouble *target_f, const cdouble *a_inc, qpms_ss_LU ludata) + const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) + const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) cdef extern from "ewald.h": struct qpms_csf_result: From b708b74292aeceafa549dbf5c2535695abe33684 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 10 Jan 2020 17:11:55 +0200 Subject: [PATCH 15/29] WIP Updating the cython scatsystem etc. Former-commit-id: ead1919c099cb2a0953310953685df69b7e1cbfb --- qpms/cytmatrices.pxd | 11 +++++++++++ qpms/cytmatrices.pyx | 4 ---- qpms/qpms_c.pyx | 4 +++- 3 files changed, 14 insertions(+), 5 deletions(-) diff --git a/qpms/cytmatrices.pxd b/qpms/cytmatrices.pxd index 51641d9..183ec29 100644 --- a/qpms/cytmatrices.pxd +++ b/qpms/cytmatrices.pxd @@ -14,6 +14,17 @@ cdef class TMatrixInterpolator: cdef inline qpms_tmatrix_interpolator_t *rawpointer(self): return self.interp +cdef class TMatrixGenerator: + cdef qpms_tmatrix_generator_t g + cdef object holder + cdef inline qpms_tmatrix_generator_t raw(self): + return self.g + cdef inline qpms_tmatrix_generator_t *rawpointer(self): + return &(self.g) + +cdef class TMatrixGeneratorTransformed: + pass + cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py! cdef readonly np.ndarray m # Numpy array holding the matrix data cdef readonly BaseSpec spec # Here we hold the base spec for the correct reference counting; TODO check if it gets copied diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx index 29c8d7c..6ca8afc 100644 --- a/qpms/cytmatrices.pyx +++ b/qpms/cytmatrices.pyx @@ -231,10 +231,6 @@ cdef class __AxialSymParams: cdef class TMatrixGenerator: - cdef qpms_tmatrix_generator_t g - cdef object holder - cdef qpms_tmatrix_generator_t raw(self): - return self.g def __init__(self, what): if isinstance(what, __MieParams): self.holder = what diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index c78e764..7f27b17 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -315,8 +315,10 @@ cdef class ScatteringSystem: Wrapper over the C qpms_scatsys_t structure. ''' cdef list basespecs # Here we keep the references to occuring basespecs + cdef list tmgens # here we keep the references to occuring TMatrixGenerators #cdef list Tmatrices # Here we keep the references to occuring T-matrices cdef qpms_scatsys_t *s + cdef qpms_tmatrix_function_t *tmg # this will ultimately contain pointers to stuff in basespecs and tmgens. def __cinit__(self, particles, FinitePointGroup sym): '''TODO doc. @@ -341,7 +343,7 @@ cdef class ScatteringSystem: for tm in tmobjs: # create references to BaseSpec objects self.basespecs.append(tm.spec) try: - orig.tm = malloc(orig.tm_count * sizeof(orig.tm[0])) + orig.tm = malloc(orig.tm_count * sizeof(orig.tm[0])) if not orig.tm: raise MemoryError orig.p = malloc(orig.p_count * sizeof(orig.p[0])) if not orig.p: raise MemoryError From 8f90842b24967307d034b89dc93cf5c19a3077c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 14 Jan 2020 09:17:13 +0200 Subject: [PATCH 16/29] Minor docs update in scatsystem.h Former-commit-id: 78caf1608c4ade295d47c17c20fb5743396cc8e8 --- qpms/scatsystem.h | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index b701dc8..61d2218 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -232,10 +232,15 @@ typedef struct qpms_scatsys_at_omega_t { * by their values with given tolerances. The T-matrix generators are expected * to preserve the point group symmetries for all frequencies. * + * This particular function will likely change in the future. + * * \returns An instance \a sso of qpms_scatsys_omega_t. Note that \a sso->ss * must be saved by the caller before destroying \a sso * (with qpms_scatsys_at_omega_free(), and destroyed only afterwards with * qpms_scatsys_free() when not needed anymore. + * \a sso->ss can be reused for different frequency by a + * qpms_scatsys_at_omega() call. + * */ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym, complex double omega, const struct qpms_tolerance_spec_t *tol); @@ -243,9 +248,15 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load. void qpms_scatsys_free(qpms_scatsys_t *s); -/// Destroys the result of qpms_scatsys_eval_at_omega() +/// Destroys a qpms_scatsys_at_omega_t. +/** Used on results of qpms_scatsys_apply_symmetry() and qpms_scatsys_at_omega(). */ void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw); +/// Evaluates scattering system T-matrices at a given frequency. +/** Free the result using qpms_scatsys_at_omega_free() when done. */ +qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, + complex double omega); + /// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis. /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient. * From c86ff698275ece6e140e818cea10190eb1358d4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 14 Jan 2020 10:09:17 +0200 Subject: [PATCH 17/29] WIP cython scatsystem Former-commit-id: f869e671148d3a75bbf34afe24aad02fd0d32611 --- qpms/cytmatrices.pxd | 2 +- qpms/qpms_c.pyx | 12 +++++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/qpms/cytmatrices.pxd b/qpms/cytmatrices.pxd index 183ec29..6e462b3 100644 --- a/qpms/cytmatrices.pxd +++ b/qpms/cytmatrices.pxd @@ -1,5 +1,5 @@ cimport numpy as np -from .qpms_cdefs cimport qpms_tmatrix_t, cdouble, qpms_tmatrix_interpolator_t +from .qpms_cdefs cimport qpms_tmatrix_t, cdouble, qpms_tmatrix_interpolator_t, qpms_tmatrix_generator_t from .cybspec cimport BaseSpec cdef class TMatrixInterpolator: diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 7f27b17..ef5554c 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -320,7 +320,7 @@ cdef class ScatteringSystem: cdef qpms_scatsys_t *s cdef qpms_tmatrix_function_t *tmg # this will ultimately contain pointers to stuff in basespecs and tmgens. - def __cinit__(self, particles, FinitePointGroup sym): + def __cinit__(self, particles, FinitePointGroup sym, omega): '''TODO doc. Takes the particles (which have to be a sequence of instances of Particle), fills them together with their t-matrices to the "proto-qpms_scatsys_t" @@ -540,6 +540,16 @@ cdef class ScatteringSystem: def scatter_solver(self, double k, iri=None): return ScatteringMatrix(self, k, iri) +cdef class ScatteringSystemAtOmega: + ''' + Wrapper over the C qpms_scatsys_at_omega_t structure + that keeps the T-matrix and background data evaluated + at specific frequency. + ''' + cdef qpms_scatsys_at_omega_t ssw + + pass #TODO + cdef class ScatteringMatrix: ''' Wrapper over the C qpms_ss_LU structure that keeps the factorised mode problem matrix. From b578f305ac8f6d2a919347acea19ff23347c7ed9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 14 Jan 2020 19:09:26 +0200 Subject: [PATCH 18/29] Support for constant T-matrix generator in cython. Former-commit-id: 184f88b0acf453d09e67f03cf41db06b4b4015bb --- qpms/cytmatrices.pyx | 4 ++++ qpms/qpms_cdefs.pxd | 1 + 2 files changed, 5 insertions(+) diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx index 6ca8afc..6cdfd4b 100644 --- a/qpms/cytmatrices.pyx +++ b/qpms/cytmatrices.pyx @@ -240,6 +240,10 @@ cdef class TMatrixGenerator: self.holder = what self.g.function = qpms_tmatrix_generator_axialsym self.g.params = (<__AxialSymParams?>self.holder).rawpointer() + elif isinstance(what, CTMatrix): + self.holder = what + self.g.function = qpms_tmatrix_generator_constant + self.g.params = (self.holder).rawpointer() # TODO INTERPOLATOR else: raise TypeError("Can't construct TMatrixGenerator from that") diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index bf083ad..9adf354 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -414,6 +414,7 @@ cdef extern from "tmatrices.h": qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, cdouble omega, const void *params) qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, cdouble omega, const void *params) qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, cdouble omega, const void *params) + qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, cdouble omega, const void *params) struct qpms_tmatrix_generator_sphere_param_t: qpms_epsmu_generator_t outside qpms_epsmu_generator_t inside From 6d83e26aa74bb27aeeddaf04c92647d8b50c1d65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 14 Jan 2020 22:03:39 +0200 Subject: [PATCH 19/29] Cython wrapper over qpms_tmatrix_function_t Former-commit-id: 85b601b7b70bc664d0348619d46fef3bac98fd17 --- qpms/cytmatrices.pxd | 10 ++++++++++ qpms/cytmatrices.pyx | 22 ++++++++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/qpms/cytmatrices.pxd b/qpms/cytmatrices.pxd index 6e462b3..9fcc1e4 100644 --- a/qpms/cytmatrices.pxd +++ b/qpms/cytmatrices.pxd @@ -22,6 +22,16 @@ cdef class TMatrixGenerator: cdef inline qpms_tmatrix_generator_t *rawpointer(self): return &(self.g) +cdef class TMatrixFunction: + cdef readonly qpms_tmatrix_function_t f + cdef readonly TMatrixGenerator generator # reference holder + cdef readonly BaseSpec spec # reference holder + cdef inline qpms_tmatrix_function_t raw(self): + return self.f + cdef inline qpms_tmatrix_function_t *rawpointer(self): + return &self.f + + cdef class TMatrixGeneratorTransformed: pass diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx index 6cdfd4b..99ca7be 100644 --- a/qpms/cytmatrices.pyx +++ b/qpms/cytmatrices.pyx @@ -228,6 +228,28 @@ cdef class __AxialSymParams: qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_BESSEL_REGULAR) return arr +cdef class TMatrixFunction: + ''' + Wrapper over qpms_tmatrix_function_t. The main functional difference between this + and TMatrixGenerator is that this remembers a specific BaseSpec + and its __call__ method takes only one mandatory argument (in addition to self). + ''' + def __init__(self, TMatrixGenerator tmg, BaseSpec spec): + self.generator = tmg + self.spec = spec + self.f.spec = self.generator.rawpointer() + self.f.gen = self.spec.rawpointer() + + def __call__(self, cdouble omega, fill = None): + cdef CTMatrix tm + if fill is None: # make a new CTMatrix + tm = CTMatrix(self.spec, None) + else: # TODO check whether fill has the same bspec as self? + tm = fill + if self.g.function(tm.rawpointer(), omega, self.f.gen.params) != 0: + raise ValueError("Something went wrong") + else: + return tm cdef class TMatrixGenerator: From 355bc523258a55ebb0d247f78788627db3b61197 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 15 Jan 2020 03:18:21 +0200 Subject: [PATCH 20/29] Rewrite ScatteringSystem. Compiles, not tested. Former-commit-id: 513741a41cd9b65348a8e91c367cd105592a0d68 --- qpms/cytmatrices.pxd | 4 +- qpms/cytmatrices.pyx | 6 +- qpms/qpms_c.pyx | 259 ++++++++++++++++++++++++++++++------------- qpms/qpms_cdefs.pxd | 11 +- qpms/scatsystem.h | 8 +- 5 files changed, 201 insertions(+), 87 deletions(-) diff --git a/qpms/cytmatrices.pxd b/qpms/cytmatrices.pxd index 9fcc1e4..3d39d65 100644 --- a/qpms/cytmatrices.pxd +++ b/qpms/cytmatrices.pxd @@ -1,5 +1,5 @@ cimport numpy as np -from .qpms_cdefs cimport qpms_tmatrix_t, cdouble, qpms_tmatrix_interpolator_t, qpms_tmatrix_generator_t +from .qpms_cdefs cimport qpms_tmatrix_t, cdouble, qpms_tmatrix_interpolator_t, qpms_tmatrix_generator_t, qpms_tmatrix_function_t from .cybspec cimport BaseSpec cdef class TMatrixInterpolator: @@ -23,7 +23,7 @@ cdef class TMatrixGenerator: return &(self.g) cdef class TMatrixFunction: - cdef readonly qpms_tmatrix_function_t f + cdef qpms_tmatrix_function_t f cdef readonly TMatrixGenerator generator # reference holder cdef readonly BaseSpec spec # reference holder cdef inline qpms_tmatrix_function_t raw(self): diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx index 99ca7be..940cd9e 100644 --- a/qpms/cytmatrices.pyx +++ b/qpms/cytmatrices.pyx @@ -237,8 +237,8 @@ cdef class TMatrixFunction: def __init__(self, TMatrixGenerator tmg, BaseSpec spec): self.generator = tmg self.spec = spec - self.f.spec = self.generator.rawpointer() - self.f.gen = self.spec.rawpointer() + self.f.gen = self.generator.rawpointer() + self.f.spec = self.spec.rawpointer() def __call__(self, cdouble omega, fill = None): cdef CTMatrix tm @@ -246,7 +246,7 @@ cdef class TMatrixFunction: tm = CTMatrix(self.spec, None) else: # TODO check whether fill has the same bspec as self? tm = fill - if self.g.function(tm.rawpointer(), omega, self.f.gen.params) != 0: + if self.f.gen.function(tm.rawpointer(), omega, self.f.gen.params) != 0: raise ValueError("Something went wrong") else: return tm diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index ef5554c..6f99173 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -12,7 +12,7 @@ from .cyquaternions cimport IRot3, CQuat from .cybspec cimport BaseSpec from .cycommon cimport make_c_string from .cycommon import string_c2py, PointGroupClass -from .cytmatrices cimport CTMatrix +from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator from libc.stdlib cimport malloc, free, calloc import warnings @@ -256,17 +256,37 @@ cdef class Particle: Wrapper over the qpms_particle_t structure. ''' cdef qpms_particle_t p - cdef readonly CTMatrix t # We hold the reference to the T-matrix to ensure correct reference counting + cdef readonly TMatrixFunction f # Reference to ensure correct reference counting - def __cinit__(Particle self, pos, CTMatrix t): + + def __cinit__(Particle self, pos, t, bspec = None): + cdef TMatrixGenerator tgen + cdef BaseSpec spec if(len(pos)>=2 and len(pos) < 4): self.p.pos.x = pos[0] self.p.pos.y = pos[1] self.p.pos.z = pos[2] if len(pos)==3 else 0 else: raise ValueError("Position argument has to contain 3 or 2 cartesian coordinates") - self.t = t - self.p.tmatrix = self.t.rawpointer() + if isinstance(t, CTMatrix): + tgen = TMatrixGenerator(t) + elif isinstance(t, TMatrixGenerator): + tgen = t + else: raise TypeError('t must be either CTMatrix or TMatrixGenerator, was %s' % str(type(t))) + if bspec is not None: + spec = bspec + else: + if isinstance(tgen.holder, CTMatrix): + spec = (tgen.holder).spec + else: + raise ValueError("bspec argument must be specified separately for str(type(t))") + self.f = TMatrixFunction(tgen, spec) + self.p.tmg = self.f.rawpointer() + # TODO non-trivial transformations later; if modified, do not forget to update ScatteringSystem constructor + self.p.op = qpms_tmatrix_operation_noop + + def __dealloc__(self): + qpms_tmatrix_operation_clear(&self.p.op) cdef qpms_particle_t *rawpointer(Particle self): '''Pointer to the qpms_particle_p structure. @@ -310,58 +330,105 @@ cpdef void scatsystem_set_nthreads(long n): qpms_scatsystem_set_nthreads(n) return + cdef class ScatteringSystem: ''' Wrapper over the C qpms_scatsys_t structure. ''' - cdef list basespecs # Here we keep the references to occuring basespecs - cdef list tmgens # here we keep the references to occuring TMatrixGenerators + cdef list tmgobjs # here we keep the references to occuring TMatrixFunctions (and hence BaseSpecs and TMatrixGenerators) #cdef list Tmatrices # Here we keep the references to occuring T-matrices cdef qpms_scatsys_t *s - cdef qpms_tmatrix_function_t *tmg # this will ultimately contain pointers to stuff in basespecs and tmgens. - def __cinit__(self, particles, FinitePointGroup sym, omega): - '''TODO doc. - Takes the particles (which have to be a sequence of instances of Particle), - fills them together with their t-matrices to the "proto-qpms_scatsys_t" - orig and calls qpms_scatsys_apply_symmetry - (and then cleans orig) - ''' + def check_s(self): # cdef instead? + if self.s == NULL: + raise ValueError("ScatteringSystem's s-pointer not set. You must not use the default constructor; use the create() method instead") + #TODO is there a way to disable the constructor outside this module? + + @staticmethod # We don't have any "standard" constructor for this right now + def create(particles, FinitePointGroup sym, cdouble omega): # TODO tolerances + # These we are going to construct + cdef ScatteringSystem self + cdef _ScatteringSystemAtOmega pyssw + cdef qpms_scatsys_t orig # This should be automatically init'd to 0 (CHECKME) - cdef qpms_ss_pi_t p_count = len(particles) - cdef qpms_ss_tmi_t tm_count = 0 + cdef qpms_ss_pi_t pi, p_count = len(particles) + cdef qpms_ss_tmi_t tmi, tm_count = 0 + cdef qpms_ss_tmgi_t tmgi, tmg_count = 0 + + cdef qpms_scatsys_at_omega_t *ssw + cdef qpms_scatsys_t *ss + + cdef Particle p + + tmgindices = dict() + tmgobjs = list() tmindices = dict() - tmobjs = list() - self.basespecs=list() - for p in particles: # find and enumerate unique t-matrices - if id(p.t) not in tmindices: - tmindices[id(p.t)] = tm_count - tmobjs.append(p.t) + tmlist = list() + for p in particles: # find and enumerate unique t-matrix generators + if p.p.op.typ != QPMS_TMATRIX_OPERATION_NOOP: + raise NotImplementedError("currently, only no-op T-matrix operations are allowed in ScatteringSystem constructor") + tmg_key = id(p.f) + if tmg_key not in tmgindices: + tmgindices[tmg_key] = tmg_count + tmgobjs.append(p.f) # Save the references on BaseSpecs and TMatrixGenerators (via TMatrixFunctions) + tmg_count += 1 + # Following lines have to be adjusted when nontrivial operations allowed: + tm_derived_key = (tmg_key, None) # TODO unique representation of p.p.op instead of None + if tm_derived_key not in tmindices: + tmindices[tm_derived_key] = tm_count + tmlist.append(tm_derived_key) tm_count += 1 + + orig.tmg_count = tmg_count orig.tm_count = tm_count orig.p_count = p_count - for tm in tmobjs: # create references to BaseSpec objects - self.basespecs.append(tm.spec) try: - orig.tm = malloc(orig.tm_count * sizeof(orig.tm[0])) + orig.tmg = malloc(orig.tmg_count * sizeof(orig.tmg[0])) + if not orig.tmg: raise MemoryError + orig.tm = malloc(orig.tm_count * sizeof(orig.tm[0])) if not orig.tm: raise MemoryError orig.p = malloc(orig.p_count * sizeof(orig.p[0])) if not orig.p: raise MemoryError + for tmgi in range(orig.tmg_count): + orig.tmg[tmgi] = (tmgobjs[tmgi]).raw() for tmi in range(tm_count): - orig.tm[tmi] = ((tmobjs[tmi])).rawpointer() + tm_derived_key = tmlist[tmi] + tmgi = tmgindices[tmg_key[0]] + orig.tm[tmi].tmgi = tmgi + orig.tm[tmi].op = qpms_tmatrix_operation_noop # TODO adjust when notrivial operations allowed for pi in range(p_count): - orig.p[pi].pos = ((particles[pi])).cval().pos - orig.p[pi].tmatrix_id = tmindices[id(particles[pi].t)] - self.s = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer()) + p = particles[pi] + tmg_key = id(p.f) + tm_derived_key = (tmg_key, None) # TODO unique representation of p.p.op instead of None + orig.p[pi].pos = p.cval().pos + orig.p[pi].tmatrix_id = tmindices[tm_derived_key] + ssw = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer(), omega, &QPMS_TOLERANCE_DEFAULT) + ss = ssw[0].ss finally: + free(orig.tmg) free(orig.tm) free(orig.p) + self = ScatteringSystem() + self.s = ss + self.tmgobjs = tmgobjs + pyssw = _ScatteringSystemAtOmega() + pyssw.ssw = ssw + pyssw.ss_pyref = self + return self, pyssw + + def __call__(self, cdouble omega): + self.check_s() + cdef _ScatteringSystemAtOmega pyssw = _ScatteringSystemAtOmega() + pyssw.ssw = qpms_scatsys_at_omega(self.s, omega) + pyssw.ss_pyref = self def __dealloc__(self): - qpms_scatsys_free(self.s) + if(self.s): + qpms_scatsys_free(self.s) property particles_tmi: def __get__(self): + self.check_s() r = list() cdef qpms_ss_pi_t pi for pi in range(self.s[0].p_count): @@ -369,20 +436,27 @@ cdef class ScatteringSystem: return r property fecv_size: - def __get__(self): return self.s[0].fecv_size + def __get__(self): + self.check_s() + return self.s[0].fecv_size property saecv_sizes: def __get__(self): + self.check_s() return [self.s[0].saecv_sizes[i] for i in range(self.s[0].sym[0].nirreps)] property irrep_names: def __get__(self): + self.check_s() return [string_c2py(self.s[0].sym[0].irreps[iri].name) if (self.s[0].sym[0].irreps[iri].name) else None for iri in range(self.s[0].sym[0].nirreps)] property nirreps: - def __get__(self): return self.s[0].sym[0].nirreps + def __get__(self): + self.check_s() + return self.s[0].sym[0].nirreps def pack_vector(self, vect, iri): + self.check_s() if len(vect) != self.fecv_size: raise ValueError("Length of a full vector has to be %d, not %d" % (self.fecv_size, len(vect))) @@ -394,6 +468,7 @@ cdef class ScatteringSystem: qpms_scatsys_irrep_pack_vector(&target_view[0], &vect_view[0], self.s, iri) return target_np def unpack_vector(self, packed, iri): + self.check_s() if len(packed) != self.saecv_sizes[iri]: raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d" % (iri, self.saecv_sizes, len(packed))) @@ -406,6 +481,7 @@ cdef class ScatteringSystem: self.s, iri, 0) return target_np def pack_matrix(self, fullmatrix, iri): + self.check_s() cdef size_t flen = self.s[0].fecv_size cdef size_t rlen = self.saecv_sizes[iri] fullmatrix = np.array(fullmatrix, dtype=complex, copy=False, order='C') @@ -420,6 +496,7 @@ cdef class ScatteringSystem: self.s, iri) return target_np def unpack_matrix(self, packedmatrix, iri): + self.check_s() cdef size_t flen = self.s[0].fecv_size cdef size_t rlen = self.saecv_sizes[iri] packedmatrix = np.array(packedmatrix, dtype=complex, copy=False, order='C') @@ -434,29 +511,8 @@ cdef class ScatteringSystem: self.s, iri, 0) return target_np - def modeproblem_matrix_full(self, double k): - cdef size_t flen = self.s[0].fecv_size - cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( - (flen,flen),dtype=complex, order='C') - cdef cdouble[:,::1] target_view = target - qpms_scatsys_build_modeproblem_matrix_full(&target_view[0][0], self.s, k) - return target - - def modeproblem_matrix_packed(self, double k, qpms_iri_t iri, version='pR'): - cdef size_t rlen = self.saecv_sizes[iri] - cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( - (rlen,rlen),dtype=complex, order='C') - cdef cdouble[:,::1] target_view = target - if (version == 'R'): - qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(&target_view[0][0], self.s, iri, k) - elif (version == 'pR'): - with nogil: - qpms_scatsys_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.s, iri, k) - else: - qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(&target_view[0][0], self.s, iri, k) - return target - def translation_matrix_full(self, double k, J = QPMS_HANKEL_PLUS): + self.check_s() cdef size_t flen = self.s[0].fecv_size cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (flen,flen),dtype=complex, order='C') @@ -465,6 +521,7 @@ cdef class ScatteringSystem: return target def translation_matrix_packed(self, double k, qpms_iri_t iri, J = QPMS_HANKEL_PLUS): + self.check_s() cdef size_t rlen = self.saecv_sizes[iri] cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (rlen,rlen),dtype=complex, order='C') @@ -475,6 +532,7 @@ cdef class ScatteringSystem: property fullvec_psizes: def __get__(self): + self.check_s() cdef np.ndarray[int32_t, ndim=1] ar = np.empty((self.s[0].p_count,), dtype=np.int32) cdef int32_t[::1] ar_view = ar for pi in range(self.s[0].p_count): @@ -484,6 +542,7 @@ cdef class ScatteringSystem: property fullvec_poffsets: def __get__(self): + self.check_s() cdef np.ndarray[intptr_t, ndim=1] ar = np.empty((self.s[0].p_count,), dtype=np.intp) cdef intptr_t[::1] ar_view = ar cdef intptr_t offset = 0 @@ -494,6 +553,7 @@ cdef class ScatteringSystem: property positions: def __get__(self): + self.check_s() cdef np.ndarray[np.double_t, ndim=2] ar = np.empty((self.s[0].p_count, 3), dtype=float) cdef np.double_t[:,::1] ar_view = ar for pi in range(self.s[0].p_count): @@ -503,6 +563,7 @@ cdef class ScatteringSystem: return ar def planewave_full(self, k_cart, E_cart): + self.check_s() k_cart = np.array(k_cart) E_cart = np.array(E_cart) if k_cart.shape != (3,) or E_cart.shape != (3,): @@ -522,7 +583,27 @@ cdef class ScatteringSystem: self.s, qpms_incfield_planewave, &p, 0) return target_np +cdef class _ScatteringSystemAtOmega: + ''' + Wrapper over the C qpms_scatsys_at_omega_t structure + that keeps the T-matrix and background data evaluated + at specific frequency. + ''' + cdef qpms_scatsys_at_omega_t *ssw + cdef ScatteringSystem ss_pyref + + def check(self): # cdef instead? + if not self.ssw: + raise ValueError("_ScatteringSystemAtOmega's ssw-pointer not set. You must not use the default constructor; ScatteringSystem.create() instead") + self.ss_pyref.check_s() + #TODO is there a way to disable the constructor outside this module? + + def __dealloc__(self): + if (self.ssw): + qpms_scatsys_at_omega_free(self.ssw) + def apply_Tmatrices_full(self, a): + self.check() if len(a) != self.fecv_size: raise ValueError("Length of a full vector has to be %d, not %d" % (self.fecv_size, len(a))) @@ -531,41 +612,67 @@ cdef class ScatteringSystem: cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty( (self.fecv_size,), dtype=complex, order='C') cdef cdouble[::1] target_view = target_np - qpms_scatsys_apply_Tmatrices_full(&target_view[0], &a_view[0], self.s) + qpms_scatsysw_apply_Tmatrices_full(&target_view[0], &a_view[0], self.ssw) return target_np - cdef qpms_scatsys_t *rawpointer(self): - return self.s + cdef qpms_scatsys_at_omega_t *rawpointer(self): + return self.ssw def scatter_solver(self, double k, iri=None): - return ScatteringMatrix(self, k, iri) + self.check() + return ScatteringMatrix(self, iri) -cdef class ScatteringSystemAtOmega: - ''' - Wrapper over the C qpms_scatsys_at_omega_t structure - that keeps the T-matrix and background data evaluated - at specific frequency. - ''' - cdef qpms_scatsys_at_omega_t ssw + property fecv_size: + def __get__(self): return self.ss_pyref.fecv_size + property saecv_sizes: + def __get__(self): return self.ss_pyref.saecv_sizes + property irrep_names: + def __get__(self): return self.ss_pyref.irrep_names + property nirreps: + def __get__(self): return self.ss_pyref.nirreps + + def modeproblem_matrix_full(self): + self.check() + cdef size_t flen = self.s[0].fecv_size + cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( + (flen,flen),dtype=complex, order='C') + cdef cdouble[:,::1] target_view = target + qpms_scatsysw_build_modeproblem_matrix_full(&target_view[0][0], self.ssw) + return target + + def modeproblem_matrix_packed(self, qpms_iri_t iri, version='pR'): + self.check() + cdef size_t rlen = self.saecv_sizes[iri] + cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( + (rlen,rlen),dtype=complex, order='C') + cdef cdouble[:,::1] target_view = target + if (version == 'R'): + qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(&target_view[0][0], self.ssw, iri) + elif (version == 'pR'): + with nogil: + qpms_scatsysw_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.ssw, iri) + else: + qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(&target_view[0][0], self.ssw, iri) + return target - pass #TODO cdef class ScatteringMatrix: ''' Wrapper over the C qpms_ss_LU structure that keeps the factorised mode problem matrix. ''' - cdef ScatteringSystem ss # Here we keep the reference to the parent scattering system + cdef _ScatteringSystemAtOmega ssw # Here we keep the reference to the parent scattering system cdef qpms_ss_LU lu - def __cinit__(self, ScatteringSystem ss, double k, iri=None): - self.ss = ss + def __cinit__(self, _ScatteringSystemAtOmega ssw, iri=None): + ssw.check() + self.ssw = ssw # TODO? pre-allocate the matrix with numpy to make it transparent? if iri is None: - self.lu = qpms_scatsys_build_modeproblem_matrix_full_LU( - NULL, NULL, ss.rawpointer(), k) + self.lu = qpms_scatsysw_build_modeproblem_matrix_full_LU( + NULL, NULL, ssw.rawpointer()) else: - self.lu = qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( - NULL, NULL, ss.rawpointer(), iri, k) + self.lu = qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( + NULL, NULL, ssw.rawpointer(), iri) def __dealloc__(self): qpms_ss_LU_free(self.lu) @@ -578,13 +685,13 @@ cdef class ScatteringMatrix: cdef size_t vlen cdef qpms_iri_t iri = -1; if self.lu.full: - vlen = self.lu.ss[0].fecv_size + vlen = self.lu.ssw[0].ss[0].fecv_size if len(a_inc) != vlen: raise ValueError("Length of a full coefficient vector has to be %d, not %d" % (vlen, len(a_inc))) else: iri = self.lu.iri - vlen = self.lu.ss[0].saecv_sizes[iri] + vlen = self.lu.ssw[0].ss[0].saecv_sizes[iri] if len(a_inc) != vlen: raise ValueError("Length of a %d. irrep packed coefficient vector has to be %d, not %d" % (iri, vlen, len(a_inc))) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 9adf354..6df0435 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -151,10 +151,10 @@ cdef extern from "qpms_error.h": qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types) qpms_dbgmsg_flags qpms_dbgmsg_disable(qpms_dbgmsg_flags types) -cdef extern from "qpms/tolerances.h": +cdef extern from "tolerances.h": struct qpms_tolerance_spec_t: pass # TODO - qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT + const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT # This is due to the fact that cython apparently cannot nest the unnamed struct/unions in an obvious way @@ -499,6 +499,7 @@ cdef extern from "tmatrices.h": qpms_tmatrix_operation_kind_t typ pass # TODO add the op union later if needed const qpms_tmatrix_operation_t qpms_tmatrix_operation_noop + void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *) cdef extern from "pointgroups.h": bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls) @@ -519,7 +520,8 @@ cdef extern from "scatsystem.h": void qpms_scatsystem_set_nthreads(long n) struct qpms_particle_t: cart3_t pos - const qpms_tmatrix_t *tmatrix + const qpms_tmatrix_function_t *tmg + qpms_tmatrix_operation_t op struct qpms_particle_tid_t: cart3_t pos qpms_ss_tmi_t tmatrix_id @@ -530,7 +532,7 @@ cdef extern from "scatsystem.h": qpms_epsmu_generator_t medium qpms_tmatrix_function_t *tmg qpms_ss_tmgi_t tmg_count - qpms_ss_derived_tmatrix_t **tm + qpms_ss_derived_tmatrix_t *tm qpms_ss_tmi_t tm_count qpms_particle_tid_t *p qpms_ss_pi_t p_count @@ -549,6 +551,7 @@ cdef extern from "scatsystem.h": cdouble wavenumber qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym, cdouble omega, const qpms_tolerance_spec_t *tol) + qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, cdouble omega) void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw) cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed, const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 61d2218..db9f807 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -20,10 +20,13 @@ void qpms_scatsystem_set_nthreads(long n); /// A particle, defined by its T-matrix and position. +/** This is rather only an auxillary intermediate structure to ultimately + * build an qpms_scatsys_t instance */ typedef struct qpms_particle_t { // Does it make sense to ever use other than cartesian coords for this? cart3_t pos; ///< Particle position in cartesian coordinates. - const qpms_tmatrix_t *tmatrix; ///< T-matrix; not owned by qpms_particle_t. + const qpms_tmatrix_function_t *tmg; ///< T-matrix function; not owned by qpms_particle_t. + qpms_tmatrix_operation_t op; ///< T-matrix transformation operation w.r.t. \a tmg. } qpms_particle_t; struct qpms_finite_group_t; @@ -123,7 +126,7 @@ typedef struct qpms_ss_particle_orbitinfo { } qpms_ss_particle_orbitinfo_t; /// Auxillary type used in qpms_scatsys_t: A recepy to create another T-matrices by symmetry operations. -typedef struct qpms_ss_derived_tmatrix { +typedef struct qpms_ss_derived_tmatrix_t { 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. } qpms_ss_derived_tmatrix_t; @@ -216,6 +219,7 @@ typedef struct qpms_scatsys_at_omega_t { * The following fields must be filled in the "proto- scattering system" \a orig: * * orig->medium – The pointer is copied to the new qpms_scatsys_t instance; * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting + * qpms_scatsys_t instances are properly destroyed. * * orig->tmg – The pointers are copied to the new qpms_scatsys_t instance; * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting * qpms_scatsys_t instances are properly destroyed. The pointers from orig->tmg, however, are copied. From c445b83593bab4b6317ade014e260c9f32513f51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 16 Jan 2020 07:52:50 +0200 Subject: [PATCH 21/29] Implement previously forgotten functions. Former-commit-id: df5215ad0349eb75bd2c7bee60f1fac50d23eb99 --- qpms/scatsystem.c | 53 ++++++++++++++++++++++++++++++++++++++++++++++- qpms/tmatrices.c | 18 ++++++++++++++++ qpms/tmatrices.h | 11 ++++++++++ 3 files changed, 81 insertions(+), 1 deletion(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index ef29643..1af486f 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -476,7 +476,58 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) { free(ss); } +void qpms_scatsys_at_omega_refill(qpms_scatsys_at_omega_t *ssw, + complex double omega) { + const qpms_scatsys_t * const ss = ssw->ss; + ssw->omega = omega; + ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega); + ssw->wavenumber = qpms_wavenumber(omega, ssw->medium); + qpms_tmatrix_t **tmatrices_preop; + QPMS_CRASHING_CALLOC(tmatrices_preop, ss->tmg_count, sizeof(*tmatrices_preop)); + for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi) + tmatrices_preop[tmgi] = qpms_tmatrix_init_from_function(ss->tmg[tmgi], omega); + for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) + qpms_tmatrix_apply_operation_replace(ssw->tm[tmi], &ss->tm[tmi].op, + tmatrices_preop[ss->tm[tmi].tmgi]); + for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi) + qpms_tmatrix_free(tmatrices_preop[tmgi]); + free(tmatrices_preop); +} +qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, + complex double omega) { + // TODO + qpms_scatsys_at_omega_t *ssw; + QPMS_CRASHING_MALLOC(ssw, sizeof(*ssw)); + ssw->omega = omega; + ssw->ss = ss; + ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega); + ssw->wavenumber = qpms_wavenumber(omega, ssw->medium); + qpms_tmatrix_t **tmatrices_preop; + QPMS_CRASHING_CALLOC(tmatrices_preop, ss->tmg_count, sizeof(*tmatrices_preop)); + for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi) + tmatrices_preop[tmgi] = qpms_tmatrix_init_from_function(ss->tmg[tmgi], omega); + for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) { + ssw->tm[tmi] = qpms_tmatrix_apply_operation(&ss->tm[tmi].op, + tmatrices_preop[ss->tm[tmi].tmgi]); //<- main difference to .._refill() + QPMS_ENSURE(ssw->tm[tmi], + "Got NULL pointer from qpms_tmatrix_apply_operation"); + } + for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi) + qpms_tmatrix_free(tmatrices_preop[tmgi]); + free(tmatrices_preop); + return ssw; +} + +void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw) { + if (ssw) { + if(ssw->tm) + for(qpms_ss_tmi_t i = 0; i < ssw->ss->tm_count; ++i) + qpms_tmatrix_free(ssw->tm[i]); + free(ssw->tm); + } + free(ssw); +} // (copypasta from symmetries.c) // TODO at some point, maybe support also other norms. @@ -1050,7 +1101,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( -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 diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 96516ba..c60f146 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1140,7 +1140,25 @@ void qpms_tmatrix_operation_compose_chain_init(qpms_tmatrix_operation_t *dest, o->opmem_size = opmem_size; } +qpms_tmatrix_t *qpms_tmatrix_mv(qpms_tmatrix_t *dest, + const qpms_tmatrix_t *orig) { + QPMS_ENSURE(qpms_vswf_set_spec_isidentical(dest->spec, orig->spec), + "Basis specifications must be identical!"); + memcpy(dest->m, orig->m, SQ(orig->spec->n)*sizeof(complex double)); + return dest; +} +qpms_tmatrix_t *qpms_tmatrix_apply_operation_replace(qpms_tmatrix_t *dest, + const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig) { + //QPMS_ENSURE(qpms_vswf_set_spec_isidentical(dest->spec, orig->spec), + // "Basis specifications must be identical!"); + // TODO should I check also for dest->owns_m? + // FIXME this is highly inoptimal; the hierarchy should be such + // that _operation() and operation_inplace() call this, and not the + // other way around + qpms_tmatrix_mv(dest, orig); + return qpms_tmatrix_apply_operation_inplace(op, dest); +} 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 8b4d2e6..117e3f9 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -27,6 +27,13 @@ qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec); /// Copies a T-matrix, allocating new array for the T-matrix data. qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T); +/// Copies a T-matrix contents between two pre-allocated T-matrices. +/** orig->spec and dest->spec must be identical. + * + * \returns \a dest + */ +qpms_tmatrix_t *qpms_tmatrix_mv(qpms_tmatrix_t *dest, const qpms_tmatrix_t *orig); + /// Destroys a T-matrix. void qpms_tmatrix_free(qpms_tmatrix_t *t); @@ -623,6 +630,10 @@ qpms_tmatrix_t *qpms_tmatrix_apply_operation(const qpms_tmatrix_operation_t *op, /// Apply an operation on a T-matrix and replace it with the result. qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(const qpms_tmatrix_operation_t *op, qpms_tmatrix_t *orig); +/// Apply an operation on a T-matrix and replace another one it with the result. +qpms_tmatrix_t *qpms_tmatrix_apply_operation_replace(qpms_tmatrix_t *dest, + const qpms_tmatrix_operation_t *op, const qpms_tmatrix_t *orig); + /// (Recursively) deallocates internal data of qpms_tmatrix_operation_t. /** It does NOT clear the memory pointed to by op it self, but only * heap-allocated data of certain operations, if labeled as owned by it. From a16cea4dca892b5803774f86eb060a6676443acd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 16 Jan 2020 08:51:05 +0200 Subject: [PATCH 22/29] Set medium generators Former-commit-id: 4eceb154349587fdb363a3e46073f065aa70fe61 --- qpms/qpms_c.pyx | 10 +++++++--- qpms/scatsystem.c | 1 + qpms/scatsystem.h | 2 +- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 6f99173..ec88dbe 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -13,6 +13,7 @@ from .cybspec cimport BaseSpec from .cycommon cimport make_c_string from .cycommon import string_c2py, PointGroupClass from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator +from .cymaterials cimport EpsMuGenerator from libc.stdlib cimport malloc, free, calloc import warnings @@ -337,6 +338,7 @@ cdef class ScatteringSystem: ''' cdef list tmgobjs # here we keep the references to occuring TMatrixFunctions (and hence BaseSpecs and TMatrixGenerators) #cdef list Tmatrices # Here we keep the references to occuring T-matrices + cdef EpsMuGenerator medium_holder # Here we keep the reference to medium generator cdef qpms_scatsys_t *s def check_s(self): # cdef instead? @@ -345,7 +347,7 @@ cdef class ScatteringSystem: #TODO is there a way to disable the constructor outside this module? @staticmethod # We don't have any "standard" constructor for this right now - def create(particles, FinitePointGroup sym, cdouble omega): # TODO tolerances + def create(particles, medium, FinitePointGroup sym, cdouble omega): # TODO tolerances # These we are going to construct cdef ScatteringSystem self cdef _ScatteringSystemAtOmega pyssw @@ -378,7 +380,8 @@ cdef class ScatteringSystem: tmindices[tm_derived_key] = tm_count tmlist.append(tm_derived_key) tm_count += 1 - + cdef EpsMuGenerator mediumgen = EpsMuGenerator(medium) + orig.medium = mediumgen.g orig.tmg_count = tmg_count orig.tm_count = tm_count orig.p_count = p_count @@ -393,7 +396,7 @@ cdef class ScatteringSystem: orig.tmg[tmgi] = (tmgobjs[tmgi]).raw() for tmi in range(tm_count): tm_derived_key = tmlist[tmi] - tmgi = tmgindices[tmg_key[0]] + tmgi = tmgindices[tm_derived_key[0]] orig.tm[tmi].tmgi = tmgi orig.tm[tmi].op = qpms_tmatrix_operation_noop # TODO adjust when notrivial operations allowed for pi in range(p_count): @@ -409,6 +412,7 @@ cdef class ScatteringSystem: free(orig.tm) free(orig.p) self = ScatteringSystem() + self.medium_holder = mediumgen self.s = ss self.tmgobjs = tmgobjs pyssw = _ScatteringSystemAtOmega() diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 1af486f..805c5ea 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -180,6 +180,7 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, QPMS_CRASHING_MALLOC(ss, sizeof(qpms_scatsys_t)); ss->lenscale = lenscale; ss->sym = sym; + ss->medium = orig->medium; // Copy the qpms_tmatrix_fuction_t from orig ss->tmg_count = orig->tmg_count; diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index db9f807..721fb69 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -217,7 +217,7 @@ typedef struct qpms_scatsys_at_omega_t { * so keep them alive until scatsys is destroyed. * * The following fields must be filled in the "proto- scattering system" \a orig: - * * orig->medium – The pointer is copied to the new qpms_scatsys_t instance; + * * orig->medium – The pointers are copied to the new qpms_scatsys_t instance; * the target qpms_abstract_tmatrix_t objects must be kept alive before all the resulting * qpms_scatsys_t instances are properly destroyed. * * orig->tmg – The pointers are copied to the new qpms_scatsys_t instance; From 4674fa58441b388298466b7d5ffb1de225389ec6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 16 Jan 2020 10:01:06 +0200 Subject: [PATCH 23/29] fix qpms_tmatrix_copy and pointer for comparison Former-commit-id: 98d91011109689512ccda2f8aab593909d0555e0 --- qpms/scatsystem.c | 2 +- qpms/tmatrices.c | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 805c5ea..2d48fd7 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -224,7 +224,7 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, 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)) { + if (qpms_tmatrix_isclose(ti, 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 diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index c60f146..25146cc 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -66,7 +66,7 @@ qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) { qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec); size_t n = T->spec->n; for(size_t i = 0; i < n*n; ++i) - t->m = T->m; + t->m[i] = T->m[i]; return t; } From 80a9f8703afecc0f240f60bb2e7624717ffe6c5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 17 Jan 2020 16:06:51 +0200 Subject: [PATCH 24/29] Fix index mismatch Former-commit-id: b1077e37785539b055c9d182b35273e321c0eda8 --- qpms/scatsystem.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 2d48fd7..0265472 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -228,8 +228,8 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, 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); - ss->tm[j].tmgi = orig->tm[j].tmgi; // T-matrix functions are preserved. + qpms_tmatrix_operation_copy(&ss->tm[j].op, &orig->tm[i].op); + ss->tm[j].tmgi = orig->tm[i].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, ssw->tm[j]->spec->lMax); From 3da4ec6b071b01e4df6c9cd000da6d836891191b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 20 Jan 2020 12:22:45 +0200 Subject: [PATCH 25/29] Disable useless compiler warnings. (To be enabled again later during code cleanup.) Former-commit-id: b89a16fad5d13d2f39f550b3e5e2e9b991908821 --- qpms/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index d90239c..66036a8 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -29,7 +29,7 @@ target_link_libraries (qpms ) target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) -target_compile_options(qpms PRIVATE -Wall -Wno-return-type) +target_compile_options(qpms PRIVATE -Wall -Wno-return-type -Wno-unused-variable -Wno-unused-function -Wno-unused-but-set-variable -Wno-unused-label) target_compile_definitions(qpms PRIVATE LATTICESUMS32 QPMS_VECTORS_NICE_TRANSFORMATIONS EWALD_AUTO_CUTOFF ) From 937757cf4763ceca6dc2a6c7a4939a2e2bf8767c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 20 Jan 2020 12:24:10 +0200 Subject: [PATCH 26/29] Fix unitialised values, assertion. Former-commit-id: ad712e7b88996e636e78350dcd23cfaf611bf0ec --- qpms/scatsystem.c | 3 ++- qpms/tmatrices.c | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 0265472..de1b6d0 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -277,6 +277,7 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, qpms_tmatrix_free(transformed); } else { // MISS, save the matrix (also the "abstract" one) ssw->tm[ss->tm_count] = transformed; + ss->tm[ss->tm_count].tmgi = ss->tm[tmi].tmgi; 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.op.compose_chain); o->ops[0] = & ss->tm[tmj].op; // Let's just borrow this @@ -1091,7 +1092,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( } fullvec_offsetC += bspecC->n; } - assert(fullvec_offsetC = full_len); + assert(fullvec_offsetC == full_len); fullvec_offsetR += bspecR->n; } assert(fullvec_offsetR == full_len); diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 25146cc..7d21a5c 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1113,6 +1113,7 @@ void qpms_tmatrix_operation_copy(qpms_tmatrix_operation_t *dest, const qpms_tmat default: QPMS_WTF; } + dest->typ = src->typ; } void qpms_tmatrix_operation_compose_chain_init(qpms_tmatrix_operation_t *dest, From 71852aa0174a5ee1492b5b26b61df4b7aa5761ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 20 Jan 2020 16:29:10 +0200 Subject: [PATCH 27/29] Fix function name in header. Former-commit-id: d9171a27990855ba0bdce741929b445b9688b444 --- qpms/scatsystem.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 721fb69..49bc9b8 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -363,7 +363,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( qpms_iri_t iri ); /// Alternative (serial reference) implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed(). -complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorder_serial( +complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( /// 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, From a2a51d0de6f566782a424577ae972e8e988fb3d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 20 Jan 2020 17:30:36 +0200 Subject: [PATCH 28/29] Fix k-argument related bugs etc. Former-commit-id: 72c955f31bcd1bfd9cd714c5b19d038f9c7ec6e3 --- qpms/qpms_c.pyx | 4 ++-- qpms/scatsystem.c | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index ec88dbe..47d3b52 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -622,7 +622,7 @@ cdef class _ScatteringSystemAtOmega: cdef qpms_scatsys_at_omega_t *rawpointer(self): return self.ssw - def scatter_solver(self, double k, iri=None): + def scatter_solver(self, iri=None): self.check() return ScatteringMatrix(self, iri) @@ -637,7 +637,7 @@ cdef class _ScatteringSystemAtOmega: def modeproblem_matrix_full(self): self.check() - cdef size_t flen = self.s[0].fecv_size + cdef size_t flen = self.ss_pyref.s[0].fecv_size cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (flen,flen),dtype=complex, order='C') cdef cdouble[:,::1] target_view = target diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index de1b6d0..b6b1149 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1675,7 +1675,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 = {ss, &opistartR, &opistartR_mutex, iri, target_packed, J}; + arg = {ss, &opistartR, &opistartR_mutex, iri, target_packed, k, J}; // FIXME THIS IS NOT PORTABLE: long nthreads; From bc5a024e86e687f9405ca10d0f50767c8d2240c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 21 Jan 2020 11:52:45 +0200 Subject: [PATCH 29/29] Fix another stupid bugs Former-commit-id: 093b5d5f09ac6a6f9be35fb7e20f73b2ba48f1d6 --- qpms/qpms_c.pyx | 1 + qpms/scatsystem.c | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 47d3b52..5a6e203 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -425,6 +425,7 @@ cdef class ScatteringSystem: cdef _ScatteringSystemAtOmega pyssw = _ScatteringSystemAtOmega() pyssw.ssw = qpms_scatsys_at_omega(self.s, omega) pyssw.ss_pyref = self + return pyssw def __dealloc__(self): if(self.s): diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index b6b1149..febef46 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -280,7 +280,7 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, ss->tm[ss->tm_count].tmgi = ss->tm[tmi].tmgi; 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.op.compose_chain); - o->ops[0] = & ss->tm[tmj].op; // Let's just borrow this + o->ops[0] = & ss->tm[tmi].op; // Let's just borrow this o->ops_owned[0] = false; o->opmem[0].typ = QPMS_TMATRIX_OPERATION_LRMATRIX; o->opmem[0].op.lrmatrix.m = m; @@ -505,6 +505,7 @@ qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, ssw->ss = ss; ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega); ssw->wavenumber = qpms_wavenumber(omega, ssw->medium); + QPMS_CRASHING_CALLOC(ssw->tm, ss->tm_count, sizeof(*ssw->tm)); qpms_tmatrix_t **tmatrices_preop; QPMS_CRASHING_CALLOC(tmatrices_preop, ss->tmg_count, sizeof(*tmatrices_preop)); for (qpms_ss_tmgi_t tmgi = 0; tmgi < ss->tmg_count; ++tmgi)