scatsystem in-place T-matrix symmetrisations

Former-commit-id: b6d744319168f0df3117d066ca3143dd3c58efc2

qpms/scatsystem.c

@@ -1,6 +1,7 @@
-#include "scatsystem.h"
-
+#include <stdlib.h>
 #include <cblas.h>
+#include "scatsystem.h"
+#include "indexing.h"
 
 qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
   qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t));
@@ -10,12 +11,21 @@ qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
     size_t n = bspec->n;
     t->m = calloc(n*n, sizeof(complex double));
     if (!t->m) abort();
+    t->owns_m = true;
   }
   return t;
 }
 
+qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) {
+  qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
+  size_t n = T->spec->n;
+  for(size_t i = 0; i < n*n; ++i)
+    t->m = T->m;
+  return t;
+}
+
 void qpms_tmatrix_free(qpms_tmatrix_t *t){
-  if(t) free(t->m);
+  if(t && t->owns_m) free(t->m);
   free(t);
 }
 
@@ -24,7 +34,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
                 const complex double *M 
                 )
 {
-  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
   const size_t n = T->spec->n;
   complex double tmp[n][n];
   // tmp = M T
@@ -41,12 +51,25 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
   return t;
 }
 
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
+                qpms_tmatrix_t *T, 
+                const complex double *M 
+                )
+{
+  qpms_tmatrix_t *t = qpms_tmatrix_apply_symop(T, M);
+  const size_t n = T->spec->n;
+  for(size_t i = 0; i < n*n; ++i)
+    T->m[i] = 0.5 * (t->m[i] + T->m[i]);
+  qpms_tmatrix_free(t);
+  return T;
+}
+
 qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
                 const qpms_tmatrix_t *T, 
                 const complex double *M 
                 )
 {
-  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
   const size_t n = T->spec->n;
   complex double tmp[n][n];
   // tmp = M T
@@ -66,35 +89,43 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
 }
 
 qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(const qpms_tmatrix_t *T) {
-  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
+  return qpms_tmatrix_symmetrise_C_inf_inplace(t);
+}
+
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(qpms_tmatrix_t *T) {
   const size_t n = T->spec->n;
   for (size_t row = 0; row < n; row++) {
     qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
     for (size_t col = 0; col < n; col++) {
       qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
       if (rm == cm)
-        t->m[n*row + col] = T->m[n*row + col];
+        ;// No-op // t->m[n*row + col] = T->m[n*row + col];
       else
-        t->m[n*row + col] = 0;
+        T->m[n*row + col] = 0;
     }
   }
-  return t;
+  return T;
 }
 
 qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(const qpms_tmatrix_t *T, int N) {
-  qpms_tmatrix_t t = qpms_tmatrix_init(T->spec);
+  qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
+  return qpms_tmatrix_symmetrise_C_N_inplace(t, N);
+}
+
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) {
   const size_t n = T->spec->n;
   for (size_t row = 0; row < n; row++) {
     qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
     for (size_t col = 0; col < n; col++) {
       qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
       if (((rm - cm) % N) == 0)
-        t->m[n*row + col] = T->m[n*row + col];
+        ; // T->m[n*row + col] = T->m[n*row + col];
       else
-        t->m[n*row + col] = 0;
+        T->m[n*row + col] = 0;
     }
   }
-  return t;
+  return T;
 }
 
 

qpms/scatsystem.h

@@ -4,6 +4,7 @@
 #ifndef QPMS_SCATSYSTEM_H
 #define QPMS_SCATSYSTEM_H
 #include "vswf.h"
+#include <stdbool.h>
 
 /// A T-matrix.
 /** In the future, I might rather use a more abstract approach in which T-matrix
@@ -21,6 +22,7 @@ typedef struct qpms_tmatrix_t {
 	 */
 	const qpms_vswf_set_spec_t *spec; 
 	complex double *m; ///< Matrix elements in row-major order.
+	bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
 } qpms_tmatrix_t;
 
 /// Returns a pointer to the beginning of the T-matrix row \a rowno.
@@ -31,6 +33,9 @@ static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
 /// Initialises a zero T-matrix.
 qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec);
 
+/// Copies a T-matrix, allocating new array for the T-matrix data.
+qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T);
+
 /// Destroys a T-matrix.
 void qpms_tmatrix_free(qpms_tmatrix_t *t);
 
@@ -84,4 +89,50 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
 		int N ///< number of z-axis rotations in the group
 		);
 
+/// Symmetrizes a T-matrix with an involution symmetry operation, rewriting the original one.
+/** The symmetry operation is expected to be a unitary (square) 
+ *  matrix \a M with the same
+ *  VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
+ *  correspond to CHECKME \f[ T' = \frac{T + MTM^\dagger}{2} \f] 
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
+		qpms_tmatrix_t *T, ///< the original T-matrix
+		const complex double *M ///< the symmetry op matrix in row-major format
+		);
+
+/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one.
+/**
+ *  \f[ {T'}_{tlm}^{\tau\lambda\mu} = T_{tlm}^{\tau\lambda\mu} \delta_{m\mu} \f]
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(
+		qpms_tmatrix_t *T ///< the original T-matrix
+		);
+/// Creates a \f$ C_N \f$ -symmetrized version of a T-matrix, rewriting the original one.
+/**
+ *  \f[ {T'}_{tlm}^{\tau\lambda\mu} = \begin{cases}
+ *  T{}_{lm}^{\lambda\mu} & (m-\mu)\mod N=0\\
+ *  0 & (m-\mu)\mod N\ne0
+ *  \end{cases} . \f]
+ */
+qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(
+		qpms_tmatrix_t *T, ///< the original T-matrix
+		int N ///< number of z-axis rotations in the group
+		);
+
+/// NOT IMPLEMENTED Loads scuff-tmatrix generated files.
+/** 
+ * freqs, freqs_su, tmatrices_array and tmdata arrays are allocated by this function
+ * and have to be freed by the caller after use.
+ * The contents of tmatrices_array is NOT supposed to be freed element per element.
+ */
+qpms_errno_t qpms_load_scuff_tmatrix(
+		const char *path, ///< Path to the TMatrix file
+		const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec
+		size_t *n, ///< Number of successfully loaded t-matrices
+		double **freqs, ///< Frequencies in SI units
+		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
+		);
+
 #endif //QPMS_SCATSYSTEM_H