Merge branch 'fixprojs'

Former-commit-id: 9e2357e6cc6844360b51a55a8af984dec82167a7

qpms/qpms_c.pyx

@@ -700,7 +700,7 @@ cdef class BaseSpec:
         if 'norm' in kwargs.keys():
             self.s.norm = kwargs['norm']
         else:
-            self.s.norm = QPMS_NORMALISATION_POWER
+            self.s.norm = QPMS_NORMALISATION_POWER_CS
         # set the other metadata
         cdef qpms_l_t l
         self.s.lMax_L = -1

qpms/qpms_error.h

@@ -19,9 +19,9 @@ 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_CRASHING_MALLOC(pointer, size) {(pointer) = malloc(size); if(!pointer && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t) (size));}
 
-#define QPMS_CRASHING_REALLOC(pointer, size) {(pointer) = realloc(pointer, size); if(!pointer) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Rellocation of %zd bytes for " #pointer " failed.", (size_t) (size));}
+#define QPMS_CRASHING_REALLOC(pointer, size) {(pointer) = realloc(pointer, size); if(!pointer && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Rellocation of %zd bytes for " #pointer " failed.", (size_t) (size));}
 
 #define QPMS_WTF qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error.")
 

qpms/scatsystem.c

@@ -561,6 +561,10 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
       current_orbit[0] = pi;
       ot_current.size = 1; 
       ot_current.tmatrices[0] = ss->p[pi].tmatrix_id;
+#ifdef DUMP_PARTICLE_POSITIONS
+      cart3_t pos = ss->p[pi].pos;
+      fprintf(stderr, "An orbit [%.4g, %.4g, %.4g] => ", pos.x, pos.y, pos.z);
+#endif
     }
 
     for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi) {
@@ -583,6 +587,10 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
         qpms_particle_tid_t newparticle = {newpoint, new_tmi};
         ss->p[ss->p_count] = newparticle;
         ++(ss->p_count);
+#ifdef DUMP_PARTICLE_POSITIONS
+        if(new_orbit)
+          fprintf(stderr, "[%.4g, %.4g, %.4g] ", newpoint.x, newpoint.y, newpoint.z);
+#endif
       }
       ss->p_sym_map[gmi + pi * sym->order] = pj;
 
@@ -600,6 +608,9 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
       }
     }
     if (new_orbit) { // Now compare if the new orbit corresponds to some of the existing types.
+#ifdef DUMP_PARTICLE_POSITIONS
+      fputc('\n', stderr);
+#endif
       qpms_ss_oti_t oti;
       for(oti = 0; oti < ss->orbit_type_count; ++oti)
         if (orbit_types_equal(&ot_current, &(ss->orbit_types[oti]))) break; // HIT, orbit type already exists
@@ -719,7 +730,7 @@ complex double *qpms_orbit_action_matrix(complex double *target,
     const qpms_gmi_t Row = ot->action[sym->order * Col + g];
     for(size_t row = 0; row < bspec->n; ++row)
       for(size_t col = 0; col < bspec->n; ++col)
-        target[n*n*N*Row + n*Col + n*N*row + col] = tmp[row][col]; //CHECKCONJ
+        target[n*n*N*Row + n*Col + n*N*row + col] = conj(tmp[row][col]); //CHECKCONJ
   }
 #ifdef DUMP_ACTIONMATRIX
   fprintf(stderr,"%d: %s\n", 
@@ -855,6 +866,119 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
   return target;
 }
 
+complex double *qpms_scatsys_irrep_transform_matrix(complex double *U,
+    const qpms_scatsys_t *ss, qpms_iri_t iri) {
+  const size_t packedlen = ss->saecv_sizes[iri];
+  const size_t full_len = ss->fecv_size;
+  if (U == NULL)
+    QPMS_CRASHING_MALLOC(U,full_len * packedlen * sizeof(complex double));
+  memset(U, 0, full_len * packedlen * sizeof(complex double));
+
+  size_t fullvec_offset = 0;
+
+  for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
+    const qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t;
+    const qpms_ss_orbit_type_t *const ot = ss->orbit_types + oti;
+    const qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn;
+    const qpms_ss_orbit_pi_t opi = ss->p_orbitinfo[pi].p;
+    // This is where the particle's orbit starts in the "packed" vector:
+    const size_t packed_orbit_offset = ss->saecv_ot_offsets[iri*ss->orbit_type_count + oti]
+      + osn * ot->irbase_sizes[iri];
+    // Orbit coeff vector's full size:
+    const size_t orbit_fullsize = ot->size * ot->bspecn;
+    const size_t particle_fullsize = ot->bspecn;
+    const size_t orbit_packedsize = ot->irbase_sizes[iri];
+    // This is the orbit-level matrix projecting the whole orbit onto the irrep.
+    const complex double *om = ot->irbases + ot->irbase_offsets[iri];
+
+    for (size_t prow = 0; prow < orbit_packedsize; ++prow) 
+      for (size_t pcol = 0; pcol < ot->bspecn; ++pcol) 
+        U[full_len * (packed_orbit_offset + prow) + (fullvec_offset + pcol)]
+          = om[orbit_fullsize * prow + (opi * ot->bspecn + pcol)];
+    fullvec_offset += ot->bspecn;
+  }
+
+  return U;
+}
+
+complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed,
+		const complex double *orig_full, const qpms_scatsys_t *ss,
+		qpms_iri_t iri){
+  const size_t packedlen = ss->saecv_sizes[iri];
+  if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
+    return target_packed; 
+  const size_t full_len = ss->fecv_size;
+  if (target_packed == NULL)
+    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
+  if (target_packed == NULL) abort();
+  memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
+
+  // Workspace for the intermediate matrix
+  complex double *tmp; 
+  QPMS_CRASHING_MALLOC(tmp, full_len * packedlen * sizeof(complex double));
+
+  complex double *U = qpms_scatsys_irrep_transform_matrix(NULL, ss, iri);
+
+  const complex double one = 1, zero = 0;
+  
+  // tmp = F U*
+  cblas_zgemm(
+      CblasRowMajor, CblasNoTrans, CblasConjTrans,
+      full_len /*M*/, packedlen /*N*/, full_len /*K*/,
+      &one /*alpha*/, orig_full/*A*/, full_len/*ldA*/, 
+      U /*B*/, full_len/*ldB*/, 
+      &zero /*beta*/, tmp /*C*/, packedlen /*LDC*/);
+  // target = U tmp
+  cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+      packedlen /*M*/, packedlen /*N*/, full_len /*K*/,
+      &one /*alpha*/, U/*A*/, full_len/*ldA*/,
+      tmp /*B*/, packedlen /*ldB*/, &zero /*beta*/,
+      target_packed /*C*/, packedlen /*ldC*/);
+  
+  free(tmp);
+  free(U);
+  return target_packed;
+}
+  
+/// Transforms a big "packed" matrix into the full basis (trivial implementation).
+complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full,
+		const complex double *orig_packed, const qpms_scatsys_t *ss,
+		qpms_iri_t iri, bool add){
+  const size_t packedlen = ss->saecv_sizes[iri];
+  const size_t full_len = ss->fecv_size;
+  if (target_full == NULL)
+    target_full = malloc(SQ(full_len)*sizeof(complex double));
+  if (target_full == NULL) abort();
+  if(!add) memset(target_full, 0, SQ(full_len)*sizeof(complex double));
+
+  if(!packedlen) return target_full; // Empty irrep, do nothing.
+
+  // Workspace for the intermediate matrix result
+  complex double *tmp; 
+  QPMS_CRASHING_MALLOC(tmp, full_len * packedlen * sizeof(complex double));
+
+  complex double *U = qpms_scatsys_irrep_transform_matrix(NULL, ss, iri);
+
+  const complex double one = 1, zero = 0;
+  
+  // tmp = P U
+  cblas_zgemm(
+      CblasRowMajor, CblasNoTrans, CblasNoTrans,
+      packedlen /*M*/, full_len /*N*/, packedlen /*K*/,
+      &one /*alpha*/, orig_packed/*A*/, packedlen/*ldA*/, 
+      U /*B*/, full_len/*ldB*/, 
+      &zero /*beta*/, tmp /*C*/, full_len /*LDC*/);
+  // target += U* tmp
+  cblas_zgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans,
+      full_len /*M*/, full_len /*N*/, packedlen /*K*/,
+      &one /*alpha*/, U/*A*/, full_len/*ldA*/,
+      tmp /*B*/, full_len /*ldB*/, &one /*beta*/,
+      target_full /*C*/, full_len /*ldC*/);
+  free(tmp);
+  free(U);
+  return target_full;
+}
+
 complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
 		const complex double *orig_full, const qpms_scatsys_t *ss,
 		qpms_iri_t iri){
@@ -929,13 +1053,14 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
         }
         fullvec_offsetC += otC->bspecn;
       }
-      fullvec_offsetR += otR->bspecn;
     }
+    fullvec_offsetR += otR->bspecn;
   }
   free(tmp);
   return target_packed;
 }
 
+
 /// Transforms a big "packed" matrix into the full basis.
 /** TODO doc */
 complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,

qpms/scatsystem.h

@@ -426,6 +426,26 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const st
 /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
 void qpms_scatsys_free(qpms_scatsys_t *s);
 
+/// 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.
+ * 
+ * TODO doc about shape etc.
+ */
+complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U,
+		const qpms_scatsys_t *ss, qpms_iri_t iri);
+
+/// Projects a "big" matrix onto an irrep (slow reference implementation).
+/** TODO doc */
+complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed,
+		const complex double *orig_full, const qpms_scatsys_t *ss,
+		qpms_iri_t iri);
+
+/// Transforms a big "packed" matrix into the full basis (slow reference implementation).
+/** TODO doc */
+complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full,
+		const complex double *orig_packed, const qpms_scatsys_t *ss,
+		qpms_iri_t iri, bool add);
+
 /// Projects a "big" matrix onto an irrep.
 /** TODO doc */
 complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,

qpms/symmetries.py

@@ -53,7 +53,7 @@ class SVWFPointGroupInfo: # only for point groups, coz in svwf_rep() I use I_tyt
         self.svwf_rep_gen_func = svwf_rep_gen_func
         self.irreps = dict()
         for irrepname, irrepgens in irrepgens_dict.items():
-            is1d = isinstance(irrepgens[0], int)
+            is1d = isinstance(irrepgens[0], (int,float,complex))
             irrepdim = 1 if is1d else irrepgens[0].shape[0]
             self.irreps[irrepname] = generate_grouprep(self.permgroup, 
                                                        1 if is1d else np.eye(irrepdim),
@@ -495,6 +495,24 @@ point_group_info = { # representation info of some useful point groups
                                    qpms.IRot3.identity(),
                                 )
     ),
+    'C2' : SVWFPointGroupInfo('C2',
+                 # permutation group generators
+                               (Permutation(0,1), # 180 deg rotation around z axis
+                               ), 
+                 # dictionary with irrep generators
+                               {
+                                    # Bradley, Cracknell p. 57;
+                                    'A': (1,),
+                                    'B': (-1,),    
+                               },
+                 # function that generates a tuple with svwf representation generators
+                               lambda lMax : (qpms.zrotN_tyty(2, lMax),),
+                 # quaternion rep generators
+                                rep3d_gens = (
+                                    qpms.IRot3.zrotN(2),
+                                )
+
+    ),
     'C2v' : SVWFPointGroupInfo('C2v',
                  # permutation group generators
                                (Permutation(0,1, size=4)(2,3), # x -> - x mirror operation (i.e. yz mirror plane)
@@ -546,6 +564,24 @@ point_group_info = { # representation info of some useful point groups
                                    qpms.IRot3.zflip(),
                                 )
     ),
+    'C4' : SVWFPointGroupInfo('C4',
+                 # permutation group generators
+                               (Permutation(0,1,2,3, size=4),), #C4 rotation
+                 # dictionary with irrep generators
+                               {
+                                    # Bradley, Cracknell p. 58
+                                    'A': (1,),
+                                    'B': (-1,),
+                                    '1E': (-1j,),
+                                    '2E': (1j,),
+                               },
+                 # function that generates a tuple with svwf representation generators
+                               lambda lMax : (qpms.zrotN_tyty(4, lMax), ),
+                 # quaternion rep generators
+                               rep3d_gens = (
+                                   qpms.IRot3.zrotN(4),
+                                )
+    ),
     'C4v' : SVWFPointGroupInfo('C4v',
                  # permutation group generators
                                (Permutation(0,1,2,3, size=4), #C4 rotation
@@ -621,4 +657,41 @@ point_group_info = { # representation info of some useful point groups
                                    qpms.IRot3.zflip(),
                                 )
     ),
+    'x_and_z_flip': SVWFPointGroupInfo(
+	'x_and_z_flip',
+	(
+	    Permutation(0,1, size=4), # x -> -x mirror op
+	    Permutation(2,3, size=4), # z -> -z mirror op
+	),
+	{
+	    "P'": (1, 1),
+	    "R'": (-1, 1),
+	    "P''": (-1,-1),
+	    "R''": (1, -1),
+	},
+	lambda lMax : (qpms.xflip_tyty(lMax), qpms.zflip_tyty(lMax)),
+        rep3d_gens = (
+            qpms.IRot3.xflip(),
+            qpms.IRot3.zflip(),
+        )
+
+    ),
+    'y_and_z_flip': SVWFPointGroupInfo(
+	'y_and_z_flip',
+	(
+	    Permutation(0,1, size=4), # y -> -y mirror op
+	    Permutation(2,3, size=4), # z -> -z mirror op
+	),
+	{
+	    "P'": (1, 1),
+	    "R'": (-1, 1),
+	    "P''": (-1,-1),
+	    "R''": (1, -1),
+	},
+	lambda lMax : (qpms.yflip_tyty(lMax), qpms.zflip_tyty(lMax)),
+        rep3d_gens = (
+            qpms.IRot3.yflip(),
+            qpms.IRot3.zflip(),
+        )
+    ),
 }

qpms/tmatrix_io.c

@@ -33,9 +33,9 @@ qpms_errno_t qpms_read_scuff_tmatrix(
   if (!(freqs && n && tmdata)) 
     qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
       "freqs, n, and tmdata are mandatory arguments and must not be NULL.");
-  if (bs->norm != QPMS_NORMALISATION_POWER) // CHECKME CORRECT?
+  if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT?
     qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
-        "Not implemented; only QPMS_NORMALISATION_POWER (CHECKME)"
+        "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)"
         " norm supported right now.");
   int n_alloc = 128; // First chunk to allocate
   *n = 0;

qpms/translations.c

@@ -229,8 +229,14 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm,
 #elif defined AN3
       * ipow (nu - n)
 #endif
+#ifdef AM1
+      * ipow(-m+mu)
+#endif //NNU
 #ifdef AM2
       * min1pow(-m+mu)
+#endif //NNU
+#ifdef AM3
+      * ipow(m-mu)
 #endif //NNU
   ;
   return presum * sum;
@@ -409,8 +415,14 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm,
 #elif defined AN3
       * ipow (nu - n)
 #endif
+#ifdef AM1
+      * ipow(-m+mu)
+#endif //NNU
 #ifdef AM2
       * min1pow(-m+mu)
+#endif //NNU
+#ifdef AM3
+      * ipow(m-mu)
 #endif //NNU
   ;
 
@@ -559,8 +571,14 @@ static void qpms_trans_calculator_multipliers_A_general(
 #elif defined AN3
       * ipow (nu - n)
 #endif
+#ifdef AM1
+      * ipow(-m+mu)
+#endif //NNU
 #ifdef AM2
       * min1pow(-m+mu)
+#endif //NNU
+#ifdef AM3
+      * ipow(m-mu)
 #endif //NNU
       ;
     // FIXME I might not need complex here
@@ -615,9 +633,15 @@ void qpms_trans_calculator_multipliers_B_general(
 #elif defined BN3
       * ipow (nu - n)
 #endif
+#ifdef BM1
+      * ipow(-m+mu)
+#endif  
 #ifdef BM2
       * min1pow(-m+mu)
 #endif  
+#ifdef BM3
+      * ipow(m-mu)
+#endif  
 #ifdef BF1
       * I
 #elif defined BF2

qpms/wigner.h

@@ -152,8 +152,8 @@ static inline cart3_t qpms_quat_rot_cart3(qpms_quat_t q, const cart3_t v) {
 	//const qpms_quat_t qc = qpms_quat_normalise(qpms_quat_pow(q, -1)); // implementation of _pow wrong!
 	const qpms_quat_t qc = qpms_quat_conj(q);
 	const qpms_quat_t vv = qpms_quat_from_cart3(v);
-	return qpms_quat_to_cart3(qpms_quat_mult(qc,
+	return qpms_quat_to_cart3(qpms_quat_mult(q,
-				qpms_quat_mult(vv, q)));
+				qpms_quat_mult(vv, qc)));
 }
 
 /// Wigner D matrix element from a rotator quaternion for integer \a l.

setup.py

@@ -31,11 +31,14 @@ qpms_c = Extension('qpms_c',
             'qpms/error.c'],
         extra_compile_args=['-std=c99','-ggdb', '-O0',
             '-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required
+            #'-DAN2',
+            #'-DBN2',
             #'-DQPMS_USE_OMP',
             '-DDISABLE_NDEBUG', # uncomment to enable assertions in the modules
             #'-fopenmp',
             ],
-        libraries=['gsl', 'lapacke', 'blas', 'gslcblas', #'omp'
+        libraries=['gsl', 'lapacke', 'cblas_LINUX', 'blas_LINUX', 
+            'gslcblas', #'omp'
             # TODO resolve the problem with openblas (missing gotoblas symbol) and preferable use other blas library
 	],
         runtime_library_dirs=os.environ['LD_LIBRARY_PATH'].split(':') if 'LD_LIBRARY_PATH' in os.environ else []

tests/fixproj-notes

@@ -0,0 +1,4 @@
+c99 -g -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss1.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c  ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c  ../qpms/wigner.c -lm -lgsl  -llapacke ~/repo/CBLAS/lib/cblas_LINUX.a ~/repo/BLAS-3.8.0/blas_LINUX.a
+
+
+LD_LIBRARY_PATH=$HOME/repo/CBLAS/lib/ LDFLAGS=-L../CBLAS/lib/ python3 setup.py install --user

tests/generate_groupcode.py

@@ -2,7 +2,8 @@
 from qpms.symmetries import point_group_info
 
 codestring = "#include <qpms/groups.h>\n"
-for name, info in point_group_info.items():
+for name in sorted(point_group_info.keys()):
+    info = point_group_info[name]
     codestring += 'const qpms_finite_group_t QPMS_FINITE_GROUP_%s = ' %name
     codestring += info.generate_c_source()
     codestring += ";\n\n"

tests/quaternion_O3.c

@@ -3,14 +3,18 @@
 
 int main() {
   cart3_t v = {1, 2, 3};
-  qpms_quat4d_t q4[4] = {
+  qpms_quat4d_t q4[8] = {
     {1, 0, 0, 0},
     {0, 1, 0, 0},
     {0, 0, 1, 0},
     {0, 0, 0, 1},
+    {-1, 0, 0, 0},
+    {0, -1, 0, 0},
+    {0, 0, -1, 0},
+    {0, 0, 0, -1},
   };
   printf("original: (%g, %g, %g)\n", v.x, v.y, v.z);
-  for(int i = 0; i < 4; ++i) {
+  for(int i = 0; i < 8; ++i) {
     for (int det = -1; det < 2; det += 2) {
       const qpms_quat4d_t qr = q4[i];
       const qpms_quat_t qc = qpms_quat_2c_from_4d(qr);

tests/scatsys.py

@@ -2,6 +2,8 @@ from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSyste
 from qpms.symmetries import point_group_info
 import numpy as np
 
+np.random.seed(444)
+
 sym = FinitePointGroup(point_group_info['D3h'])
 bspec2 = BaseSpec(lMax=2)
 bspec1 = BaseSpec(lMax=1)

tests/ss_syms_packs.c

@@ -35,20 +35,27 @@ int main()
   for(qpms_l_t l = 1; l <= 1; ++l)
     for (qpms_m_t m = -l; m <= l; ++m)
       qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
-  for(qpms_l_t l = 1; l <= 2; ++l)
+  for(qpms_l_t l = 1; l <= 1; ++l)
     for (qpms_m_t m = -l; m <= l; ++m)
       qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
 #endif
   qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1);
   qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2);
 
+#if 0
   // Random diagonal T-matrices
   for(size_t i = 0; i < b1->n; ++i)
     t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
   for(size_t i = 0; i < b2->n; ++i)
     t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+#else
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = 1;
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = 1;
+#endif
 
-  const cart3_t pp1 = {1, 2, 0}, pp2 = {1, 2, 3}, pp3 = {0.1, 2, 0};
+  const cart3_t pp1 = {0, 0, 1}, pp2 = {0,0, 2}, pp3 = {0,0 , 0};
   qpms_tmatrix_t * tmlist[] = {t1, t2};
   qpms_particle_tid_t plist[] = {{pp1, 0}, {pp2, 0}, {pp3, 1}};
 
@@ -58,7 +65,7 @@ int main()
   protoss.p = plist;
   protoss.p_count=3;
 
-  qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D4h);
+  qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D3h);
 
   printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n",
       (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps,

tests/sss1.c

@@ -0,0 +1,105 @@
+//  c99 -g -I.. ss_syms_packs.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c  ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c  ../qpms/wigner.c -lm -lgsl -lblas -llapacke
+typedef int qpms_gmi_t;// There is something wrong in the includes, apparently.
+#include <qpms/qpms_types.h>
+#include <qpms/scatsystem.h>
+#include <stdlib.h>
+#include <qpms/vswf.h>
+#include <qpms/indexing.h>
+#include <stdio.h>
+#include "staticgroups.h"
+
+const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h;
+const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v;
+const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g;
+const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v;
+const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h;
+const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h;
+
+double uniform_random(double min, double max) {
+  double random_value = min + (max-min)*(double)rand()/RAND_MAX;
+  return random_value;
+}
+
+int main()
+{
+  srand(666);
+#if 0
+  qpms_vswf_set_spec_t 
+    *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS),
+    *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS);
+#else
+  // Only electric waves
+  qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), 
+                       *b2 = qpms_vswf_set_spec_init();
+  b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS;
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -l; m <= l; ++m)
+      qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -l; m <= l; ++m)
+      qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+#endif
+  qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1);
+  qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2);
+
+#if 0
+  // Random diagonal T-matrices
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+#else
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = 1;
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = 1;
+#endif
+
+  const cart3_t pp1 = {1, 1, 0};
+  qpms_tmatrix_t * tmlist[] = {t1, t2};
+  qpms_particle_tid_t plist[] = {{pp1, 0}};
+
+  qpms_scatsys_t protoss;
+  protoss.tm = tmlist;
+  protoss.tm_count=2;
+  protoss.p = plist;
+  protoss.p_count=1;
+
+  qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D2h);
+
+  printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n",
+      (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps,
+      (int)ss->orbit_type_count);
+
+  const double k = 1.7;
+
+  complex double *S_full = qpms_scatsys_build_translation_matrix_full(
+      NULL, ss, k); 
+  complex double *S_packed[ss->sym->nirreps];
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri)
+    S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL,
+        S_full, ss, iri);
+
+  complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL,
+      S_packed[0], ss, 0, false);
+  for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri)
+    qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri],
+        ss, iri, true);
+
+  double maxerr = 0;
+  for (size_t i = 0; i < ss->fecv_size; ++i) {
+    double err = cabs(S_full[i] - S_recfull[i]);
+    maxerr = (err > maxerr) ? err : maxerr;
+  }
+
+  printf("maxerr: %lg\n", maxerr);
+  
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]);
+  free(S_full);
+  qpms_scatsys_free(ss);
+  qpms_tmatrix_free(t1);
+  qpms_tmatrix_free(t2);
+  qpms_vswf_set_spec_free(b1);
+  qpms_vswf_set_spec_free(b2);
+  return 0;
+}

tests/sss2.c

@@ -0,0 +1,180 @@
+//  c99 -g -DZLINE -DDAGRUP=D3h -DDUMP_PARTICLE_POSITIONS -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX  -I.. sss2.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c  ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c  ../qpms/wigner.c -lm -lgsl -lblas -llapacke
+typedef int qpms_gmi_t;// There is something wrong in the includes, apparently.
+#include <qpms/qpms_types.h>
+#include <qpms/scatsystem.h>
+#include <stdlib.h>
+#include <qpms/vswf.h>
+#include <qpms/indexing.h>
+#include <stdio.h>
+#include "staticgroups.h"
+
+const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h;
+const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v;
+const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g;
+const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v;
+const qpms_finite_group_t *C2 = &QPMS_FINITE_GROUP_C2;
+const qpms_finite_group_t *C4 = &QPMS_FINITE_GROUP_C4;
+const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h;
+const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h;
+const qpms_finite_group_t *x_and_z_flip = &QPMS_FINITE_GROUP_x_and_z_flip;
+const qpms_finite_group_t *y_and_z_flip = &QPMS_FINITE_GROUP_y_and_z_flip;
+
+#ifndef DAGRUP
+#define DAGRUP D4h
+#endif
+
+double uniform_random(double min, double max) {
+  double random_value = min + (max-min)*(double)rand()/RAND_MAX;
+  return random_value;
+}
+
+int main()
+{
+  srand(666);
+#if 0
+  qpms_vswf_set_spec_t 
+    *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS),
+    *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS);
+#else
+  // Only electric waves
+  qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), 
+                       *b2 = qpms_vswf_set_spec_init();
+  b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS;
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -0l; m <= l; m += 2)
+      qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -0l; m <= l; m += 2)
+      qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+#endif
+  qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1);
+  qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2);
+
+#if 0
+  // Random diagonal T-matrices
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+#else
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = 1;
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = 1;
+#endif
+
+#ifdef YLINE
+  const cart3_t pp1 = {0, 1.1, 0};
+  const cart3_t pp2 = {0, 1.4, 0};
+#elif defined XLINE
+  const cart3_t pp1 = {1.1, 0, 0};
+  const cart3_t pp2 = {1.4, 0, 0};
+#elif defined ZLINE
+  const cart3_t pp1 = {0, 0, 1.1};
+  const cart3_t pp2 = {0, 0, 1.4};
+#else
+  const cart3_t pp1 = {1.1, 1, 0};
+  const cart3_t pp2 = {0, 1.4, 0};
+#endif
+  const cart3_t pp3 = {0, 0, 1};
+  qpms_tmatrix_t * tmlist[] = {t1, t2};
+  qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0},  {pp3, 1},
+  };
+
+  qpms_scatsys_t protoss;
+  protoss.tm = tmlist;
+  protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *);
+  protoss.p = plist;
+  protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t);
+
+  qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, DAGRUP);
+
+  printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n",
+      (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps,
+      (int)ss->orbit_type_count);
+
+  const double k = 1.7;
+
+  complex double *S_full = qpms_scatsys_build_translation_matrix_full(
+      NULL, ss, k); 
+  {
+    const size_t full_len = ss->fecv_size;
+    size_t fullvec_offset_dest = 0;
+    for (qpms_ss_pi_t pdest = 0; pdest < ss->p_count; pdest++) {
+      size_t fullvec_offset_src = 0;
+      const size_t bspecn_dest = ss->tm[ss->p[pdest].tmatrix_id]->spec->n;
+      for (qpms_ss_pi_t psrc = 0; psrc < ss->p_count; psrc++) {
+        const size_t bspecn_src = ss->tm[ss->p[psrc].tmatrix_id]->spec->n;
+        fprintf(stderr, "Translation matrix element %d<-%d; (%g %g %g)<-(%g %g %g):\n",
+            (int)pdest, (int)psrc, ss->p[pdest].pos.x, ss->p[pdest].pos.y, ss->p[pdest].pos.z,
+            ss->p[psrc].pos.x, ss->p[psrc].pos.y, ss->p[psrc].pos.z);
+
+        for(size_t row = 0; row < bspecn_dest; ++row) {
+          for(size_t col = 0; col < bspecn_src; ++col)
+            fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]),
+                cimag(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]));
+          fputc('\n', stderr);
+        }
+        fullvec_offset_src += bspecn_src;
+      }
+      fullvec_offset_dest += bspecn_dest;
+    }
+  }
+  {
+    fputs("\n\n", stderr);
+    const size_t full_len = ss->fecv_size;
+    for (size_t row = 0 ; row < full_len; ++row) {
+      for (size_t col = 0 ; col < full_len; ++col)
+        fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * row + col]), cimag(S_full[full_len * row + col]));
+      fputc('\n', stderr);
+    }
+  }
+
+  complex double *S_packed[ss->sym->nirreps];
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri)
+    S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL,
+        S_full, ss, iri);
+
+  complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL,
+      S_packed[0], ss, 0, false);
+  for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri)
+    qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri],
+        ss, iri, true);
+  {
+    fputs("\n\n", stderr);
+    const size_t full_len = ss->fecv_size;
+    for (size_t row = 0 ; row < full_len; ++row) {
+      for (size_t col = 0 ; col < full_len; ++col)
+        fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_recfull[full_len * row + col]), cimag(S_recfull[full_len * row + col]));
+      fputc('\n', stderr);
+    }
+  }
+
+  double maxerr = 0;
+  for (size_t i = 0; i < ss->fecv_size; ++i) {
+    double err = cabs(S_full[i] - S_recfull[i]);
+    maxerr = (err > maxerr) ? err : maxerr;
+  }
+
+
+  printf("maxerr: %lg\n", maxerr);
+  
+  fprintf(stderr, "pi\tpos\toti\tosn\tp\n");
+  for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
+    cart3_t pos = ss->p[pi].pos;
+    qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t;
+    qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn;
+    qpms_ss_orbit_pi_t p = ss->p_orbitinfo[pi].p;
+    fprintf(stderr, "%d\t(%.3g,%.3g,%.3g)\t%d\t%d\t%d\n",
+        (int)pi, pos.x, pos.y, pos.z, (int)oti, (int)osn, (int)p);
+  }
+
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]);
+  free(S_full);
+  qpms_scatsys_free(ss);
+  qpms_tmatrix_free(t1);
+  qpms_tmatrix_free(t2);
+  qpms_vswf_set_spec_free(b1);
+  qpms_vswf_set_spec_free(b2);
+  return 0;
+}

tests/sss3.c

@@ -0,0 +1,236 @@
+// c99 -g -DNLINE -DDAGRUP=C4v -DDUMP_PARTICLE_POSITIONS -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss3.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c  ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c  ../qpms/wigner.c -lm -lgsl  -llapacke ~/repo/CBLAS/lib/cblas_LINUX.a ~/repo/BLAS-3.8.0/blas_LINUX.a
+typedef int qpms_gmi_t;// There is something wrong in the includes, apparently.
+#include <qpms/qpms_types.h>
+#include <qpms/scatsystem.h>
+#include <stdlib.h>
+#include <qpms/vswf.h>
+#include <qpms/indexing.h>
+#include <stdio.h>
+#include "staticgroups.h"
+
+const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h;
+const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v;
+const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g;
+const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v;
+const qpms_finite_group_t *C2 = &QPMS_FINITE_GROUP_C2;
+const qpms_finite_group_t *C4 = &QPMS_FINITE_GROUP_C4;
+const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h;
+const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h;
+const qpms_finite_group_t *x_and_z_flip = &QPMS_FINITE_GROUP_x_and_z_flip;
+const qpms_finite_group_t *y_and_z_flip = &QPMS_FINITE_GROUP_y_and_z_flip;
+
+#ifndef DAGRUP
+#define DAGRUP D4h
+#endif
+
+double uniform_random(double min, double max) {
+  double random_value = min + (max-min)*(double)rand()/RAND_MAX;
+  return random_value;
+}
+
+int main()
+{
+  srand(666);
+#if 0
+  qpms_vswf_set_spec_t 
+    *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS),
+    *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS);
+#else
+  // Only electric waves
+  qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), 
+                       *b2 = qpms_vswf_set_spec_init();
+  b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS;
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -0l; m <= l; m += 2)
+      qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -0l; m <= l; m += 2)
+      qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+#endif
+  qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1);
+  qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2);
+
+#if 0
+  // Random diagonal T-matrices
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+#else
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = 1;
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = 1;
+#endif
+
+#ifdef YLINE
+  const cart3_t pp1 = {0, 1.1, 0};
+  const cart3_t pp2 = {0, 1.4, 0};
+#elif defined XLINE
+  const cart3_t pp1 = {1.1, 0, 0};
+  const cart3_t pp2 = {1.4, 0, 0};
+#elif defined ZLINE
+  const cart3_t pp1 = {0, 0, 1.1};
+  const cart3_t pp2 = {0, 0, 1.4};
+#else
+  const cart3_t pp1 = {1.1, 1, 0};
+  const cart3_t pp2 = {0, 1.4, 0};
+#endif
+  const cart3_t pp3 = {0, 0, 1};
+  qpms_tmatrix_t * tmlist[] = {t1, t2};
+  qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0},  {pp3, 1},
+  };
+
+  qpms_scatsys_t protoss;
+  protoss.tm = tmlist;
+  protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *);
+  protoss.p = plist;
+  protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t);
+
+  qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, DAGRUP);
+
+  printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n",
+      (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps,
+      (int)ss->orbit_type_count);
+
+  fputs("Orbit projection matrices:\n", stderr);
+  for (qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) {
+    fprintf(stderr, "Orbit type %d:\n", (int)oti);
+    const qpms_ss_orbit_type_t *ot = &(ss->orbit_types[oti]);
+    size_t row = 0;
+    for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
+      assert(row*ot->size*ot->bspecn == ot->irbase_offsets[iri]);
+      fprintf(stderr, "---------- IR %d (%s) -------------\n", (int)iri, ss->sym->irreps[iri].name);
+      for (size_t irbi = 0; irbi < ot->irbase_sizes[iri]; ++irbi) {
+        for(qpms_ss_orbit_pi_t opi = 0; opi < ot->size; ++opi){
+          fputs("| ", stderr);
+          for(size_t i = 0; i < ot->bspecn; ++i) {
+            const complex double elem = ot->irbases[row * ot->size * ot->bspecn
+            + opi * ot->bspecn + i];
+            fprintf(stderr, "%+.3f%+.3fj ", creal(elem), cimag(elem));
+          }
+        }
+        fputs("|\n", stderr);
+        ++row;
+      }
+    }
+    fputs("------------------------\n\n",stderr);
+  }
+
+  const double k = 1.7;
+
+  complex double *S_full = qpms_scatsys_build_translation_matrix_full(
+      NULL, ss, k); 
+  {
+    const size_t full_len = ss->fecv_size;
+    size_t fullvec_offset_dest = 0;
+    for (qpms_ss_pi_t pdest = 0; pdest < ss->p_count; pdest++) {
+      size_t fullvec_offset_src = 0;
+      const size_t bspecn_dest = ss->tm[ss->p[pdest].tmatrix_id]->spec->n;
+      for (qpms_ss_pi_t psrc = 0; psrc < ss->p_count; psrc++) {
+        const size_t bspecn_src = ss->tm[ss->p[psrc].tmatrix_id]->spec->n;
+        fprintf(stderr, "Translation matrix element %d<-%d; (%g %g %g)<-(%g %g %g):\n",
+            (int)pdest, (int)psrc, ss->p[pdest].pos.x, ss->p[pdest].pos.y, ss->p[pdest].pos.z,
+            ss->p[psrc].pos.x, ss->p[psrc].pos.y, ss->p[psrc].pos.z);
+
+        for(size_t row = 0; row < bspecn_dest; ++row) {
+          for(size_t col = 0; col < bspecn_src; ++col)
+            fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]),
+                cimag(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]));
+          fputc('\n', stderr);
+        }
+        fullvec_offset_src += bspecn_src;
+      }
+      fullvec_offset_dest += bspecn_dest;
+    }
+  }
+  {
+    fputs("\n\n", stderr);
+    const size_t full_len = ss->fecv_size;
+    for (size_t row = 0 ; row < full_len; ++row) {
+      for (size_t col = 0 ; col < full_len; ++col)
+        fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * row + col]), cimag(S_full[full_len * row + col]));
+      fputc('\n', stderr);
+    }
+  }
+
+  complex double *S_packed[ss->sym->nirreps];
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
+    S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL,
+        S_full, ss, iri);
+    fprintf(stderr, "--- Packed matrix for irrep %d (%s):\n", (int) iri, ss->sym->irreps[iri].name);
+    for (size_t row = 0; row < ss->saecv_sizes[iri]; ++row) {
+      for (size_t col = 0; col < ss->saecv_sizes[iri]; ++col) {
+        complex double elem = S_packed[iri][row * ss->saecv_sizes[iri] + col];
+        fprintf(stderr, "%+.3f+%.3fj ", creal(elem), cimag(elem));
+      }
+      fputc('\n', stderr);
+    }
+  }
+  {
+    complex double *S_partrecfull = qpms_scatsys_irrep_unpack_matrix(NULL,
+        S_packed[0], ss, 0, false);
+    for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
+      qpms_scatsys_irrep_unpack_matrix(S_partrecfull, S_packed[iri],
+          ss, iri, false);
+      fprintf(stderr, "\nPartial reconstruction %d (%s):\n", (int)iri, ss->sym->irreps[iri].name);
+      const size_t full_len = ss->fecv_size;
+      for (size_t row = 0 ; row < full_len; ++row) {
+        for (size_t col = 0 ; col < full_len; ++col)
+          fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_partrecfull[full_len * row + col]), cimag(S_partrecfull[full_len * row + col]));
+        fputc('\n', stderr);
+      }
+    }
+
+    double maxerr = 0;
+    for (size_t i = 0; i < ss->fecv_size; ++i) {
+      double err = cabs(S_full[i] - S_partrecfull[i]);
+      maxerr = (err > maxerr) ? err : maxerr;
+    }
+    free(S_partrecfull);
+  }
+
+
+  complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL,
+      S_packed[0], ss, 0, false);
+  for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri)
+    qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri],
+        ss, iri, true);
+  {
+    fputs("\n\n", stderr);
+    const size_t full_len = ss->fecv_size;
+    for (size_t row = 0 ; row < full_len; ++row) {
+      for (size_t col = 0 ; col < full_len; ++col)
+        fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_recfull[full_len * row + col]), cimag(S_recfull[full_len * row + col]));
+      fputc('\n', stderr);
+    }
+  }
+
+  double maxerr = 0;
+  for (size_t i = 0; i < ss->fecv_size; ++i) {
+    double err = cabs(S_full[i] - S_recfull[i]);
+    maxerr = (err > maxerr) ? err : maxerr;
+  }
+
+
+  printf("maxerr: %lg\n", maxerr);
+  
+  fprintf(stderr, "pi\tpos\toti\tosn\tp\n");
+  for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
+    cart3_t pos = ss->p[pi].pos;
+    qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t;
+    qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn;
+    qpms_ss_orbit_pi_t p = ss->p_orbitinfo[pi].p;
+    fprintf(stderr, "%d\t(%.3g,%.3g,%.3g)\t%d\t%d\t%d\n",
+        (int)pi, pos.x, pos.y, pos.z, (int)oti, (int)osn, (int)p);
+  }
+
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]);
+  free(S_full);
+  qpms_scatsys_free(ss);
+  qpms_tmatrix_free(t1);
+  qpms_tmatrix_free(t2);
+  qpms_vswf_set_spec_free(b1);
+  qpms_vswf_set_spec_free(b2);
+  return 0;
+}

tests/sss4.c

@@ -0,0 +1,237 @@
+// c99 -g -DNLINE -DDAGRUP=C4v -DDUMP_PARTICLE_POSITIONS -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss3.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c  ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c  ../qpms/wigner.c -lm -lgsl  -llapacke ~/repo/CBLAS/lib/cblas_LINUX.a ~/repo/BLAS-3.8.0/blas_LINUX.a
+typedef int qpms_gmi_t;// There is something wrong in the includes, apparently.
+#include <qpms/qpms_types.h>
+#include <qpms/scatsystem.h>
+#include <stdlib.h>
+#include <qpms/vswf.h>
+#include <qpms/indexing.h>
+#include <stdio.h>
+#include "staticgroups.h"
+
+const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h;
+const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v;
+const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g;
+const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v;
+const qpms_finite_group_t *C2 = &QPMS_FINITE_GROUP_C2;
+const qpms_finite_group_t *C4 = &QPMS_FINITE_GROUP_C4;
+const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h;
+const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h;
+const qpms_finite_group_t *x_and_z_flip = &QPMS_FINITE_GROUP_x_and_z_flip;
+const qpms_finite_group_t *y_and_z_flip = &QPMS_FINITE_GROUP_y_and_z_flip;
+
+#ifndef DAGRUP
+#define DAGRUP D4h
+#endif
+
+double uniform_random(double min, double max) {
+  double random_value = min + (max-min)*(double)rand()/RAND_MAX;
+  return random_value;
+}
+
+int main()
+{
+  srand(666);
+#if 1
+  qpms_vswf_set_spec_t 
+    *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS),
+    *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS);
+#else
+  // Only electric waves
+  qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), 
+                       *b2 = qpms_vswf_set_spec_init();
+  b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS;
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -l; m <= l; m += 1)
+      qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+  for(qpms_l_t l = 1; l <= 1; ++l)
+    for (qpms_m_t m = -l; m <= l; m += 1)
+      qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
+#endif
+  qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1);
+  qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2);
+
+#if 0
+  // Random diagonal T-matrices
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
+#else
+  for(size_t i = 0; i < b1->n; ++i)
+    t1->m[i + i*b1->n] = 1;
+  for(size_t i = 0; i < b2->n; ++i)
+    t2->m[i + i*b2->n] = 1;
+#endif
+
+#ifdef YLINE
+  const cart3_t pp1 = {0, 1.1, 0};
+  const cart3_t pp2 = {0, 1.4, 0};
+#elif defined XLINE
+  const cart3_t pp1 = {1.1, 0, 0};
+  const cart3_t pp2 = {1.4, 0, 0};
+#elif defined ZLINE
+  const cart3_t pp1 = {0, 0, 1.1};
+  const cart3_t pp2 = {0, 0, 1.4};
+#else
+  const cart3_t pp1 = {1.1, 1, 0};
+  const cart3_t pp2 = {0, 1.4, 0};
+  const cart3_t pp4 = {2, 1.4, 0};
+#endif
+  const cart3_t pp3 = {0, 0, 0};
+  qpms_tmatrix_t * tmlist[] = {t1, t2};
+  qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0},  {pp3, 1}, {pp4,1}
+  };
+
+  qpms_scatsys_t protoss;
+  protoss.tm = tmlist;
+  protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *);
+  protoss.p = plist;
+  protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t);
+
+  qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, DAGRUP);
+
+  printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n",
+      (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps,
+      (int)ss->orbit_type_count);
+
+  fputs("Orbit projection matrices:\n", stderr);
+  for (qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) {
+    fprintf(stderr, "Orbit type %d:\n", (int)oti);
+    const qpms_ss_orbit_type_t *ot = &(ss->orbit_types[oti]);
+    size_t row = 0;
+    for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
+      assert(row*ot->size*ot->bspecn == ot->irbase_offsets[iri]);
+      fprintf(stderr, "---------- IR %d (%s) -------------\n", (int)iri, ss->sym->irreps[iri].name);
+      for (size_t irbi = 0; irbi < ot->irbase_sizes[iri]; ++irbi) {
+        for(qpms_ss_orbit_pi_t opi = 0; opi < ot->size; ++opi){
+          fputs("| ", stderr);
+          for(size_t i = 0; i < ot->bspecn; ++i) {
+            const complex double elem = ot->irbases[row * ot->size * ot->bspecn
+            + opi * ot->bspecn + i];
+            fprintf(stderr, "%+.3f%+.3fj ", creal(elem), cimag(elem));
+          }
+        }
+        fputs("|\n", stderr);
+        ++row;
+      }
+    }
+    fputs("------------------------\n\n",stderr);
+  }
+
+  const double k = 1.7;
+
+  complex double *S_full = qpms_scatsys_build_translation_matrix_full(
+      NULL, ss, k); 
+  {
+    const size_t full_len = ss->fecv_size;
+    size_t fullvec_offset_dest = 0;
+    for (qpms_ss_pi_t pdest = 0; pdest < ss->p_count; pdest++) {
+      size_t fullvec_offset_src = 0;
+      const size_t bspecn_dest = ss->tm[ss->p[pdest].tmatrix_id]->spec->n;
+      for (qpms_ss_pi_t psrc = 0; psrc < ss->p_count; psrc++) {
+        const size_t bspecn_src = ss->tm[ss->p[psrc].tmatrix_id]->spec->n;
+        fprintf(stderr, "Translation matrix element %d<-%d; (%g %g %g)<-(%g %g %g):\n",
+            (int)pdest, (int)psrc, ss->p[pdest].pos.x, ss->p[pdest].pos.y, ss->p[pdest].pos.z,
+            ss->p[psrc].pos.x, ss->p[psrc].pos.y, ss->p[psrc].pos.z);
+
+        for(size_t row = 0; row < bspecn_dest; ++row) {
+          for(size_t col = 0; col < bspecn_src; ++col)
+            fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]),
+                cimag(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]));
+          fputc('\n', stderr);
+        }
+        fullvec_offset_src += bspecn_src;
+      }
+      fullvec_offset_dest += bspecn_dest;
+    }
+  }
+  {
+    fputs("\n\n", stderr);
+    const size_t full_len = ss->fecv_size;
+    for (size_t row = 0 ; row < full_len; ++row) {
+      for (size_t col = 0 ; col < full_len; ++col)
+        fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * row + col]), cimag(S_full[full_len * row + col]));
+      fputc('\n', stderr);
+    }
+  }
+
+  complex double *S_packed[ss->sym->nirreps];
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
+    S_packed[iri] = qpms_scatsys_irrep_pack_matrix_stupid(NULL,
+        S_full, ss, iri);
+    fprintf(stderr, "--- Packed matrix for irrep %d (%s):\n", (int) iri, ss->sym->irreps[iri].name);
+    for (size_t row = 0; row < ss->saecv_sizes[iri]; ++row) {
+      for (size_t col = 0; col < ss->saecv_sizes[iri]; ++col) {
+        complex double elem = S_packed[iri][row * ss->saecv_sizes[iri] + col];
+        fprintf(stderr, "%+.3f+%.3fj ", creal(elem), cimag(elem));
+      }
+      fputc('\n', stderr);
+    }
+  }
+  {
+    complex double *S_partrecfull = qpms_scatsys_irrep_unpack_matrix_stupid(NULL,
+        S_packed[0], ss, 0, false);
+    for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
+      qpms_scatsys_irrep_unpack_matrix_stupid(S_partrecfull, S_packed[iri],
+          ss, iri, false);
+      fprintf(stderr, "\nPartial reconstruction %d (%s):\n", (int)iri, ss->sym->irreps[iri].name);
+      const size_t full_len = ss->fecv_size;
+      for (size_t row = 0 ; row < full_len; ++row) {
+        for (size_t col = 0 ; col < full_len; ++col)
+          fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_partrecfull[full_len * row + col]), cimag(S_partrecfull[full_len * row + col]));
+        fputc('\n', stderr);
+      }
+    }
+
+    double maxerr = 0;
+    for (size_t i = 0; i < ss->fecv_size; ++i) {
+      double err = cabs(S_full[i] - S_partrecfull[i]);
+      maxerr = (err > maxerr) ? err : maxerr;
+    }
+    free(S_partrecfull);
+  }
+
+
+  complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix_stupid(NULL,
+      S_packed[0], ss, 0, false);
+  for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri)
+    qpms_scatsys_irrep_unpack_matrix_stupid(S_recfull, S_packed[iri],
+        ss, iri, true);
+  {
+    fputs("\n\n", stderr);
+    const size_t full_len = ss->fecv_size;
+    for (size_t row = 0 ; row < full_len; ++row) {
+      for (size_t col = 0 ; col < full_len; ++col)
+        fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_recfull[full_len * row + col]), cimag(S_recfull[full_len * row + col]));
+      fputc('\n', stderr);
+    }
+  }
+
+  double maxerr = 0;
+  for (size_t i = 0; i < ss->fecv_size; ++i) {
+    double err = cabs(S_full[i] - S_recfull[i]);
+    maxerr = (err > maxerr) ? err : maxerr;
+  }
+
+
+  printf("maxerr: %lg\n", maxerr);
+  
+  fprintf(stderr, "pi\tpos\toti\tosn\tp\n");
+  for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
+    cart3_t pos = ss->p[pi].pos;
+    qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t;
+    qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn;
+    qpms_ss_orbit_pi_t p = ss->p_orbitinfo[pi].p;
+    fprintf(stderr, "%d\t(%.3g,%.3g,%.3g)\t%d\t%d\t%d\n",
+        (int)pi, pos.x, pos.y, pos.z, (int)oti, (int)osn, (int)p);
+  }
+
+  for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]);
+  free(S_full);
+  qpms_scatsys_free(ss);
+  qpms_tmatrix_free(t1);
+  qpms_tmatrix_free(t2);
+  qpms_vswf_set_spec_free(b1);
+  qpms_vswf_set_spec_free(b2);
+  return 0;
+}

tests/staticgroups.c

@@ -1,4 +1,42 @@
-#include "staticgroups.h"
+#include <qpms/groups.h>
+const qpms_finite_group_t QPMS_FINITE_GROUP_C2 = {
+  "C2", // name
+  2, // order
+  0, // idi
+  (qpms_gmi_t[]) { // mt
+    0, 1,
+    1, 0,
+  },
+  (qpms_gmi_t[]) { // invi
+    0, 1
+  },
+  (qpms_gmi_t[]) {1}, // gens
+  1, // ngens
+  (qpms_permutation_t[]){ // permrep
+    "(1)",
+    "(0 1)",
+  },
+  NULL, // elemlabels
+  2, // permrep_nelem
+  (qpms_irot3_t[]) { // rep3d
+   {{1.0+0.0*I, 0.0+0.0*I}, 1},
+   {{6.123233995736766e-17+1.0*I, 0.0+0.0*I}, 1},
+  },
+  2, // nirreps
+  (struct qpms_finite_group_irrep_t[]) { // irreps
+    {
+      1, // dim
+      "B", //name
+      (complex double []) {1, -1} // m
+    },
+    {
+      1, // dim
+      "A", //name
+      (complex double []) {1, 1} // m
+    },
+  } // end of irreps
+};
+
 const qpms_finite_group_t QPMS_FINITE_GROUP_C2v = {
   "C2v", // name
   4, // order
@@ -35,6 +73,11 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_C2v = {
       "A1", //name
       (complex double []) {1, 1, 1, 1} // m
     },
+    {
+      1, // dim
+      "A2", //name
+      (complex double []) {1, -1, 1, -1} // m
+    },
     {
       1, // dim
       "B2", //name
@@ -45,10 +88,241 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_C2v = {
       "B1", //name
       (complex double []) {1, 1, -1, -1} // m
     },
+  } // end of irreps
+};
+
+const qpms_finite_group_t QPMS_FINITE_GROUP_C4 = {
+  "C4", // name
+  4, // order
+  0, // idi
+  (qpms_gmi_t[]) { // mt
+    0, 1, 2, 3,
+    1, 2, 3, 0,
+    2, 3, 0, 1,
+    3, 0, 1, 2,
+  },
+  (qpms_gmi_t[]) { // invi
+    0, 3, 2, 1
+  },
+  (qpms_gmi_t[]) {1}, // gens
+  1, // ngens
+  (qpms_permutation_t[]){ // permrep
+    "(3)",
+    "(0 1 2 3)",
+    "(0 2)(1 3)",
+    "(0 3 2 1)",
+  },
+  NULL, // elemlabels
+  4, // permrep_nelem
+  (qpms_irot3_t[]) { // rep3d
+   {{1.0+0.0*I, 0.0+0.0*I}, 1},
+   {{0.7071067811865476+0.7071067811865475*I, 0.0+0.0*I}, 1},
+   {{2.220446049250313e-16+1.0*I, 0.0+0.0*I}, 1},
+   {{-0.7071067811865474+0.7071067811865477*I, 0.0+0.0*I}, 1},
+  },
+  4, // nirreps
+  (struct qpms_finite_group_irrep_t[]) { // irreps
+    {
+      1, // dim
+      "B", //name
+      (complex double []) {1, -1, 1, -1} // m
+    },
+    {
+      1, // dim
+      "2E", //name
+      (complex double []) {1, 1j, (-1+0j), (-0-1j)} // m
+    },
+    {
+      1, // dim
+      "A", //name
+      (complex double []) {1, 1, 1, 1} // m
+    },
+    {
+      1, // dim
+      "1E", //name
+      (complex double []) {1, -1j, (-1+0j), 1j} // m
+    },
+  } // end of irreps
+};
+
+const qpms_finite_group_t QPMS_FINITE_GROUP_C4v = {
+  "C4v", // name
+  8, // order
+  0, // idi
+  (qpms_gmi_t[]) { // mt
+    0, 1, 2, 3, 4, 5, 6, 7,
+    1, 2, 3, 0, 7, 4, 5, 6,
+    2, 3, 0, 1, 6, 7, 4, 5,
+    3, 0, 1, 2, 5, 6, 7, 4,
+    4, 5, 6, 7, 0, 1, 2, 3,
+    5, 6, 7, 4, 3, 0, 1, 2,
+    6, 7, 4, 5, 2, 3, 0, 1,
+    7, 4, 5, 6, 1, 2, 3, 0,
+  },
+  (qpms_gmi_t[]) { // invi
+    0, 3, 2, 1, 4, 5, 6, 7
+  },
+  (qpms_gmi_t[]) {1, 7}, // gens
+  2, // ngens
+  (qpms_permutation_t[]){ // permrep
+    "(3)",
+    "(0 1 2 3)",
+    "(0 2)(1 3)",
+    "(0 3 2 1)",
+    "(3)(0 2)",
+    "(0 3)(1 2)",
+    "(1 3)",
+    "(0 1)(2 3)",
+  },
+  NULL, // elemlabels
+  4, // permrep_nelem
+  (qpms_irot3_t[]) { // rep3d
+   {{1.0+0.0*I, 0.0+0.0*I}, 1},
+   {{0.7071067811865476+0.7071067811865475*I, 0.0+0.0*I}, 1},
+   {{2.220446049250313e-16+1.0*I, 0.0+0.0*I}, 1},
+   {{-0.7071067811865474+0.7071067811865477*I, 0.0+0.0*I}, 1},
+   {{0.0+0.0*I, 0.7071067811865477-0.7071067811865474*I}, -1},
+   {{0.0+0.0*I, 1.0+2.220446049250313e-16*I}, -1},
+   {{0.0+0.0*I, 0.7071067811865475+0.7071067811865476*I}, -1},
+   {{0.0+0.0*I, 0.0+1.0*I}, -1},
+  },
+  5, // nirreps
+  (struct qpms_finite_group_irrep_t[]) { // irreps
+    {
+      1, // dim
+      "A1", //name
+      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m
+    },
     {
       1, // dim
       "A2", //name
-      (complex double []) {1, -1, 1, -1} // m
+      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m
+    },
+    {
+      2, // dim
+      "E", //name
+      (complex double []) {
+        // (3)
+        1.0, 0.0, 
+        0.0, 1.0, 
+        // (0 1 2 3)
+        0.0, -1.0, 
+        1.0, 0.0, 
+        // (0 2)(1 3)
+        -1.0, 0.0, 
+        0.0, -1.0, 
+        // (0 3 2 1)
+        0.0, 1.0, 
+        -1.0, 0.0, 
+        // (3)(0 2)
+        0.0, 1.0, 
+        1.0, 0.0, 
+        // (0 3)(1 2)
+        1.0, 0.0, 
+        0.0, -1.0, 
+        // (1 3)
+        0.0, -1.0, 
+        -1.0, 0.0, 
+        // (0 1)(2 3)
+        -1.0, 0.0, 
+        0.0, 1.0, 
+      }
+    },
+    {
+      1, // dim
+      "B2", //name
+      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m
+    },
+    {
+      1, // dim
+      "B1", //name
+      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m
+    },
+  } // end of irreps
+};
+
+const qpms_finite_group_t QPMS_FINITE_GROUP_D2h = {
+  "D2h", // name
+  8, // order
+  0, // idi
+  (qpms_gmi_t[]) { // mt
+    0, 1, 2, 3, 4, 5, 6, 7,
+    1, 0, 3, 2, 5, 4, 7, 6,
+    2, 3, 0, 1, 6, 7, 4, 5,
+    3, 2, 1, 0, 7, 6, 5, 4,
+    4, 5, 6, 7, 0, 1, 2, 3,
+    5, 4, 7, 6, 1, 0, 3, 2,
+    6, 7, 4, 5, 2, 3, 0, 1,
+    7, 6, 5, 4, 3, 2, 1, 0,
+  },
+  (qpms_gmi_t[]) { // invi
+    0, 1, 2, 3, 4, 5, 6, 7
+  },
+  (qpms_gmi_t[]) {1, 3, 7}, // gens
+  3, // ngens
+  (qpms_permutation_t[]){ // permrep
+    "(5)",
+    "(5)(0 1)(2 3)",
+    "(5)(0 2)(1 3)",
+    "(5)(0 3)(1 2)",
+    "(0 3)(1 2)(4 5)",
+    "(0 2)(1 3)(4 5)",
+    "(0 1)(2 3)(4 5)",
+    "(4 5)",
+  },
+  NULL, // elemlabels
+  6, // permrep_nelem
+  (qpms_irot3_t[]) { // rep3d
+   {{1.0+0.0*I, 0.0+0.0*I}, 1},
+   {{0.0+0.0*I, 0.0+1.0*I}, -1},
+   {{0.0+1.0*I, 0.0+0.0*I}, 1},
+   {{0.0+0.0*I, 1.0+0.0*I}, -1},
+   {{0.0+0.0*I, 0.0+1.0*I}, 1},
+   {{-1.0+0.0*I, 0.0+0.0*I}, -1},
+   {{0.0+0.0*I, -1.0+0.0*I}, 1},
+   {{0.0+1.0*I, 0.0+0.0*I}, -1},
+  },
+  8, // nirreps
+  (struct qpms_finite_group_irrep_t[]) { // irreps
+    {
+      1, // dim
+      "A2\'", //name
+      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m
+    },
+    {
+      1, // dim
+      "B1\'", //name
+      (complex double []) {1, 1, -1, -1, -1, -1, 1, 1} // m
+    },
+    {
+      1, // dim
+      "A2\'\'", //name
+      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m
+    },
+    {
+      1, // dim
+      "B2\'", //name
+      (complex double []) {1, -1, -1, 1, 1, -1, -1, 1} // m
+    },
+    {
+      1, // dim
+      "A1\'", //name
+      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m
+    },
+    {
+      1, // dim
+      "A1\'\'", //name
+      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m
+    },
+    {
+      1, // dim
+      "B2\'\'", //name
+      (complex double []) {1, 1, -1, -1, 1, 1, -1, -1} // m
+    },
+    {
+      1, // dim
+      "B1\'\'", //name
+      (complex double []) {1, -1, -1, 1, -1, 1, 1, -1} // m
     },
   } // end of irreps
 };
@@ -108,6 +382,11 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D3h = {
   },
   6, // nirreps
   (struct qpms_finite_group_irrep_t[]) { // irreps
+    {
+      1, // dim
+      "A2\'", //name
+      (complex double []) {1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1} // m
+    },
     {
       2, // dim
       "E\'", //name
@@ -150,6 +429,21 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D3h = {
         0.0, 0.9999999999999996, 
       }
     },
+    {
+      1, // dim
+      "A2\'\'", //name
+      (complex double []) {1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1} // m
+    },
+    {
+      1, // dim
+      "A1\'", //name
+      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m
+    },
+    {
+      1, // dim
+      "A1\'\'", //name
+      (complex double []) {1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1} // m
+    },
     {
       2, // dim
       "E\'\'", //name
@@ -192,26 +486,6 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D3h = {
         0.0, -0.9999999999999996, 
       }
     },
-    {
-      1, // dim
-      "A2\'", //name
-      (complex double []) {1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1} // m
-    },
-    {
-      1, // dim
-      "A1\'\'", //name
-      (complex double []) {1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1} // m
-    },
-    {
-      1, // dim
-      "A1\'", //name
-      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m
-    },
-    {
-      1, // dim
-      "A2\'\'", //name
-      (complex double []) {1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1} // m
-    },
   } // end of irreps
 };
 
@@ -284,8 +558,13 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D4h = {
   (struct qpms_finite_group_irrep_t[]) { // irreps
     {
       1, // dim
-      "B1\'\'", //name
+      "A2\'", //name
-      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1} // m
+      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1} // m
+    },
+    {
+      1, // dim
+      "B1\'", //name
+      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1} // m
     },
     {
       2, // dim
@@ -341,6 +620,36 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D4h = {
         0.0, 1.0, 
       }
     },
+    {
+      1, // dim
+      "A2\'\'", //name
+      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1} // m
+    },
+    {
+      1, // dim
+      "B2\'", //name
+      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1} // m
+    },
+    {
+      1, // dim
+      "A1\'", //name
+      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m
+    },
+    {
+      1, // dim
+      "A1\'\'", //name
+      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1} // m
+    },
+    {
+      1, // dim
+      "B2\'\'", //name
+      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1} // m
+    },
+    {
+      1, // dim
+      "B1\'\'", //name
+      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1} // m
+    },
     {
       2, // dim
       "E\'\'", //name
@@ -395,41 +704,6 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D4h = {
         0.0, -1.0, 
       }
     },
-    {
-      1, // dim
-      "B1\'", //name
-      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1} // m
-    },
-    {
-      1, // dim
-      "A2\'\'", //name
-      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1} // m
-    },
-    {
-      1, // dim
-      "B2\'", //name
-      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1} // m
-    },
-    {
-      1, // dim
-      "A1\'\'", //name
-      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1} // m
-    },
-    {
-      1, // dim
-      "A2\'", //name
-      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1} // m
-    },
-    {
-      1, // dim
-      "B2\'\'", //name
-      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1} // m
-    },
-    {
-      1, // dim
-      "A1\'", //name
-      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m
-    },
   } // end of irreps
 };
 
@@ -463,184 +737,110 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_trivial_g = {
   } // end of irreps
 };
 
-const qpms_finite_group_t QPMS_FINITE_GROUP_C4v = {
+const qpms_finite_group_t QPMS_FINITE_GROUP_x_and_z_flip = {
-  "C4v", // name
+  "x_and_z_flip", // name
-  8, // order
+  4, // order
   0, // idi
   (qpms_gmi_t[]) { // mt
-    0, 1, 2, 3, 4, 5, 6, 7,
+    0, 1, 2, 3,
-    1, 2, 3, 0, 7, 4, 5, 6,
+    1, 0, 3, 2,
-    2, 3, 0, 1, 6, 7, 4, 5,
+    2, 3, 0, 1,
-    3, 0, 1, 2, 5, 6, 7, 4,
+    3, 2, 1, 0,
-    4, 5, 6, 7, 0, 1, 2, 3,
-    5, 6, 7, 4, 3, 0, 1, 2,
-    6, 7, 4, 5, 2, 3, 0, 1,
-    7, 4, 5, 6, 1, 2, 3, 0,
   },
   (qpms_gmi_t[]) { // invi
-    0, 3, 2, 1, 4, 5, 6, 7
+    0, 1, 2, 3
   },
-  (qpms_gmi_t[]) {1, 7}, // gens
+  (qpms_gmi_t[]) {1, 3}, // gens
   2, // ngens
   (qpms_permutation_t[]){ // permrep
     "(3)",
-    "(0 1 2 3)",
+    "(3)(0 1)",
-    "(0 2)(1 3)",
-    "(0 3 2 1)",
-    "(3)(0 2)",
-    "(0 3)(1 2)",
-    "(1 3)",
     "(0 1)(2 3)",
+    "(2 3)",
   },
   NULL, // elemlabels
   4, // permrep_nelem
   (qpms_irot3_t[]) { // rep3d
    {{1.0+0.0*I, 0.0+0.0*I}, 1},
-   {{0.7071067811865476+0.7071067811865475*I, 0.0+0.0*I}, 1},
-   {{2.220446049250313e-16+1.0*I, 0.0+0.0*I}, 1},
-   {{-0.7071067811865474+0.7071067811865477*I, 0.0+0.0*I}, 1},
-   {{0.0+0.0*I, 0.7071067811865477-0.7071067811865474*I}, -1},
-   {{0.0+0.0*I, 1.0+2.220446049250313e-16*I}, -1},
-   {{0.0+0.0*I, 0.7071067811865475+0.7071067811865476*I}, -1},
    {{0.0+0.0*I, 0.0+1.0*I}, -1},
-  },
-  5, // nirreps
-  (struct qpms_finite_group_irrep_t[]) { // irreps
-    {
-      2, // dim
-      "E", //name
-      (complex double []) {
-        // (3)
-        1.0, 0.0, 
-        0.0, 1.0, 
-        // (0 1 2 3)
-        0.0, -1.0, 
-        1.0, 0.0, 
-        // (0 2)(1 3)
-        -1.0, 0.0, 
-        0.0, -1.0, 
-        // (0 3 2 1)
-        0.0, 1.0, 
-        -1.0, 0.0, 
-        // (3)(0 2)
-        0.0, 1.0, 
-        1.0, 0.0, 
-        // (0 3)(1 2)
-        1.0, 0.0, 
-        0.0, -1.0, 
-        // (1 3)
-        0.0, -1.0, 
-        -1.0, 0.0, 
-        // (0 1)(2 3)
-        -1.0, 0.0, 
-        0.0, 1.0, 
-      }
-    },
-    {
-      1, // dim
-      "A1", //name
-      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m
-    },
-    {
-      1, // dim
-      "B2", //name
-      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m
-    },
-    {
-      1, // dim
-      "B1", //name
-      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m
-    },
-    {
-      1, // dim
-      "A2", //name
-      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m
-    },
-  } // end of irreps
-};
-
-const qpms_finite_group_t QPMS_FINITE_GROUP_D2h = {
-  "D2h", // name
-  8, // order
-  0, // idi
-  (qpms_gmi_t[]) { // mt
-    0, 1, 2, 3, 4, 5, 6, 7,
-    1, 0, 3, 2, 5, 4, 7, 6,
-    2, 3, 0, 1, 6, 7, 4, 5,
-    3, 2, 1, 0, 7, 6, 5, 4,
-    4, 5, 6, 7, 0, 1, 2, 3,
-    5, 4, 7, 6, 1, 0, 3, 2,
-    6, 7, 4, 5, 2, 3, 0, 1,
-    7, 6, 5, 4, 3, 2, 1, 0,
-  },
-  (qpms_gmi_t[]) { // invi
-    0, 1, 2, 3, 4, 5, 6, 7
-  },
-  (qpms_gmi_t[]) {1, 3, 7}, // gens
-  3, // ngens
-  (qpms_permutation_t[]){ // permrep
-    "(5)",
-    "(5)(0 1)(2 3)",
-    "(5)(0 2)(1 3)",
-    "(5)(0 3)(1 2)",
-    "(0 3)(1 2)(4 5)",
-    "(0 2)(1 3)(4 5)",
-    "(0 1)(2 3)(4 5)",
-    "(4 5)",
-  },
-  NULL, // elemlabels
-  6, // permrep_nelem
-  (qpms_irot3_t[]) { // rep3d
-   {{1.0+0.0*I, 0.0+0.0*I}, 1},
-   {{0.0+0.0*I, 0.0+1.0*I}, -1},
-   {{0.0+1.0*I, 0.0+0.0*I}, 1},
-   {{0.0+0.0*I, 1.0+0.0*I}, -1},
-   {{0.0+0.0*I, 0.0+1.0*I}, 1},
-   {{-1.0+0.0*I, 0.0+0.0*I}, -1},
    {{0.0+0.0*I, -1.0+0.0*I}, 1},
    {{0.0+1.0*I, 0.0+0.0*I}, -1},
   },
-  8, // nirreps
+  4, // nirreps
   (struct qpms_finite_group_irrep_t[]) { // irreps
     {
       1, // dim
-      "B1\'\'", //name
+      "P\'\'", //name
-      (complex double []) {1, -1, -1, 1, -1, 1, 1, -1} // m
+      (complex double []) {1, -1, 1, -1} // m
     },
     {
       1, // dim
-      "A2\'\'", //name
+      "P\'", //name
-      (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m
+      (complex double []) {1, 1, 1, 1} // m
     },
     {
       1, // dim
-      "B1\'", //name
+      "R\'\'", //name
-      (complex double []) {1, 1, -1, -1, -1, -1, 1, 1} // m
+      (complex double []) {1, 1, -1, -1} // m
     },
     {
       1, // dim
-      "A2\'", //name
+      "R\'", //name
-      (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m
+      (complex double []) {1, -1, -1, 1} // m
+    },
+  } // end of irreps
+};
+
+const qpms_finite_group_t QPMS_FINITE_GROUP_y_and_z_flip = {
+  "y_and_z_flip", // name
+  4, // order
+  0, // idi
+  (qpms_gmi_t[]) { // mt
+    0, 1, 2, 3,
+    1, 0, 3, 2,
+    2, 3, 0, 1,
+    3, 2, 1, 0,
+  },
+  (qpms_gmi_t[]) { // invi
+    0, 1, 2, 3
+  },
+  (qpms_gmi_t[]) {1, 3}, // gens
+  2, // ngens
+  (qpms_permutation_t[]){ // permrep
+    "(3)",
+    "(3)(0 1)",
+    "(0 1)(2 3)",
+    "(2 3)",
+  },
+  NULL, // elemlabels
+  4, // permrep_nelem
+  (qpms_irot3_t[]) { // rep3d
+   {{1.0+0.0*I, 0.0+0.0*I}, 1},
+   {{0.0+0.0*I, 1.0+0.0*I}, -1},
+   {{0.0+0.0*I, 0.0+1.0*I}, 1},
+   {{0.0+1.0*I, 0.0+0.0*I}, -1},
+  },
+  4, // nirreps
+  (struct qpms_finite_group_irrep_t[]) { // irreps
+    {
+      1, // dim
+      "P\'\'", //name
+      (complex double []) {1, -1, 1, -1} // m
     },
     {
       1, // dim
-      "A1\'\'", //name
+      "P\'", //name
-      (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m
+      (complex double []) {1, 1, 1, 1} // m
     },
     {
       1, // dim
-      "B2\'", //name
+      "R\'\'", //name
-      (complex double []) {1, -1, -1, 1, 1, -1, -1, 1} // m
+      (complex double []) {1, 1, -1, -1} // m
     },
     {
       1, // dim
-      "B2\'\'", //name
+      "R\'", //name
-      (complex double []) {1, 1, -1, -1, 1, 1, -1, -1} // m
+      (complex double []) {1, -1, -1, 1} // m
-    },
-    {
-      1, // dim
-      "A1\'", //name
-      (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m
     },
   } // end of irreps
 };

tests/staticgroups.h

@@ -1,10 +1,14 @@
 #ifndef STATICGROUPS_H
 #define STATICGROUPS_H
 #include <qpms/groups.h>
+extern const qpms_finite_group_t QPMS_FINITE_GROUP_C2;
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_C2v;
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_trivial_g;
+extern const qpms_finite_group_t QPMS_FINITE_GROUP_C4;
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_C4v;
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_D3h;
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_D2h;
 extern const qpms_finite_group_t QPMS_FINITE_GROUP_D4h;
+extern const qpms_finite_group_t QPMS_FINITE_GROUP_x_and_z_flip;
+extern const qpms_finite_group_t QPMS_FINITE_GROUP_y_and_z_flip;
 #endif