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

scatsystem: mode problem matrix for xy-plane periodic systems 

Former-commit-id: 755bb86f966b031861a2338eb97e68111694b1d8
This commit is contained in:
Marek Nečada 2020-02-28 13:53:53 +02:00
parent 2b628736f0
commit 1fb4e5760a
2 changed files with 97 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

135
qpms/scatsystem.c
View File

@ -43,19 +43,19 @@ void qpms_scatsystem_set_nthreads(long n) {
qpms_scatsystem_nthreads_override = n; qpms_scatsystem_nthreads_override = n;
} }
static inline void qpms_ss_ensure_periodic(qpms_scatsys_t *ss) { static inline void qpms_ss_ensure_periodic(const qpms_scatsys_t *ss) {
QPMS_ENSURE(ss->lattice_dimension > 0, "This method is applicable only to periodic systems."); QPMS_ENSURE(ss->lattice_dimension > 0, "This method is applicable only to periodic systems.");
} }
static inline void qpms_ss_ensure_periodic_a(qpms_scatsys_t *ss, const char *s) { static inline void qpms_ss_ensure_periodic_a(const qpms_scatsys_t *ss, const char *s) {
QPMS_ENSURE(ss->lattice_dimension > 0, "This method is applicable only to periodic systems. Use %s instead.", s); QPMS_ENSURE(ss->lattice_dimension > 0, "This method is applicable only to periodic systems. Use %s instead.", s);
} }
static inline void qpms_ss_ensure_nonperiodic(qpms_scatsys_t *ss) { static inline void qpms_ss_ensure_nonperiodic(const qpms_scatsys_t *ss) {
QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems."); QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems.");
} }
static inline void qpms_ss_ensure_nonperiodic_a(qpms_scatsys_t *ss, const char *s) { static inline void qpms_ss_ensure_nonperiodic_a(const qpms_scatsys_t *ss, const char *s) {
QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems. Use %s instead.", s); QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems. Use %s instead.", s);
} }
@ -1145,6 +1145,20 @@ complex double *qpms_scatsys_build_translation_matrix_full(
target, ss, k, QPMS_HANKEL_PLUS); target, ss, k, QPMS_HANKEL_PLUS);
} }
complex double *qpms_scatsyswk_build_translation_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
const qpms_scatsys_at_omega_k_t *sswk
)
{
const qpms_scatsys_at_omega_t *ssw = sswk->ssw;
const complex double wavenumber = ssw->wavenumber;
const qpms_scatsys_t *ss = ssw->ss;
qpms_ss_ensure_periodic(ss);
const cart3_t k_cart3 = cart3_from_double_array(sswk->k);
return qpms_scatsys_periodic_build_translation_matrix_full(target, ss, wavenumber, &k_cart3);
}
complex double *qpms_scatsys_build_translation_matrix_e_full( complex double *qpms_scatsys_build_translation_matrix_e_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,
@ -1184,9 +1198,8 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
return target; return target;
} }
static inline int qpms_ss_ppair_W32xy(qpms_scatsys_t *ss, static inline int qpms_ss_ppair_W32xy(const qpms_scatsys_t *ss,
qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, const cart2_t kvector,
const cart2_t kvector,
complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride, complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride,
qpms_ewald_part parts) { qpms_ewald_part parts) {
const qpms_vswf_set_spec_t *srcspec = qpms_ss_bspec_pi(ss, psrc); const qpms_vswf_set_spec_t *srcspec = qpms_ss_bspec_pi(ss, psrc);
@ -1198,11 +1211,25 @@ static inline int qpms_ss_ppair_W32xy(qpms_scatsys_t *ss,
return qpms_trans_calculator_get_trans_array_e32_e(ss->c, return qpms_trans_calculator_get_trans_array_e32_e(ss->c,
target, NULL /*err*/, destspec, deststride, srcspec, srcstride, target, NULL /*err*/, destspec, deststride, srcspec, srcstride,
ss->eta, wavenumber, ss->per->eta, wavenumber,
cart3xy2cart2(ss->per->lattice_basis[0]), cart3xy2cart2(ss->per->lattice_basis[1]), cart3xy2cart2(ss->per->lattice_basis[0]), cart3xy2cart2(ss->per->lattice_basis[1]),
kvector, kvector,
cart2_substract(cart3xy2cart2(ss->p[pdest].pos), cart3xy2cart2(ss->p[psrc].pos)) cart2_substract(cart3xy2cart2(ss->p[pdest].pos), cart3xy2cart2(ss->p[psrc].pos)),
u->maxR, u->maxK, parts); maxR, maxK, parts);
}
static inline int qpms_ss_ppair_W(const qpms_scatsys_t *ss,
qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, const double wavevector[],
complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride,
qpms_ewald_part parts) {
if(ss->lattice_dimension == 2 && // Currently, we can only the xy-plane
!ss->per->lattice_basis[0].z && !ss->per->lattice_basis[1].z &&
!wavevector[2])
return qpms_ss_ppair_W32xy(ss, pdest, psrc, wavenumber, cart2_from_double_array(wavevector),
target + deststride * ss->fecv_pstarts[pdest] + srcstride * ss->fecv_pstarts[psrc],
deststride, srcstride, parts);
else
QPMS_NOT_IMPLEMENTED("Only 2D xy-lattices currently supported");
} }
complex double *qpms_scatsys_periodic_build_translation_matrix_full( complex double *qpms_scatsys_periodic_build_translation_matrix_full(
@ -1216,7 +1243,7 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full(
const ptrdiff_t deststride = ss->fecv_size, srcstride = 1; const ptrdiff_t deststride = ss->fecv_size, srcstride = 1;
// We have some limitations in the current implementation // We have some limitations in the current implementation
if(ss->lattice_dimension == 2 && // Currently, we can only the xy-plane if(ss->lattice_dimension == 2 && // Currently, we can only the xy-plane
!ss->per->lattice_basis[0].z && !ss->per_lattice_basis[1].z && !ss->per->lattice_basis[0].z && !ss->per->lattice_basis[1].z &&
!wavevector->z) { !wavevector->z) {
for (qpms_ss_pi_t pd = 0; pd < ss->p_count; ++pd) for (qpms_ss_pi_t pd = 0; pd < ss->p_count; ++pd)
for (qpms_ss_pi_t ps = 0; ps < ss->p_count; ++ps) { for (qpms_ss_pi_t ps = 0; ps < ss->p_count; ++ps) {
@ -1229,21 +1256,23 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full(
return target; return target;
} }
complex double *qpms_scatsysw_build_modeproblem_matrix_full( // Common implementation of qpms_scatsysw[k]_build_modeproblem_matrix_full
static inline complex double *qpms_scatsysw_scatsyswk_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_at_omega_t *ssw const qpms_scatsys_at_omega_t *ssw,
const double k[] // NULL if non-periodic
) )
{ {
const complex double k = ssw->wavenumber; const complex double wavenumber = ssw->wavenumber;
const qpms_scatsys_t *ss = ssw->ss; const qpms_scatsys_t *ss = ssw->ss;
qpms_ss_ensure_nonperiodic(ss); qpms_ss_ensure_nonperiodic(ss);
const size_t full_len = ss->fecv_size; const size_t full_len = ss->fecv_size;
if(!target) if(!target)
QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
complex double *tmp; 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)); //unnecessary? 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;
@ -1255,49 +1284,55 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full(
const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m; const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m;
for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) {
const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec; const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec;
if(piC != piR) { // The diagonal will be dealt with later. if (k == NULL) { // non-periodic case
const cart3_t posC = ss->p[piC].pos; if(piC != piR) { // No "self-interaction" in non-periodic case
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, const cart3_t posC = ss->p[piC].pos;
tmp, // tmp is S(piR<-piC) QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c,
bspecR, bspecC->n, bspecC, 1, tmp, // tmp is S(piR<-piC)
k, posR, posC, QPMS_HANKEL_PLUS)); bspecR, bspecC->n, bspecC, 1,
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, wavenumber, posR, posC, QPMS_HANKEL_PLUS));
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, }
&minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, } else { // periodic case
tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, QPMS_ENSURE_SUCCESS(qpms_ss_ppair_W(ss, piR, piC, wavenumber, k,
target + fullvec_offsetR*full_len + fullvec_offsetC /*c*/, tmp /*target*/, bspecC->n /*deststride*/, 1 /*srcstride*/,
full_len /*ldc*/); QPMS_EWALD_FULL));
} }
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/,
&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*/);
fullvec_offsetC += bspecC->n; fullvec_offsetC += bspecC->n;
} }
fullvec_offsetR += bspecR->n; 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; // Add the identity, diagonal part M[pi,pi] += 1
for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] += 1;
free(tmp); free(tmp);
return target; return target;
} }
complex double *qpms_scatsyswk_build_modeproblem_matrix_full( complex double *qpms_scatsysw_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_at_omega_t *ssw) {
complex double *target, qpms_ss_ensure_nonperiodic_a(ssw->ss, "qpms_scatsyswk_build_modeproblem_matrix_full()");
const qpms_scatsys_at_omega_k_t *sswk return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full(
) target, ssw, NULL);
{
const qpms_scatsys_at_omega_t *ssw = sswk->ssw;
const complex double k = ssw->wavenumber;
const qpms_scatsys_t *ss = ssw->ss;
qpms_ss_ensure_periodic(ss);
const size_t full_len = ss->fecv_size;
if(!target)
QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
//////////// TODO ////////////////
QPMS_NOT_IMPLEMENTED("TODO");
return target;
} }
complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
complex double *target, const qpms_scatsys_at_omega_k_t *sswk)
{
qpms_ss_ensure_periodic_a(sswk->ssw->ss, "qpms_scatsysw_build_modeproblem_matrix_full()");
return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full(target, sswk->ssw, sswk->k);
}
// Serial reference implementation. // Serial reference implementation.
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
/// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated.
@ -1990,7 +2025,7 @@ void qpms_ss_LU_free(qpms_ss_LU lu) {
} }
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmatrix_full, qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmatrix_full,
int *target_piv, const qpms_scatsys_at_omega_t *ssw, const qpms_scatsys_at_k_omega_t *sswk) { int *target_piv, const qpms_scatsys_at_omega_t *ssw, const qpms_scatsys_at_omega_k_t *sswk) {
if (sswk) { if (sswk) {
QPMS_ASSERT(sswk->ssw == ssw || !ssw); QPMS_ASSERT(sswk->ssw == ssw || !ssw);
ssw = sswk->ssw; ssw = sswk->ssw;
@ -2016,7 +2051,7 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmat
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed, qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed,
int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) { int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) {
QPMS_ENSURE(mpmatrix_packed, "A non-NULL pointer to the pre-calculated mode matrix is required"); QPMS_ENSURE(mpmatrix_packed, "A non-NULL pointer to the pre-calculated mode matrix is required");
qpms_scatsys_ensure_nonperiodic(ssw->ss); qpms_ss_ensure_nonperiodic(ssw->ss);
size_t n = ssw->ss->saecv_sizes[iri]; size_t n = ssw->ss->saecv_sizes[iri];
if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, n * sizeof(int)); if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, n * sizeof(int));
QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n,
@ -2033,7 +2068,7 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(complex doubl
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
complex double *target, int *target_piv, complex double *target, int *target_piv,
const qpms_scatsys_at_omega_t *ssw){ const qpms_scatsys_at_omega_t *ssw){
qpms_scatsys_ensure_nonperiodic_a(ssw->ss, "qpms_scatsyswk_build_modeproblem_matrix_full_LU()"); qpms_ss_ensure_nonperiodic_a(ssw->ss, "qpms_scatsyswk_build_modeproblem_matrix_full_LU()");
target = qpms_scatsysw_build_modeproblem_matrix_full(target, ssw); target = qpms_scatsysw_build_modeproblem_matrix_full(target, ssw);
return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, ssw, NULL); return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, ssw, NULL);
} }

20
qpms/scatsystem.h
View File

@ -356,7 +356,7 @@ complex double *qpms_scatsys_build_translation_matrix_full(
complex double k ///< Wave number to use in the translation matrix. complex double k ///< Wave number to use in the translation matrix.
); );
/// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system. NI /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
/** /**
* \returns \a target on success, NULL on error. * \returns \a target on success, NULL on error.
*/ */
@ -422,10 +422,11 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
); );
struct qpms_scatsys_at_omega_k_t; // Defined below.
/// LU factorisation (LAPACKE_zgetrf) result holder. /// LU factorisation (LAPACKE_zgetrf) result holder.
typedef struct qpms_ss_LU { typedef struct qpms_ss_LU {
const qpms_scatsys_at_omega_t *ssw; const qpms_scatsys_at_omega_t *ssw;
const qpms_scatsys_at_omega_k_t *sswk; ///< Only for periodic systems, otherwise NULL. const struct qpms_scatsys_at_omega_k_t *sswk; ///< Only for periodic systems, otherwise NULL.
bool full; ///< true if full matrix; false if irrep-packed. bool full; ///< true if full matrix; false if irrep-packed.
qpms_iri_t iri; ///< Irrep index if `full == false`. qpms_iri_t iri; ///< Irrep index if `full == false`.
/// LU decomposition array. /// LU decomposition array.
@ -455,7 +456,7 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory. complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems. const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems.
const qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL. const struct qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL.
); );
/// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents. /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents.
@ -476,12 +477,15 @@ complex double *qpms_scatsys_scatter_solve(
// ======================= Periodic system -only related stuff ============================= // ======================= Periodic system -only related stuff =============================
/// Scattering system at a given frequency and k-vector. Used only with periodic systems. /// Scattering system at a given frequency and k-vector. Used only with periodic systems.
/**
* N.B. use as a stack variable now, but this might become heap-allocated in the future (with own con- and destructor)
*/
typedef struct qpms_scatsys_at_omega_k_t { typedef struct qpms_scatsys_at_omega_k_t {
const qpms_scatsys_at_omega_t *ssw; const qpms_scatsys_at_omega_t *ssw;
double k[3]; ///< The k-vector's cartesian coordinates. double k[3]; ///< The k-vector's cartesian coordinates.
} qpms_scatsys_at_omega_k_t; } qpms_scatsys_at_omega_k_t;
/// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system. NI /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
/** /**
* \returns \a target on success, NULL on error. * \returns \a target on success, NULL on error.
*/ */
@ -497,19 +501,19 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full(
complex double *target, complex double *target,
const qpms_scatsys_t *ss, const qpms_scatsys_t *ss,
complex double wavenumber, ///< Wave number to use in the translation matrix. complex double wavenumber, ///< Wave number to use in the translation matrix.
const double wavevector[] ///< Wavevector / pseudomomentum in cartesian coordinates. const cart3_t *wavevector ///< Wavevector / pseudomomentum in cartesian coordinates.
); );
/// Global translation matrix. NI /// Global translation matrix.
complex double *qpms_scatsyswk_build_translation_matrix_full( complex double *qpms_scatsyswk_build_translation_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_omega_k_t *sswk const qpms_scatsys_at_omega_k_t *sswk
); );
/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only. /// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only.
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_k_t *sswk const qpms_scatsys_at_omega_k_t *sswk