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

Merge branch 'abstract_scatsystem' 

"Abstract" scattering system for the finite case.


Former-commit-id: 1be9cb6196f660beaca04e8bd998b225cca30e94
This commit is contained in:
Marek Nečada 2020-01-21 15:07:46 +02:00
parent 3b6fb71f2e bc5a024e86
commit 76171179e7
13 changed files with 1063 additions and 274 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
qpms/CMakeLists.txt
View File

@ -3,6 +3,9 @@ find_package(GSL 2.0 REQUIRED)
find_package(BLAS REQUIRED) find_package(BLAS REQUIRED)
find_package(LAPACK 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 #includes
set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS}) set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS})
@ -26,7 +29,7 @@ target_link_libraries (qpms
) )
target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) 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 target_compile_definitions(qpms PRIVATE LATTICESUMS32 QPMS_VECTORS_NICE_TRANSFORMATIONS
EWALD_AUTO_CUTOFF EWALD_AUTO_CUTOFF
) )

23
qpms/cytmatrices.pxd
View File

@ -1,5 +1,5 @@
cimport numpy as np 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, qpms_tmatrix_function_t
from .cybspec cimport BaseSpec from .cybspec cimport BaseSpec
cdef class TMatrixInterpolator: cdef class TMatrixInterpolator:
@ -14,6 +14,27 @@ cdef class TMatrixInterpolator:
cdef inline qpms_tmatrix_interpolator_t *rawpointer(self): cdef inline qpms_tmatrix_interpolator_t *rawpointer(self):
return self.interp 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 TMatrixFunction:
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):
return self.f
cdef inline qpms_tmatrix_function_t *rawpointer(self):
return &self.f
cdef class TMatrixGeneratorTransformed:
pass
cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py! 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 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 cdef readonly BaseSpec spec # Here we hold the base spec for the correct reference counting; TODO check if it gets copied

26
qpms/cytmatrices.pyx
View File

@ -228,13 +228,31 @@ cdef class __AxialSymParams:
qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_BESSEL_REGULAR) qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_BESSEL_REGULAR)
return arr 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.gen = self.generator.rawpointer()
self.f.spec = 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.f.gen.function(tm.rawpointer(), omega, self.f.gen.params) != 0:
raise ValueError("Something went wrong")
else:
return tm
cdef class TMatrixGenerator: 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): def __init__(self, what):
if isinstance(what, __MieParams): if isinstance(what, __MieParams):
self.holder = what self.holder = what

6
qpms/materials.h
View File

@ -24,6 +24,12 @@ typedef struct qpms_epsmu_generator_t {
const void *params; const void *params;
} qpms_epsmu_generator_t; } 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. /// Constant optical property "generator" for qpms_epsmu_generator_t.
qpms_epsmu_t qpms_epsmu_const_g(complex double omega, ///< Frequency ignored. qpms_epsmu_t qpms_epsmu_const_g(complex double omega, ///< Frequency ignored.
const void *epsmu ///< Points to the qpms_epsmu_t to be returned. const void *epsmu ///< Points to the qpms_epsmu_t to be returned.

258
qpms/qpms_c.pyx
View File

@ -12,7 +12,8 @@ from .cyquaternions cimport IRot3, CQuat
from .cybspec cimport BaseSpec from .cybspec cimport BaseSpec
from .cycommon cimport make_c_string from .cycommon cimport make_c_string
from .cycommon import string_c2py, PointGroupClass from .cycommon import string_c2py, PointGroupClass
from .cytmatrices cimport CTMatrix from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator
from .cymaterials cimport EpsMuGenerator
from libc.stdlib cimport malloc, free, calloc from libc.stdlib cimport malloc, free, calloc
import warnings import warnings
@ -256,17 +257,37 @@ cdef class Particle:
Wrapper over the qpms_particle_t structure. Wrapper over the qpms_particle_t structure.
''' '''
cdef qpms_particle_t p 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): if(len(pos)>=2 and len(pos) < 4):
self.p.pos.x = pos[0] self.p.pos.x = pos[0]
self.p.pos.y = pos[1] self.p.pos.y = pos[1]
self.p.pos.z = pos[2] if len(pos)==3 else 0 self.p.pos.z = pos[2] if len(pos)==3 else 0
else: else:
raise ValueError("Position argument has to contain 3 or 2 cartesian coordinates") raise ValueError("Position argument has to contain 3 or 2 cartesian coordinates")
self.t = t if isinstance(t, CTMatrix):
self.p.tmatrix = self.t.rawpointer() tgen = TMatrixGenerator(t)
elif isinstance(t, TMatrixGenerator):
tgen = <TMatrixGenerator>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 = (<CTMatrix>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): cdef qpms_particle_t *rawpointer(Particle self):
'''Pointer to the qpms_particle_p structure. '''Pointer to the qpms_particle_p structure.
@ -310,56 +331,109 @@ cpdef void scatsystem_set_nthreads(long n):
qpms_scatsystem_set_nthreads(n) qpms_scatsystem_set_nthreads(n)
return return
cdef class ScatteringSystem: cdef class ScatteringSystem:
''' '''
Wrapper over the C qpms_scatsys_t structure. Wrapper over the C qpms_scatsys_t structure.
''' '''
cdef list basespecs # Here we keep the references to occuring basespecs 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 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 cdef qpms_scatsys_t *s
def __cinit__(self, particles, FinitePointGroup sym): def check_s(self): # cdef instead?
'''TODO doc. if self.s == <qpms_scatsys_t *>NULL:
Takes the particles (which have to be a sequence of instances of Particle), raise ValueError("ScatteringSystem's s-pointer not set. You must not use the default constructor; use the create() method instead")
fills them together with their t-matrices to the "proto-qpms_scatsys_t" #TODO is there a way to disable the constructor outside this module?
orig and calls qpms_scatsys_apply_symmetry
(and then cleans orig) @staticmethod # We don't have any "standard" constructor for this right now
''' def create(particles, medium, 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_scatsys_t orig # This should be automatically init'd to 0 (CHECKME)
cdef qpms_ss_pi_t p_count = len(particles) cdef qpms_ss_pi_t pi, p_count = len(particles)
cdef qpms_ss_tmi_t tm_count = 0 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() tmindices = dict()
tmobjs = list() tmlist = list()
self.basespecs=list() for p in particles: # find and enumerate unique t-matrix generators
for p in particles: # find and enumerate unique t-matrices if p.p.op.typ != QPMS_TMATRIX_OPERATION_NOOP:
if id(p.t) not in tmindices: raise NotImplementedError("currently, only no-op T-matrix operations are allowed in ScatteringSystem constructor")
tmindices[id(p.t)] = tm_count tmg_key = id(p.f)
tmobjs.append(p.t) 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 tm_count += 1
cdef EpsMuGenerator mediumgen = EpsMuGenerator(medium)
orig.medium = mediumgen.g
orig.tmg_count = tmg_count
orig.tm_count = tm_count orig.tm_count = tm_count
orig.p_count = p_count orig.p_count = p_count
for tm in tmobjs: # create references to BaseSpec objects
self.basespecs.append(tm.spec)
try: try:
orig.tm = <qpms_tmatrix_t **>malloc(orig.tm_count * sizeof(orig.tm[0])) orig.tmg = <qpms_tmatrix_function_t *>malloc(orig.tmg_count * sizeof(orig.tmg[0]))
if not orig.tmg: raise MemoryError
orig.tm = <qpms_ss_derived_tmatrix_t *>malloc(orig.tm_count * sizeof(orig.tm[0]))
if not orig.tm: raise MemoryError if not orig.tm: raise MemoryError
orig.p = <qpms_particle_tid_t *>malloc(orig.p_count * sizeof(orig.p[0])) orig.p = <qpms_particle_tid_t *>malloc(orig.p_count * sizeof(orig.p[0]))
if not orig.p: raise MemoryError if not orig.p: raise MemoryError
for tmgi in range(orig.tmg_count):
orig.tmg[tmgi] = (<TMatrixFunction?>tmgobjs[tmgi]).raw()
for tmi in range(tm_count): for tmi in range(tm_count):
orig.tm[tmi] = (<CTMatrix?>(tmobjs[tmi])).rawpointer() tm_derived_key = tmlist[tmi]
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): for pi in range(p_count):
orig.p[pi].pos = (<Particle?>(particles[pi])).cval().pos p = particles[pi]
orig.p[pi].tmatrix_id = tmindices[id(particles[pi].t)] tmg_key = id(p.f)
self.s = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer()) 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: finally:
free(orig.tmg)
free(orig.tm) free(orig.tm)
free(orig.p) free(orig.p)
self = ScatteringSystem()
self.medium_holder = mediumgen
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
return pyssw
def __dealloc__(self): def __dealloc__(self):
qpms_scatsys_free(self.s) if(self.s):
qpms_scatsys_free(self.s)
property particles_tmi: property particles_tmi:
def __get__(self): def __get__(self):
self.check_s()
r = list() r = list()
cdef qpms_ss_pi_t pi cdef qpms_ss_pi_t pi
for pi in range(self.s[0].p_count): for pi in range(self.s[0].p_count):
@ -367,20 +441,27 @@ cdef class ScatteringSystem:
return r return r
property fecv_size: 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: property saecv_sizes:
def __get__(self): def __get__(self):
self.check_s()
return [self.s[0].saecv_sizes[i] return [self.s[0].saecv_sizes[i]
for i in range(self.s[0].sym[0].nirreps)] for i in range(self.s[0].sym[0].nirreps)]
property irrep_names: property irrep_names:
def __get__(self): def __get__(self):
self.check_s()
return [string_c2py(self.s[0].sym[0].irreps[iri].name) return [string_c2py(self.s[0].sym[0].irreps[iri].name)
if (self.s[0].sym[0].irreps[iri].name) else None if (self.s[0].sym[0].irreps[iri].name) else None
for iri in range(self.s[0].sym[0].nirreps)] for iri in range(self.s[0].sym[0].nirreps)]
property 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): def pack_vector(self, vect, iri):
self.check_s()
if len(vect) != self.fecv_size: if len(vect) != self.fecv_size:
raise ValueError("Length of a full vector has to be %d, not %d" raise ValueError("Length of a full vector has to be %d, not %d"
% (self.fecv_size, len(vect))) % (self.fecv_size, len(vect)))
@ -392,6 +473,7 @@ cdef class ScatteringSystem:
qpms_scatsys_irrep_pack_vector(&target_view[0], &vect_view[0], self.s, iri) qpms_scatsys_irrep_pack_vector(&target_view[0], &vect_view[0], self.s, iri)
return target_np return target_np
def unpack_vector(self, packed, iri): def unpack_vector(self, packed, iri):
self.check_s()
if len(packed) != self.saecv_sizes[iri]: if len(packed) != self.saecv_sizes[iri]:
raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d" raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d"
% (iri, self.saecv_sizes, len(packed))) % (iri, self.saecv_sizes, len(packed)))
@ -404,6 +486,7 @@ cdef class ScatteringSystem:
self.s, iri, 0) self.s, iri, 0)
return target_np return target_np
def pack_matrix(self, fullmatrix, iri): def pack_matrix(self, fullmatrix, iri):
self.check_s()
cdef size_t flen = self.s[0].fecv_size cdef size_t flen = self.s[0].fecv_size
cdef size_t rlen = self.saecv_sizes[iri] cdef size_t rlen = self.saecv_sizes[iri]
fullmatrix = np.array(fullmatrix, dtype=complex, copy=False, order='C') fullmatrix = np.array(fullmatrix, dtype=complex, copy=False, order='C')
@ -418,6 +501,7 @@ cdef class ScatteringSystem:
self.s, iri) self.s, iri)
return target_np return target_np
def unpack_matrix(self, packedmatrix, iri): def unpack_matrix(self, packedmatrix, iri):
self.check_s()
cdef size_t flen = self.s[0].fecv_size cdef size_t flen = self.s[0].fecv_size
cdef size_t rlen = self.saecv_sizes[iri] cdef size_t rlen = self.saecv_sizes[iri]
packedmatrix = np.array(packedmatrix, dtype=complex, copy=False, order='C') packedmatrix = np.array(packedmatrix, dtype=complex, copy=False, order='C')
@ -432,29 +516,8 @@ cdef class ScatteringSystem:
self.s, iri, 0) self.s, iri, 0)
return target_np 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): def translation_matrix_full(self, double k, J = QPMS_HANKEL_PLUS):
self.check_s()
cdef size_t flen = self.s[0].fecv_size cdef size_t flen = self.s[0].fecv_size
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
(flen,flen),dtype=complex, order='C') (flen,flen),dtype=complex, order='C')
@ -463,6 +526,7 @@ cdef class ScatteringSystem:
return target return target
def translation_matrix_packed(self, double k, qpms_iri_t iri, J = QPMS_HANKEL_PLUS): 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 size_t rlen = self.saecv_sizes[iri]
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
(rlen,rlen),dtype=complex, order='C') (rlen,rlen),dtype=complex, order='C')
@ -473,6 +537,7 @@ cdef class ScatteringSystem:
property fullvec_psizes: property fullvec_psizes:
def __get__(self): 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 np.ndarray[int32_t, ndim=1] ar = np.empty((self.s[0].p_count,), dtype=np.int32)
cdef int32_t[::1] ar_view = ar cdef int32_t[::1] ar_view = ar
for pi in range(self.s[0].p_count): for pi in range(self.s[0].p_count):
@ -482,6 +547,7 @@ cdef class ScatteringSystem:
property fullvec_poffsets: property fullvec_poffsets:
def __get__(self): 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 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[::1] ar_view = ar
cdef intptr_t offset = 0 cdef intptr_t offset = 0
@ -492,6 +558,7 @@ cdef class ScatteringSystem:
property positions: property positions:
def __get__(self): 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.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 cdef np.double_t[:,::1] ar_view = ar
for pi in range(self.s[0].p_count): for pi in range(self.s[0].p_count):
@ -501,6 +568,7 @@ cdef class ScatteringSystem:
return ar return ar
def planewave_full(self, k_cart, E_cart): def planewave_full(self, k_cart, E_cart):
self.check_s()
k_cart = np.array(k_cart) k_cart = np.array(k_cart)
E_cart = np.array(E_cart) E_cart = np.array(E_cart)
if k_cart.shape != (3,) or E_cart.shape != (3,): if k_cart.shape != (3,) or E_cart.shape != (3,):
@ -520,7 +588,27 @@ cdef class ScatteringSystem:
self.s, qpms_incfield_planewave, <void *>&p, 0) self.s, qpms_incfield_planewave, <void *>&p, 0)
return target_np 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): def apply_Tmatrices_full(self, a):
self.check()
if len(a) != self.fecv_size: if len(a) != self.fecv_size:
raise ValueError("Length of a full vector has to be %d, not %d" raise ValueError("Length of a full vector has to be %d, not %d"
% (self.fecv_size, len(a))) % (self.fecv_size, len(a)))
@ -529,31 +617,67 @@ cdef class ScatteringSystem:
cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty( cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
(self.fecv_size,), dtype=complex, order='C') (self.fecv_size,), dtype=complex, order='C')
cdef cdouble[::1] target_view = target_np 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 return target_np
cdef qpms_scatsys_t *rawpointer(self): cdef qpms_scatsys_at_omega_t *rawpointer(self):
return self.s return self.ssw
def scatter_solver(self, iri=None):
self.check()
return ScatteringMatrix(self, iri)
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.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
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
def scatter_solver(self, double k, iri=None):
return ScatteringMatrix(self, k, iri)
cdef class ScatteringMatrix: cdef class ScatteringMatrix:
''' '''
Wrapper over the C qpms_ss_LU structure that keeps the factorised mode problem matrix. 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 cdef qpms_ss_LU lu
def __cinit__(self, ScatteringSystem ss, double k, iri=None): def __cinit__(self, _ScatteringSystemAtOmega ssw, iri=None):
self.ss = ss ssw.check()
self.ssw = ssw
# TODO? pre-allocate the matrix with numpy to make it transparent? # TODO? pre-allocate the matrix with numpy to make it transparent?
if iri is None: if iri is None:
self.lu = qpms_scatsys_build_modeproblem_matrix_full_LU( self.lu = qpms_scatsysw_build_modeproblem_matrix_full_LU(
NULL, NULL, ss.rawpointer(), k) NULL, NULL, ssw.rawpointer())
else: else:
self.lu = qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU( self.lu = qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
NULL, NULL, ss.rawpointer(), iri, k) NULL, NULL, ssw.rawpointer(), iri)
def __dealloc__(self): def __dealloc__(self):
qpms_ss_LU_free(self.lu) qpms_ss_LU_free(self.lu)
@ -566,13 +690,13 @@ cdef class ScatteringMatrix:
cdef size_t vlen cdef size_t vlen
cdef qpms_iri_t iri = -1; cdef qpms_iri_t iri = -1;
if self.lu.full: 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: if len(a_inc) != vlen:
raise ValueError("Length of a full coefficient vector has to be %d, not %d" raise ValueError("Length of a full coefficient vector has to be %d, not %d"
% (vlen, len(a_inc))) % (vlen, len(a_inc)))
else: else:
iri = self.lu.iri 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: if len(a_inc) != vlen:
raise ValueError("Length of a %d. irrep packed coefficient vector has to be %d, not %d" raise ValueError("Length of a %d. irrep packed coefficient vector has to be %d, not %d"
% (iri, vlen, len(a_inc))) % (iri, vlen, len(a_inc)))

86
qpms/qpms_cdefs.pxd
View File

@ -102,6 +102,7 @@ cdef extern from "qpms_types.h":
QPMS_VSWF_ELECTRIC QPMS_VSWF_ELECTRIC
QPMS_VSWF_MAGNETIC QPMS_VSWF_MAGNETIC
QPMS_VSWF_LONGITUDINAL QPMS_VSWF_LONGITUDINAL
ctypedef int32_t qpms_ss_tmgi_t
ctypedef int32_t qpms_ss_tmi_t ctypedef int32_t qpms_ss_tmi_t
ctypedef int32_t qpms_ss_pi_t ctypedef int32_t qpms_ss_pi_t
ctypedef int qpms_gmi_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_enable(qpms_dbgmsg_flags types)
qpms_dbgmsg_flags qpms_dbgmsg_disable(qpms_dbgmsg_flags types) qpms_dbgmsg_flags qpms_dbgmsg_disable(qpms_dbgmsg_flags types)
cdef extern from "tolerances.h":
struct qpms_tolerance_spec_t:
pass # TODO
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 # 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: ctypedef union qpms_incfield_planewave_params_k:
@ -476,6 +482,24 @@ cdef extern from "tmatrices.h":
double omega, cdouble epsilon_fg, cdouble epsilon_bg) double omega, cdouble epsilon_fg, cdouble epsilon_bg)
qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(cdouble *target, cdouble omega, 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) 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
void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *)
cdef extern from "pointgroups.h": cdef extern from "pointgroups.h":
bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls) bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls)
@ -496,12 +520,19 @@ cdef extern from "scatsystem.h":
void qpms_scatsystem_set_nthreads(long n) void qpms_scatsystem_set_nthreads(long n)
struct qpms_particle_t: struct qpms_particle_t:
cart3_t pos cart3_t pos
const qpms_tmatrix_t *tmatrix const qpms_tmatrix_function_t *tmg
qpms_tmatrix_operation_t op
struct qpms_particle_tid_t: struct qpms_particle_tid_t:
cart3_t pos cart3_t pos
qpms_ss_tmi_t tmatrix_id 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: 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_ss_tmi_t tm_count
qpms_particle_tid_t *p qpms_particle_tid_t *p
qpms_ss_pi_t p_count qpms_ss_pi_t p_count
@ -509,10 +540,19 @@ cdef extern from "scatsystem.h":
size_t fecv_size size_t fecv_size
size_t *saecv_sizes size_t *saecv_sizes
const qpms_finite_group_t *sym 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) void qpms_scatsys_free(qpms_scatsys_t *s)
qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI
qpms_scatsys_t *qpms_scatsys_load(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)
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, cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed,
const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full, cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full,
@ -521,41 +561,43 @@ cdef extern from "scatsystem.h":
const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full, cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full,
const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add) const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
cdouble *qpms_scatsys_build_modeproblem_matrix_full(cdouble *target, cdouble *qpms_scatsysw_build_modeproblem_matrix_full(cdouble *target,
const qpms_scatsys_t *ss, cdouble k) const qpms_scatsys_at_omega_t *ssw)
cdouble *qpms_scatsys_build_translation_matrix_full(cdouble *target, cdouble *qpms_scatsys_build_translation_matrix_full(cdouble *target,
const qpms_scatsys_t *ss, cdouble k) const qpms_scatsys_t *ss, cdouble k)
cdouble *qpms_scatsys_build_translation_matrix_e_full(cdouble *target, cdouble *qpms_scatsys_build_translation_matrix_e_full(cdouble *target,
const qpms_scatsys_t *ss, cdouble k, qpms_bessel_t J) const qpms_scatsys_t *ss, cdouble k, qpms_bessel_t J)
cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed(cdouble *target, cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(cdouble *target,
const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil
cdouble *qpms_scatsys_build_translation_matrix_e_irrep_packed(cdouble *target, 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 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 *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil cdouble *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil
cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial( cdouble *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) nogil cdouble *target, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) nogil
cdouble *qpms_scatsys_incident_field_vector_full(cdouble *target_full, cdouble *qpms_scatsys_incident_field_vector_full(cdouble *target_full,
const qpms_scatsys_t *ss, qpms_incfield_t field_at_point, const qpms_scatsys_t *ss, qpms_incfield_t field_at_point,
const void *args, bint add) const void *args, bint add)
cdouble *qpms_scatsys_apply_Tmatrices_full(cdouble *target_full, const cdouble *inc_full, cdouble *qpms_scatsysw_apply_Tmatrices_full(cdouble *target_full, const cdouble *inc_full,
const qpms_scatsys_t *ss) const qpms_scatsys_at_omega_t *ssw)
struct qpms_ss_LU: struct qpms_ss_LU:
const qpms_scatsys_t *ss const qpms_scatsys_at_omega_t *ssw
bint full bint full
qpms_iri_t iri qpms_iri_t iri
cdouble *a cdouble *a
int *ipiv int *ipiv
void qpms_ss_LU_free(qpms_ss_LU lu) void qpms_ss_LU_free(qpms_ss_LU lu)
qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_full_LU(cdouble *target, qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(cdouble *target,
int *target_piv, const qpms_scatsys_t *ss, cdouble k) int *target_piv, const qpms_scatsys_at_omega_t *ssw)
qpms_ss_LU qpms_scatsys_build_modeproblem_matrix_irrep_packed_LU(cdouble *target, qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(cdouble *target,
int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri, cdouble k) int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri)
qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full, qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(cdouble *modeproblem_matrix_full,
int *target_piv, const qpms_scatsys_t *ss) int *target_piv, const qpms_scatsys_at_omega_t *ssw)
qpms_ss_LU qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full, qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(cdouble *modeproblem_matrix_full,
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)
cdouble *qpms_scatsys_scatter_solve(cdouble *target_f, const cdouble *a_inc, qpms_ss_LU ludata) 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": cdef extern from "ewald.h":
struct qpms_csf_result: struct qpms_csf_result:

9
qpms/qpms_error.h
View File

@ -56,6 +56,15 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
(size_t) (size));\ (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) {\ #define QPMS_CRASHING_CALLOC(pointer, nmemb, size) {\
(pointer) = calloc((nmemb), (size));\ (pointer) = calloc((nmemb), (size));\
if(QPMS_UNLIKELY(!(pointer) && (nmemb) && (size)))\ if(QPMS_UNLIKELY(!(pointer) && (nmemb) && (size)))\

4
qpms/qpms_types.h
View File

@ -297,6 +297,9 @@ typedef struct qpms_vswf_set_spec_t {
/// T-matrix index used in qpms_scatsys_t and related structures. /// T-matrix index used in qpms_scatsys_t and related structures.
typedef int32_t qpms_ss_tmi_t; 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. /// Particle index used in qpms_scatsys_t and related structures.
typedef int32_t qpms_ss_pi_t; typedef int32_t qpms_ss_pi_t;
@ -407,6 +410,7 @@ typedef struct qpms_epsmu_t {
complex double mu; ///< Relative permeability. complex double mu; ///< Relative permeability.
} qpms_epsmu_t; } qpms_epsmu_t;
struct qpms_tolerance_spec_t; // See tolerances.h
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l)) #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
#endif // QPMS_TYPES #endif // QPMS_TYPES

361
qpms/scatsystem.c
View File

@ -21,6 +21,7 @@
#include "tmatrices.h" #include "tmatrices.h"
#include <pthread.h> #include <pthread.h>
#include "kahansum.h" #include "kahansum.h"
#include "tolerances.h"
#ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
#include "qpmsblas.h" #include "qpmsblas.h"
@ -80,7 +81,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]); qpms_ss_orbit_type_t * const ot_new = & (ss->orbit_types[ss->orbit_type_count]);
ot_new->size = ot_current->size; 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; const size_t bspecn = bspec->n;
ot_new->bspecn = bspecn; ot_new->bspecn = bspecn;
@ -143,7 +144,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! // 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_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 // TODO check data sanity
qpms_l_t lMax = 0; // the overall lMax of all base specs. qpms_l_t lMax = 0; // the overall lMax of all base specs.
@ -171,49 +174,82 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
for (qpms_ss_pi_t j = 0; j < i; ++j) 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)); 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 // Allocate T-matrix, particle and particle orbit info arrays
qpms_scatsys_t *ss = malloc(sizeof(qpms_scatsys_t)); qpms_scatsys_t *ss;
QPMS_CRASHING_MALLOC(ss, sizeof(qpms_scatsys_t));
ss->lenscale = lenscale; ss->lenscale = lenscale;
ss->sym = sym; ss->sym = sym;
ss->medium = orig->medium;
// 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)));
ss->tm_capacity = sym->order * orig->tm_count; 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_capacity = sym->order * orig->p_count;
ss->p = malloc(ss->p_capacity * sizeof(qpms_particle_tid_t)); QPMS_CRASHING_MALLOC(ss->p, 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_orbitinfo, ss->p_capacity * sizeof(qpms_ss_particle_orbitinfo_t));
for (qpms_ss_pi_t pi = 0; pi < ss->p_capacity; ++pi) { 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].t = QPMS_SS_P_ORBITINFO_UNDEF;
ss->p_orbitinfo[pi].p = QPMS_SS_P_ORBITINFO_UNDEF; ss->p_orbitinfo[pi].p = QPMS_SS_P_ORBITINFO_UNDEF;
} }
// Copy T-matrices; checking for duplicities // 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));
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
// Evaluate T-matrices at omega; checking for duplicities
ss->max_bspecn = 0; // We'll need it later.for memory alloc estimates. ss->max_bspecn = 0; // We'll need it later.for memory alloc estimates.
qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities; VLA!
ss->tm_count = 0; ss->tm_count = 0;
for (qpms_ss_tmi_t i = 0; i < orig->tm_count; ++i) { for (qpms_ss_tmi_t i = 0; i < orig->tm_count; ++i) {
qpms_tmatrix_t *ti = qpms_tmatrix_apply_operation(&(orig->tm[i].op), tm_orig_omega[orig->tm[i].tmgi]);
qpms_ss_tmi_t j; qpms_ss_tmi_t j;
for (j = 0; j < ss->tm_count; ++j) for (j = 0; j < ss->tm_count; ++j)
if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) { if (qpms_tmatrix_isclose(ti, ssw->tm[j], tol->rtol, tol->atol)) {
break; break;
} }
if (j == ss->tm_count) { // duplicity not found, copy the t-matrix if (j == ss->tm_count) { // duplicity not found, save both the "abstract" and "at omega" T-matrices
ss->tm[j] = qpms_tmatrix_copy(orig->tm[i]); qpms_tmatrix_operation_copy(&ss->tm[j].op, &orig->tm[i].op);
ss->max_bspecn = MAX(ss->tm[j]->spec->n, ss->max_bspecn); ss->tm[j].tmgi = orig->tm[i].tmgi; // T-matrix functions are preserved.
lMax = MAX(lMax, ss->tm[j]->spec->lMax); ssw->tm[j] = ti;
ss->max_bspecn = MAX(ssw->tm[j]->spec->n, ss->max_bspecn);
lMax = MAX(lMax, ssw->tm[j]->spec->lMax);
++(ss->tm_count); ++(ss->tm_count);
} }
else qpms_tmatrix_free(ti);
tm_dupl_remap[i] = j; tm_dupl_remap[i] = j;
if (normalisation == QPMS_NORMALISATION_UNDEF) if (normalisation == QPMS_NORMALISATION_UNDEF)
normalisation = ss->tm[i]->spec->norm; normalisation = ssw->tm[i]->spec->norm;
// We expect all bspec norms to be the same. // We expect all bspec norms to be the same.
else QPMS_ENSURE(normalisation == ss->tm[j]->spec->norm, else QPMS_ENSURE(normalisation == ssw->tm[j]->spec->norm,
"Normalisation convention must be the same for all T-matrices." "Normalisation convention must be the same for all T-matrices."
" %d != %d\n", normalisation, ss->tm[j]->spec->norm); " %d != %d\n", normalisation, ssw->tm[j]->spec->norm);
} }
// Free the original T-matrices at omega
for(qpms_ss_tmgi_t i = 0; i < orig->tmg_count; ++i)
qpms_tmatrix_free(tm_orig_omega[i]);
free(tm_orig_omega);
// Copy particles, remapping the t-matrix indices // Copy particles, remapping the t-matrix indices
for (qpms_ss_pi_t i = 0; i < orig->p_count; ++(i)) { for (qpms_ss_pi_t i = 0; i < orig->p_count; ++(i)) {
ss->p[i] = orig->p[i]; ss->p[i] = orig->p[i];
@ -222,44 +258,57 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
ss->p_count = orig->p_count; ss->p_count = orig->p_count;
// allocate t-matrix symmetry map // 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 // Extend the T-matrices list by the symmetry operations
for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi) for (qpms_ss_tmi_t tmi = 0; tmi < ss->tm_count; ++tmi)
for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){ for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){
const size_t d = ss->tm[tmi]->spec->n; const size_t d = ssw->tm[tmi]->spec->n;
complex double M[d][d]; // transformation matrix complex double *m;
qpms_irot3_uvswfi_dense(M[0], ss->tm[tmi]->spec, sym->rep3d[gmi]); QPMS_CRASHING_MALLOC(m, d*d*sizeof(complex double)); // ownership passed to ss->tm[ss->tm_count].op
qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ss->tm[tmi], M[0]); qpms_irot3_uvswfi_dense(m, ssw->tm[tmi]->spec, sym->rep3d[gmi]);
qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ssw->tm[tmi], m);
qpms_ss_tmi_t tmj; qpms_ss_tmi_t tmj;
for (tmj = 0; tmj < ss->tm_count; ++tmj) for (tmj = 0; tmj < ss->tm_count; ++tmj)
if (qpms_tmatrix_isclose(transformed, ss->tm[tmj], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) if (qpms_tmatrix_isclose(transformed, ssw->tm[tmj], tol->rtol, tol->atol))
break; break;
if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists
//TODO some "rounding error cleanup" symmetrisation could be performed here?
qpms_tmatrix_free(transformed); qpms_tmatrix_free(transformed);
} else { // MISS, save the matrix and increment the count } else { // MISS, save the matrix (also the "abstract" one)
ss->tm[ss->tm_count] = transformed; ssw->tm[ss->tm_count] = transformed;
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[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;
o->opmem[0].op.lrmatrix.m_size = d * d;
o->ops[1] = o->opmem;
o->ops_owned[1] = true;
++(ss->tm_count); ++(ss->tm_count);
} }
ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix
} }
// Possibly free some space using the new ss->tm_count instead of (old) ss->tm_count*sym->order // 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. // tm could be realloc'd as well, but those are just pointers, not likely many.
// allocate particle symmetry map // 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) // allocate orbit type array (TODO realloc later if too long)
ss->orbit_type_count = 0; 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_orbit_pi_t) * sym->order * sym->order
+ sizeof(qpms_ss_tmi_t) * sym->order + sizeof(qpms_ss_tmi_t) * sym->order
+ 3*sizeof(size_t) * sym->nirreps + 3*sizeof(size_t) * sym->nirreps
+ sizeof(complex double) * SQ(sym->order * ss->max_bspecn)) * ss->p_count + sizeof(complex double) * SQ(sym->order * ss->max_bspecn)) * ss->p_count
); );
ss->otspace_end = ss->otspace;
// Workspace for the orbit type determination // Workspace for the orbit type determination
qpms_ss_orbit_type_t ot_current; qpms_ss_orbit_type_t ot_current;
@ -348,15 +397,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 // 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); QPMS_CRASHING_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); QPMS_CRASHING_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_orbitinfo, sizeof(qpms_ss_particle_orbitinfo_t)*ss->p_count);
ss->p_capacity = ss->p_count; ss->p_capacity = ss->p_count;
{ // Reallocate the orbit type data space and update the pointers if needed. { // Reallocate the orbit type data space and update the pointers if needed.
size_t otspace_sz = ss->otspace_end - ss->otspace; size_t otspace_sz = ss->otspace_end - ss->otspace;
char *old_otspace = 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; ptrdiff_t shift = ss->otspace - old_otspace;
if(shift) { if(shift) {
for (size_t oi = 0; oi < ss->orbit_type_count; ++oi) { for (size_t oi = 0; oi < ss->orbit_type_count; ++oi) {
@ -373,15 +422,14 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
// Set ss->fecv_size and ss->fecv_pstarts // Set ss->fecv_size and ss->fecv_pstarts
ss->fecv_size = 0; 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) { for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
ss->fecv_pstarts[pi] = ss->fecv_size; ss->fecv_pstarts[pi] = ss->fecv_size;
ss->fecv_size += ss->tm[ss->p[pi].tmatrix_id]->spec->n; // That's a lot of dereferencing! ss->fecv_size += ssw->tm[ss->p[pi].tmatrix_id]->spec->n; // That's a lot of dereferencing!
} }
ss->saecv_sizes = malloc(sizeof(size_t) * sym->nirreps); if (!ss->saecv_sizes) abort(); QPMS_CRASHING_MALLOC(ss->saecv_sizes, sizeof(size_t) * sym->nirreps);
ss->saecv_ot_offsets = malloc(sizeof(size_t) * sym->nirreps * ss->orbit_type_count); QPMS_CRASHING_MALLOC(ss->saecv_ot_offsets, sizeof(size_t) * sym->nirreps * ss->orbit_type_count);
if (!ss->saecv_ot_offsets) abort();
for(qpms_iri_t iri = 0; iri < sym->nirreps; ++iri) { for(qpms_iri_t iri = 0; iri < sym->nirreps; ++iri) {
ss->saecv_sizes[iri] = 0; ss->saecv_sizes[iri] = 0;
for(qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) { for(qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) {
@ -408,13 +456,14 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
} }
ss->c = qpms_trans_calculator_init(lMax, normalisation); ss->c = qpms_trans_calculator_init(lMax, normalisation);
return ss; return ssw;
} }
void qpms_scatsys_free(qpms_scatsys_t *ss) { void qpms_scatsys_free(qpms_scatsys_t *ss) {
if(ss) { if(ss) {
free(ss->tm); free(ss->tm);
free(ss->tmg);
free(ss->p); free(ss->p);
free(ss->fecv_pstarts); free(ss->fecv_pstarts);
free(ss->tm_sym_map); free(ss->tm_sym_map);
@ -429,7 +478,59 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) {
free(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_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)
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) // (copypasta from symmetries.c)
// TODO at some point, maybe support also other norms. // TODO at some point, maybe support also other norms.
@ -442,7 +543,7 @@ static inline void check_norm_compat(const qpms_vswf_set_spec_t *s)
case QPMS_NORMALISATION_NORM_SPHARM: case QPMS_NORMALISATION_NORM_SPHARM:
break; break;
default: default:
abort(); // Only SPHARM and POWER norms are supported right now. QPMS_WTF; // Only SPHARM and POWER norms are supported right now.
} }
} }
@ -455,9 +556,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 // check_norm_compat(bspec); not needed here, the qpms_irot3_uvswfi_dense should complain if necessary
const size_t n = bspec->n; const size_t n = bspec->n;
const qpms_gmi_t N = ot->size; const qpms_gmi_t N = ot->size;
if (target == NULL) if (target == NULL)
target = malloc(n*n*N*N*sizeof(complex double)); 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)); memset(target, 0, n*n*N*N*sizeof(complex double));
complex double tmp[n][n]; // this is the 'per-particle' action complex double tmp[n][n]; // this is the 'per-particle' action
qpms_irot3_uvswfi_dense(tmp[0], bspec, sym->rep3d[g]); qpms_irot3_uvswfi_dense(tmp[0], bspec, sym->rep3d[g]);
@ -500,8 +600,7 @@ complex double *qpms_orbit_irrep_projector_matrix(complex double *target,
const size_t n = bspec->n; const size_t n = bspec->n;
const qpms_gmi_t N = ot->size; const qpms_gmi_t N = ot->size;
if (target == NULL) if (target == NULL)
target = malloc(n*n*N*N*sizeof(complex double)); 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)); memset(target, 0, n*n*N*N*sizeof(complex double));
// Workspace for the orbit group action matrices // Workspace for the orbit group action matrices
complex double *tmp = malloc(n*n*N*N*sizeof(complex double)); complex double *tmp = malloc(n*n*N*N*sizeof(complex double));
@ -556,7 +655,6 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
const bool newtarget = (target == NULL); const bool newtarget = (target == NULL);
if (newtarget) if (newtarget)
QPMS_CRASHING_MALLOC(target,n*n*N*N*sizeof(complex double)); 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)); memset(target, 0, n*n*N*N*sizeof(complex double));
// Get the projector (also workspace for right sg. vect.) // Get the projector (also workspace for right sg. vect.)
@ -564,12 +662,10 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
ot, bspec, sym, iri); ot, bspec, sym, iri);
if(!projector) abort(); if(!projector) abort();
// Workspace for the right singular vectors. // Workspace for the right singular vectors.
complex double *V_H = malloc(n*n*N*N*sizeof(complex double)); complex double *V_H; QPMS_CRASHING_MALLOC(V_H, n*n*N*N*sizeof(complex double));
if(!V_H) abort();
// THIS SHOULD NOT BE NECESSARY // THIS SHOULD NOT BE NECESSARY
complex double *U = malloc(n*n*N*N*sizeof(complex double)); complex double *U; QPMS_CRASHING_MALLOC(U, n*n*N*N*sizeof(complex double));
if(!U) abort(); double *s; QPMS_CRASHING_MALLOC(s, n*N*sizeof(double));
double *s = malloc(n*N*sizeof(double)); if(!s) abort();
int info = LAPACKE_zgesdd(LAPACK_ROW_MAJOR, int info = LAPACKE_zgesdd(LAPACK_ROW_MAJOR,
'A', // jobz; overwrite projector with left sg.vec. and write right into V_H 'A', // jobz; overwrite projector with left sg.vec. and write right into V_H
@ -646,8 +742,7 @@ complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_pac
return target_packed; return target_packed;
const size_t full_len = ss->fecv_size; const size_t full_len = ss->fecv_size;
if (target_packed == NULL) if (target_packed == NULL)
target_packed = malloc(SQ(packedlen)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
if (target_packed == NULL) abort();
memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
// Workspace for the intermediate matrix // Workspace for the intermediate matrix
@ -684,8 +779,7 @@ complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_f
const size_t packedlen = ss->saecv_sizes[iri]; const size_t packedlen = ss->saecv_sizes[iri];
const size_t full_len = ss->fecv_size; const size_t full_len = ss->fecv_size;
if (target_full == NULL) if (target_full == NULL)
target_full = malloc(SQ(full_len)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_full, SQ(full_len)*sizeof(complex double));
if (target_full == NULL) abort();
if(!add) memset(target_full, 0, 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. if(!packedlen) return target_full; // Empty irrep, do nothing.
@ -724,13 +818,12 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
return target_packed; return target_packed;
const size_t full_len = ss->fecv_size; const size_t full_len = ss->fecv_size;
if (target_packed == NULL) if (target_packed == NULL)
target_packed = malloc(SQ(packedlen)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
if (target_packed == NULL) abort();
memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
// Workspace for the intermediate particle-orbit matrix result // Workspace for the intermediate particle-orbit matrix result
complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn) complex double *tmp;
* ss->sym->order); if (!tmp) abort(); QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order);
const complex double one = 1, zero = 0; const complex double one = 1, zero = 0;
@ -806,15 +899,14 @@ complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
const size_t packedlen = ss->saecv_sizes[iri]; const size_t packedlen = ss->saecv_sizes[iri];
const size_t full_len = ss->fecv_size; const size_t full_len = ss->fecv_size;
if (target_full == NULL) if (target_full == NULL)
target_full = malloc(SQ(full_len)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_full, SQ(full_len)*sizeof(complex double));
if (target_full == NULL) abort();
if(!add) memset(target_full, 0, 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. if(!packedlen) return target_full; // Empty irrep, do nothing.
// Workspace for the intermediate particle-orbit matrix result // Workspace for the intermediate particle-orbit matrix result
complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn) complex double *tmp;
* ss->sym->order); if (!tmp) abort(); QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order);
const complex double one = 1, zero = 0; const complex double one = 1, zero = 0;
@ -889,8 +981,7 @@ complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
const size_t packedlen = ss->saecv_sizes[iri]; const size_t packedlen = ss->saecv_sizes[iri];
if (!packedlen) return target_packed; // Empty irrep if (!packedlen) return target_packed; // Empty irrep
if (target_packed == NULL) if (target_packed == NULL)
target_packed = malloc(packedlen*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_packed, packedlen*sizeof(complex double));
if (target_packed == NULL) abort();
memset(target_packed, 0, packedlen*sizeof(complex double)); memset(target_packed, 0, packedlen*sizeof(complex double));
const complex double one = 1; const complex double one = 1;
@ -928,8 +1019,7 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
const qpms_iri_t iri, bool add) { const qpms_iri_t iri, bool add) {
const size_t full_len = ss->fecv_size; const size_t full_len = ss->fecv_size;
if (target_full == NULL) if (target_full == NULL)
target_full = malloc(full_len*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_full, full_len*sizeof(complex double));
if (target_full == NULL) abort();
if (!add) memset(target_full, 0, full_len*sizeof(complex double)); if (!add) memset(target_full, 0, full_len*sizeof(complex double));
const complex double one = 1; const complex double one = 1;
@ -989,11 +1079,11 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
{ // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC) { // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC)
size_t fullvec_offsetR = 0; size_t fullvec_offsetR = 0;
for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) { 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; const cart3_t posR = ss->p[piR].pos;
size_t fullvec_offsetC = 0; size_t fullvec_offsetC = 0;
for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { 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. if(piC != piR) { // The diagonal will be dealt with later.
const cart3_t posC = ss->p[piC].pos; const cart3_t posC = ss->p[piC].pos;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
@ -1003,7 +1093,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
} }
fullvec_offsetC += bspecC->n; fullvec_offsetC += bspecC->n;
} }
assert(fullvec_offsetC = full_len); assert(fullvec_offsetC == full_len);
fullvec_offsetR += bspecR->n; fullvec_offsetR += bspecR->n;
} }
assert(fullvec_offsetR == full_len); assert(fullvec_offsetR == full_len);
@ -1014,13 +1104,14 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
complex double *qpms_scatsys_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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, const qpms_scatsys_at_omega_t *ssw
complex double k ///< Wave number to use in the translation matrix.
) )
{ {
const complex double k = ssw->wavenumber;
const qpms_scatsys_t *ss = ssw->ss;
const size_t full_len = ss->fecv_size; const size_t full_len = ss->fecv_size;
if(!target) if(!target)
QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
@ -1031,13 +1122,13 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
{ // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC) { // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC)
size_t fullvec_offsetR = 0; size_t fullvec_offsetR = 0;
for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) { 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; const cart3_t posR = ss->p[piR].pos;
size_t fullvec_offsetC = 0; size_t fullvec_offsetC = 0;
// dest particle T-matrix // 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) { 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. if(piC != piR) { // The diagonal will be dealt with later.
const cart3_t posC = ss->p[piC].pos; const cart3_t posC = ss->p[piC].pos;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
@ -1064,19 +1155,20 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
} }
// Serial reference implementation. // 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. /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, qpms_iri_t iri, const qpms_scatsys_at_omega_t *ssw,
complex double k ///< Wave number to use in the translation matrix. 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]; const size_t packedlen = ss->saecv_sizes[iri];
if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
return target_packed; return target_packed;
if (target_packed == NULL) if (target_packed == NULL)
target_packed = malloc(SQ(packedlen)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
if (target_packed == NULL) abort();
memset(target_packed, 0, 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. // some of the following workspaces are probably redundant; TODO optimize later.
@ -1102,7 +1194,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
const size_t packed_orbit_offsetR = const size_t packed_orbit_offsetR =
ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiR] ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiR]
+ osnR * otR->irbase_sizes[iri]; + 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: // Orbit coeff vector's full size:
const size_t orbit_fullsizeR = otR->size * otR->bspecn; const size_t orbit_fullsizeR = otR->size * otR->bspecn;
const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n
@ -1112,7 +1204,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
const cart3_t posR = ss->p[piR].pos; const cart3_t posR = ss->p[piR].pos;
if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
// dest particle T-matrix // 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 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_oti_t otiC = ss->p_orbitinfo[piC].t;
const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC; const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC;
@ -1122,7 +1214,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
const size_t packed_orbit_offsetC = const size_t packed_orbit_offsetC =
ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC] ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
+ osnC * otC->irbase_sizes[iri]; + 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: // Orbit coeff vector's full size:
const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t orbit_fullsizeC = otC->size * otC->bspecn;
const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@ -1174,19 +1266,19 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_serial(
return target_packed; 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. /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, qpms_iri_t iri, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri
complex double k ///< Wave number to use in the translation matrix.
) )
{ {
const qpms_scatsys_t *ss = ssw->ss;
const complex double k = ssw->wavenumber;
const size_t packedlen = ss->saecv_sizes[iri]; const size_t packedlen = ss->saecv_sizes[iri];
if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
return target_packed; return target_packed;
if (target_packed == NULL) if (target_packed == NULL)
target_packed = malloc(SQ(packedlen)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
if (target_packed == NULL) abort();
memset(target_packed, 0, 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. // some of the following workspaces are probably redundant; TODO optimize later.
@ -1215,7 +1307,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n 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. // This is the orbit-level matrix projecting the whole orbit onto the irrep.
const complex double *omR = otR->irbases + otR->irbase_offsets[iri]; const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
// Orbit coeff vector's full size: // Orbit coeff vector's full size:
@ -1231,7 +1323,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
assert(ss->p_orbitinfo[piR].osn == osnR); assert(ss->p_orbitinfo[piR].osn == osnR);
const cart3_t posR = ss->p[piR].pos; const cart3_t posR = ss->p[piR].pos;
// dest particle T-matrix // 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 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_oti_t otiC = ss->p_orbitinfo[piC].t;
const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC; const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC;
@ -1241,7 +1333,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
const size_t packed_orbit_offsetC = const size_t packed_orbit_offsetC =
ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC] ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
+ osnC * otC->irbase_sizes[iri]; + 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: // Orbit coeff vector's full size:
const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t orbit_fullsizeC = otC->size * otC->bspecn;
const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@ -1294,20 +1386,21 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
return target_packed; 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_t *ss; const qpms_scatsys_at_omega_t *ssw;
qpms_ss_pi_t *opistartR_ptr; qpms_ss_pi_t *opistartR_ptr;
pthread_mutex_t *opistartR_mutex; pthread_mutex_t *opistartR_mutex;
qpms_iri_t iri; qpms_iri_t iri;
complex double *target_packed; complex double *target_packed;
complex double k;
}; };
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; *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 qpms_iri_t iri = a->iri;
const size_t packedlen = ss->saecv_sizes[iri]; const size_t packedlen = ss->saecv_sizes[iri];
@ -1347,7 +1440,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n 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. // This is the orbit-level matrix projecting the whole orbit onto the irrep.
const complex double *omR = otR->irbases + otR->irbase_offsets[iri]; const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
// Orbit coeff vector's full size: // Orbit coeff vector's full size:
@ -1363,7 +1456,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
assert(ss->p_orbitinfo[piR].osn == osnR); assert(ss->p_orbitinfo[piR].osn == osnR);
const cart3_t posR = ss->p[piR].pos; const cart3_t posR = ss->p[piR].pos;
// dest particle T-matrix // 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 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_oti_t otiC = ss->p_orbitinfo[piC].t;
const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC; const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC;
@ -1373,7 +1466,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
const size_t packed_orbit_offsetC = const size_t packed_orbit_offsetC =
ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC] ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
+ osnC * otC->irbase_sizes[iri]; + 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: // Orbit coeff vector's full size:
const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t orbit_fullsizeC = otC->size * otC->bspecn;
const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@ -1387,7 +1480,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
Sblock, // Sblock is S(piR->piC) Sblock, // Sblock is S(piR->piC)
bspecR, bspecC->n, bspecC, 1, bspecR, bspecC->n, bspecC, 1,
a->k, posR, posC, QPMS_HANKEL_PLUS)); k, posR, posC, QPMS_HANKEL_PLUS));
SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans, SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans,
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
@ -1484,7 +1577,7 @@ static void *qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thre
if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n 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 = qpms_ss_bspec_pi(ss, orbitstartpiR);
// This is the orbit-level matrix projecting the whole orbit onto the irrep. // This is the orbit-level matrix projecting the whole orbit onto the irrep.
const complex double *omR = otR->irbases + otR->irbase_offsets[iri]; const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
// Orbit coeff vector's full size: // Orbit coeff vector's full size:
@ -1508,7 +1601,7 @@ static void *qpms_scatsys_build_translation_matrix_e_irrep_packed_parallelR_thre
const size_t packed_orbit_offsetC = const size_t packed_orbit_offsetC =
ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC] ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
+ osnC * otC->irbase_sizes[iri]; + 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 = qpms_ss_bspec_pi(ss, piC);
// Orbit coeff vector's full size: // Orbit coeff vector's full size:
const size_t orbit_fullsizeC = otC->size * otC->bspecn; const size_t orbit_fullsizeC = otC->size * otC->bspecn;
const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
@ -1567,7 +1660,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
qpms_iri_t iri, qpms_iri_t iri,
complex double k, ///< Wave number to use in the translation matrix. const complex double k,
qpms_bessel_t J qpms_bessel_t J
) )
{ {
@ -1576,8 +1669,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? if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
return target_packed; return target_packed;
if (target_packed == NULL) if (target_packed == NULL)
target_packed = malloc(SQ(packedlen)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_packed, SQ(packedlen)*sizeof(complex double));
if (target_packed == NULL) abort();
memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
qpms_ss_pi_t opistartR = 0; qpms_ss_pi_t opistartR = 0;
@ -1618,26 +1710,24 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
// Parallel implementation, now default // 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. /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
complex double *target_packed, complex double *target_packed,
const qpms_scatsys_t *ss, qpms_iri_t iri, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri
complex double k ///< Wave number to use in the translation matrix.
) )
{ {
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? if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
return target_packed; return target_packed;
if (target_packed == NULL) if (target_packed == NULL)
target_packed = malloc(SQ(packedlen)*sizeof(complex double)); QPMS_CRASHING_MALLOC(target_packed,SQ(packedlen)*sizeof(complex double));
if (target_packed == NULL) abort();
memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
qpms_ss_pi_t opistartR = 0; qpms_ss_pi_t opistartR = 0;
pthread_mutex_t opistartR_mutex; pthread_mutex_t opistartR_mutex;
QPMS_ENSURE_SUCCESS(pthread_mutex_init(&opistartR_mutex, NULL)); 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 = {ss, &opistartR, &opistartR_mutex, iri, target_packed, k}; arg = {ssw, &opistartR, &opistartR_mutex, iri, target_packed};
// FIXME THIS IS NOT PORTABLE: // FIXME THIS IS NOT PORTABLE:
long nthreads; long nthreads;
@ -1658,7 +1748,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
pthread_t thread_ids[nthreads]; pthread_t thread_ids[nthreads];
for(long thi = 0; thi < nthreads; ++thi) for(long thi = 0; thi < nthreads; ++thi)
QPMS_ENSURE_SUCCESS(pthread_create(thread_ids + thi, NULL, 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)); (void *) &arg));
for(long thi = 0; thi < nthreads; ++thi) { for(long thi = 0; thi < nthreads; ++thi) {
void *retval; void *retval;
@ -1696,18 +1786,18 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
#endif #endif
complex double *qpms_scatsys_apply_Tmatrices_full( complex double *qpms_scatsysw_apply_Tmatrices_full(
complex double *target_full, const complex double *inc_full, complex double *target_full, const complex double *inc_full,
const qpms_scatsys_t *ss) { const qpms_scatsys_at_omega_t *ssw) {
QPMS_UNTESTED; QPMS_UNTESTED;
const qpms_scatsys_t *ss = ssw->ss;
if (!target_full) QPMS_CRASHING_CALLOC(target_full, ss->fecv_size, if (!target_full) QPMS_CRASHING_CALLOC(target_full, ss->fecv_size,
sizeof(complex double)); sizeof(complex double));
for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
complex double *ptarget = target_full + ss->fecv_pstarts[pi]; complex double *ptarget = target_full + ss->fecv_pstarts[pi];
const complex double *psrc = inc_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. // 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); qpms_apply_tmatrix(ptarget, psrc, T);
} }
return target_full; return target_full;
@ -1751,8 +1841,9 @@ void qpms_ss_LU_free(qpms_ss_LU lu) {
free(lu.ipiv); 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_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"); 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)); 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, QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, ss->fecv_size, ss->fecv_size,
@ -1760,45 +1851,45 @@ qpms_ss_LU qpms_scatsys_modeproblem_matrix_full_factorise(complex double *mpmatr
qpms_ss_LU lu; qpms_ss_LU lu;
lu.a = mpmatrix_full; lu.a = mpmatrix_full;
lu.ipiv = target_piv; lu.ipiv = target_piv;
lu.ss = ss; lu.ssw = ssw;
lu.full = true; lu.full = true;
lu.iri = -1; lu.iri = -1;
return lu; 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_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"); 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)); if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, n * sizeof(int));
QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n,
mpmatrix_packed, n, target_piv)); mpmatrix_packed, n, target_piv));
qpms_ss_LU lu; qpms_ss_LU lu;
lu.a = mpmatrix_packed; lu.a = mpmatrix_packed;
lu.ipiv = target_piv; lu.ipiv = target_piv;
lu.ss = ss; lu.ssw = ssw;
lu.full = false; lu.full = false;
lu.iri = iri; lu.iri = iri;
return lu; 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, complex double *target, int *target_piv,
const qpms_scatsys_t *ss, complex double k){ const qpms_scatsys_at_omega_t *ssw){
target = qpms_scatsys_build_modeproblem_matrix_full(target, ss, k); target = qpms_scatsysw_build_modeproblem_matrix_full(target, ssw);
return qpms_scatsys_modeproblem_matrix_full_factorise(target, target_piv, ss); 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, complex double *target, int *target_piv,
const qpms_scatsys_t *ss, qpms_iri_t iri, complex double k){ const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri){
target = qpms_scatsys_build_modeproblem_matrix_irrep_packed(target, ss, iri, k); target = qpms_scatsysw_build_modeproblem_matrix_irrep_packed(target, ssw, iri);
return qpms_scatsys_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ss, iri); return qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(target, target_piv, ssw, iri);
} }
complex double *qpms_scatsys_scatter_solve( complex double *qpms_scatsys_scatter_solve(
complex double *f, const complex double *a_inc, qpms_ss_LU lu) { 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)); if (!f) QPMS_CRASHING_MALLOC(f, n * sizeof(complex double));
memcpy(f, a_inc, n*sizeof(complex double)); // It will be rewritten by zgetrs 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*/, QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/, n /*n*/, 1 /*nrhs number of right hand sides*/,

160
qpms/scatsystem.h
View File

@ -12,6 +12,7 @@
#define QPMS_SCATSYSTEM_H #define QPMS_SCATSYSTEM_H
#include "qpms_types.h" #include "qpms_types.h"
#include "vswf.h" #include "vswf.h"
#include "tmatrices.h"
#include <stdbool.h> #include <stdbool.h>
/// Overrides the number of threads spawned by the paralellized functions. /// Overrides the number of threads spawned by the paralellized functions.
@ -19,10 +20,13 @@
void qpms_scatsystem_set_nthreads(long n); void qpms_scatsystem_set_nthreads(long n);
/// A particle, defined by its T-matrix and position. /// 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 { typedef struct qpms_particle_t {
// Does it make sense to ever use other than cartesian coords for this? // Does it make sense to ever use other than cartesian coords for this?
cart3_t pos; ///< Particle position in cartesian coordinates. 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; } qpms_particle_t;
struct qpms_finite_group_t; struct qpms_finite_group_t;
@ -121,13 +125,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_orbit_pi_t p; ///< Order (sija, ei rankki) of the particle inside that orbit type.
} qpms_ss_particle_orbitinfo_t; } 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_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;
struct qpms_trans_calculator; struct qpms_trans_calculator;
struct qpms_epsmu_generator_t;
typedef struct qpms_scatsys_t { typedef struct qpms_scatsys_t {
// TODO does bspec belong here? struct qpms_epsmu_generator_t medium; ///< Optical properties of the background medium.
qpms_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_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.
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_ss_tmi_t tm_capacity; ///< Capacity of tm[].
qpms_particle_tid_t *p; ///< Particles. qpms_particle_tid_t *p; ///< Particles.
qpms_ss_pi_t p_count; ///< Size of particles array. qpms_ss_pi_t p_count; ///< Size of particles array.
qpms_ss_pi_t p_capacity; ///< Capacity of p[]. qpms_ss_pi_t p_capacity; ///< Capacity of p[].
@ -167,23 +190,77 @@ typedef struct qpms_scatsys_t {
struct qpms_trans_calculator *c; struct qpms_trans_calculator *c;
} qpms_scatsys_t; } qpms_scatsys_t;
/// Convenience function to access pi'th particle's bspec. /// Retrieve the bspec of \a tmi'th element of \a ss->tm.
static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) { static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi) {
return ss->tm[ss->p[pi].tmatrix_id]->spec; return ss->tmg[ss->tm[tmi].tmgi].spec;
} }
/// Creates a new scatsys by applying a symmetry group, copying particles if needed. /// Retrieve the bspec of \a pi'th particle in \a ss->p.
/** In fact, it copies everything except the vswf set specs, so keep them alive until scatsys is destroyed. static inline const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi) {
* The following fields must be filled: return ss->tmg[ss->tm[ss->p[pi].tmatrix_id].tmgi].spec;
* orig->tm }
* orig->tm_count
* orig->p typedef struct qpms_scatsys_at_omega_t {
* orig->p_count 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 \a ss->tm.
*/
qpms_tmatrix_t **tm;
complex double omega; ///< Angular frequency
qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency
complex double wavenumber; ///< Background medium wavenumber
} 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.
*
* The following fields must be filled in the "proto- scattering system" \a orig:
* * 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;
* 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
*
* 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.
*
* 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_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,
complex double omega, const struct qpms_tolerance_spec_t *tol);
/// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load. /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
void qpms_scatsys_free(qpms_scatsys_t *s); void qpms_scatsys_free(qpms_scatsys_t *s);
/// 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. /// 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. /** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient.
* *
@ -266,37 +343,36 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
); );
/// Creates the full \f$ (I - TS) \f$ matrix of the scattering system. /// Creates the full \f$ (I - TS) \f$ matrix of the scattering system.
complex double *qpms_scatsys_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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, const qpms_scatsys_at_omega_t *ssw
complex double k ///< Wave number to use in the translation matrix.
); );
/// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form. /// 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, qpms_iri_t iri, const qpms_scatsys_at_omega_t *ssw,
complex double k ///< Wave number to use in the translation matrix. qpms_iri_t iri
); );
/// Alternative implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). /// Alternative implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_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->fecv_size**2 elements. If NULL, new will be allocated. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, qpms_iri_t iri, const qpms_scatsys_at_omega_t *ssw,
complex double k ///< Wave number to use in the translation matrix. qpms_iri_t iri
); );
/// Alternative (serial reference) implementation of qpms_scatsys_build_modeproblem_matrix_irrep_packed(). /// Alternative (serial reference) implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed().
complex double *qpms_scatsys_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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, qpms_iri_t iri, const qpms_scatsys_at_omega_t *ssw,
complex double k ///< Wave number to use in the translation matrix. qpms_iri_t iri
); );
/// LU factorisation (LAPACKE_zgetrf) result holder. /// LU factorisation (LAPACKE_zgetrf) result holder.
typedef struct qpms_ss_LU { 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. bool full; ///< true if full matrix; false if irrep-packed.
qpms_iri_t iri; ///< Irrep index if `full == false`. qpms_iri_t iri; ///< Irrep index if `full == false`.
/// LU decomposition array. /// LU decomposition array.
@ -307,33 +383,33 @@ typedef struct qpms_ss_LU {
void qpms_ss_LU_free(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. /// 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). 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). 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
complex double k ///< Wave number to use in the translation matrix.
); );
/// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch. /// 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). 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). 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,
complex double k ///< Wave number to use in the translation matrix. qpms_iri_t iri
); );
/// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents. /// 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. 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). 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. /// 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. 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). 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$. /// Solves a (possibly partial, irrep-packed) scattering problem \f$ (I-TS)f = Ta_\mathrm{inc} \f$ using a pre-factorised \f$ (I-TS) \f$.
@ -422,10 +498,10 @@ complex double *qpms_scatsys_incident_field_vector_full(
); );
/// Applies T-matrices onto an incident field vector in the full basis. /// 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. 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 complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
const qpms_scatsys_t *ss const qpms_scatsys_at_omega_t *ssw
); );

230
qpms/tmatrices.c
View File

@ -53,6 +53,15 @@ qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
return t; 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 *qpms_tmatrix_copy(const qpms_tmatrix_t *T) {
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec); qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
size_t n = T->spec->n; size_t n = T->spec->n;
@ -991,3 +1000,224 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
return QPMS_SUCCESS; return QPMS_SUCCESS;
} }
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);
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_SYMMETRISE: // 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->opmem[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)
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);
}
free(o->ops);
}
break;
default:
QPMS_WTF;
}
}
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->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;
default:
QPMS_WTF;
}
dest->typ = src->typ;
}
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;
}
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) {
// 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_NOOP:
return T;
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;
}
}

154
qpms/tmatrices.h
View File

@ -8,6 +8,8 @@
#include "materials.h" #include "materials.h"
#include <stdio.h> #include <stdio.h>
struct qpms_finite_group_t; struct qpms_finite_group_t;
typedef struct qpms_finite_group_t qpms_finite_group_t; typedef struct qpms_finite_group_t qpms_finite_group_t;
@ -17,11 +19,21 @@ static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
} }
/// Initialises a zero T-matrix. /// 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); 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. /// Copies a T-matrix, allocating new array for the T-matrix data.
qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T); 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. /// Destroys a T-matrix.
void qpms_tmatrix_free(qpms_tmatrix_t *t); void qpms_tmatrix_free(qpms_tmatrix_t *t);
@ -233,7 +245,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
complex double *qpms_apply_tmatrix( 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. 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 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. /// Generic T-matrix generator function that fills a pre-initialised qpms_tmatrix_t given a frequency.
@ -251,6 +263,13 @@ typedef struct qpms_tmatrix_generator_t {
const void *params; ///< Parameter pointer passed to the function. const void *params; ///< Parameter pointer passed to the function.
} qpms_tmatrix_generator_t; } qpms_tmatrix_generator_t;
/// Initialises and evaluates a new T-matrix using a generator.
qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
const qpms_vswf_set_spec_t *bspec,
qpms_tmatrix_generator_t gen,
complex double omega);
/// Implementation of qpms_matrix_generator_t that just copies a constant matrix. /// 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 /** N.B. this does almost no checks at all, so use it only for t-matrices with
* the same base spec. * the same base spec.
@ -427,7 +446,7 @@ qpms_arc_function_retval_t qpms_arc_cylinder(double theta,
); );
/// Arc parametrisation of spherical particle; for qpms_arc_function_t. /// 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, qpms_arc_function_retval_t qpms_arc_sphere(double theta,
const void *R; ///< Points to double containing particle's radius. const void *R; ///< Points to double containing particle's radius.
); );
@ -503,6 +522,137 @@ 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;
/// 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, *f.gen, omega);
}
/// 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.
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! */
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;
};
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.
*/
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.)
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.
};
/// 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! */
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.
};
/// 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;
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_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 *finite_group;
} 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);
/// 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.
* 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);
/// (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);
/// Inits a new "chain" of composed operations, some of which might be owned.
void qpms_tmatrix_operation_compose_chain_init(
qpms_tmatrix_operation_t *target, ///< The operation structure that will be set to the chain.
size_t nops, ///< Number of chained operations (length of the \a ops array)
size_t opmem_size ///< Size of the own operations buffer (length of the \a opmem array)
);
#if 0 #if 0
// Abstract types that describe T-matrix/particle/scatsystem symmetries // Abstract types that describe T-matrix/particle/scatsystem symmetries
// To be implemented later. See also the thoughts in the beginning of groups.h. // To be implemented later. See also the thoughts in the beginning of groups.h.

15
qpms/tolerances.h Normal file
View File

@ -0,0 +1,15 @@
/*! \file tolerances.h */
#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