Fix projections - work in progress

Former-commit-id: f2f6d3d762fc7c2d728af0cf31c1229ed8db1b1c

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

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;

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 []

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

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);

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)

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,

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 <qpms/qpms_types.h>
+#include <qpms/scatsystem.h>
+#include <stdlib.h>
+#include <qpms/vswf.h>
+#include <qpms/indexing.h>
+#include <stdio.h>
+#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;
+}

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 <qpms/qpms_types.h>
+#include <qpms/scatsystem.h>
+#include <stdlib.h>
+#include <qpms/vswf.h>
+#include <qpms/indexing.h>
+#include <stdio.h>
+#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;
+}