"Direct" generation of irrep-packed mode problem matrices.

Former-commit-id: e0833046f2f6052e5e0523dce6047e89d23509ff

qpms/qpms_c.pyx

@@ -1464,6 +1464,14 @@ cdef class ScatteringSystem:
         qpms_scatsys_build_modeproblem_matrix_full(&target_view[0][0], self.s, k)
         return target
 
+    def modeproblem_matrix_packed(self, double k, qpms_iri_t iri):
+        cdef size_t rlen = self.saecv_sizes[iri]
+        cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
+                (rlen,rlen),dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target
+        qpms_scatsys_build_modeproblem_matrix_irrep_packed(&target_view[0][0], self.s, iri, k)
+        return target
+
     def translation_matrix_full(self, double k):
         cdef size_t flen = self.s[0].fecv_size
         cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(

qpms/scatsystem.c

@@ -1314,13 +1314,115 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
   return target;
 }
 
-
-
-
 complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
     /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
-    complex double *target,
+    complex double *target_packed,
     const qpms_scatsys_t *ss, qpms_iri_t iri,
     double k ///< Wave number to use in the translation matrix.
-    );
+    )
+{
+  const size_t packedlen = ss->saecv_sizes[iri];
+  if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps?
+    return target_packed; 
+  const size_t full_len = ss->fecv_size;
+  if (target_packed == NULL)
+    target_packed = malloc(SQ(packedlen)*sizeof(complex double));
+  if (target_packed == NULL) abort();
+  memset(target_packed, 0, SQ(packedlen)*sizeof(complex double));
+
+  // 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
+  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));
+
+  // 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;
+
+  for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) { //Row loop
+    const qpms_ss_oti_t otiR = ss->p_orbitinfo[piR].t;
+    const qpms_ss_orbit_type_t *const otR = ss->orbit_types + otiR;
+    const qpms_ss_osn_t osnR = ss->p_orbitinfo[piR].osn;
+    const qpms_ss_orbit_pi_t opiR = ss->p_orbitinfo[piR].p;
+    // This is where the particle's 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];
+    const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec;
+    // Orbit coeff vector's full size:
+    const size_t orbit_fullsizeR = otR->size * otR->bspecn;
+    const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n
+    const size_t orbit_packedsizeR = otR->irbase_sizes[iri];
+    // This is the orbit-level matrix projecting the whole orbit onto the irrep.
+    const complex double *omR = otR->irbases + otR->irbase_offsets[iri];
+    const cart3_t posR = ss->p[piR].pos;
+    if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep
+      // dest particle T-matrix
+      const complex double *tmmR = ss->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 = ss->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;
+            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));
+
+						cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+								bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
+								&one/*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]
+          cblas_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[...]
+          cblas_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*/,
+              target_packed + packedlen*packed_orbit_offsetR + packed_orbit_offsetC /*C*/,
+              packedlen /*ldC*/);
+        }
+      }
+    }
+  }
+  free(tmp);
+  free(Sblock);
+  free(TSblock);
+  return target_packed;
+}
+
 

qpms/scatsystem.h

@@ -482,7 +482,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
 		const qpms_scatsys_t *ss,
 		double k ///< Wave number to use in the translation matrix.
 		);
-complex double *qpms_scatsys_build_modeproblem_irrep_packed(
+complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,
 		const qpms_scatsys_t *ss, qpms_iri_t iri,