diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 1e9a85b..21b32aa 100644 --- a/qpms/scatsystem.c +++ b/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; }