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

Reverse modeproblem matrix sign so it applies directly to scattering 

Former-commit-id: d05f8207bd86263ecdc93d3838e948c4fe90359a
This commit is contained in:
Marek Nečada 2019-07-23 08:43:25 +03:00
parent e05a79483d
commit 5758c5d587
2 changed files with 34 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
qpms/scatsystem.c
View File

@ -1022,8 +1022,8 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
complex double *tmp; complex double *tmp;
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)); //unnecessary? memset(target, 0, SQ(full_len) * sizeof(complex double)); //unnecessary?
complex double one = 1, zero = 0; 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;
for(qpms_ss_pi_t piR = 0; piR < ss->p_count; ++piR) { 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; 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)); k, posR, posC, QPMS_HANKEL_PLUS));
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 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*/, tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,
target + fullvec_offsetR*full_len + fullvec_offsetC /*c*/, target + fullvec_offsetR*full_len + fullvec_offsetC /*c*/,
full_len /*ldc*/); full_len /*ldc*/);
@ -1051,8 +1051,8 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full(
fullvec_offsetR += bspecR->n; fullvec_offsetR += bspecR->n;
} }
} }
// diagonal part M[pi,pi] = -1 // diagonal part M[pi,pi] = +1
for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] = -1; for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] = +1;
free(tmp); free(tmp);
return target; 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. // some of the following workspaces are probably redundant; TODO optimize later.
// workspaces for the uncompressed particle<-particle tranlation matrix block // 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; complex double *Sblock, *TSblock;
QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn)); QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn));
QPMS_CRASHING_MALLOC(TSblock, 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; complex double *tmp;
QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); 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 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; 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, cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 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*/, Sblock/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,
TSblock /*c*/, bspecC->n /*ldc*/); 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 row = 0; row < bspecR->n; ++row)
for (size_t col = 0; col < bspecC->n; ++col) 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] // 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. // some of the following workspaces are probably redundant; TODO optimize later.
// workspaces for the uncompressed particle<-particle tranlation matrix block // 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; complex double *Sblock, *TSblock;
QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn)); QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn));
QPMS_CRASHING_MALLOC(TSblock, 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; complex double *tmp;
QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); 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; 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 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, cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 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*/, Sblock/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,
TSblock /*c*/, bspecC->n /*ldc*/); 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 row = 0; row < bspecR->n; ++row)
for (size_t col = 0; col < bspecC->n; ++col) 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] // 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. // some of the following workspaces are probably redundant; TODO optimize later.
// workspaces for the uncompressed particle<-particle tranlation matrix block // 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; complex double *Sblock, *TSblock;
QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn)); QPMS_CRASHING_MALLOC(Sblock, sizeof(complex double)*SQ(ss->max_bspecn));
QPMS_CRASHING_MALLOC(TSblock, 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; complex double *tmp;
QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); 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) { while(1) {
// In the beginning, pick a target (row) orbit for this thread // 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, SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans,
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 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*/, Sblock/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,
TSblock /*c*/, bspecC->n /*ldc*/); 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 row = 0; row < bspecR->n; ++row)
for (size_t col = 0; col < bspecC->n; ++col) 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] // tmp[oiR|piR,piC] = ∑_K M[piR,K] U*[K,piC]

16
qpms/scatsystem.h
View File

@ -111,6 +111,7 @@ typedef struct qpms_ss_orbit_type_t {
qpms_ss_pi_t p_offset; qpms_ss_pi_t p_offset;
} qpms_ss_orbit_type_t; } qpms_ss_orbit_type_t;
typedef ptrdiff_t qpms_ss_osn_t; ///< "serial number" of av orbit in a given type. 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. /// 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 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( complex double *qpms_scatsys_build_modeproblem_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, complex double *target,
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
/*COMPLEXIFY*/double k ///< Wave number to use in the translation matrix. /*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( 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. /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target, 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. /*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. /// NOT IMPLEMENTED Dumps a qpms_scatsys_t structure to a file.
qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path); 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. complex double k ///< Wave number.
); );
#endif #endif
#endif //QPMS_SCATSYSTEM_H #endif //QPMS_SCATSYSTEM_H