Implementation of qpms_tmatrix_apply_operation().

Former-commit-id: 6773f0e1d02d5f929c2039f99338f08c25d0ccab

qpms/CMakeLists.txt

@@ -3,6 +3,9 @@ find_package(GSL 2.0 REQUIRED)
 find_package(BLAS REQUIRED)
 find_package(LAPACK REQUIRED)
 
+# disable an annoying warning that gives false positives probably due to a bug in gcc
+# and other not very relevant warnings
+set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wno-int-in-bool-context -Wno-comment")
 
 #includes
 set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS})

qpms/scatsystem.h

@@ -122,9 +122,10 @@ typedef struct qpms_ss_particle_orbitinfo {
 } qpms_ss_particle_orbitinfo_t;
 
 struct qpms_trans_calculator;
+struct qpms_epsmu_generator_t;
 
 typedef struct qpms_scatsys_t {
-	qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium.
+	struct qpms_qpms_epsmu_generator_t *medium; ///< Optical properties of the background medium.
 	qpms_abstract_tmatrix_t **tm; ///< T-matrices in the system
 	qpms_ss_tmi_t tm_count; ///< Number of all different T-matrices
 	qpms_ss_tmi_t tm_capacity; ///< Capacity of tm[].

qpms/tmatrices.c

@@ -1007,7 +1007,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
       if(f->op.irot3arr.owns_ops)
         free(f->op.irot3arr.ops);
       break;
-    case QPMS_TMATRIX_OPERATION_FINITE_GROUP: // Group never owned
+    case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE: // Group never owned
       break;
     case QPMS_TMATRIX_OPERATION_COMPOSE_SUM:
       {
@@ -1015,7 +1015,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
           &(f->op.compose_sum);
         if(o->opmem) {
           for(size_t i = 0; i < o->n; ++i) 
-            qpms_tmatrix_operation_clear(o->ops[i]);
+            qpms_tmatrix_operation_clear(&(o->opmem[i]));
           free(o->opmem);
         }
         free(o->ops);
@@ -1027,7 +1027,7 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
           &(f->op.compose_chain);
         if(o->opmem) {
           for(size_t i = 0; i < o->n; ++i) 
-            qpms_tmatrix_operation_clear(o->ops[i]);
+            qpms_tmatrix_operation_clear(&(o->opmem[i]));
           free(o->opmem);
         }
         free(o->ops);
@@ -1037,3 +1037,62 @@ void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
       QPMS_WTF;
   }
 }
+
+qpms_tmatrix_t *qpms_tmatrix_apply_operation(
+    const qpms_tmatrix_operation_t *f, const qpms_tmatrix_t *orig) {
+  // Certain operations could be optimized, but the effect would be marginal.
+  qpms_tmatrix_t *res = qpms_tmatrix_copy(orig); 
+  return qpms_tmatrix_apply_operation_inplace(f, res);
+}
+
+static qpms_tmatrix_t *qtao_compose_sum_inplace(qpms_tmatrix_t *T, 
+    const struct qpms_tmatrix_operation_compose_sum *cs) {
+  qpms_tmatrix_t *tmp_target = qpms_tmatrix_init(T->spec);
+  qpms_tmatrix_t *sum = qpms_tmatrix_init(T->spec);
+  for (size_t i = 0; i < cs->n; ++i) {
+    memcpy(tmp_target->m, T->m, SQ(T->spec->n) * sizeof(complex double));
+    QPMS_ENSURE(qpms_tmatrix_apply_operation_inplace(cs->ops[i] , tmp_target),
+        "Got NULL pointer from qpms_tmatrix_apply_operation_inplace, hupsis!");
+    for (size_t j = 0; j < SQ(T->spec->n); ++j)
+      sum->m[j] += tmp_target->m[j];
+  }
+  for(size_t j = 0; j < SQ(T->spec->n); ++j)
+    T->m[j] = sum->m[j] * cs->factor;
+  qpms_tmatrix_free(sum);
+  qpms_tmatrix_free(tmp_target);
+  return T;
+}
+
+static qpms_tmatrix_t *qtao_compose_chain_inplace(qpms_tmatrix_t *T,
+    const struct qpms_tmatrix_operation_compose_chain *cc) {
+  for(size_t i = 0; i < cc->n; ++i)
+    qpms_tmatrix_apply_operation_inplace(cc->ops[i], T);
+  return T;
+}
+
+static qpms_tmatrix_t *qtao_scmulz_inplace(qpms_tmatrix_t *T,
+    const struct qpms_tmatrix_operation_scmulz *s) {
+  for(size_t i = 0; i < SQ(T->spec->n); ++i)
+    T->m[i] *= s->m[i];
+  return T;
+}
+
+qpms_tmatrix_t *qpms_tmatrix_apply_operation_inplace(
+    const qpms_tmatrix_operation_t *f, qpms_tmatrix_t *T) {
+  switch(f->typ) {
+    case QPMS_TMATRIX_OPERATION_LRMATRIX:
+      return qpms_tmatrix_apply_symop_inplace(T, f->op.lrmatrix.m);
+    case QPMS_TMATRIX_OPERATION_IROT3:
+      return qpms_tmatrix_symmetrise_irot3arr_inplace(T, 1, &(f->op.irot3));
+    case QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE:
+      return qpms_tmatrix_symmetrise_finite_group_inplace(T, f->op.finite_group);
+    case QPMS_TMATRIX_OPERATION_COMPOSE_SUM:
+      return qtao_compose_sum_inplace(T, &(f->op.compose_sum));
+    case QPMS_TMATRIX_OPERATION_COMPOSE_CHAIN:
+      return qtao_compose_chain_inplace(T, &(f->op.compose_chain));
+    case QPMS_TMATRIX_OPERATION_SCMULZ:
+      return qtao_scmulz_inplace(T, &(f->op.scmulz));
+    default:
+      QPMS_WTF;
+  }
+}

qpms/tmatrices.h

@@ -235,7 +235,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
 complex double *qpms_apply_tmatrix(
 		complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
 		const complex double *a, ///< Incident field coefficient array of size T->spec->n.
-		const qpms_tmatrix_t *T
+		const qpms_tmatrix_t *T ///< T-matrix \a T to apply.
 		);
 
 /// Generic T-matrix generator function that fills a pre-initialised qpms_tmatrix_t given a frequency.
@@ -429,7 +429,7 @@ qpms_arc_function_retval_t qpms_arc_cylinder(double theta,
 		);
 
 /// Arc parametrisation of spherical particle; for qpms_arc_function_t.
-/** Useful mostly only for benchmarks, as one can use the Mie-Lorentz solution. */
+/** Useful mostly only for benchmarks or debugging, as one can use the Mie-Lorentz solution. */
 qpms_arc_function_retval_t qpms_arc_sphere(double theta,
 		const void *R; ///< Points to double containing particle's radius.
 		);
@@ -520,7 +520,7 @@ typedef struct qpms_tmatrix_function_t {
 
 /// Specifies different kinds of operations done on T-matrices
 typedef enum {
-	QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f[ T' = MTM^\dagger \f]; see @ref qpms_tmatrix_operation_lrmatrix.
+	QPMS_TMATRIX_OPERATION_LRMATRIX, ///< General matrix transformation \f$ T' = MTM^\dagger \f$; see @ref qpms_tmatrix_operation_lrmatrix.
 	QPMS_TMATRIX_OPERATION_IROT3, ///< Single rotoreflection specified by a qpms_irot3_t.
 	QPMS_TMATRIX_OPERATION_IROT3ARR, ///< Symmetrise using an array of rotoreflection; see @ref qpms_tmatrix_operation_irot3arr.
 	QPMS_TMATRIX_OPERATION_COMPOSE_SUM, ///< Apply several transformations and sum the results, see @ref qpms_tmatrix_operation_compose_sum.
@@ -534,13 +534,13 @@ typedef enum {
 struct qpms_tmatrix_operation_lrmatrix {
 	/// Raw matrix data of \a M in row-major order.
 	/** The matrix must be taylored for the given bspec! */
-	const complex double *m; 
+	complex double *m; 
 	bool owns_m; ///< Whether \a m is owned by this;
 };
 
 struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations.
 
-/// Specifies a composed operation of type \f$ T' = c\Sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t.
+/// Specifies a composed operation of type \f$ T' = c\sum_i f_i'(T) \f$ for qpms_tmatrix_operation_t.
 struct qpms_tmatrix_operation_compose_sum {
 	size_t n; ///< Number of operations in ops;
 	const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers array \a ops[] always owned by this.)
@@ -550,21 +550,21 @@ struct qpms_tmatrix_operation_compose_sum {
 	 * the memory held by \a opmem and to be properly initialised.
 	 * Each \a ops member has to point to _different_ elements of \a opmem.
 	 */ 
-	qpms_tmatrix_operation_t *opmem; 
+	struct qpms_tmatrix_operation_t *opmem; 
 };
 
 /// Specifies a composed operation of type \f$ T' =  f_{n-1}(f_{n-2}(\dots f_0(T)\dots))) \f$ for qpms_tmatrix_operation_t.
 struct qpms_tmatrix_operation_compose_chain {
 	size_t n; ///< Number of operations in ops;
 	const struct qpms_tmatrix_operation_t **ops; ///< Operations array. (Pointers owned by this.)
-	qpms_tmatrix_operation_t *opmem; ///< (Optional) operations buffer into which elements of \a ops point. (Owned by this or NULL.)
+	struct qpms_tmatrix_operation_t *opmem; ///< (Optional) operations buffer into which elements of \a ops point. (Owned by this or NULL.)
 };
 
 /// Specifies an elementwise complex multiplication of type \f$ T'_{ij} = M_{ij}T_{ij} \f$ for qpms_tmatrix_operation_t.
 struct qpms_tmatrix_operation_scmulz {
 	/// Raw matrix data of \a M in row-major order.
 	/** The matrix must be taylored for the given bspec! */
-	const complex double *m; 
+	complex double *m; 
 	bool owns_m; ///< Whether \a m is owned by this.
 };
 
@@ -573,7 +573,7 @@ struct qpms_tmatrix_operation_scmulz {
  * or qpms_symmetrise_tmdata_irot3arr_inplace(). */
 struct qpms_tmatrix_operation_irot3arr {
 	size_t n; ///< Number of rotoreflections;
-	const qpms_irot3_t *ops; ///< Rotoreflection array of size \a n.
+	qpms_irot3_t *ops; ///< Rotoreflection array of size \a n.
 	bool owns_ops; ///< Whether \a ops array is owned by this.
 };
 
@@ -584,14 +584,13 @@ typedef struct qpms_tmatrix_operation_t {
 	union {
 		struct qpms_tmatrix_operation_lrmatrix lrmatrix;
 		struct qpms_tmatrix_operation_scmulz scmulz;
-		qpms_irot3 irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3
+		qpms_irot3_t irot3; ///< Single rotoreflection; \a typ = QPMS_TMATRIX_OPERATION_IROT3
-		struct qpms_tmatrix_operation_finite_group finite_group;
 		struct qpms_tmatrix_operation_irot3arr irot3arr;
 		struct qpms_tmatrix_operation_compose_sum compose_sum;
 		struct qpms_tmatrix_operation_compose_chain compose_chain;
 		/// Finite group for QPMS_TMATRIX_OPERATION_FINITE_GROUP_SYMMETRISE.
 		/** Not owned by this; \a rep3d must be filled. */
-		const qpms_finite_group_t *finitegroup; 
+		const qpms_finite_group_t *finite_group; 
 	} op; ///< Operation data; actual type is determined by \a typ.
 } qpms_tmatrix_operation_t;