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

WIP trying to fix the half-ellipse contour. 

Former-commit-id: cbba6fffacbb38322a590478899f38c7d7dafe8e
This commit is contained in:
Marek Nečada 2019-09-26 11:36:41 +03:00
parent b3c3eeb2a2
commit 2e7b925be2
3 changed files with 16 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
qpms/beyn.c
View File

@ -125,7 +125,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe,
} }
for(size_t i = 0; i < nline; ++i) { for(size_t i = 0; i < nline; ++i) {
double t = 0.5 * (1 - (double) nline) + i; double t = 0.5 * (1 - (double) nline) + i;
complex double z = centre + faktor * t * l; complex double z = centre + faktor * t * 2 * l / nline;
switch(or) { // ensure correct zero signs; CHECKME!!! switch(or) { // ensure correct zero signs; CHECKME!!!
case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS: case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
if(creal(z) == 0 && signbit(creal(z))) if(creal(z) == 0 && signbit(creal(z)))
@ -146,7 +146,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe,
default: QPMS_WTF; default: QPMS_WTF;
} }
c->z_dz[narc + i][0] = z; c->z_dz[narc + i][0] = z;
c->z_dz[narc + i][1] = faktor * l / nline; c->z_dz[narc + i][1] = faktor * 2 * l / nline;
} }
// We hide the half-axes metadata after the discretisation points. // We hide the half-axes metadata after the discretisation points.
c->z_dz[n][0] = rRe; c->z_dz[n][0] = rRe;

5
tests/CMakeLists.txt
View File

@ -28,6 +28,11 @@ target_link_libraries( tbeyn3a qpms gsl lapacke amos m )
target_include_directories( tbeyn3a PRIVATE .. ) target_include_directories( tbeyn3a PRIVATE .. )
target_compile_definitions( tbeyn3a PRIVATE VARIANTA ) target_compile_definitions( tbeyn3a PRIVATE VARIANTA )
add_executable( tbeyn3a_implus tbeyn3.c )
target_link_libraries( tbeyn3a_implus qpms gsl lapacke amos m )
target_include_directories( tbeyn3a_implus PRIVATE .. )
target_compile_definitions( tbeyn3a_implus PRIVATE VARIANTA IMPLUS )
add_executable( tbeyn3b tbeyn3.c ) add_executable( tbeyn3b tbeyn3.c )
target_link_libraries( tbeyn3b qpms gsl lapacke amos m ) target_link_libraries( tbeyn3b qpms gsl lapacke amos m )
target_include_directories( tbeyn3b PRIVATE .. ) target_include_directories( tbeyn3b PRIVATE .. )

11
tests/tbeyn3.c
View File

@ -1,8 +1,10 @@
#define _GNU_SOURCE
#include <qpms/beyn.h> #include <qpms/beyn.h>
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include <stdlib.h> #include <stdlib.h>
#include <math.h> #include <math.h>
#include <fenv.h>
static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); } static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); }
// Normal distribution via Box-Muller transform // Normal distribution via Box-Muller transform
@ -45,13 +47,14 @@ int M_function(complex double *target, const size_t m, const complex double z, v
} }
int main(int argc, char **argv) { int main(int argc, char **argv) {
complex double z0 = 0; feenableexcept(FE_INVALID | FE_OVERFLOW);
complex double z0 = 0+3e-1*I;
#ifdef RXSMALL #ifdef RXSMALL
double Rx = .1; double Rx = .1;
#else #else
double Rx = .3; // Variant B will fail in this case due to large number of eigenvalues (>30) double Rx = .3; // Variant B will fail in this case due to large number of eigenvalues (>30)
#endif #endif
double Ry = .3; double Ry = .25;
#ifdef VARIANTF #ifdef VARIANTF
int L = 10, N = 150, dim = 10; int L = 10, N = 150, dim = 10;
#else #else
@ -59,7 +62,11 @@ int main(int argc, char **argv) {
#endif #endif
if (argc > 1) N = atoi(argv[1]); if (argc > 1) N = atoi(argv[1]);
if (argc > 2) L = atoi(argv[2]); if (argc > 2) L = atoi(argv[2]);
#ifdef IMPLUS
beyn_contour_t *contour = beyn_contour_halfellipse(z0, Rx, Ry, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
#else
beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
#endif
struct param p; struct param p;
p.T0 = malloc(dim*dim*sizeof(double)); p.T0 = malloc(dim*dim*sizeof(double));
p.T1 = malloc(dim*dim*sizeof(double)); p.T1 = malloc(dim*dim*sizeof(double));