Scatsystem matrix packing unfinished

Former-commit-id: 9f9806529550b3e109d284395ce53c61df648d7e

qpms/scatsystem.c

@@ -783,10 +783,60 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
   return target;
 }
 
-
+#if 0
 complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
 		const complex double *orig_full, const qpms_scatsys_t *ss,
-		qpms_iri_t iri);
+		qpms_iri_t iri){
+  const size_t packedlen = ss->saecv_sizes[iri];
+  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));
+
+  const complex double one = 1;
+
+  size_t fullvec_offsetR = 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[pi].osnR;
+    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];
+    // Orbit coeff vector's full size:
+    const size_t orbit_fullsizeR = otR->size * otR->bspecn;
+    const size_t particle_fullsizeR = otR->bspecn;
+    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];
+
+    size_t fullvec_offsetC = 0;
+    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[pi].osnC;
+      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];
+      // Orbit coeff vector's full size:
+      const size_t orbit_fullsizeC = otC->size * otC->bspecn;
+      const size_t particle_fullsizeC = otC->bspecn;
+      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];
+
+    cblas_zgemv(CblasRowMajor, CblasNoTrans,
+        orbit_packedsize, particle_fullsize, &one, om + opi*particle_fullsize, orbit_fullsize,
+        orig_full+fullvec_offset, 1,
+        &one, target_packed+packed_orbit_offset, 1);
+
+    fullvec_offsetR += ot->bspecnR;
+  }
+  return target_packed;
+#endif
 
 /// Transforms a big "packed" matrix into the full basis.
 /** TODO doc */
@@ -809,7 +859,6 @@ complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
 
   size_t fullvec_offset = 0;
   for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
-    const size_t bspecn = ss->tm[ss->p[pi].tmatrix_id]->spec->n;
     const qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t;
     const qpms_ss_orbit_type_t *const ot = ss->orbit_types + oti;
     const qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn;
@@ -829,7 +878,7 @@ complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
         orig_full+fullvec_offset, 1,
         &one, target_packed+packed_orbit_offset, 1);
 
-    fullvec_offset += bspecn;
+    fullvec_offset += ot->bspecn;
   }
   return target_packed;
 }
@@ -849,7 +898,6 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
  
   size_t fullvec_offset = 0;
   for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
-    const size_t bspecn = ss->tm[ss->p[pi].tmatrix_id]->spec->n;
     const qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t;
     const qpms_ss_orbit_type_t *const ot = ss->orbit_types + oti;
     const qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn;
@@ -869,7 +917,7 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
         orig_packed+packed_orbit_offset, 1,
         &one, target_full+fullvec_offset, 1);
 
-    fullvec_offset += bspecn;
+    fullvec_offset += ot->bspecn;
   }
   return target_full; 
 }