From 08ef7ea39336e50e8e766a852a40cde9c656ab2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 21 Feb 2019 05:18:50 +0000 Subject: [PATCH] scatsystem in-place T-matrix symmetrisations Former-commit-id: b6d744319168f0df3117d066ca3143dd3c58efc2 --- qpms/scatsystem.c | 57 ++++++++++++++++++++++++++++++++++++----------- qpms/scatsystem.h | 51 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+), 13 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index f6a3881..4021eed 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1,6 +1,7 @@ -#include "scatsystem.h" - +#include #include +#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; } diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 6680b77..cbaf479 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -4,6 +4,7 @@ #ifndef QPMS_SCATSYSTEM_H #define QPMS_SCATSYSTEM_H #include "vswf.h" +#include /// 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