From 5758c5d587548953b43a900e40d163b46b3009be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 23 Jul 2019 08:43:25 +0300 Subject: [PATCH] Reverse modeproblem matrix sign so it applies directly to scattering Former-commit-id: d05f8207bd86263ecdc93d3838e948c4fe90359a --- qpms/scatsystem.c | 40 ++++++++++++++++++++-------------------- qpms/scatsystem.h | 16 ++++++++++++++-- 2 files changed, 34 insertions(+), 22 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 02e04da..630f6bf 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1022,8 +1022,8 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full( complex double *tmp; QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double)); memset(target, 0, SQ(full_len) * sizeof(complex double)); //unnecessary? - complex double one = 1, zero = 0; - { // Non-diagonal part; M[piR, piC] = T[piR] S(piR<-piC) + const complex double zero = 0, minusone = -1; + { // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC) size_t fullvec_offsetR = 0; for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) { const qpms_vswf_set_spec_t *bspecR = ss->tm[ss->p[piR].tmatrix_id]->spec; @@ -1041,7 +1041,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full( k, posR, posC, QPMS_HANKEL_PLUS)); cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, - &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, + &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, target + fullvec_offsetR*full_len + fullvec_offsetC /*c*/, full_len /*ldc*/); @@ -1051,8 +1051,8 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full( fullvec_offsetR += bspecR->n; } } - // diagonal part M[pi,pi] = -1 - for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] = -1; + // diagonal part M[pi,pi] = +1 + for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] = +1; free(tmp); return target; @@ -1076,7 +1076,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( // some of the following workspaces are probably redundant; TODO optimize later. // workspaces for the uncompressed particle<-particle tranlation matrix block - // and the result of multiplying with a T-matrix + // and the result of multiplying with a T-matrix (times -1) complex double *Sblock, *TSblock; QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn)); QPMS_CRASHING_MALLOC(TSblock, sizeof(complex double)*SQ(ss->max_bspecn)); @@ -1085,7 +1085,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( complex double *tmp; QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); - const complex double one = 1, zero = 0; + const complex double one = 1, zero = 0, minusone = -1; 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; @@ -1134,13 +1134,13 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, - &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, + &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, Sblock/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, TSblock /*c*/, bspecC->n /*ldc*/); - } else { // diagonal, fill with diagonal -1 + } else { // diagonal, fill with diagonal +1 for (size_t row = 0; row < bspecR->n; ++row) for (size_t col = 0; col < bspecC->n; ++col) - TSblock[row * bspecC->n + col] = (row == col)? -1 : 0; + TSblock[row * bspecC->n + col] = (row == col)? +1 : 0; } // tmp[oiR|piR,piC] = ∑_K M[piR,K] U*[K,piC] @@ -1186,7 +1186,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( // some of the following workspaces are probably redundant; TODO optimize later. // workspaces for the uncompressed particle<-particle tranlation matrix block - // and the result of multiplying with a T-matrix + // and the result of multiplying with a T-matrix (times -1) complex double *Sblock, *TSblock; QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn)); QPMS_CRASHING_MALLOC(TSblock, sizeof(complex double)*SQ(ss->max_bspecn)); @@ -1195,7 +1195,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( complex double *tmp; QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); - const complex double one = 1, zero = 0; + const complex double one = 1, zero = 0, minusone = -1; for(qpms_ss_pi_t opistartR = 0; opistartR < ss->p_count; opistartR += ss->orbit_types[ss->p_orbitinfo[ss->p_by_orbit[opistartR]].t].size //orbit_p_countR; might write a while() instead @@ -1253,13 +1253,13 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR( cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, - &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, + &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, Sblock/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, TSblock /*c*/, bspecC->n /*ldc*/); - } else { // diagonal, fill with diagonal -1 + } else { // diagonal, fill with diagonal +1 for (size_t row = 0; row < bspecR->n; ++row) for (size_t col = 0; col < bspecC->n; ++col) - TSblock[row * bspecC->n + col] = (row == col)? -1 : 0; + TSblock[row * bspecC->n + col] = (row == col)? +1 : 0; } // tmp[oiR|piR,piC] = ∑_K M[piR,K] U*[K,piC] @@ -1308,7 +1308,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread // some of the following workspaces are probably redundant; TODO optimize later. // workspaces for the uncompressed particle<-particle tranlation matrix block - // and the result of multiplying with a T-matrix + // and the result of multiplying with a T-matrix (times -1) complex double *Sblock, *TSblock; QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn)); QPMS_CRASHING_MALLOC(TSblock, sizeof(complex double)*SQ(ss->max_bspecn)); @@ -1317,7 +1317,7 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread complex double *tmp; QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); - const complex double one = 1, zero = 0; + const complex double one = 1, zero = 0, minusone = -1; while(1) { // In the beginning, pick a target (row) orbit for this thread @@ -1385,13 +1385,13 @@ static void *qpms_scatsys_build_modeproblem_matrix_irrep_packed_parallelR_thread SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, - &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, + &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, Sblock/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, TSblock /*c*/, bspecC->n /*ldc*/); - } else { // diagonal, fill with diagonal -1 + } else { // diagonal, fill with diagonal +1 for (size_t row = 0; row < bspecR->n; ++row) for (size_t col = 0; col < bspecC->n; ++col) - TSblock[row * bspecC->n + col] = (row == col)? -1 : 0; + TSblock[row * bspecC->n + col] = (row == col)? +1 : 0; } // tmp[oiR|piR,piC] = ∑_K M[piR,K] U*[K,piC] diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index a3a88f1..2daddb0 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -111,6 +111,7 @@ typedef struct qpms_ss_orbit_type_t { qpms_ss_pi_t p_offset; } 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. @@ -265,14 +266,14 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( qpms_bessel_t J ); -/// Creates the full \f$ (-I + TS) \f$ matrix of the scattering system. +/// Creates the full \f$ (I - TS) \f$ matrix of the scattering system. complex double *qpms_scatsys_build_modeproblem_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_t *ss, /*COMPLEXIFY*/double k ///< Wave number to use in the translation matrix. ); -/// Creates the mode problem matrix \f$ (-I + TS) \f$ directly in the irrep-packed form. +/// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form. complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, @@ -294,6 +295,16 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorder_pa /*COMPLEXIFY*/double k ///< Wave number to use in the translation matrix. ); +/// LU factorisation (LAPACKE_zgetrf) result holder. +typedef struct qpms_ss_LU { + /// LU decomposition array. + complex double *a; + /// Pivot index array, size at least max(1,min(m, n)). + lapack_int *ipiv; +} qpms_ss_LU; + +void qpms_ss_LU_free(qpms_ss_LU *); + /// NOT IMPLEMENTED Dumps a qpms_scatsys_t structure to a file. qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path); @@ -421,4 +432,5 @@ ccart3_t qpms_scatsys_eval_E_irrep(const qpms_scatsys_t *ss, complex double k ///< Wave number. ); #endif + #endif //QPMS_SCATSYSTEM_H