Scatsystem: excitation coefficient vector irrep packing and unpacking.
Former-commit-id: 5a4b983f43a0e429a8d92eff435d763358a3ad99
This commit is contained in:
Marek Nečada
parent
0053eeb953
commit
519e60a501
|
|
@ -362,6 +362,7 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu
|
|||
|
||||
const qpms_vswf_set_spec_t *bspec = ss->tm[ot_current->tmatrices[0]]->spec;
|
||||
const size_t bspecn = bspec->n;
|
||||
ot_new->bspecn = bspecn;
|
||||
|
||||
const size_t actionsiz = sizeof(ot_current->action[0]) * ot_current->size
|
||||
* ss->sym->order;
|
||||
|
|
@ -381,6 +382,8 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu
|
|||
ss->otspace_end += irbase_sizes_siz;
|
||||
ot_new->irbase_cumsizes = (void *) (ss->otspace_end);
|
||||
ss->otspace_end += irbase_sizes_siz;
|
||||
ot_new->irbase_offsets = (void *) (ss->otspace_end);
|
||||
ss->otspace_end += irbase_sizes_siz;
|
||||
|
||||
const size_t irbases_siz = sizeof(ot_new->irbases[0]) * SQ(ot_new->size * bspecn);
|
||||
ot_new->irbases = (void *) (ss->otspace_end);
|
||||
|
|
@ -388,6 +391,7 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu
|
|||
|
||||
size_t lastbs, bs_cumsum = 0;
|
||||
for(qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
|
||||
ot_new->irbase_offsets[iri] = bs_cumsum * bspecn * ot_new->size;
|
||||
qpms_orbit_irrep_basis(&lastbs,
|
||||
ot_new->irbases + bs_cumsum*ot_new->size*bspecn,
|
||||
ot_new, bspec, ss->sym, iri);
|
||||
|
|
@ -447,7 +451,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
|
|||
|
||||
// Copy T-matrices; checking for duplicities
|
||||
|
||||
size_t max_bspecn = 0; // We'll need it later.for memory alloc estimates.
|
||||
ss->max_bspecn = 0; // We'll need it later.for memory alloc estimates.
|
||||
|
||||
qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities
|
||||
ss->tm_count = 0;
|
||||
|
|
@ -462,7 +466,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
|
|||
++(ss->tm_count);
|
||||
}
|
||||
tm_dupl_remap[i] = j;
|
||||
max_bspecn = MAX(ss->tm[i]->spec->n, max_bspecn);
|
||||
ss->max_bspecn = MAX(ss->tm[i]->spec->n, ss->max_bspecn);
|
||||
}
|
||||
|
||||
// Copy particles, remapping the t-matrix indices
|
||||
|
|
@ -508,8 +512,8 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
|
|||
ss->otspace_end = ss->otspace = malloc( // reallocate later
|
||||
(sizeof(qpms_ss_orbit_pi_t) * sym->order * sym->order
|
||||
+ sizeof(qpms_ss_tmi_t) * sym->order
|
||||
+ 2*sizeof(size_t) * sym->nirreps
|
||||
+ sizeof(complex double) * SQ(sym->order * max_bspecn)) * ss->p_count
|
||||
+ 3*sizeof(size_t) * sym->nirreps
|
||||
+ sizeof(complex double) * SQ(sym->order * ss->max_bspecn)) * ss->p_count
|
||||
);
|
||||
|
||||
// Workspace for the orbit type determination
|
||||
|
|
@ -614,6 +618,18 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp
|
|||
ss->fecv_size += ss->tm[ss->p[pi].tmatrix_id]->spec->n; // That's a lot of dereferencing!
|
||||
}
|
||||
|
||||
ss->saecv_sizes = malloc(sizeof(size_t) * sym->nirreps); if (!ss->saecv_sizes) abort();
|
||||
ss->saecv_ot_offsets = malloc(sizeof(size_t) * sym->nirreps * ss->orbit_type_count);
|
||||
if (!ss->saecv_ot_offsets) abort();
|
||||
for(qpms_iri_t iri = 0; iri < sym->nirreps; ++iri) {
|
||||
ss->saecv_sizes[iri] = 0;
|
||||
for(qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) {
|
||||
ss->saecv_ot_offsets[iri * ss->orbit_type_count + oti] = ss->saecv_sizes[iri];
|
||||
const qpms_ss_orbit_type_t *ot = &(ss->orbit_types[oti]);
|
||||
ss->saecv_sizes[iri] += ot->instance_count * ot->irbase_sizes[iri];
|
||||
}
|
||||
}
|
||||
|
||||
return ss;
|
||||
}
|
||||
|
||||
|
|
@ -628,6 +644,7 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) {
|
|||
free(ss->otspace);
|
||||
free(ss->p_orbitinfo);
|
||||
free(ss->orbit_types);
|
||||
free(ss->saecv_sizes);
|
||||
}
|
||||
free(ss);
|
||||
}
|
||||
|
|
@ -767,4 +784,94 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size,
|
|||
}
|
||||
|
||||
|
||||
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);
|
||||
|
||||
/// Transforms a big "packed" matrix into the full basis.
|
||||
/** TODO doc */
|
||||
complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
|
||||
const complex double *orig_packed, const qpms_scatsys_t *ss,
|
||||
qpms_iri_t iri, bool add);
|
||||
|
||||
/// Projects a "big" vector onto an irrep.
|
||||
/** TODO doc */
|
||||
complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
|
||||
const complex double *orig_full, const qpms_scatsys_t *ss,
|
||||
const qpms_iri_t iri) {
|
||||
const size_t packedlen = ss->saecv_sizes[iri];
|
||||
if (target_packed == NULL)
|
||||
target_packed = malloc(packedlen*sizeof(complex double));
|
||||
if (target_packed == NULL) abort();
|
||||
memset(target_packed, 0, packedlen*sizeof(complex double));
|
||||
|
||||
const complex double one = 1;
|
||||
|
||||
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;
|
||||
const qpms_ss_orbit_pi_t opi = ss->p_orbitinfo[pi].p;
|
||||
// This is where the particle's orbit starts in the "packed" vector:
|
||||
const size_t packed_orbit_offset = ss->saecv_ot_offsets[iri*ss->orbit_type_count + oti]
|
||||
+ osn * ot->irbase_sizes[iri];
|
||||
// Orbit coeff vector's full size:
|
||||
const size_t orbit_fullsize = ot->size * ot->bspecn;
|
||||
const size_t particle_fullsize = ot->bspecn;
|
||||
const size_t orbit_packedsize = ot->irbase_sizes[iri];
|
||||
// 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,
|
||||
orbit_packedsize, particle_fullsize, &one, om + opi*particle_fullsize, orbit_fullsize,
|
||||
orig_full+fullvec_offset, 1,
|
||||
&one, target_packed+packed_orbit_offset, 1);
|
||||
|
||||
fullvec_offset += bspecn;
|
||||
}
|
||||
return target_packed;
|
||||
}
|
||||
|
||||
/// Transforms a big "packed" vector into the full basis.
|
||||
/** TODO doc */
|
||||
complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
|
||||
const complex double *orig_packed, const qpms_scatsys_t *ss,
|
||||
const qpms_iri_t iri, bool add) {
|
||||
const size_t full_len = ss->fecv_size;
|
||||
if (target_full == NULL)
|
||||
target_full = malloc(full_len*sizeof(complex double));
|
||||
if (target_full == NULL) abort();
|
||||
if (!add) memset(target_full, 0, full_len*sizeof(complex double));
|
||||
|
||||
const complex double one = 1;
|
||||
|
||||
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;
|
||||
const qpms_ss_orbit_pi_t opi = ss->p_orbitinfo[pi].p;
|
||||
// This is where the particle's orbit starts in the "packed" vector:
|
||||
const size_t packed_orbit_offset = ss->saecv_ot_offsets[iri*ss->orbit_type_count + oti]
|
||||
+ osn * ot->irbase_sizes[iri];
|
||||
// Orbit coeff vector's full size:
|
||||
const size_t orbit_fullsize = ot->size * ot->bspecn;
|
||||
const size_t particle_fullsize = ot->bspecn;
|
||||
const size_t orbit_packedsize = ot->irbase_sizes[iri];
|
||||
// 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,
|
||||
orbit_packedsize, particle_fullsize, &one, om + opi*particle_fullsize, orbit_fullsize,
|
||||
orig_packed+packed_orbit_offset, 1,
|
||||
&one, target_full+fullvec_offset, 1);
|
||||
|
||||
fullvec_offset += bspecn;
|
||||
}
|
||||
return target_full;
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -314,7 +314,8 @@ typedef qpms_ss_tmi_t qpms_ss_oti_t; ///< Auxilliary type used for labeling orbi
|
|||
*
|
||||
*/
|
||||
typedef struct qpms_ss_orbit_type_t {
|
||||
qpms_ss_orbit_pi_t size; ///< Size of the orbit (a divisor of the group order).
|
||||
qpms_ss_orbit_pi_t size; ///< Size of the orbit (a divisor of the group order), i.e. number of particles on the orbit.
|
||||
size_t bspecn; ///< Individual particle's coefficient vector length. The same as ss->tm[this.tmatrices[0]]->spec->n.
|
||||
/// Action of the group elements onto the elements in this orbit.
|
||||
/** Its size is sym->order * this.size
|
||||
* and its values lie between 0 and \a this.size − 1.
|
||||
|
|
@ -341,8 +342,11 @@ typedef struct qpms_ss_orbit_type_t {
|
|||
* TODO doc.
|
||||
*/
|
||||
size_t *irbase_sizes;
|
||||
//The following are pretty redundant, TODO reduce them at some point.
|
||||
/// Cumulative sums of irbase_sizes.
|
||||
size_t *irbase_cumsizes;
|
||||
/// TODO doc.
|
||||
size_t *irbase_offsets;
|
||||
/// Per-orbit irreducible representation orthonormal bases.
|
||||
/** This also defines the unitary operator that transforms the orbital excitation coefficients
|
||||
* in the symmetry-adapted basis.
|
||||
|
|
@ -356,11 +360,13 @@ typedef struct qpms_ss_orbit_type_t {
|
|||
size_t instance_count;
|
||||
} qpms_ss_orbit_type_t;
|
||||
|
||||
typedef ptrdiff_t qpms_ss_osn_t; ///< "serial number" of av orbit in a given type.
|
||||
|
||||
/// Auxillary type used in qpms_scatsys_t that identifies the particle's orbit and its id inside that orbit.
|
||||
typedef struct qpms_ss_particle_orbitinfo {
|
||||
qpms_ss_oti_t t; ///< Orbit type.
|
||||
#define QPMS_SS_P_ORBITINFO_UNDEF (-1) ///< This labels that the particle has not yet been assigned to an orbit.
|
||||
ptrdiff_t osn; ///< "Serial number" of the orbit in the given type. TODO type and more doc.
|
||||
qpms_ss_osn_t osn; ///< "Serial number" of the orbit in the given type. TODO type and more doc.
|
||||
qpms_ss_orbit_pi_t p; ///< Order (sija, ei rankki) of the particle inside that orbit type.
|
||||
} qpms_ss_particle_orbitinfo_t;
|
||||
|
||||
|
|
@ -386,15 +392,19 @@ typedef struct qpms_scatsys_t {
|
|||
qpms_ss_particle_orbitinfo_t *p_orbitinfo; ///< Orbit type identification of each particle. (Array length is \a p_count.)
|
||||
|
||||
size_t fecv_size; ///< Number of elements of a full excitation coefficient vector size.
|
||||
//size_t *saecv_sizes; ///< NI. Number of elements of symmetry-adjusted coefficient vector sizes (order as in sym->irreps).
|
||||
size_t *saecv_sizes; ///< Number of elements of symmetry-adjusted coefficient vector sizes (order as in sym->irreps).
|
||||
|
||||
size_t *fecv_pstarts; ///< Indices of where pi'th particle's excitation coeffs start in a full excitation coefficient vector.
|
||||
size_t *saecv_ot_offsets; ///< TODO DOC. In the packed vector, saecv_ot_offsets[iri * orbit_type_count + oti] indicates start of ot
|
||||
/**< TODO maybe move it to qpms_ss_orbit_type_t, ffs. */
|
||||
//size_t **saecv_pstarts; ///< NI. Indices of where pi'th particle's excitation coeff start in a symmetry-adjusted e.c.v.
|
||||
///**< First index is irrep index as in sym->irreps, second index is particle index. */
|
||||
|
||||
// TODO shifted origin of the symmetry group etc.
|
||||
// TODO some indices for fast operations here.
|
||||
// private
|
||||
size_t max_bspecn; ///< Maximum tm[...]->spec->n. Mainly for workspace allocation.
|
||||
|
||||
|
||||
// We keep the p_orbitinfo arrays in this chunk in order to avoid memory fragmentation
|
||||
char *otspace;
|
||||
|
|
@ -416,27 +426,27 @@ void qpms_scatsys_free(qpms_scatsys_t *s);
|
|||
|
||||
/// Projects a "big" matrix onto an irrep.
|
||||
/** TODO doc */
|
||||
qpms_errno_t qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
|
||||
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);
|
||||
|
||||
/// Transforms a big "packed" matrix into the full basis.
|
||||
/** TODO doc */
|
||||
qpms_errno_t qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
|
||||
complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
|
||||
const complex double *orig_packed, const qpms_scatsys_t *ss,
|
||||
bool add);
|
||||
qpms_iri_t iri, bool add);
|
||||
|
||||
/// Projects a "big" vector onto an irrep.
|
||||
/** TODO doc */
|
||||
qpms_errno_t qpms_scatsys_irrep_pack_vector(complex double *target_packed,
|
||||
complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
|
||||
const complex double *orig_full, const qpms_scatsys_t *ss,
|
||||
qpms_iri_t iri);
|
||||
|
||||
/// Transforms a big "packed" vector into the full basis.
|
||||
/** TODO doc */
|
||||
qpms_errno_t qpms_scatsys_irrep_unpack_vector(complex double *target_full,
|
||||
complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
|
||||
const complex double *orig_packed, const qpms_scatsys_t *ss,
|
||||
bool add);
|
||||
qpms_iri_t iri, bool add);
|
||||
|
||||
/// NOT IMPLEMENTED Dumps a qpms_scatsys_t structure to a file.
|
||||
qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path);
|
||||
|
|
|
|||
Loading…
Reference in New Issue