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

Fix unitialised diagonal blocks in full modeproblem matrix. 

The bug in qpms_scatsysw_build_modeproblem_matrix_full()
via qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full()
affected finite-size arrays: instead of filling the diagonal
interaction blocks with zeros, they were filled with undefined
values from a superfluous zgemm() call.
This commit is contained in:
Marek Nečada 2021-12-21 11:29:44 +02:00
parent 77a3566aba
commit aea0afee15
1 changed files with 8 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
qpms/scatsystem.c
View File

@ -1295,7 +1295,6 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
complex double *tmp; // translation matrix, S or W complex double *tmp; // translation matrix, S or W
QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double)); QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double));
memset(target, 0, SQ(full_len) * sizeof(complex double));
const complex double zero = 0, minusone = -1; const complex double zero = 0, minusone = -1;
{ // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC) { // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC)
size_t fullvec_offsetR = 0; size_t fullvec_offsetR = 0;
@ -1321,7 +1320,14 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
QPMS_EWALD_FULL, eta)); QPMS_EWALD_FULL, eta));
} }
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, if (k == NULL && piC == piR) {
// non-periodic case diagonal block: no "self-interaction", just
// fill with zeroes (the ones on the diagonal are added in the very end)
QPMS_ENSURE_SUCCESS(LAPACKE_zlaset(CblasRowMajor, 'x',
bspecR->n /*m*/, bspecC->n/*n*/, 0 /*alfa: offdiag*/, 0 /*beta: diag*/,
target + fullvec_offsetR*full_len + fullvec_offsetC /*a*/,
full_len /*lda*/));
} else cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
&minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/,
tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,