oldtests/sss4.c

   1 // c99 -g -DNLINE -DDAGRUP=C4v -DDUMP_PARTICLE_POSITIONS -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss3.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
   2 typedef int qpms_gmi_t;// There is something wrong in the includes, apparently.
   3 #include <qpms/qpms_types.h>
   4 #include <qpms/scatsystem.h>
   5 #include <stdlib.h>
   6 #include <qpms/vswf.h>
   7 #include <qpms/indexing.h>
   8 #include <stdio.h>
   9 #include "staticgroups.h"
  10
  11 const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h;
  12 const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v;
  13 const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g;
  14 const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v;
  15 const qpms_finite_group_t *C2 = &QPMS_FINITE_GROUP_C2;
  16 const qpms_finite_group_t *C4 = &QPMS_FINITE_GROUP_C4;
  17 const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h;
  18 const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h;
  19 const qpms_finite_group_t *x_and_z_flip = &QPMS_FINITE_GROUP_x_and_z_flip;
  20 const qpms_finite_group_t *y_and_z_flip = &QPMS_FINITE_GROUP_y_and_z_flip;
  21
  22 #ifndef DAGRUP
  23 #define DAGRUP D4h
  24 #endif
  25
  26 double uniform_random(double min, double max) {
  27   double random_value = min + (max-min)*(double)rand()/RAND_MAX;
  28   return random_value;
  29 }
  30
  31 int main()
  32 {
  33   srand(666);
  34 #if 1
  35   qpms_vswf_set_spec_t
  36     *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS),
  37     *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS);
  38 #else
  39   // Only electric waves
  40   qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(),
  41                        *b2 = qpms_vswf_set_spec_init();
  42   b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS;
  43   for(qpms_l_t l = 1; l <= 1; ++l)
  44     for (qpms_m_t m = -l; m <= l; m += 1)
  45       qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
  46   for(qpms_l_t l = 1; l <= 1; ++l)
  47     for (qpms_m_t m = -l; m <= l; m += 1)
  48       qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l));
  49 #endif
  50   qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1);
  51   qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2);
  52
  53 #if 0
  54   // Random diagonal T-matrices
  55   for(size_t i = 0; i < b1->n; ++i)
  56     t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
  57   for(size_t i = 0; i < b2->n; ++i)
  58     t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1);
  59 #else
  60   for(size_t i = 0; i < b1->n; ++i)
  61     t1->m[i + i*b1->n] = 1;
  62   for(size_t i = 0; i < b2->n; ++i)
  63     t2->m[i + i*b2->n] = 1;
  64 #endif
  65
  66 #ifdef YLINE
  67   const cart3_t pp1 = {0, 1.1, 0};
  68   const cart3_t pp2 = {0, 1.4, 0};
  69 #elif defined XLINE
  70   const cart3_t pp1 = {1.1, 0, 0};
  71   const cart3_t pp2 = {1.4, 0, 0};
  72 #elif defined ZLINE
  73   const cart3_t pp1 = {0, 0, 1.1};
  74   const cart3_t pp2 = {0, 0, 1.4};
  75 #else
  76   const cart3_t pp1 = {1.1, 1, 0};
  77   const cart3_t pp2 = {0, 1.4, 0};
  78   const cart3_t pp4 = {2, 1.4, 0};
  79 #endif
  80   const cart3_t pp3 = {0, 0, 0};
  81   qpms_tmatrix_t * tmlist[] = {t1, t2};
  82   qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0},  {pp3, 1}, {pp4,1}
  83   };
  84
  85   qpms_scatsys_t protoss;
  86   protoss.tm = tmlist;
  87   protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *);
  88   protoss.p = plist;
  89   protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t);
  90
  91   qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, DAGRUP);
  92
  93   printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n",
  94       (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps,
  95       (int)ss->orbit_type_count);
  96
  97   fputs("Orbit projection matrices:\n", stderr);
  98   for (qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) {
  99     fprintf(stderr, "Orbit type %d:\n", (int)oti);
 100     const qpms_ss_orbit_type_t *ot = &(ss->orbit_types[oti]);
 101     size_t row = 0;
 102     for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
 103       assert(row*ot->size*ot->bspecn == ot->irbase_offsets[iri]);
 104       fprintf(stderr, "---------- IR %d (%s) -------------\n", (int)iri, ss->sym->irreps[iri].name);
 105       for (size_t irbi = 0; irbi < ot->irbase_sizes[iri]; ++irbi) {
 106         for(qpms_ss_orbit_pi_t opi = 0; opi < ot->size; ++opi){
 107           fputs("| ", stderr);
 108           for(size_t i = 0; i < ot->bspecn; ++i) {
 109             const complex double elem = ot->irbases[row * ot->size * ot->bspecn
 110             + opi * ot->bspecn + i];
 111             fprintf(stderr, "%+.3f%+.3fj ", creal(elem), cimag(elem));
 112           }
 113         }
 114         fputs("|\n", stderr);
 115         ++row;
 116       }
 117     }
 118     fputs("------------------------\n\n",stderr);
 119   }
 120
 121   const double k = 1.7;
 122
 123   complex double *S_full = qpms_scatsys_build_translation_matrix_full(
 124       NULL, ss, k);
 125   {
 126     const size_t full_len = ss->fecv_size;
 127     size_t fullvec_offset_dest = 0;
 128     for (qpms_ss_pi_t pdest = 0; pdest < ss->p_count; pdest++) {
 129       size_t fullvec_offset_src = 0;
 130       const size_t bspecn_dest = ss->tm[ss->p[pdest].tmatrix_id]->spec->n;
 131       for (qpms_ss_pi_t psrc = 0; psrc < ss->p_count; psrc++) {
 132         const size_t bspecn_src = ss->tm[ss->p[psrc].tmatrix_id]->spec->n;
 133         fprintf(stderr, "Translation matrix element %d<-%d; (%g %g %g)<-(%g %g %g):\n",
 134             (int)pdest, (int)psrc, ss->p[pdest].pos.x, ss->p[pdest].pos.y, ss->p[pdest].pos.z,
 135             ss->p[psrc].pos.x, ss->p[psrc].pos.y, ss->p[psrc].pos.z);
 136
 137         for(size_t row = 0; row < bspecn_dest; ++row) {
 138           for(size_t col = 0; col < bspecn_src; ++col)
 139             fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]),
 140                 cimag(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]));
 141           fputc('\n', stderr);
 142         }
 143         fullvec_offset_src += bspecn_src;
 144       }
 145       fullvec_offset_dest += bspecn_dest;
 146     }
 147   }
 148   {
 149     fputs("\n\n", stderr);
 150     const size_t full_len = ss->fecv_size;
 151     for (size_t row = 0 ; row < full_len; ++row) {
 152       for (size_t col = 0 ; col < full_len; ++col)
 153         fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * row + col]), cimag(S_full[full_len * row + col]));
 154       fputc('\n', stderr);
 155     }
 156   }
 157
 158   complex double *S_packed[ss->sym->nirreps];
 159   for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
 160     S_packed[iri] = qpms_scatsys_irrep_pack_matrix_stupid(NULL,
 161         S_full, ss, iri);
 162     fprintf(stderr, "--- Packed matrix for irrep %d (%s):\n", (int) iri, ss->sym->irreps[iri].name);
 163     for (size_t row = 0; row < ss->saecv_sizes[iri]; ++row) {
 164       for (size_t col = 0; col < ss->saecv_sizes[iri]; ++col) {
 165         complex double elem = S_packed[iri][row * ss->saecv_sizes[iri] + col];
 166         fprintf(stderr, "%+.3f+%.3fj ", creal(elem), cimag(elem));
 167       }
 168       fputc('\n', stderr);
 169     }
 170   }
 171   {
 172     complex double *S_partrecfull = qpms_scatsys_irrep_unpack_matrix_stupid(NULL,
 173         S_packed[0], ss, 0, false);
 174     for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) {
 175       qpms_scatsys_irrep_unpack_matrix_stupid(S_partrecfull, S_packed[iri],
 176           ss, iri, false);
 177       fprintf(stderr, "\nPartial reconstruction %d (%s):\n", (int)iri, ss->sym->irreps[iri].name);
 178       const size_t full_len = ss->fecv_size;
 179       for (size_t row = 0 ; row < full_len; ++row) {
 180         for (size_t col = 0 ; col < full_len; ++col)
 181           fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_partrecfull[full_len * row + col]), cimag(S_partrecfull[full_len * row + col]));
 182         fputc('\n', stderr);
 183       }
 184     }
 185
 186     double maxerr = 0;
 187     for (size_t i = 0; i < ss->fecv_size; ++i) {
 188       double err = cabs(S_full[i] - S_partrecfull[i]);
 189       maxerr = (err > maxerr) ? err : maxerr;
 190     }
 191     free(S_partrecfull);
 192   }
 193
 194
 195   complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix_stupid(NULL,
 196       S_packed[0], ss, 0, false);
 197   for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri)
 198     qpms_scatsys_irrep_unpack_matrix_stupid(S_recfull, S_packed[iri],
 199         ss, iri, true);
 200   {
 201     fputs("\n\n", stderr);
 202     const size_t full_len = ss->fecv_size;
 203     for (size_t row = 0 ; row < full_len; ++row) {
 204       for (size_t col = 0 ; col < full_len; ++col)
 205         fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_recfull[full_len * row + col]), cimag(S_recfull[full_len * row + col]));
 206       fputc('\n', stderr);
 207     }
 208   }
 209
 210   double maxerr = 0;
 211   for (size_t i = 0; i < ss->fecv_size; ++i) {
 212     double err = cabs(S_full[i] - S_recfull[i]);
 213     maxerr = (err > maxerr) ? err : maxerr;
 214   }
 215
 216
 217   printf("maxerr: %lg\n", maxerr);
 218
 219   fprintf(stderr, "pi\tpos\toti\tosn\tp\n");
 220   for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
 221     cart3_t pos = ss->p[pi].pos;
 222     qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t;
 223     qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn;
 224     qpms_ss_orbit_pi_t p = ss->p_orbitinfo[pi].p;
 225     fprintf(stderr, "%d\t(%.3g,%.3g,%.3g)\t%d\t%d\t%d\n",
 226         (int)pi, pos.x, pos.y, pos.z, (int)oti, (int)osn, (int)p);
 227   }
 228
 229   for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]);
 230   free(S_full);
 231   qpms_scatsys_free(ss);
 232   qpms_tmatrix_free(t1);
 233   qpms_tmatrix_free(t2);
 234   qpms_vswf_set_spec_free(b1);
 235   qpms_vswf_set_spec_free(b2);
 236   return 0;
 237 }