QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Matrix symmetry-adapted base packing/unpacking fixed. 

Former-commit-id: 4993a41082c5f14e6b8cc3adb41ceab902c16820
This commit is contained in:
Marek Nečada 2019-03-09 00:15:48 +00:00
parent 5f3759bdb3
commit 4dfb358de6
2 changed files with 20 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
qpms/scatsystem.c
View File

@ -839,16 +839,16 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
particle_fullsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeC /*K*/, particle_fullsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeC /*K*/,
&one /*alpha*/, orig_full + full_len*fullvec_offsetR + fullvec_offsetC/*A*/, &one /*alpha*/, orig_full + full_len*fullvec_offsetR + fullvec_offsetC/*A*/,
full_len/*ldA*/, full_len/*ldA*/,
omC + full_len*packed_orbit_offsetC + fullvec_offsetC /*B*/, omC + opiC*particle_fullsizeC /*B*/,
full_len/*ldB*/, &zero /*beta*/, orbit_fullsizeC/*ldB*/, &zero /*beta*/,
tmp /*C*/, orbit_packedsizeC /*LDC*/); tmp /*C*/, orbit_packedsizeC /*LDC*/);
// target[oiR|piR,oiC|piC] += U[...] tmp[...] // target[oiR|piR,oiC|piC] += U[...] tmp[...]
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
orbit_packedsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeR /*K*/, orbit_packedsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeR /*K*/,
&one /*alpha*/, omR + full_len*packed_orbit_offsetR + fullvec_offsetR/*A*/, &one /*alpha*/, omR + opiR*particle_fullsizeR/*A*/,
full_len/*ldA*/, orbit_fullsizeR/*ldA*/,
tmp /*B*/, particle_fullsizeR /*ldB*/, &one /*beta*/, tmp /*B*/, orbit_packedsizeC /*ldB*/, &one /*beta*/,
target_packed + packedlen*packed_orbit_offsetR + packed_orbit_offsetC /*C*/, target_packed + packedlen*packed_orbit_offsetR + packed_orbit_offsetC /*C*/,
packedlen /*ldC*/); packedlen /*ldC*/);
@ -913,21 +913,21 @@ complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
// This is the orbit-level matrix projecting the whole orbit onto the irrep. // This is the orbit-level matrix projecting the whole orbit onto the irrep.
const complex double *omC = otC->irbases + otC->irbase_offsets[iri]; const complex double *omC = otC->irbases + otC->irbase_offsets[iri];
// tmp = P U // tmp = P U
// tmp[oiR|piR,piC] = ∑_K M[piR,K] U*[K,piC] // tmp[oiR|piR,piC] = ∑_K M[piR,K] U[K,piC]
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
orbit_packedsizeR /*M*/, particle_fullsizeC /*N*/, orbit_packedsizeC /*K*/, orbit_packedsizeR /*M*/, particle_fullsizeC /*N*/, orbit_packedsizeC /*K*/,
&one /*alpha*/, orig_packed + packedlen*packed_orbit_offsetR + packed_orbit_offsetC/*A*/, &one /*alpha*/, orig_packed + packedlen*packed_orbit_offsetR + packed_orbit_offsetC/*A*/,
packedlen/*ldA*/, packedlen/*ldA*/,
omC + full_len*packed_orbit_offsetC + fullvec_offsetC /*B*/, omC + opiC*particle_fullsizeC /*B*/,
full_len/*ldB*/, &zero /*beta*/, orbit_fullsizeC/*ldB*/, &zero /*beta*/,
tmp /*C*/, particle_fullsizeC /*LDC*/); tmp /*C*/, particle_fullsizeC /*LDC*/);
// target[oiR|piR,oiC|piC] += U[...] tmp[...] // target[oiR|piR,oiC|piC] += U*[...] tmp[...]
cblas_zgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans, cblas_zgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans,
particle_fullsizeR /*M*/, particle_fullsizeC /*N*/, orbit_packedsizeR /*K*/, particle_fullsizeR /*M*/, particle_fullsizeC /*N*/, orbit_packedsizeR /*K*/,
&one /*alpha*/, omR + full_len*packed_orbit_offsetR + fullvec_offsetR/*A*/, &one /*alpha*/, omR + opiR*particle_fullsizeR/*A*/,
packed_len/*ldA*/, orbit_fullsizeR/*ldA*/,
tmp /*B*/, orbit_packedsizeR /*ldB*/, &one /*beta*/, tmp /*B*/, particle_fullsizeC /*ldB*/, &one /*beta*/,
target_full + full_len*fullvec_offsetR + fullvec_offsetC /*C*/, target_full + full_len*fullvec_offsetR + fullvec_offsetC /*C*/,
full_len /*ldC*/); full_len /*ldC*/);
@ -952,6 +952,7 @@ complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
memset(target_packed, 0, packedlen*sizeof(complex double)); memset(target_packed, 0, packedlen*sizeof(complex double));
const complex double one = 1; const complex double one = 1;
const size_t full_len = ss->fecv_size;
size_t fullvec_offset = 0; size_t fullvec_offset = 0;
for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
@ -992,6 +993,7 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
if (!add) memset(target_full, 0, full_len*sizeof(complex double)); if (!add) memset(target_full, 0, full_len*sizeof(complex double));
const complex double one = 1; const complex double one = 1;
const size_t packedlen = ss->saecv_sizes[iri];
size_t fullvec_offset = 0; size_t fullvec_offset = 0;
for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {

6
tests/scatsys.py
View File

@ -30,3 +30,9 @@ for iri in range(ss.nirreps):
fullmatrix = np.zeros((ss.fecv_size, ss.fecv_size), dtype=complex) fullmatrix = np.zeros((ss.fecv_size, ss.fecv_size), dtype=complex)
for iri, m in packedmatrices: for iri, m in packedmatrices:
fullmatrix += ss.unpack_matrix(m, iri) fullmatrix += ss.unpack_matrix(m, iri)
for iri, m in packedmatrices:
print (m.shape)
repackedmatrix = ss.pack_matrix(fullmatrix,iri)
print(np.amax(abs(repackedmatrix-m)))