ss matrix projections and unpacking (Broken)

Former-commit-id: 58ed4c05f8b042feb896d5dd49f5d3a5616e979f

qpms/qpms_c.pyx

@@ -1392,31 +1392,73 @@ cdef class ScatteringSystem:
     property fecv_size: 
         def __get__(self): return self.s[0].fecv_size
     property saecv_sizes: 
-        def __get__(self): return [self.s[0].saecv_sizes[i] for i in range(self.s[0].sym[0].nirreps)]
+        def __get__(self): 
+            return [self.s[0].saecv_sizes[i] 
+                for i in range(self.s[0].sym[0].nirreps)]
     property irrep_names: 
         def __get__(self): 
-            return [string_c2py(self.s[0].sym[0].irreps[iri].name) if (self.s[0].sym[0].irreps[iri].name) else None
+            return [string_c2py(self.s[0].sym[0].irreps[iri].name) 
+                    if (self.s[0].sym[0].irreps[iri].name) else None
                 for iri in range(self.s[0].sym[0].nirreps)]
     property nirreps: 
         def __get__(self): return self.s[0].sym[0].nirreps
 
     def pack_vector(self, vect, iri):
-        if len(vect) != self.fecv_size: raise ValueError("Length of a full vector has to be %d, not %d" % (self.fecv_size, len(vect)))
+        if len(vect) != self.fecv_size: 
+            raise ValueError("Length of a full vector has to be %d, not %d" 
+                    % (self.fecv_size, len(vect)))
         vect = np.array(vect, dtype=complex, copy=False, order='C')
         cdef cdouble[::1] vect_view = vect;
-        cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty((self.saecv_sizes[iri],), dtype=complex)
+        cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
+                (self.saecv_sizes[iri],), dtype=complex, order='C')
         cdef cdouble[::1] target_view = target_np
         qpms_scatsys_irrep_pack_vector(&target_view[0], &vect_view[0], self.s, iri)
         return target_np
     def unpack_vector(self, packed, iri):
-        if len(packed) != self.saecv_sizes[iri]: raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d"
+        if len(packed) != self.saecv_sizes[iri]: 
-                                                                    % (iri, self.saecv_sizes, len(packed)))
+            raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d"
+                    % (iri, self.saecv_sizes, len(packed)))
         packed = np.array(packed, dtype=complex, copy=False, order='C')
         cdef cdouble[::1] packed_view = packed
-        cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty((self.fecv_size,), dtype=complex)
+        cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty(
+                (self.fecv_size,), dtype=complex)
         cdef cdouble[::1] target_view = target_np
-        qpms_scatsys_irrep_unpack_vector(&target_view[0], &packed_view[0], self.s, iri, 0)
+        qpms_scatsys_irrep_unpack_vector(&target_view[0], &packed_view[0], 
+                self.s, iri, 0)
         return target_np
+    def pack_matrix(self, fullmatrix, iri):
+        (iri < self.nirreps)
+        cdef size_t flen = self.s[0].fecv_size
+        cdef size_t rlen = self.saecv_sizes[iri]
+        fullmatrix = np.array(fullmatrix, dtype=complex, copy=False, order='C')
+        if fullmatrix.shape != (flen, flen):
+            raise ValueError("Full matrix shape should be (%d,%d), is %s."
+                    % (flen, flen, repr(fullmatrix.shape)))
+        cdef cdouble[:,::1] fullmatrix_view = fullmatrix
+        cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
+                (rlen, rlen), dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target_np
+        qpms_scatsys_irrep_pack_matrix(&target_view[0][0], &fullmatrix_view[0][0],
+                self.s, iri)
+        return target_np
+    def unpack_matrix(self, packedmatrix, iri):
+        (iri < self.nirreps)
+        cdef size_t flen = self.s[0].fecv_size
+        cdef size_t rlen = self.saecv_sizes[iri]
+        packedmatrix = np.array(packedmatrix, dtype=complex, copy=False, order='C')
+        if packedmatrix.shape != (rlen, rlen):
+            raise ValueError("Packed matrix shape should be (%d,%d), is %s."
+                    % (rlen, rlen, repr(packedmatrix.shape)))
+        cdef cdouble[:,::1] packedmatrix_view = packedmatrix
+        cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
+                (flen, flen), dtype=complex, order='C')
+        cdef cdouble[:,::1] target_view = target_np
+        qpms_scatsys_irrep_unpack_matrix(&target_view[0][0], &packedmatrix_view[0][0],
+                self.s, iri, 0)
+        return target_np
+
+
+
 
 def tlm2uvswfi(t, l, m):
     ''' TODO doc

qpms/scatsystem.c

@@ -926,7 +926,7 @@ complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
       cblas_zgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans,
           particle_fullsizeR /*M*/, particle_fullsizeC /*N*/, orbit_packedsizeR /*K*/,
           &one /*alpha*/, omR + full_len*packed_orbit_offsetR + fullvec_offsetR/*A*/,
-          full_len/*ldA*/,
+          packed_len/*ldA*/,
           tmp /*B*/, orbit_packedsizeR /*ldB*/, &one /*beta*/,
           target_full + full_len*fullvec_offsetR + fullvec_offsetC /*C*/,
           full_len /*ldC*/);
@@ -969,10 +969,11 @@ complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
     // This is the orbit-level matrix projecting the whole orbit onto the irrep.
     const complex double *om = ot->irbases + ot->irbase_offsets[iri];
 
-    cblas_zgemv(CblasRowMajor, CblasNoTrans,
+    cblas_zgemv(CblasRowMajor/*order*/, CblasNoTrans/*transA*/,
-        orbit_packedsize, particle_fullsize, &one, om + opi*particle_fullsize, orbit_fullsize,
+        orbit_packedsize/*M*/, particle_fullsize/*N*/, &one/*alpha*/, 
-        orig_full+fullvec_offset, 1,
+        om + opi*particle_fullsize/*A*/, orbit_fullsize/*lda*/,
-        &one, target_packed+packed_orbit_offset, 1);
+        orig_full+fullvec_offset/*X*/, 1/*incX*/,
+        &one/*beta*/, target_packed+packed_orbit_offset/*Y*/, 1/*incY*/);
 
     fullvec_offset += ot->bspecn;
   }
@@ -1008,10 +1009,10 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
     // This is the orbit-level matrix projecting the whole orbit onto the irrep.
     const complex double *om = ot->irbases + ot->irbase_offsets[iri];
 
-    cblas_zgemv(CblasRowMajor, CblasConjTrans,
+    cblas_zgemv(CblasRowMajor/*order*/, CblasConjTrans/*transA*/,
-        orbit_packedsize, particle_fullsize, &one, om + opi*particle_fullsize, orbit_fullsize,
+        orbit_packedsize/*M*/, particle_fullsize/*N*/, &one/*alpha*/, om + opi*particle_fullsize/*A*/,
-        orig_packed+packed_orbit_offset, 1,
+        orbit_fullsize/*lda*/, orig_packed+packed_orbit_offset /*X*/, 1/*incX*/, &one/*beta*/, 
-        &one, target_full+fullvec_offset, 1);
+        target_full+fullvec_offset/*Y*/, 1/*incY*/);
 
     fullvec_offset += ot->bspecn;
   }

tests/scatsys.py

@@ -20,3 +20,13 @@ unpackedvectors = np.array([ss.unpack_vector(v[1], v[0]) for v in packedvectors]
 rec_fullvector = np.sum(unpackedvectors, axis=0)
 thediff = np.amax(abs(rec_fullvector-fullvector))
 assert(thediff < 1e-14)
+
+packedmatrices = list()
+for iri in range(ss.nirreps):
+    d = ss.saecv_sizes[iri]
+    m = np.random.rand(d,d)+1j*np.random.rand(d,d)
+    packedmatrices.append((iri,m))
+
+fullmatrix = np.zeros((ss.fecv_size, ss.fecv_size), dtype=complex)
+for iri, m in packedmatrices:
+    fullmatrix += ss.unpack_matrix(m, iri)