diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index bffd0a4..eae1737 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -700,7 +700,7 @@ cdef class BaseSpec: if 'norm' in kwargs.keys(): self.s.norm = kwargs['norm'] else: - self.s.norm = QPMS_NORMALISATION_POWER + self.s.norm = QPMS_NORMALISATION_POWER_CS # set the other metadata cdef qpms_l_t l self.s.lMax_L = -1 diff --git a/qpms/tmatrix_io.c b/qpms/tmatrix_io.c index 2108de9..ac77e26 100644 --- a/qpms/tmatrix_io.c +++ b/qpms/tmatrix_io.c @@ -33,9 +33,9 @@ qpms_errno_t qpms_read_scuff_tmatrix( if (!(freqs && n && tmdata)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, "freqs, n, and tmdata are mandatory arguments and must not be NULL."); - if (bs->norm != QPMS_NORMALISATION_POWER) // CHECKME CORRECT? + if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT? qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "Not implemented; only QPMS_NORMALISATION_POWER (CHECKME)" + "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)" " norm supported right now."); int n_alloc = 128; // First chunk to allocate *n = 0; diff --git a/setup.py b/setup.py index 9993d9d..51696a0 100755 --- a/setup.py +++ b/setup.py @@ -31,11 +31,14 @@ qpms_c = Extension('qpms_c', 'qpms/error.c'], extra_compile_args=['-std=c99','-ggdb', '-O0', '-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required + #'-DAN2', + #'-DBN2', #'-DQPMS_USE_OMP', '-DDISABLE_NDEBUG', # uncomment to enable assertions in the modules #'-fopenmp', ], - libraries=['gsl', 'lapacke', 'blas', 'gslcblas', #'omp' + libraries=['gsl', 'lapacke', 'cblas_LINUX', 'blas_LINUX', + 'gslcblas', #'omp' # TODO resolve the problem with openblas (missing gotoblas symbol) and preferable use other blas library ], runtime_library_dirs=os.environ['LD_LIBRARY_PATH'].split(':') if 'LD_LIBRARY_PATH' in os.environ else [] diff --git a/tests/fixproj-notes b/tests/fixproj-notes new file mode 100644 index 0000000..6b14a7e --- /dev/null +++ b/tests/fixproj-notes @@ -0,0 +1,4 @@ +c99 -g -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss1.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -llapacke ~/repo/CBLAS/lib/cblas_LINUX.a ~/repo/BLAS-3.8.0/blas_LINUX.a + + +LD_LIBRARY_PATH=$HOME/repo/CBLAS/lib/ LDFLAGS=-L../CBLAS/lib/ python3 setup.py install --user diff --git a/tests/quaternion_O3.c b/tests/quaternion_O3.c index 7b880b3..da234eb 100644 --- a/tests/quaternion_O3.c +++ b/tests/quaternion_O3.c @@ -3,14 +3,18 @@ int main() { cart3_t v = {1, 2, 3}; - qpms_quat4d_t q4[4] = { + qpms_quat4d_t q4[8] = { {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, + {-1, 0, 0, 0}, + {0, -1, 0, 0}, + {0, 0, -1, 0}, + {0, 0, 0, -1}, }; printf("original: (%g, %g, %g)\n", v.x, v.y, v.z); - for(int i = 0; i < 4; ++i) { + for(int i = 0; i < 8; ++i) { for (int det = -1; det < 2; det += 2) { const qpms_quat4d_t qr = q4[i]; const qpms_quat_t qc = qpms_quat_2c_from_4d(qr); diff --git a/tests/scatsys.py b/tests/scatsys.py index 7206cd2..eedb12e 100644 --- a/tests/scatsys.py +++ b/tests/scatsys.py @@ -2,6 +2,8 @@ from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSyste from qpms.symmetries import point_group_info import numpy as np +np.random.seed(444) + sym = FinitePointGroup(point_group_info['D3h']) bspec2 = BaseSpec(lMax=2) bspec1 = BaseSpec(lMax=1) diff --git a/tests/ss_syms_packs.c b/tests/ss_syms_packs.c index 470fc5f..1af2c09 100644 --- a/tests/ss_syms_packs.c +++ b/tests/ss_syms_packs.c @@ -35,20 +35,27 @@ int main() for(qpms_l_t l = 1; l <= 1; ++l) for (qpms_m_t m = -l; m <= l; ++m) qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); - for(qpms_l_t l = 1; l <= 2; ++l) + for(qpms_l_t l = 1; l <= 1; ++l) for (qpms_m_t m = -l; m <= l; ++m) qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); #endif qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); +#if 0 // Random diagonal T-matrices for(size_t i = 0; i < b1->n; ++i) t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); for(size_t i = 0; i < b2->n; ++i) t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif - const cart3_t pp1 = {1, 2, 0}, pp2 = {1, 2, 3}, pp3 = {0.1, 2, 0}; + const cart3_t pp1 = {0, 0, 1}, pp2 = {0,0, 2}, pp3 = {0,0 , 0}; qpms_tmatrix_t * tmlist[] = {t1, t2}; qpms_particle_tid_t plist[] = {{pp1, 0}, {pp2, 0}, {pp3, 1}}; @@ -58,7 +65,7 @@ int main() protoss.p = plist; protoss.p_count=3; - qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D4h); + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D3h); printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, diff --git a/tests/sss1.c b/tests/sss1.c new file mode 100644 index 0000000..1e94b95 --- /dev/null +++ b/tests/sss1.c @@ -0,0 +1,105 @@ +// c99 -g -I.. ss_syms_packs.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -lblas -llapacke +typedef int qpms_gmi_t;// There is something wrong in the includes, apparently. +#include +#include +#include +#include +#include +#include +#include "staticgroups.h" + +const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h; +const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v; +const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g; +const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v; +const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h; +const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h; + +double uniform_random(double min, double max) { + double random_value = min + (max-min)*(double)rand()/RAND_MAX; + return random_value; +} + +int main() +{ + srand(666); +#if 0 + qpms_vswf_set_spec_t + *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS), + *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS); +#else + // Only electric waves + qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), + *b2 = qpms_vswf_set_spec_init(); + b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS; + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; ++m) + qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; ++m) + qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); +#endif + qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); + qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); + +#if 0 + // Random diagonal T-matrices + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif + + const cart3_t pp1 = {1, 1, 0}; + qpms_tmatrix_t * tmlist[] = {t1, t2}; + qpms_particle_tid_t plist[] = {{pp1, 0}}; + + qpms_scatsys_t protoss; + protoss.tm = tmlist; + protoss.tm_count=2; + protoss.p = plist; + protoss.p_count=1; + + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D2h); + + printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", + (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, + (int)ss->orbit_type_count); + + const double k = 1.7; + + complex double *S_full = qpms_scatsys_build_translation_matrix_full( + NULL, ss, k); + complex double *S_packed[ss->sym->nirreps]; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) + S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL, + S_full, ss, iri); + + complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri) + qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri], + ss, iri, true); + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_recfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + + printf("maxerr: %lg\n", maxerr); + + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]); + free(S_full); + qpms_scatsys_free(ss); + qpms_tmatrix_free(t1); + qpms_tmatrix_free(t2); + qpms_vswf_set_spec_free(b1); + qpms_vswf_set_spec_free(b2); + return 0; +} diff --git a/tests/sss2.c b/tests/sss2.c new file mode 100644 index 0000000..0e9b458 --- /dev/null +++ b/tests/sss2.c @@ -0,0 +1,105 @@ +// c99 -g -I.. ss_syms_packs.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -lblas -llapacke +typedef int qpms_gmi_t;// There is something wrong in the includes, apparently. +#include +#include +#include +#include +#include +#include +#include "staticgroups.h" + +const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h; +const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v; +const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g; +const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v; +const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h; +const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h; + +double uniform_random(double min, double max) { + double random_value = min + (max-min)*(double)rand()/RAND_MAX; + return random_value; +} + +int main() +{ + srand(666); +#if 0 + qpms_vswf_set_spec_t + *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS), + *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS); +#else + // Only electric waves + qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), + *b2 = qpms_vswf_set_spec_init(); + b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS; + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; ++m) + qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; ++m) + qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); +#endif + qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); + qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); + +#if 0 + // Random diagonal T-matrices + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif + + const cart3_t pp1 = {0, 0, 1}; + qpms_tmatrix_t * tmlist[] = {t1, t2}; + qpms_particle_tid_t plist[] = {{pp1, 0}}; + + qpms_scatsys_t protoss; + protoss.tm = tmlist; + protoss.tm_count=2; + protoss.p = plist; + protoss.p_count=1; + + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D2h); + + printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", + (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, + (int)ss->orbit_type_count); + + const double k = 1.7; + + complex double *S_full = qpms_scatsys_build_translation_matrix_full( + NULL, ss, k); + complex double *S_packed[ss->sym->nirreps]; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) + S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL, + S_full, ss, iri); + + complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri) + qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri], + ss, iri, true); + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_recfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + + printf("maxerr: %lg\n", maxerr); + + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]); + free(S_full); + qpms_scatsys_free(ss); + qpms_tmatrix_free(t1); + qpms_tmatrix_free(t2); + qpms_vswf_set_spec_free(b1); + qpms_vswf_set_spec_free(b2); + return 0; +}