Translation-cached version of modeproblem matrix.

Former-commit-id: c76945d915138870e1a4f150038705a5fa82ce48

qpms/scatsys_private.h

@@ -1,11 +1,19 @@
 #ifndef QPMS_SCATSYS_PRIVATE_H
 #define QPMS_SCATSYS_PRIVATE_H
+/* 
+ * This file includes some definitions shared between scatsystem.c
+ * and scatsys_translation_booster.c that are not needed anywhere 
+ * else.
+ */
 
 #include "scatsystem.h"
 
+#include <pthread.h>
+
 complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted(
 		complex double *target, const qpms_scatsys_at_omega_t *ssw);
-
+complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_boosted(
+		complex double *target_packed, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri);
 /// "private" destructor, called by qpms_scatsys_free()
 void qpms_scatsys_translation_booster_free(struct qpms_scatsys_translation_booster *);
 /// "private" constructor, use qpms_ss_create_translation_cache() instead.
@@ -18,4 +26,16 @@ qpms_scatsysw_translation_booster_create(const qpms_scatsys_at_omega_t *ssw);
 
 /// "private" destructor, called by qpms_scatsys_at_omega_free()
 void qpms_scatsysw_translation_booster_free(struct qpms_scatsysw_translation_booster *);
-#endif 
+
+
+struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{
+  const qpms_scatsys_at_omega_t *ssw;
+  qpms_ss_pi_t *opistartR_ptr;
+  pthread_mutex_t *opistartR_mutex;
+  qpms_iri_t iri;
+  complex double *target_packed;
+};
+
+void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_boosted(void *arg);
+
+#endif //QPMS_SCATSYS_PRIVATE_H

qpms/scatsys_translation_booster.c

@@ -10,10 +10,18 @@
 #include "vectors.h"
 #include "qpms_error.h"
 #include "qpms_specfunc.h"
+#include "groups.h"
 
 #define SQ(x) ((x)*(x))
 #define MAX(x, y) (((x) > (y)) ? (x) : (y))
 
+#ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
+#include "qpmsblas.h"
+#define SERIAL_ZGEMM qpms_zgemm
+#else
+#define SERIAL_ZGEMM cblas_zgemm
+#endif
+
 typedef size_t ppid_t;
 typedef size_t uoppid_t;
 
@@ -140,7 +148,7 @@ void qpms_scatsys_translation_booster_free(booster_t *b) {
 }
 
 int qpms_ss_create_translation_cache(qpms_scatsys_t *ss, qpms_ss_caching_mode_t m) {
-QPMS_ASSERT(ss);
+  QPMS_ASSERT(ss);
   if (ss->tbooster) {
     QPMS_WARN("Translation cache already created?");
     return 0;
@@ -154,8 +162,8 @@ QPMS_ASSERT(ss);
       ss->tbooster = qpms_scatsys_translation_booster_create(ss);
       if (ss->tbooster) return 0;
       else {
-          QPMS_WARN("Failed to create tranlation operator cache");
+        QPMS_WARN("Failed to create tranlation operator cache");
-          return -1;
+        return -1;
       }
     default:
       QPMS_WTF;
@@ -225,6 +233,10 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted(
     QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
   complex double *tmp;
   QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double));
+  // Workspaces for the translation operator A and B matrices
+  complex double *A, *B;
+  QPMS_CRASHING_MALLOC(A, SQ(c->nelem) * sizeof(*A));
+  QPMS_CRASHING_MALLOC(B, SQ(c->nelem) * sizeof(*B));
   memset(target, 0, SQ(full_len) * sizeof(complex double)); //unnecessary?
   double legendre_buf[gsl_sf_legendre_array_n(2*c->lMax + 1)]; //VLA, workspace for legendre arrays
   const complex double zero = 0, minusone = -1;
@@ -251,18 +263,16 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted(
             QPMS_PARANOID_ASSERT(c->normalisation == bspecC->norm && c->normalisation == bspecR->norm);
             QPMS_PARANOID_ASSERT(c->lMax >= bspecC->lMax && c->lMax >= bspecR->lMax);
             QPMS_PARANOID_ASSERT(bspecC->lMax_L < 0 && bspecR->lMax_L < 0);
-            complex double A[pair_nelem][pair_nelem]; // VLA, TODO flatten
-            complex double B[pair_nelem][pair_nelem]; // VLA, TODO flatten
             { // this replaces qpms_trans_calculator_get_AB_arrays() and ..._buf()
               const double costheta = cos(dlj.theta);
               QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,
-                    2*c->lMax+1, costheta, -1, legendre_buf));
+                    2*pair_lMax+1, costheta, -1, legendre_buf));
               const double * const legendres = legendre_buf;
               const complex double * const bessels = bw->bessels + b->bessel_offsets_r[ri];
-              qpms_trans_calculator_get_AB_arrays_precalcbuf(c, pair_lMax, A[0], B[0],
+              qpms_trans_calculator_get_AB_arrays_precalcbuf(c, pair_lMax, A, B,
                   /*deststride*/ pair_nelem, /*srcstride*/ 1, dlj.phi, bessels, legendres);
               qpms_trans_array_from_AB(tmp,// tmp is S(piR<-piC) 
-                  bspecR, bspecC->n, bspecC, 1, A[0], B[0], pair_lMax);
+                  bspecR, bspecC->n, bspecC, 1, A, B, pair_lMax);
             }
           }
           cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
@@ -281,5 +291,170 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted(
   for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] = +1;
 
   free(tmp);
+  free(A);
+  free(B);
   return target;
 }
+
+void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_boosted(void *arg)
+{
+  const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg
+    *a = arg;
+  const qpms_scatsys_at_omega_t *ssw = a->ssw;
+  const complex double k = ssw->wavenumber;
+  const qpms_scatsys_t *ss = ssw->ss;
+  const qpms_iri_t iri = a->iri;
+  const size_t packedlen = ss->saecv_sizes[iri];
+
+  QPMS_ASSERT(ssw->translation_cache && ssw->ss->tbooster);
+  const booster_t *const b = ss->tbooster;
+  const boosterw_t *const bw = ssw->translation_cache;
+  const qpms_trans_calculator *const c = ss->c;
+
+  // some of the following workspaces are probably redundant; TODO optimize later.
+
+  // workspaces for the uncompressed particle<-particle tranlation matrix block
+  // and the result of multiplying with a T-matrix (times -1)
+  complex double *Sblock, *TSblock;
+  QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn));
+  QPMS_CRASHING_MALLOC(TSblock, sizeof(complex double)*SQ(ss->max_bspecn));
+  // Workspaces for the translation operator A and B matrices
+  complex double *A, *B;
+  QPMS_CRASHING_MALLOC(A, SQ(c->nelem) * sizeof(*A));
+  QPMS_CRASHING_MALLOC(B, SQ(c->nelem) * sizeof(*B));
+  double legendre_buf[gsl_sf_legendre_array_n(2*c->lMax + 1)]; //VLA, workspace for legendre arrays
+
+  // Workspace for the intermediate particle-orbit matrix result
+  complex double *tmp;
+  QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order);
+
+  const complex double one = 1, zero = 0, minusone = -1;
+
+  while(1) {
+    // In the beginning, pick a target (row) orbit for this thread
+    QPMS_ENSURE_SUCCESS(pthread_mutex_lock(a->opistartR_mutex));
+    if(*(a->opistartR_ptr) >= ss->p_count) {// Everything is already done, end
+      QPMS_ENSURE_SUCCESS(pthread_mutex_unlock(a->opistartR_mutex));
+      break;
+    }
+    const qpms_ss_pi_t opistartR = *(a->opistartR_ptr);
+    // Now increment it for another thread:
+    *(a->opistartR_ptr) += ss->orbit_types[ss->p_orbitinfo[ss->p_by_orbit[opistartR]].t].size;
+    QPMS_ENSURE_SUCCESS(pthread_mutex_unlock(a->opistartR_mutex));
+
+    // Orbit picked (defined by opistartR), do the work.
+    const qpms_ss_pi_t orbitstartpiR = ss->p_by_orbit[opistartR];
+    const qpms_ss_oti_t otiR = ss->p_orbitinfo[orbitstartpiR].t;
+    const qpms_ss_osn_t osnR = ss->p_orbitinfo[orbitstartpiR].osn;
+    const qpms_ss_orbit_type_t *const otR = ss->orbit_types + otiR;
+    const qpms_ss_orbit_pi_t orbit_p_countR = otR->size;
+    const size_t orbit_packedsizeR = otR->irbase_sizes[iri];
+
+    if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
+      const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n
+      const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec;
+      // This is the orbit-level matrix projecting the whole orbit onto the irrep.
+      const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
+      // Orbit coeff vector's full size:
+      const size_t orbit_fullsizeR = otR->size * otR->bspecn;
+      // This is where the orbit starts in the "packed" vector:
+      const size_t packed_orbit_offsetR =
+        ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiR]
+        + osnR * otR->irbase_sizes[iri];
+      for(qpms_ss_orbit_pi_t opiR = 0; opiR < orbit_p_countR; ++opiR) {
+        for(qpms_ss_orbit_pi_t opiR = 0; opiR < orbit_p_countR; ++opiR) {
+          qpms_ss_pi_t piR = ss->p_by_orbit[opistartR + opiR];
+          assert(opiR == ss->p_orbitinfo[piR].p);
+          assert(otiR == ss->p_orbitinfo[piR].t);
+          assert(ss->p_orbitinfo[piR].osn == osnR);
+          const cart3_t posR = ss->p[piR].pos;
+          // dest particle T-matrix
+          const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m;
+          for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop
+            const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t;
+            const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC;
+            const qpms_ss_osn_t osnC = ss->p_orbitinfo[piC].osn;
+            const qpms_ss_orbit_pi_t opiC = ss->p_orbitinfo[piC].p;
+            // This is where the particle's orbit starts in the "packed" vector:
+            const size_t packed_orbit_offsetC =
+              ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC]
+              + osnC * otC->irbase_sizes[iri];
+            const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec;
+            // Orbit coeff vector's full size:
+            const size_t orbit_fullsizeC = otC->size * otC->bspecn;
+            const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n
+            const size_t orbit_packedsizeC = otC->irbase_sizes[iri];
+            // This is the orbit-level matrix projecting the whole orbit onto the irrep.
+            const complex double *omC = otC->irbases + otC->irbase_offsets[iri];
+
+            if(orbit_packedsizeC) { // avoid zgemm crash on empty irrep
+              if(piC != piR) { // non-diagonal, calculate TS
+                const cart3_t posC = ss->p[piC].pos;
+#if 0
+                QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
+                      Sblock, // Sblock is S(piR->piC)
+                      bspecR, bspecC->n, bspecC, 1,
+                      k, posR, posC, QPMS_HANKEL_PLUS));
+#endif
+                { // this block replaces qpms_trans_calculator_get_trans_array():
+                  // R is dest, C is src
+                  const sph_t dlj = cart2sph(cart3_substract(posR, posC));
+                  const uoppid_t pid = uopairid(ss->p_count, piC, piR);
+                  const size_t ri  = b->r_map[pid];
+                  QPMS_PARANOID_ASSERT(dlj.r == b->r[ri]);
+                  const qpms_l_t pair_lMax = b->lMax_r[ri];
+                  const qpms_y_t pair_nelem = qpms_lMax2nelem(pair_lMax);
+
+                  { // this replaces qpms_trans_calculator_get_AB_arrays() and _buf()
+                    const double costheta = cos(dlj.theta);
+                    QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,
+                          2*pair_lMax+1, costheta, -1, legendre_buf));
+                    const double * const legendres = legendre_buf;
+                    const complex double * const bessels = bw->bessels + b->bessel_offsets_r[ri];
+                    qpms_trans_calculator_get_AB_arrays_precalcbuf(c, pair_lMax, A, B,
+                        /*deststride*/ pair_nelem, /*srcstride*/ 1, dlj.phi, bessels, legendres);
+                  }
+                  qpms_trans_array_from_AB(Sblock, // Sblock is S(piR->piC)
+                      bspecR, bspecC->n, bspecC, 1, A, B, pair_lMax);
+                }
+
+                SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+                    bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
+                    &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/,
+                    Sblock/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,
+                    TSblock /*c*/, bspecC->n /*ldc*/);
+              } else { // diagonal, fill with diagonal +1
+                for (size_t row = 0; row < bspecR->n; ++row)
+                  for (size_t col = 0; col < bspecC->n; ++col)
+                    TSblock[row * bspecC->n + col] = (row == col)? +1 : 0;
+              }
+
+              // tmp[oiR|piR,piC] = ∑_K M[piR,K] U*[K,piC]
+              SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasConjTrans,
+                  particle_fullsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeC /*K*/,
+                  &one /*alpha*/, TSblock/*A*/, particle_fullsizeC/*ldA*/,
+                  omC + opiC*particle_fullsizeC /*B*/,
+                  orbit_fullsizeC/*ldB*/, &zero /*beta*/,
+                  tmp /*C*/, orbit_packedsizeC /*LDC*/);
+
+              // target[oiR|piR,oiC|piC] += U[...] tmp[...]
+              SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+                  orbit_packedsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeR /*K*/,
+                  &one /*alpha*/, omR + opiR*particle_fullsizeR/*A*/, orbit_fullsizeR/*ldA*/,
+                  tmp /*B*/, orbit_packedsizeC /*ldB*/, &one /*beta*/,
+                  a->target_packed + packedlen*packed_orbit_offsetR + packed_orbit_offsetC /*C*/,
+                  packedlen /*ldC*/);
+            }
+          }
+        }
+      }
+    }
+  }
+  free(tmp);
+  free(A);
+  free(B);
+  free(Sblock);
+  free(TSblock);
+  return NULL;
+}
+

qpms/scatsystem.c

@@ -1403,14 +1403,6 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
   return target_packed;
 }
 
-struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{
-  const qpms_scatsys_at_omega_t *ssw;
-  qpms_ss_pi_t *opistartR_ptr;
-  pthread_mutex_t *opistartR_mutex;
-  qpms_iri_t iri;
-  complex double *target_packed;
-};
-
 static void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg)
 {
   const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg 
@@ -1529,9 +1521,6 @@ static void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_threa
         }
       }
     }
-
-
-
   }
   free(tmp);
   free(Sblock);
@@ -1733,6 +1722,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
     const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri
     )
 {
+  const _Bool use_translation_cache = (ssw->translation_cache != NULL);
   const size_t packedlen = ssw->ss->saecv_sizes[iri];
   if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
     return target_packed; 
@@ -1765,7 +1755,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
   pthread_t thread_ids[nthreads];
   for(long thi = 0; thi < nthreads; ++thi)
     QPMS_ENSURE_SUCCESS(pthread_create(thread_ids + thi, NULL,
-      qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread,
+      use_translation_cache 
+        ? qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_boosted
+        : qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread,
       (void *) &arg));
   for(long thi = 0; thi < nthreads; ++thi) {
     void *retval;