Scatsys full mode problem matrix; wrong result.

Former-commit-id: d09e4210fa543068e0e6a7df80a1bef134e60ca7

qpms/qpms_c.pyx

@@ -1427,7 +1427,6 @@ cdef class ScatteringSystem:
                 self.s, iri, 0)
         return target_np
     def pack_matrix(self, fullmatrix, iri):
-        (iri < self.nirreps)
         cdef size_t flen = self.s[0].fecv_size
         cdef size_t rlen = self.saecv_sizes[iri]
         fullmatrix = np.array(fullmatrix, dtype=complex, copy=False, order='C')
@@ -1442,7 +1441,6 @@ cdef class ScatteringSystem:
                 self.s, iri)
         return target_np
     def unpack_matrix(self, packedmatrix, iri):
-        (iri < self.nirreps)
         cdef size_t flen = self.s[0].fecv_size
         cdef size_t rlen = self.saecv_sizes[iri]
         packedmatrix = np.array(packedmatrix, dtype=complex, copy=False, order='C')
@@ -1457,7 +1455,13 @@ cdef class ScatteringSystem:
                 self.s, iri, 0)
         return target_np
 
-
+    def modeproblem_matrix_full(self, double k):
+        cdef size_t flen = self.s[0].fecv_size
+        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
+                (flen,flen),dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target
+        qpms_scatsys_build_modeproblem_matrix_full(&target_view[0][0], self.s, k)
+        return target
 
 
 def tlm2uvswfi(t, l, m):

qpms/qpms_cdefs.pxd

@@ -285,6 +285,10 @@ cdef extern from "scatsystem.h":
             const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
     cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full,
             const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
+    cdouble *qpms_scatsys_build_modeproblem_matrix_full(cdouble *target,
+            const qpms_scatsys_t *ss, double k)
+    cdouble *qpms_scatsys_build_modeproblem_matrix_irrep_packed(cdouble *target,
+            const qpms_scatsys_t *ss, qpms_iri_t iri, double k)
 
 
 

qpms/qpms_error.h

@@ -19,4 +19,10 @@ void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum,
 //void qpms_error_at_line(const char *filename, unsigned int linenum,
 //		const char *fmt, ...);
 
+#define QPMS_CRASHING_MALLOC(pointer, size) {(pointer) = malloc(size); if(!pointer) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t) (size));}
+
+#define QPMS_WTF qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error.")
+
+#define QPMS_ENSURE_SUCCESS(x) {if(x) QPMS_WTF;}
+
 #endif 

qpms/scatsystem.c

@@ -13,6 +13,7 @@
 #include "wigner.h"
 #include <string.h>
 #include "qpms_error.h"
+#include "translations.h"
 
 #define SQ(x) ((x)*(x))
 #define QPMS_SCATSYS_LEN_RTOL 1e-13
@@ -413,6 +414,9 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu
 qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) {
   // TODO check data sanity
 
+  qpms_l_t lMax = 0; // the overall lMax of all base specs.
+  qpms_normalisation_t normalisation = QPMS_NORMALISATION_UNDEF;
+
   // First, determine the rough radius of the array
   double lenscale = 0;
   {
@@ -468,6 +472,11 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
     } 
     tm_dupl_remap[i] = j;
     ss->max_bspecn = MAX(ss->tm[i]->spec->n, ss->max_bspecn);
+    if (normalisation == QPMS_NORMALISATION_UNDEF)
+      normalisation = ss->tm[i]->spec->norm;
+    // We expect all bspec norms to be the same.
+    else assert(normalisation == ss->tm[i]->spec->norm); 
+    lMax = MAX(lMax, ss->tm[i]->spec->lMax);
   }
 
   // Copy particles, remapping the t-matrix indices
@@ -631,6 +640,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
     }
   }
   
+  ss->c = qpms_trans_calculator_init(lMax, normalisation);
   return ss;
 }
 
@@ -646,6 +656,7 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) {
     free(ss->p_orbitinfo);
     free(ss->orbit_types);
     free(ss->saecv_sizes);
+    qpms_trans_calculator_free(ss->c);
   }
   free(ss);
 }
@@ -796,7 +807,7 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
 
   // Workspace for the intermediate particle-orbit matrix result
   complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn)
-      * ss->sym->order); if (!tmp) abort;
+      * ss->sym->order); if (!tmp) abort();
 
   const complex double one = 1, zero = 0;
 
@@ -875,7 +886,7 @@ complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
 
   // Workspace for the intermediate particle-orbit matrix result
   complex double *tmp = malloc(sizeof(complex double) * SQ(ss->max_bspecn)
-      * ss->sym->order); if (!tmp) abort;
+      * ss->sym->order); if (!tmp) abort();
 
   const complex double one = 1, zero = 0;
 
@@ -1021,4 +1032,61 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
   return target_full; 
 }
 
+complex double *qpms_scatsys_build_modeproblem_matrix_full(
+    /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
+    complex double *target,
+    const qpms_scatsys_t *ss,
+    double k ///< Wave number to use in the translation matrix.
+    )
+{
+  const size_t full_len = ss->fecv_size;
+  if(!target)
+    QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
+  complex double *tmp;
+  QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double));
+  memset(target, 0, SQ(full_len) * sizeof(complex double)); //unnecessary?
+  complex double one = 1, zero = 0;
+  { // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC)
+    size_t fullvec_offsetR = 0;
+    for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) {
+      const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec;
+      const cart3_t posR = ss->p[piR].pos;
+      size_t fullvec_offsetC = 0;
+      // dest particle T-matrix
+      const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m;
+      for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) {
+        if(piC == piR) continue; // The diagonal will be dealt with later.
+        const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec;
+        const cart3_t posC = ss->p[piC].pos;
+        QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
+              tmp, // tmp is S(piR<-piC)
+              bspecR, bspecC->n, bspecC, 1,
+              k, posR, posC));      
+        cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+              bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 
+              &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/,
+              tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,
+              target + fullvec_offsetR*full_len + fullvec_offsetC /*c*/,
+              full_len /*ldc*/);
+        fullvec_offsetC += bspecC->n;
+      }
+      fullvec_offsetR += bspecR->n;
+    }
+  }
+  // diagonal part M[pi,pi] = -1
+  for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] = -1;
+  
+  free(tmp);
+  return target;
+}
+
+
+
+
+complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
+    /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
+    complex double *target,
+    const qpms_scatsys_t *ss, qpms_iri_t iri,
+    double k ///< Wave number to use in the translation matrix.
+    );
 

qpms/scatsystem.h

@@ -370,6 +370,7 @@ typedef struct qpms_ss_particle_orbitinfo {
 	qpms_ss_orbit_pi_t p; ///< Order (sija, ei rankki) of the particle inside that orbit type.
 } qpms_ss_particle_orbitinfo_t;
 
+struct qpms_trans_calculator;
 
 typedef struct qpms_scatsys_t {
 	// TODO does bspec belong here?
@@ -410,6 +411,7 @@ typedef struct qpms_scatsys_t {
 	char *otspace;
 	char *otspace_end;
 	double lenscale; // radius of the array, used as a relative tolerance measure
+	struct qpms_trans_calculator *c;
 } qpms_scatsys_t;
 
 /// Creates a new scatsys by applying a symmetry group, copying particles if needed.
@@ -448,6 +450,19 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
 		const complex double *orig_packed, const qpms_scatsys_t *ss,
 		qpms_iri_t iri, bool add);
 
+complex double *qpms_scatsys_build_modeproblem_matrix_full(
+		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
+		complex double *target,
+		const qpms_scatsys_t *ss,
+		double k ///< Wave number to use in the translation matrix.
+		);
+complex double *qpms_scatsys_build_modeproblem_irrep_packed(
+		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
+		complex double *target,
+		const qpms_scatsys_t *ss, qpms_iri_t iri,
+		double k ///< Wave number to use in the translation matrix.
+		);
+
 /// NOT IMPLEMENTED Dumps a qpms_scatsys_t structure to a file.
 qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path);
 
@@ -512,6 +527,7 @@ complex double *qpms_orbit_irrep_basis(
 		const qpms_iri_t iri);
 
 
+
 #if 0
 // Abstract types that describe T-matrix/particle/scatsystem symmetries
 // To be implemented later. See also the thoughts in the beginning of groups.h.

qpms/translations.c

@@ -1185,7 +1185,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
         = (srct == destt) ? A[desty][srcy] : B[desty][srcy];
     }
   }
-  return QPMS_SUCCESS;
+  return retval;
 }
 
 /// Version with \a k and cartesian particle positions

tests/scatsys.py

@@ -19,7 +19,7 @@ packedvectors = [(iri, ss.pack_vector(fullvector, iri)) for iri in range(ss.nirr
 unpackedvectors = np.array([ss.unpack_vector(v[1], v[0]) for v in packedvectors])
 rec_fullvector = np.sum(unpackedvectors, axis=0)
 thediff = np.amax(abs(rec_fullvector-fullvector))
-assert(thediff < 1e-14)
+assert(thediff < 1e-8)
 
 packedmatrices = list()
 for iri in range(ss.nirreps):
@@ -36,3 +36,10 @@ for iri, m in packedmatrices:
     repackedmatrix = ss.pack_matrix(fullmatrix,iri)
     print(np.amax(abs(repackedmatrix-m)))
 
+k = 1.7
+modematrix_full = ss.modeproblem_matrix_full(k)
+modematrix_packed_list = [(iri, ss.pack_matrix(modematrix_full,iri)) for iri in range(ss.nirreps)]
+modematrix_full_rec = np.empty((ss.fecv_size, ss.fecv_size), dtype=complex)
+for iri, m in modematrix_packed_list:
+    modematrix_full_rec += ss.unpack_matrix(m,iri)
+print(np.amax(abs(modematrix_full-modematrix_full_rec)))