Symmetrisation of tmdata.

Former-commit-id: 504457d0259578a8eaa8db720c0f1995c5bbc6cc

qpms/scatsystem.c

@@ -11,6 +11,7 @@
 #include "vectors.h"
 #include "wigner.h"
 #include <string.h>
+#include "qpms_error.h"
 
 
 #define QPMS_SCATSYS_LEN_RTOL 1e-13
@@ -86,6 +87,71 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
   return t;
 }
 
+qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
+    complex double *tmdata, const size_t tmcount,
+    const qpms_vswf_set_spec_t *bspec,
+    const size_t n_symops, const qpms_irot3_t *symops) {
+  const size_t n = bspec->n;
+  qpms_tmatrix_t *tmcopy = qpms_tmatrix_init(bspec);
+  complex double *symop_matrices = malloc(n*n*sizeof(complex double) * n_symops);
+  if(!symop_matrices) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
+      "malloc() failed.");
+  for (size_t i = 0; i < n_symops; ++i) 
+    qpms_irot3_uvswfi_dense(symop_matrices + i*n*n, bspec, symops[i]);
+  complex double tmp[n][n];
+  const complex double one = 1, zero = 0;
+  for (size_t tmi = 0; tmi < tmcount; ++tmi) {
+    // Move the data in tmcopy; we will then write the sum directly into tmdata.
+    memcpy(tmcopy->m, tmdata+n*n*tmi, n*n*sizeof(complex double));
+    memset(tmdata+n*n*tmi, 0, n*n*sizeof(complex double));
+    for (size_t i = 0; i < n_symops; ++i) {
+      const complex double *const M = symop_matrices + i*n*n;
+      // tmp = M T
+      cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+          n, n, n, &one, M, n, tmcopy->m, n, &zero, tmp, n);
+      // tmdata[...] += tmp M* = M T M*
+      cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans,
+          n, n, n, &one, tmp, n, M, n, &one, tmdata + tmi*n*n, n);
+    }
+    for (size_t ii = 0; ii < n*n; ++ii)
+      tmdata[n*n*tmi+ii] /= n_symops;
+  }
+  free(symop_matrices);
+  qpms_tmatrix_free(tmcopy);
+  return QPMS_SUCCESS;
+}
+
+qpms_errno_t qpms_symmetrise_tmdata_finite_group(
+    complex double *tmdata, const size_t tmcount,
+    const qpms_vswf_set_spec_t *bspec,
+    const qpms_finite_group_t *pointgroup) {
+  if (!(pointgroup->rep3d)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
+      "This function requires pointgroup->rep3d to be set correctly!");
+  return qpms_symmetrise_tmdata_irot3arr(tmdata, tmcount, bspec,
+      pointgroup->order, pointgroup->rep3d);
+}
+
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
+    qpms_tmatrix_t *T,
+    size_t n_symops,
+    const qpms_irot3_t *symops
+    ) {
+  if(qpms_symmetrise_tmdata_irot3arr(T->m, 1,
+        T->spec, n_symops, symops) != QPMS_SUCCESS)
+    return NULL;
+  else return T;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
+    qpms_tmatrix_t *T,
+    const qpms_finite_group_t *pointgroup
+    ) {
+  if(qpms_symmetrise_tmdata_finite_group(T->m, 1,
+        T->spec, pointgroup) != QPMS_SUCCESS)
+    return NULL;
+  else return T;
+}
+
 qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
                 qpms_tmatrix_t *T, 
                 const complex double *M 

qpms/scatsystem.h

@@ -1,5 +1,12 @@
 /*! \file scatsystem.h
  * \brief Modern interface for finite lattice calculations, including symmetries.
+ *
+ * N.B. Only "reasonably normalised" waves are supported now in most of the
+ * functions defined here, i.e. those that can be rotated by the usual
+ * Wigner matrices, i.e. the "power" or "spharm" -normalised ones.
+ * 
+ * TODO FIXME check whether Condon-Shortley phase can have some nasty influence
+ * here; I fear that yes.
  */
 #ifndef QPMS_SCATSYSTEM_H
 #define QPMS_SCATSYSTEM_H
@@ -27,6 +34,9 @@ typedef struct qpms_tmatrix_t {
 	bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
 } qpms_tmatrix_t;
 
+struct qpms_finite_group_t;
+typedef struct qpms_finite_group_t qpms_finite_group_t;
+
 /// Returns a pointer to the beginning of the T-matrix row \a rowno.
 static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
 	return t->m + (t->spec->n * rowno);
@@ -162,7 +172,6 @@ qpms_errno_t qpms_load_scuff_tmatrix(
 		qpms_tmatrix_t **tmatrices_array, 
 		complex double **tmdata ///< The T-matrices raw contents
 		);
-
 /// Loads a scuff-tmatrix generated file.
 /** A simple wrapper over qpms_read_scuff_tmatrix() that needs a 
  * path instead of open FILE.
@@ -185,20 +194,53 @@ qpms_errno_t qpms_read_scuff_tmatrix(
 		complex double ** tmdata 
 		);
 
-/// Loads scuff-tmatrix generated files.
+/// In-place application of point group elements on raw T-matrix data.
-/** 
+/** \a tmdata can be e.g. obtained by qpms_load_scuff_tmatrix().
- * freqs, freqs_su, tmatrices_array and tmdata arrays are allocated by this function
+ * The \a symops array should always contain all elements of a finite
- * and have to be freed by the caller after use.
+ * point (sub)group, including the identity operation.
- * The contents of tmatrices_array is NOT supposed to be freed element per element.
+ *
+ * TODO more doc.
  */
-qpms_errno_t qpms_load_scuff_tmatrix(
+qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
-		const char *path, ///< Path to the TMatrix file
+		complex double *tmdata, const size_t tmcount,
-		const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec
+		const qpms_vswf_set_spec_t *bspec,
-		size_t *n, ///< Number of successfully loaded t-matrices
+		size_t n_symops,
-		double **freqs, ///< Frequencies in SI units
+		const qpms_irot3_t *symops
-		double **freqs_su, ///< Frequencies in SCUFF units (optional)
+		);
-		qpms_tmatrix_t **tmatrices_array, ///< The resulting T-matrices.
+
-		complex double **tmdata ///< The t-matrices raw contents
+/// In-place application of a point group on raw T-matrix data.
+/** This does the same as qpms_symmetrise_tmdata_irot3arr(),
+ * but takes a valid finite point group as an argument.
+ *
+ * TODO more doc.
+ */
+qpms_errno_t qpms_symmetrise_tmdata_finite_group(
+		complex double *tmdata, const size_t tmcount,
+		const qpms_vswf_set_spec_t *bspec,
+		const qpms_finite_group_t *pointgroup
+		);
+
+/// In-place application of point group elements on a T-matrix.
+/** The \a symops array should always contain all elements of a finite
+ * point (sub)group, including the identity operation.
+ *
+ * TODO more doc.
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
+		qpms_tmatrix_t *T,
+		size_t n_symops,
+		const qpms_irot3_t *symops
+		);
+
+/// In-place application of point group elements on a T-matrix.
+/** This does the same as qpms_tmatrix_symmetrise_irot3arr(),
+ * but takes a valid finite point group as an argument.
+ *
+ * TODO more doc.
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
+		qpms_tmatrix_t *T,
+		const qpms_finite_group_t *pointgroup
 		);
 
 /* Fuck this, include the whole <gsl/gsl_spline.h>
@@ -244,7 +286,6 @@ typedef struct qpms_particle_tid_t {
 	qpms_ss_tmi_t tmatrix_id; ///< T-matrix index
 } qpms_particle_tid_t;
 
-struct qpms_finite_group_t;
 
 typedef qpms_gmi_t qpms_ss_orbit_pi_t; ///< Auxilliary type used in qpms_ss_orbit_type_t for labeling particles inside orbits.
 typedef qpms_ss_tmi_t qpms_ss_oti_t; ///< Auxilliary type used for labeling orbit types.

qpms/tmatrix_io.c

@@ -29,6 +29,10 @@ 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_CS) // CHECKME CORRECT?
+    qpms_pr_error_at_flf(__FILE__, __LINE__, __func__,
+        "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)"
+        " norm supported right now.");
   int n_alloc = 128; // First chunk to allocate
   *n = 0;
   *freqs = malloc(n_alloc * sizeof(double));