Irrep-decomposed scatsys beyn; fix missing FinitePointGroup reference

Former-commit-id: fa1032eb8fcb8ce1018b69fff5af6375b34115be

qpms/groups.h

@@ -81,6 +81,12 @@ static inline qpms_gmi_t qpms_finite_group_inv(const qpms_finite_group_t *G,
 	return G->invi[a];
 }
 
+static inline _Bool qpms_iri_is_valid(const qpms_finite_group_t *G, qpms_iri_t iri) {
+	return (iri > G->nirreps || iri < 0) ? 0 : 1;
+}
+
+
+
 /// NOT IMPLEMENTED Get irrep index by name.
 qpms_iri_t qpms_finite_group_find_irrep_by_name(qpms_finite_group_t *G, char *name);
 

qpms/qpms_c.pyx

@@ -280,7 +280,7 @@ cdef class Particle:
             if isinstance(tgen.holder, CTMatrix):
                 spec = (<CTMatrix>tgen.holder).spec
             else:
-                raise ValueError("bspec argument must be specified separately for str(type(t))")
+                raise ValueError("bspec argument must be specified separately for %s" % 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
@@ -340,6 +340,16 @@ cdef class ScatteringSystem:
     #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 FinitePointGroup sym
+
+    cdef qpms_iri_t iri_py2c(self, iri, allow_None = True):
+        if iri is None and allow_None:
+            return QPMS_NO_IRREP
+        cdef qpms_iri_t nir = self.nirreps
+        cdef qpms_iri_t ciri = iri
+        if ciri < 0 or ciri > nir:
+            raise ValueError("Invalid irrep index %s (of %d irreps)", str(iri), self.nirreps)
+        return ciri
 
     def check_s(self): # cdef instead?
         if self.s == <qpms_scatsys_t *>NULL:
@@ -416,6 +426,7 @@ cdef class ScatteringSystem:
         self.medium_holder = mediumgen
         self.s = ss
         self.tmgobjs = tmgobjs
+        self.sym = sym
         pyssw = _ScatteringSystemAtOmega()
         pyssw.ssw = ssw
         pyssw.ss_pyref = self
@@ -589,7 +600,7 @@ cdef class ScatteringSystem:
                 self.s, qpms_incfield_planewave, <void *>&p, 0)
         return target_np
 
-    def find_modes(self, cdouble omega_centre, double omega_rr, double omega_ri,
+    def find_modes(self, cdouble omega_centre, double omega_rr, double omega_ri, iri = None,
             size_t contour_points = 20, double rank_tol = 1e-4, size_t rank_min_sel=1,
             double res_tol = 0):
         """
@@ -597,6 +608,7 @@ cdef class ScatteringSystem:
         """
 
         cdef beyn_result_t *res = qpms_scatsys_finite_find_eigenmodes(self.s, 
+                self.iri_py2c(iri),
                 omega_centre, omega_rr, omega_ri, contour_points, 
                 rank_tol, rank_min_sel, res_tol)
         if res == NULL: raise RuntimeError
@@ -635,6 +647,7 @@ cdef class ScatteringSystem:
                 'residuals':residuals,
                 'eigval_err':eigval_err,
                 'ranktest_SV':ranktest_SV,
+                'iri': iri,
         }
 
         return retdict

qpms/qpms_cdefs.pxd

@@ -107,6 +107,7 @@ cdef extern from "qpms_types.h":
     ctypedef int32_t qpms_ss_pi_t
     ctypedef int qpms_gmi_t
     ctypedef int qpms_iri_t
+    qpms_iri_t QPMS_NO_IRREP
     ctypedef const char * qpms_permutation_t
     struct qpms_tmatrix_t:
         qpms_vswf_set_spec_t *spec
@@ -598,9 +599,9 @@ cdef extern from "scatsystem.h":
     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)
-    beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, cdouble omega_centre,
+    beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, qpms_iri_t iri,
-            double omega_rr, double omega_ri, size_t contour_npoints, double rank_tol, size_t rank_sel_min,
+            cdouble omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, 
-            double res_tol)
+            double rank_tol, size_t rank_sel_min, double res_tol)
 
 cdef extern from "ewald.h":
     struct qpms_csf_result:

qpms/qpms_types.h

@@ -309,6 +309,8 @@ typedef int qpms_gmi_t;
 
 /// Irreducible representation index. See also groups.h.
 typedef int qpms_iri_t;
+/// Constant labeling that no irrep decomposition is done (full system solved instead).
+#define QPMS_NO_IRREP ((qpms_iri_t) -1)
 
 /// Permutation representation of a group element.
 /** For now, it's just a string of the form "(0,1)(3,4,5)"

qpms/scatsystem.c

@@ -1895,29 +1895,50 @@ complex double *qpms_scatsys_scatter_solve(
   return f;
 }
 
+struct qpms_scatsys_finite_eval_Beyn_ImTS_param {
+  const qpms_scatsys_t *ss;
+  qpms_iri_t iri;
+};
+
 /// Wrapper for Beyn algorithm (non-periodic system)
 static int qpms_scatsys_finite_eval_Beyn_ImTS(complex double *target,
     size_t m, complex double omega, void *params) {
-  const qpms_scatsys_t *ss = params;
+  const struct qpms_scatsys_finite_eval_Beyn_ImTS_param *p = params;
-  qpms_scatsys_at_omega_t *ssw = qpms_scatsys_at_omega(ss, omega);
+  qpms_scatsys_at_omega_t *ssw = qpms_scatsys_at_omega(p->ss, omega);
   QPMS_ENSURE(ssw != NULL, "qpms_scatsys_at_omega() returned NULL");
-  QPMS_ENSURE(NULL != qpms_scatsysw_build_modeproblem_matrix_full(target, ssw),
+  if (p->iri == QPMS_NO_IRREP) {
-    "qpms_scatsysw_build_modeproblem_matrix_full() returned NULL");
+    QPMS_ASSERT(m == p->ss->fecv_size);
+    QPMS_ENSURE(NULL != qpms_scatsysw_build_modeproblem_matrix_full(
+          target, ssw),
+      "qpms_scatsysw_build_modeproblem_matrix_full() returned NULL");
+  } else {
+    QPMS_ASSERT(m == p->ss->saecv_sizes[p->iri]);
+    QPMS_ENSURE(NULL != qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
+          target, ssw, p->iri),
+      "qpms_scatsysw_build_modeproblem_matrix_irrep_packed() returned NULL");
+  }
   qpms_scatsys_at_omega_free(ssw);
   return QPMS_SUCCESS;
 }
 
-beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss,
+beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
+    const qpms_scatsys_t * const ss, const qpms_iri_t iri,
     complex double omega_centre, double omega_rr, double omega_ri, 
     size_t contour_npoints, 
     double rank_tol, size_t rank_sel_min, double res_tol) {
+  size_t n; // matrix dimension
+  if (qpms_iri_is_valid(ss->sym, iri)) {
+    n = ss->saecv_sizes[iri];
+  } else if (iri == QPMS_NO_IRREP) {
+    n = ss->fecv_size;
+  } else QPMS_WTF;
 
   beyn_contour_t *contour = beyn_contour_ellipse(omega_centre,
       omega_rr, omega_ri, contour_npoints);
-
+  
-  beyn_result_t *result = beyn_solve(ss->fecv_size, 
+  struct qpms_scatsys_finite_eval_Beyn_ImTS_param p = {ss, iri};
-      ss->fecv_size /* possibly make smaller? */,
+  beyn_result_t *result = beyn_solve(n, n /* possibly make smaller? */,
-      qpms_scatsys_finite_eval_Beyn_ImTS, NULL, (void *) ss,
+      qpms_scatsys_finite_eval_Beyn_ImTS, NULL, (void *) &p,
       contour, rank_tol, rank_sel_min, res_tol);
   QPMS_ENSURE(result != NULL, "beyn_solve() returned NULL");
 

qpms/scatsystem.h

@@ -528,7 +528,9 @@ struct beyn_result_t; // See beyn.h for full definition
  * rather than the linear algebra.
  */
 struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
-		const qpms_scatsys_t *ss,
+		const qpms_scatsys_t *ss, 
+		/// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system.
+		qpms_iri_t iri,
 		complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
 		double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
 		double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.