Extract and inline translation matrix reordering procedure.

Former-commit-id: 327e1ccf9efc26647896db5bbe1a1a4e4d034d3c

qpms/CMakeLists.txt

@@ -31,7 +31,7 @@ target_include_directories (qpms PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 
 target_compile_options(qpms PRIVATE -Wall -Wno-return-type -Wno-unused-variable -Wno-unused-function -Wno-unused-but-set-variable -Wno-unused-label)
 target_compile_definitions(qpms PRIVATE LATTICESUMS32 QPMS_VECTORS_NICE_TRANSFORMATIONS
-	EWALD_AUTO_CUTOFF
+	EWALD_AUTO_CUTOFF QPMS_EVALUATE_PARANOID_ASSERTS
 	)
 
 install(TARGETS qpms

qpms/qpms_error.h

@@ -100,6 +100,13 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 		qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error (%d)", errorcode); \
 }
 
+// Same as previous, but with message.
+#define QPMS_ENSURE_SUCCESS_M(x, msg, ...) { \
+	int errorcode = (x); \
+	if(QPMS_UNLIKELY(errorcode)) \
+		qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); \
+}
+
 
 #define QPMS_ENSURE(x, msg, ...) {if(QPMS_UNLIKELY(!(x))) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); }
 
@@ -109,6 +116,16 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
 			"Unexpected error. This is most certainly a bug.");\
 }
 
+#ifdef QPMS_EVALUATE_PARANOID_ASSERTS
+  #define QPMS_PARANOID_ASSERT(x) {\
+	if(QPMS_UNLIKELY(!(x)))\
+       		qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,\
+			"Unexpected error. This is most certainly a bug.");\
+}
+#else
+  #define QPMS_PARANOID_ASSERT(x) {;}
+#endif
+
 #define QPMS_NOT_IMPLEMENTED(msg, ...) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, \
 		"Not implemented:" msg, ##__VA_ARGS__)
 

qpms/translations.c

@@ -13,6 +13,7 @@
 #include <gsl/gsl_sf_coupling.h>
 #include "qpms_error.h"
 #include "normalisation.h"
+#include "translations_inlines.h"
 
 /*
  * Define macros with additional factors that "should not be there" according
@@ -717,20 +718,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
   qpms_errno_t retval = qpms_trans_calculator_get_AB_arrays(c,
       A[0], B[0], c->nelem, 1,
       kdlj, r_ge_d, J);
-  for (size_t desti = 0; desti < destspec->n; ++desti) {
-    qpms_y_t desty; qpms_vswf_type_t destt;
-    if(QPMS_SUCCESS != qpms_uvswfi2ty(destspec->ilist[desti], &destt, &desty))
-        qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,
-          "Invalid u. vswf index %llx.", destspec->ilist[desti]);
-    for (size_t srci = 0; srci < srcspec->n; ++srci){
-      qpms_y_t srcy; qpms_vswf_type_t srct;
-      if(QPMS_SUCCESS != qpms_uvswfi2ty(srcspec->ilist[srci], &srct, &srcy))
-          qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,
-            "Invalid u. vswf index %llx.", srcspec->ilist[srci]);
-      target[srci * srcstride + desti * deststride]
-        = (srct == destt) ? A[desty][srcy] : B[desty][srcy];
-    }
-  }
+  qpms_trans_array_from_AB(target, destspec, deststride, srcspec, srcstride,
+      A[0], B[0], c->lMax);
   return retval;
 }
 
@@ -766,7 +755,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul
   qpms_errno_t retval = qpms_trans_calculator_get_AB_arrays_e32_e(c,
       A, Aerr, B, Berr, ldAB, 1, 
       eta, k, b1, b2, beta, particle_shift, maxR, maxK, parts);
-  for (size_t desti = 0; desti < destspec->n; ++desti) {
+  for (size_t desti = 0; desti < destspec->n; ++desti) { 
+    // TODO replace with (modified) qpms_trans_array_from_AB()
     qpms_y_t desty; qpms_vswf_type_t destt;
     if(QPMS_SUCCESS != qpms_uvswfi2ty(destspec->ilist[desti], &destt, &desty))
         qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,

qpms/translations_inlines.h

@@ -0,0 +1,38 @@
+#ifndef TRANSLATIONS_INLINES_H
+#define TRANSLATIONS_INLINES_H
+#include "translations.h"
+#include "indexing.h"
+
+/// Rearranges the default-ordered "A,B" array elements into "bspec"-defined matrix.
+// TODO DOC
+static inline void qpms_trans_array_from_AB(
+		complex double *t,
+		const qpms_vswf_set_spec_t *const t_destspec,
+		const size_t t_deststride,
+		const qpms_vswf_set_spec_t *const t_srcspec,
+		const size_t t_srcstride,
+		const complex double *const A, const complex double *const B,
+		/// A and B matrices' lMax.
+		/** This also determines their size and stride: they are assumed to
+		 * be square matrices of size `nelem * nelem` where
+		 * `nelem = qpms_lMax2nelem(lMax_AB)`
+		 */
+		const qpms_l_t lMax_AB
+		) {
+	QPMS_PARANOID_ASSERT(lMax_AB >= t_srcspec->lMax && lMax_AB >= t_destspec->lMax);
+	const qpms_y_t nelem_AB = qpms_lMax2nelem(lMax_AB);
+	for (size_t desti = 0; desti < t_destspec->n; ++desti) {
+		qpms_y_t desty; qpms_vswf_type_t destt;
+		QPMS_ENSURE_SUCCESS_M(qpms_uvswfi2ty(t_destspec->ilist[desti], &destt, &desty),
+				"Invalid u. vswf index %llx.", t_destspec->ilist[desti]);
+		for (size_t srci = 0; srci < t_srcspec->n; ++srci){
+			qpms_y_t srcy; qpms_vswf_type_t srct;
+			QPMS_ENSURE_SUCCESS_M(qpms_uvswfi2ty(t_srcspec->ilist[srci], &srct, &srcy),
+					"Invalid u. vswf index %llx.", t_srcspec->ilist[srci]);                                                                t[srci * t_srcstride + desti * t_deststride]
+				= (srct == destt) ? A[desty*nelem_AB + srcy] : B[desty*nelem_AB + srcy];
+		}
+	}
+}
+
+
+#endif