diff --git a/qpms/scatsys_private.h b/qpms/scatsys_private.h index b6e2b09..fd6e68a 100644 --- a/qpms/scatsys_private.h +++ b/qpms/scatsys_private.h @@ -1,11 +1,19 @@ #ifndef QPMS_SCATSYS_PRIVATE_H #define QPMS_SCATSYS_PRIVATE_H +/* + * This file includes some definitions shared between scatsystem.c + * and scatsys_translation_booster.c that are not needed anywhere + * else. + */ #include "scatsystem.h" +#include + complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted( complex double *target, const qpms_scatsys_at_omega_t *ssw); - +complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_boosted( + complex double *target_packed, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri); /// "private" destructor, called by qpms_scatsys_free() void qpms_scatsys_translation_booster_free(struct qpms_scatsys_translation_booster *); /// "private" constructor, use qpms_ss_create_translation_cache() instead. @@ -18,4 +26,16 @@ qpms_scatsysw_translation_booster_create(const qpms_scatsys_at_omega_t *ssw); /// "private" destructor, called by qpms_scatsys_at_omega_free() void qpms_scatsysw_translation_booster_free(struct qpms_scatsysw_translation_booster *); -#endif + + +struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{ + const qpms_scatsys_at_omega_t *ssw; + qpms_ss_pi_t *opistartR_ptr; + pthread_mutex_t *opistartR_mutex; + qpms_iri_t iri; + complex double *target_packed; +}; + +void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_boosted(void *arg); + +#endif //QPMS_SCATSYS_PRIVATE_H diff --git a/qpms/scatsys_translation_booster.c b/qpms/scatsys_translation_booster.c index 1a6be26..80d4fc2 100644 --- a/qpms/scatsys_translation_booster.c +++ b/qpms/scatsys_translation_booster.c @@ -10,10 +10,18 @@ #include "vectors.h" #include "qpms_error.h" #include "qpms_specfunc.h" +#include "groups.h" #define SQ(x) ((x)*(x)) #define MAX(x, y) (((x) > (y)) ? (x) : (y)) +#ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS +#include "qpmsblas.h" +#define SERIAL_ZGEMM qpms_zgemm +#else +#define SERIAL_ZGEMM cblas_zgemm +#endif + typedef size_t ppid_t; typedef size_t uoppid_t; @@ -140,7 +148,7 @@ void qpms_scatsys_translation_booster_free(booster_t *b) { } int qpms_ss_create_translation_cache(qpms_scatsys_t *ss, qpms_ss_caching_mode_t m) { -QPMS_ASSERT(ss); + QPMS_ASSERT(ss); if (ss->tbooster) { QPMS_WARN("Translation cache already created?"); return 0; @@ -154,8 +162,8 @@ QPMS_ASSERT(ss); ss->tbooster = qpms_scatsys_translation_booster_create(ss); if (ss->tbooster) return 0; else { - QPMS_WARN("Failed to create tranlation operator cache"); - return -1; + QPMS_WARN("Failed to create tranlation operator cache"); + return -1; } default: QPMS_WTF; @@ -225,6 +233,10 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted( QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); complex double *tmp; QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double)); + // Workspaces for the translation operator A and B matrices + complex double *A, *B; + QPMS_CRASHING_MALLOC(A, SQ(c->nelem) * sizeof(*A)); + QPMS_CRASHING_MALLOC(B, SQ(c->nelem) * sizeof(*B)); memset(target, 0, SQ(full_len) * sizeof(complex double)); //unnecessary? double legendre_buf[gsl_sf_legendre_array_n(2*c->lMax + 1)]; //VLA, workspace for legendre arrays const complex double zero = 0, minusone = -1; @@ -251,18 +263,16 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted( QPMS_PARANOID_ASSERT(c->normalisation == bspecC->norm && c->normalisation == bspecR->norm); QPMS_PARANOID_ASSERT(c->lMax >= bspecC->lMax && c->lMax >= bspecR->lMax); QPMS_PARANOID_ASSERT(bspecC->lMax_L < 0 && bspecR->lMax_L < 0); - complex double A[pair_nelem][pair_nelem]; // VLA, TODO flatten - complex double B[pair_nelem][pair_nelem]; // VLA, TODO flatten { // this replaces qpms_trans_calculator_get_AB_arrays() and ..._buf() const double costheta = cos(dlj.theta); QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE, - 2*c->lMax+1, costheta, -1, legendre_buf)); + 2*pair_lMax+1, costheta, -1, legendre_buf)); const double * const legendres = legendre_buf; const complex double * const bessels = bw->bessels + b->bessel_offsets_r[ri]; - qpms_trans_calculator_get_AB_arrays_precalcbuf(c, pair_lMax, A[0], B[0], + qpms_trans_calculator_get_AB_arrays_precalcbuf(c, pair_lMax, A, B, /*deststride*/ pair_nelem, /*srcstride*/ 1, dlj.phi, bessels, legendres); qpms_trans_array_from_AB(tmp,// tmp is S(piR<-piC) - bspecR, bspecC->n, bspecC, 1, A[0], B[0], pair_lMax); + bspecR, bspecC->n, bspecC, 1, A, B, pair_lMax); } } cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, @@ -281,5 +291,170 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full_boosted( for (size_t i = 0; i < full_len; ++i) target[full_len * i + i] = +1; free(tmp); + free(A); + free(B); return target; } + +void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_boosted(void *arg) +{ + const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg + *a = arg; + const qpms_scatsys_at_omega_t *ssw = a->ssw; + const complex double k = ssw->wavenumber; + const qpms_scatsys_t *ss = ssw->ss; + const qpms_iri_t iri = a->iri; + const size_t packedlen = ss->saecv_sizes[iri]; + + QPMS_ASSERT(ssw->translation_cache && ssw->ss->tbooster); + const booster_t *const b = ss->tbooster; + const boosterw_t *const bw = ssw->translation_cache; + const qpms_trans_calculator *const c = ss->c; + + // 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 (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)); + // Workspaces for the translation operator A and B matrices + complex double *A, *B; + QPMS_CRASHING_MALLOC(A, SQ(c->nelem) * sizeof(*A)); + QPMS_CRASHING_MALLOC(B, SQ(c->nelem) * sizeof(*B)); + double legendre_buf[gsl_sf_legendre_array_n(2*c->lMax + 1)]; //VLA, workspace for legendre arrays + + // Workspace for the intermediate particle-orbit matrix result + complex double *tmp; + QPMS_CRASHING_MALLOC(tmp, sizeof(complex double) * SQ(ss->max_bspecn) * ss->sym->order); + + const complex double one = 1, zero = 0, minusone = -1; + + while(1) { + // In the beginning, pick a target (row) orbit for this thread + QPMS_ENSURE_SUCCESS(pthread_mutex_lock(a->opistartR_mutex)); + if(*(a->opistartR_ptr) >= ss->p_count) {// Everything is already done, end + QPMS_ENSURE_SUCCESS(pthread_mutex_unlock(a->opistartR_mutex)); + break; + } + const qpms_ss_pi_t opistartR = *(a->opistartR_ptr); + // Now increment it for another thread: + *(a->opistartR_ptr) += ss->orbit_types[ss->p_orbitinfo[ss->p_by_orbit[opistartR]].t].size; + QPMS_ENSURE_SUCCESS(pthread_mutex_unlock(a->opistartR_mutex)); + + // Orbit picked (defined by opistartR), do the work. + const qpms_ss_pi_t orbitstartpiR = ss->p_by_orbit[opistartR]; + const qpms_ss_oti_t otiR = ss->p_orbitinfo[orbitstartpiR].t; + const qpms_ss_osn_t osnR = ss->p_orbitinfo[orbitstartpiR].osn; + const qpms_ss_orbit_type_t *const otR = ss->orbit_types + otiR; + const qpms_ss_orbit_pi_t orbit_p_countR = otR->size; + const size_t orbit_packedsizeR = otR->irbase_sizes[iri]; + + if(orbit_packedsizeR) { // avoid zgemm crash on empty irrep + const size_t particle_fullsizeR = otR->bspecn; // == bspecR->n + const qpms_vswf_set_spec_t *bspecR = ssw->tm[ss->p[orbitstartpiR].tmatrix_id]->spec; + // This is the orbit-level matrix projecting the whole orbit onto the irrep. + const complex double *omR = otR->irbases + otR->irbase_offsets[iri]; + // Orbit coeff vector's full size: + const size_t orbit_fullsizeR = otR->size * otR->bspecn; + // This is where the orbit starts in the "packed" vector: + const size_t packed_orbit_offsetR = + ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiR] + + osnR * otR->irbase_sizes[iri]; + for(qpms_ss_orbit_pi_t opiR = 0; opiR < orbit_p_countR; ++opiR) { + for(qpms_ss_orbit_pi_t opiR = 0; opiR < orbit_p_countR; ++opiR) { + qpms_ss_pi_t piR = ss->p_by_orbit[opistartR + opiR]; + assert(opiR == ss->p_orbitinfo[piR].p); + assert(otiR == ss->p_orbitinfo[piR].t); + assert(ss->p_orbitinfo[piR].osn == osnR); + const cart3_t posR = ss->p[piR].pos; + // dest particle T-matrix + const complex double *tmmR = ssw->tm[ss->p[piR].tmatrix_id]->m; + for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { //Column loop + const qpms_ss_oti_t otiC = ss->p_orbitinfo[piC].t; + const qpms_ss_orbit_type_t *const otC = ss->orbit_types + otiC; + const qpms_ss_osn_t osnC = ss->p_orbitinfo[piC].osn; + const qpms_ss_orbit_pi_t opiC = ss->p_orbitinfo[piC].p; + // This is where the particle's orbit starts in the "packed" vector: + const size_t packed_orbit_offsetC = + ss->saecv_ot_offsets[iri*ss->orbit_type_count + otiC] + + osnC * otC->irbase_sizes[iri]; + const qpms_vswf_set_spec_t *bspecC = ssw->tm[ss->p[piC].tmatrix_id]->spec; + // Orbit coeff vector's full size: + const size_t orbit_fullsizeC = otC->size * otC->bspecn; + const size_t particle_fullsizeC = otC->bspecn; // == bspecC->n + const size_t orbit_packedsizeC = otC->irbase_sizes[iri]; + // This is the orbit-level matrix projecting the whole orbit onto the irrep. + const complex double *omC = otC->irbases + otC->irbase_offsets[iri]; + + if(orbit_packedsizeC) { // avoid zgemm crash on empty irrep + if(piC != piR) { // non-diagonal, calculate TS + const cart3_t posC = ss->p[piC].pos; +#if 0 + QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, + Sblock, // Sblock is S(piR->piC) + bspecR, bspecC->n, bspecC, 1, + k, posR, posC, QPMS_HANKEL_PLUS)); +#endif + { // this block replaces qpms_trans_calculator_get_trans_array(): + // R is dest, C is src + const sph_t dlj = cart2sph(cart3_substract(posR, posC)); + const uoppid_t pid = uopairid(ss->p_count, piC, piR); + const size_t ri = b->r_map[pid]; + QPMS_PARANOID_ASSERT(dlj.r == b->r[ri]); + const qpms_l_t pair_lMax = b->lMax_r[ri]; + const qpms_y_t pair_nelem = qpms_lMax2nelem(pair_lMax); + + { // this replaces qpms_trans_calculator_get_AB_arrays() and _buf() + const double costheta = cos(dlj.theta); + QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE, + 2*pair_lMax+1, costheta, -1, legendre_buf)); + const double * const legendres = legendre_buf; + const complex double * const bessels = bw->bessels + b->bessel_offsets_r[ri]; + qpms_trans_calculator_get_AB_arrays_precalcbuf(c, pair_lMax, A, B, + /*deststride*/ pair_nelem, /*srcstride*/ 1, dlj.phi, bessels, legendres); + } + qpms_trans_array_from_AB(Sblock, // Sblock is S(piR->piC) + bspecR, bspecC->n, bspecC, 1, A, B, pair_lMax); + } + + SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans, + bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, + &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 + 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; + } + + // tmp[oiR|piR,piC] = ∑_K M[piR,K] U*[K,piC] + SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasConjTrans, + particle_fullsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeC /*K*/, + &one /*alpha*/, TSblock/*A*/, particle_fullsizeC/*ldA*/, + omC + opiC*particle_fullsizeC /*B*/, + orbit_fullsizeC/*ldB*/, &zero /*beta*/, + tmp /*C*/, orbit_packedsizeC /*LDC*/); + + // target[oiR|piR,oiC|piC] += U[...] tmp[...] + SERIAL_ZGEMM(CblasRowMajor, CblasNoTrans, CblasNoTrans, + orbit_packedsizeR /*M*/, orbit_packedsizeC /*N*/, particle_fullsizeR /*K*/, + &one /*alpha*/, omR + opiR*particle_fullsizeR/*A*/, orbit_fullsizeR/*ldA*/, + tmp /*B*/, orbit_packedsizeC /*ldB*/, &one /*beta*/, + a->target_packed + packedlen*packed_orbit_offsetR + packed_orbit_offsetC /*C*/, + packedlen /*ldC*/); + } + } + } + } + } + } + free(tmp); + free(A); + free(B); + free(Sblock); + free(TSblock); + return NULL; +} + diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 283dca1..9d8944a 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1403,14 +1403,6 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( return target_packed; } -struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg{ - const qpms_scatsys_at_omega_t *ssw; - qpms_ss_pi_t *opistartR_ptr; - pthread_mutex_t *opistartR_mutex; - qpms_iri_t iri; - complex double *target_packed; -}; - static void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread(void *arg) { const struct qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_arg @@ -1529,9 +1521,6 @@ static void *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_threa } } } - - - } free(tmp); free(Sblock); @@ -1733,6 +1722,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ) { + const _Bool use_translation_cache = (ssw->translation_cache != NULL); const size_t packedlen = ssw->ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; @@ -1765,7 +1755,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( pthread_t thread_ids[nthreads]; for(long thi = 0; thi < nthreads; ++thi) QPMS_ENSURE_SUCCESS(pthread_create(thread_ids + thi, NULL, - qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread, + use_translation_cache + ? qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread_boosted + : qpms_scatsysw_build_modeproblem_matrix_irrep_packed_parallelR_thread, (void *) &arg)); for(long thi = 0; thi < nthreads; ++thi) { void *retval;