From 37d022e6d0f24fef901f2a2e7c28315c238eea72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 15:43:28 +0200 Subject: [PATCH] Fix offset incrementation on the diagonal. Solves the translation matrix projection discrepancy. Former-commit-id: c9681f2b8254e8dcb6dc6f05b0bd53b3e2a49730 --- qpms/scatsystem.c | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 07f70a4..7a66196 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1120,13 +1120,14 @@ complex double *qpms_scatsys_build_translation_matrix_full( const cart3_t posR = ss->p[piR].pos; size_t fullvec_offsetC = 0; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { - if(piC == piR) continue; // The diagonal will be dealt with later. const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec; - const cart3_t posC = ss->p[piC].pos; - QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, - target + fullvec_offsetR*full_len + fullvec_offsetC, - bspecR, full_len, bspecC, 1, - k, posR, posC)); + if(piC != piR) { // The diagonal will be dealt with later. + const cart3_t posC = ss->p[piC].pos; + QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, + target + fullvec_offsetR*full_len + fullvec_offsetC, + bspecR, full_len, bspecC, 1, + k, posR, posC)); + } fullvec_offsetC += bspecC->n; } assert(fullvec_offsetC = full_len); @@ -1162,19 +1163,20 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full( // dest particle T-matrix const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { - if(piC == piR) continue; // The diagonal will be dealt with later. const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec; - const cart3_t posC = ss->p[piC].pos; - QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, - tmp, // tmp is S(piR<-piC) - bspecR, bspecC->n, bspecC, 1, - k, posR, posC)); - cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + if(piC != piR) { // The diagonal will be dealt with later. + const cart3_t posC = ss->p[piC].pos; + QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, + tmp, // tmp is S(piR<-piC) + bspecR, bspecC->n, bspecC, 1, + k, posR, posC)); + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, target + fullvec_offsetR*full_len + fullvec_offsetC /*c*/, full_len /*ldc*/); + } fullvec_offsetC += bspecC->n; } fullvec_offsetR += bspecR->n;