WIP Rewriting scatsystem.
Former-commit-id: 17f0f48ab54b84c4701b17846f941dd0142eb668
This commit is contained in:
Marek Nečada
parent
5dd93235f0
commit
541af5a984
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -21,6 +21,7 @@
|
|||
#include "tmatrices.h"
|
||||
#include <pthread.h>
|
||||
#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*/,
|
||||
|
|
|
|||
|
|
@ -12,6 +12,7 @@
|
|||
#define QPMS_SCATSYSTEM_H
|
||||
#include "qpms_types.h"
|
||||
#include "vswf.h"
|
||||
#include "tmatrices.h"
|
||||
#include <stdbool.h>
|
||||
|
||||
/// 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
|
||||
);
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Reference in New Issue