From cada0093286ce3226c4709e2ae9122c080421247 Mon Sep 17 00:00:00 2001 From: Sander Wahls Date: Tue, 13 Jan 2026 12:33:43 +0100 Subject: [PATCH 01/82] fnft_kdvv_inverse: added skeleton files --- examples/mex_fnft_kdvv_inverse_example_1.m | 35 ++++ include/fnft_kdvv_inverse.h | 82 ++++++++++ matlab/mex_fnft_kdvv_inverse.c | 177 +++++++++++++++++++++ matlab/mex_fnft_kdvv_inverse.m | 49 ++++++ src/fnft_kdvv_inverse.c | 39 +++++ 5 files changed, 382 insertions(+) create mode 100644 examples/mex_fnft_kdvv_inverse_example_1.m create mode 100644 include/fnft_kdvv_inverse.h create mode 100644 matlab/mex_fnft_kdvv_inverse.c create mode 100644 matlab/mex_fnft_kdvv_inverse.m create mode 100644 src/fnft_kdvv_inverse.c diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m new file mode 100644 index 0000000..3ff7505 --- /dev/null +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -0,0 +1,35 @@ +% This file is part of FNFT. +% +% FNFT is free software; you can redistribute it and/or +% modify it under the terms of the version 2 of the GNU General +% Public License as published by the Free Software Foundation. +% +% FNFT is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% +% Contributors: +% Sander Wahls (KIT) 2026. + +clear all; +close all; + +contspec = []; +XI = [0 1]; +bound_states = [1i 2i]; +normconsts = complex([1 -1]); +D = 256; +T = [-10 10]; +% q = MEX_FNFT_KDVVV_INVERSE(contspec, XI, bound_states, ... +% normconsts_or_residuals, D, T); + +q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T) + +t = linspace(T(1), T(2), D); +plot(t, q) +xlabel('t') +ylabel('q(t)') diff --git a/include/fnft_kdvv_inverse.h b/include/fnft_kdvv_inverse.h new file mode 100644 index 0000000..1977d42 --- /dev/null +++ b/include/fnft_kdvv_inverse.h @@ -0,0 +1,82 @@ +/* + * This file is part of FNFT. + * + * FNFT is free software; you can redistribute it and/or + * modify it under the terms of the version 2 of the GNU General + * Public License as published by the Free Software Foundation. + * + * FNFT is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + * Contributors: + * Sander Wahls (KIT) 2026 + */ + +/** + * @file fnft_kdvv_inverse.h + * @brief Fast inverse nonlinear Fourier transform for the vanishing + * Korteweg-de Vries equation. + * @ingroup fnft_inverse + */ + +#ifndef FNFT_KDVV_INVERSE_H +#define FNFT_KDVV_INVERSE_H + +#include "fnft__errwarn.h" + +/** + * @brief Fast inverse nonlinear Fourier transform for the + * Korteweg-de Vries equation with vanishing boundary conditions. + * + * @param[in] M Number of samples of the continuous spectrum. + * @param[in,out] contspec Array of length M, contains samples + * \f$ \hat{q}(\xi_n) \f$, where \f$ \xi_n = XI[0] + n(XI[1]-XI[0])/(M-1) \f$ + * and \f$n=0,1,\dots,M-1\f$, of the to-be-inverted continuous spectrum in + * ascending order (i.e., + * \f$ \hat{q}(\xi_0), \hat{q}(\xi_1), \dots, \hat{q}(\xi_{M-1}) \f$). + * Please note that the routine currently might overwrite this array. + * @param[in] XI Array of length 2, contains the position of the first and the last + * sample of the continuous spectrum. + * @param[in] K Number of discrete spectrum points. + * @param[in] bound_states Complex array of length K. Complex roots of + * \f$ a(\xi) \f$ in the upper half of the complex-plane. + * @param[in] normconsts_or_residues Complex array of length K. Values of + * either the norming constants \f$ b(\xi) \f$ or the residues + * \f$ \frac{b(\xi)}{\partial{a(\xi)}/\partial{\xi}}\f$ at the values bound_states. + * @param[in] D Number of samples of the to be generated signal q. Should be a + * positive power of two. + * @param[out] q Array of length D. Is filled with samples + * \f$ q(t_n) \f$, where \f$ t_n = T[0] + n(T[1]-T[0])/(D-1) \f$ + * and \f$n=0,1,\dots,D-1\f$, of the to-be-generated signal in ascending order + * (i.e., \f$ q(t_0), q(t_1), \dots, q(t_{D-1}) \f$). + * @param[in] T Array of length 2, contains the position in time of the first and + * of the last sample of q. It should be T[0]. +* +* Contributors: +* Sander Wahls (KIT) 2026. +*/ + +#include +#include +#include "mex.h" +#ifndef SKIP_MATRIX_H +#include "matrix.h" +#endif +#include "fnft_kdvv_inverse.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + FNFT_UINT D; + FNFT_COMPLEX * q = NULL; + FNFT_REAL * T; + FNFT_UINT M, K; + FNFT_COMPLEX * contspec = NULL; + FNFT_COMPLEX * bound_states = NULL; + FNFT_COMPLEX * normconsts_or_residues = NULL; + FNFT_REAL * XI; + FNFT_UINT i; + double *re, *im; + char msg[128]; // buffer for error messages + FNFT_INT ret_code; + FNFT_UINT k; + + if (nlhs < 1) + return; + + /* Check types and dimensions of the first seven inputs: contspec, XI, + bound_states, normconsts_or_residues, D, T */ + + if ( nrhs < 6 ) + mexErrMsgTxt("At least seven inputs expected."); + if ( !mxIsEmpty(prhs[0]) && (!mxIsDouble(prhs[0]) || !mxIsComplex(prhs[0]) || mxGetM(prhs[0]) != 1) ) + mexErrMsgTxt("First input contspec should be a complex row vector (double precision) or []. Try passing complex(double(contspec(:)'))."); + if ( mxIsComplex(prhs[1]) || !mxIsDouble(prhs[1]) || mxGetM(prhs[1]) != 1 || mxGetN(prhs[1]) != 2 ) + mexErrMsgTxt("Second input XI should be a real 1x2 vector (double precision)."); + if ( !mxIsEmpty(prhs[2]) && (!mxIsDouble(prhs[2]) || !mxIsComplex(prhs[2]) || mxGetM(prhs[2]) != 1) ) + mexErrMsgTxt("Third input bound_states should be a complex row vector (double precision) or []. Try passing complex(double(bound_states(:)'))."); + if ( !mxIsEmpty(prhs[3]) && (!mxIsDouble(prhs[3]) || !mxIsComplex(prhs[3]) || mxGetM(prhs[3]) != 1) ) + mexErrMsgTxt("Fourth input normconsts_or_residues should be a complex row vector (double precision) or []. Try passing complex(double(normconsts_or_residues(:)'))."); + if ( mxIsComplex(prhs[4]) || !mxIsDouble(prhs[4]) || mxGetNumberOfElements(prhs[4]) != 1 ) + mexErrMsgTxt("Fifth input D should be a real scalar (double precision)."); + if ( mxIsComplex(prhs[5]) || !mxIsDouble(prhs[5]) || mxGetM(prhs[5]) != 1 || mxGetN(prhs[5]) != 2 ) + mexErrMsgTxt("Sixth input T should be a real 1x2 vector (double precision)."); + + M = mxGetNumberOfElements(prhs[0]); + K = mxGetNumberOfElements(prhs[2]); + T = mxGetPr(prhs[5]); + D = (unsigned int)mxGetScalar(prhs[4]); + XI = mxGetPr(prhs[1]); + + /* Check values of first four inputs */ + + if ( K != mxGetNumberOfElements(prhs[3]) ) + mexErrMsgTxt("bound_states and normconsts_or_residues should have the same lengths."); + if ( T[0] >= T[1] ) + mexErrMsgTxt("T(1) >= T(2)."); + if ( XI[0] >= XI[1] ) + mexErrMsgTxt("XI(1) >= XI(2)."); + if ( D<2 ) + mexErrMsgTxt("D < 2."); + + /* Redirect FNFT error messages and warnings to Matlabs command window */ + + fnft_errwarn_setprintf(mexPrintf); + + /* Check remaining inputs, if any */ + + for (k=7; k<(FNFT_UINT)nrhs; k++) { + + /* Check if current input is a string as desired and convert it */ + if ( !mxIsChar(prhs[k]) ) { + snprintf(msg, sizeof msg, "%uth input should be a string.", + (unsigned int)(k+1)); + goto on_error; + } + char *str = mxArrayToString(prhs[k]); + if ( str == NULL ) { + snprintf(msg, sizeof msg, "Out of memory."); + goto on_error; + } + + /* Try to interpret value of string input */ + if ( strcmp(str, "quiet") == 0 ) { + + fnft_errwarn_setprintf(NULL); + + } else { + snprintf(msg, sizeof msg, "%uth input has invalid value.", + (unsigned int)(k+1)); + goto on_error; + } + } + + /* Allocate memory */ + + q = mxMalloc(D * sizeof(FNFT_COMPLEX)); + if (M>0) + contspec = mxMalloc(M * sizeof(FNFT_COMPLEX)); + if (K>0) { + bound_states = mxMalloc(K * sizeof(FNFT_COMPLEX)); + normconsts_or_residues = mxMalloc(K * sizeof(FNFT_COMPLEX)); + } + if ( q == NULL || (M>0 && contspec == NULL) || (K>0 && bound_states == NULL) + || (K>0 && normconsts_or_residues == NULL) ) { + snprintf(msg, sizeof msg, "Out of memory."); + goto on_error; + } + + /* Convert inputs */ + + re = mxGetPr(prhs[0]); + im = mxGetPi(prhs[0]); + for (i=0; i=D, contains the samples +% of the reflection coefficient, the b-scattering +% coefficient or the inverse Fourier transform of the +% b-scattering coefficient on an equidistant grid. +% Pass [] if the continuous spectrum is zero +% (i.e., a multi-soliton is desired) +% XI Real 1x2 vector, contains the location of the first and +% the last sample in contspec +% bound_states Complex row vector, contains the desired bound states. +% Pass [] if the discrete spectrum is empty. +% normconsts_or_residues Complex row vector, same length as bound_states. +% Contains the corresponding norming constants (default) or, +% if the corresponding option is passed, residues. Pass [] +% if the discrete spectrum is empty. +% D Real scalar, number of time domain samples; must be a +% positive power of two +% T Real 1x2 vector, contains the location of the first and +% the last sample in q +% kappa +1.0 or -1.0 +% +% OUTPUTS +% q Complex row vector of length D + +% This file is part of FNFT. +% +% FNFT is free software; you can redistribute it and/or +% modify it under the terms of the version 2 of the GNU General +% Public License as published by the Free Software Foundation. +% +% FNFT is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% +% Contributors: +% Sander Wahls (KIT) 2026. \ No newline at end of file diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c new file mode 100644 index 0000000..23a93f5 --- /dev/null +++ b/src/fnft_kdvv_inverse.c @@ -0,0 +1,39 @@ +/* + * This file is part of FNFT. + * + * FNFT is free software; you can redistribute it and/or + * modify it under the terms of the version 2 of the GNU General + * Public License as published by the Free Software Foundation. + * + * FNFT is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + * Contributors: + * Sander Wahls (KIT) 2026 + */ + +#define FNFT_ENABLE_SHORT_NAMES + +#include "fnft_kdvv_inverse.h" +INT fnft_kdvv_inverse( + const UINT M, + COMPLEX * const contspec, + REAL const * const XI, + UINT const K, + COMPLEX const * const bound_states, + COMPLEX const * const normconsts_or_residues, + const UINT D, + COMPLEX * const q, + REAL const * const T, + void *opts_ptr) +{ + for (UINT n=0; n Date: Mon, 26 Jan 2026 12:49:24 +0100 Subject: [PATCH 02/82] first working version with early stage of crum transformation feature. One test case is hard implemented yet. --- examples/mex_fnft_kdvv_inverse_example_1.m | 11 +- .../mod_mex_fnft_kdvv_inverse_example_1.m | 33 ++ src/fnft_kdvv_inverse.c | 302 +++++++++++++++++- 3 files changed, 337 insertions(+), 9 deletions(-) create mode 100644 examples/mod_mex_fnft_kdvv_inverse_example_1.m diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index 3ff7505..09e5cf0 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -20,13 +20,10 @@ contspec = []; XI = [0 1]; -bound_states = [1i 2i]; -normconsts = complex([1 -1]); -D = 256; -T = [-10 10]; -% q = MEX_FNFT_KDVVV_INVERSE(contspec, XI, bound_states, ... -% normconsts_or_residuals, D, T); - +bound_states = 1i*[1, 2, 3, 4] /2; +normconsts = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); +D = 1001; +T = [-20 20]; q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T) t = linspace(T(1), T(2), D); diff --git a/examples/mod_mex_fnft_kdvv_inverse_example_1.m b/examples/mod_mex_fnft_kdvv_inverse_example_1.m new file mode 100644 index 0000000..99546f2 --- /dev/null +++ b/examples/mod_mex_fnft_kdvv_inverse_example_1.m @@ -0,0 +1,33 @@ +% This file is part of FNFT. +% +% FNFT is free software; you can redistribute it and/or +% modify it under the terms of the version 2 of the GNU General +% Public License as published by the Free Software Foundation. +% +% FNFT is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% +% Contributors: +% Sander Wahls (KIT) 2026. + +close all; + +contspec = []; +XI = [0 1]; +% bound_states = 1i*[1, 2, 3, 4] /2; +bound_states = 1i*[1] /2; +% normconsts = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); +normconsts = complex([1].*[10]); +D = 1001; +T = [-20 20]; +q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T); + +t = linspace(T(1), T(2), D); +plot(t, q) +xlabel('t') +ylabel('q(t)') diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 23a93f5..1646a33 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -20,6 +20,246 @@ #define FNFT_ENABLE_SHORT_NAMES #include "fnft_kdvv_inverse.h" +#include // for printf + + +// INT input_multisoliton_E_pairs( +// UINT const K, +// COMPLEX const * const bound_states, +// COMPLEX const * const normconsts_or_residues, +// const UINT D, +// COMPLEX * const q, +// REAL const * const T) +// { +// INT ret_code = SUCCESS; + +// // Initialize local variables +// COMPLEX tmp; +// COMPLEX * bnd_states = NULL; +// COMPLEX * norm_consts = NULL; + +// bnd_states = malloc(K * sizeof(COMPLEX)); +// CHECK_NOMEM(bnd_states, ret_code, leave_fun); +// norm_consts = malloc(K * sizeof(COMPLEX)); +// CHECK_NOMEM(norm_consts, ret_code, leave_fun); + +// //sorting bnd_states and according norm_const in descending order based on magnitude of imaginary part +// for (UINT i = 0; i < K; ++i){ +// for (UINT j = i + 1; j < K; ++j){ +// if (CIMAG(bnd_states[i]) < CIMAG(bnd_states[j])){ +// tmp = bnd_states[i]; +// bnd_states[i] = bnd_states[j]; +// bnd_states[j] = tmp; +// tmp = norm_consts[i]; +// norm_consts[i] = norm_consts[j]; +// norm_consts[j] = tmp; +// } +// } +// } + + + +// // Apply normconst +// //TODO + +// // Classic Crum transformation + +// UINT N_rest = K; + +// while(N_rest > 0){ + +// // Scale trajectories +// //TODO + +// ret_code = add_two_solitons_E(q, x_grid, ); +// CHECK_RETCODE(ret_code, leave_fun); +// } + +// leave_fun: +// free(bound_states); +// free(norm_consts); + +// return ret_code; +// } + + +INT add_one_soliton_E( + COMPLEX const * const bound_state, + REAL const * const theta_E1, + REAL const * const theta_E2, + REAL const * const x_grid, + const UINT D, + COMPLEX * const q, + REAL const * const T) +{ + INT ret_code = SUCCESS; + + COMPLEX * M_min1_11 = NULL; + M_min1_11 = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(M_min1_11, ret_code, leave_fun); + + //TODO: I weg!! + COMPLEX k1 = bound_state[0]*I; + + // printf("Marker1"); + + // Determine where x>0 + UINT i_pos; + + for (UINT i = 0; i < D; i++){ + if (x_grid[i] > 0){ + i_pos = i; + break; + } + } + + UINT n_pos=D-i_pos; + UINT n_neg=i_pos; + + // Divide x_grid into positive and negative + REAL * x_pos = NULL; + REAL * x_neg = NULL; + + x_pos = malloc(n_pos * sizeof(REAL)); + CHECK_NOMEM(x_pos, ret_code, leave_fun); + + x_neg = malloc(n_neg * sizeof(REAL)); + CHECK_NOMEM(x_neg, ret_code, leave_fun); + + for (UINT i = 0; i < D; i++){ + if (i < i_pos){ + x_neg[i] = x_grid[i]; + } else { + x_pos[i - i_pos] = x_grid[i]; + } + } + + // divide theta_E + COMPLEX * th1_pos = NULL; + COMPLEX * th2_pos = NULL; + COMPLEX * th1_neg = NULL; + COMPLEX * th2_neg = NULL; + + th1_pos = malloc(n_pos * sizeof(COMPLEX)); + CHECK_NOMEM(th1_pos, ret_code, leave_fun); + th2_pos = malloc(n_pos * sizeof(COMPLEX)); + CHECK_NOMEM(th2_pos, ret_code, leave_fun); + th1_neg = malloc(n_neg * sizeof(COMPLEX)); + CHECK_NOMEM(th1_neg, ret_code, leave_fun); + th2_neg = malloc(n_neg * sizeof(COMPLEX)); + CHECK_NOMEM(th2_neg, ret_code, leave_fun); + + for (UINT i=0; i Date: Mon, 26 Jan 2026 14:39:20 +0100 Subject: [PATCH 03/82] Refactoring, Fix scaling problem due to missing sqrt in mod_mex_fnft_kdvv_inverse_example_1.m --- .../mod_mex_fnft_kdvv_inverse_example_1.m | 2 +- src/fnft_kdvv_inverse.c | 55 ++++--------------- 2 files changed, 11 insertions(+), 46 deletions(-) diff --git a/examples/mod_mex_fnft_kdvv_inverse_example_1.m b/examples/mod_mex_fnft_kdvv_inverse_example_1.m index 99546f2..cf42565 100644 --- a/examples/mod_mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mod_mex_fnft_kdvv_inverse_example_1.m @@ -20,7 +20,7 @@ contspec = []; XI = [0 1]; % bound_states = 1i*[1, 2, 3, 4] /2; -bound_states = 1i*[1] /2; +bound_states = 1i*sqrt([1] /2); % normconsts = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); normconsts = complex([1].*[10]); D = 1001; diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 1646a33..2734702 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -89,8 +89,7 @@ INT add_one_soliton_E( REAL const * const theta_E2, REAL const * const x_grid, const UINT D, - COMPLEX * const q, - REAL const * const T) + COMPLEX * const q) { INT ret_code = SUCCESS; @@ -98,7 +97,7 @@ INT add_one_soliton_E( M_min1_11 = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(M_min1_11, ret_code, leave_fun); - //TODO: I weg!! + // Transformation into a real number. COMPLEX k1 = bound_state[0]*I; // printf("Marker1"); @@ -206,28 +205,7 @@ INT add_one_soliton_E( // Update output for (UINT i=0; i Date: Mon, 26 Jan 2026 17:23:55 +0100 Subject: [PATCH 04/82] Refactoring- Preparations for usage with multiple eigenvalues to add. Sorting for bound_states and normconsts implemented. --- src/fnft_kdvv_inverse.c | 140 ++++++++++++++++------------------------ 1 file changed, 57 insertions(+), 83 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 2734702..21a2701 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -23,70 +23,10 @@ #include // for printf -// INT input_multisoliton_E_pairs( -// UINT const K, -// COMPLEX const * const bound_states, -// COMPLEX const * const normconsts_or_residues, -// const UINT D, -// COMPLEX * const q, -// REAL const * const T) -// { -// INT ret_code = SUCCESS; - -// // Initialize local variables -// COMPLEX tmp; -// COMPLEX * bnd_states = NULL; -// COMPLEX * norm_consts = NULL; - -// bnd_states = malloc(K * sizeof(COMPLEX)); -// CHECK_NOMEM(bnd_states, ret_code, leave_fun); -// norm_consts = malloc(K * sizeof(COMPLEX)); -// CHECK_NOMEM(norm_consts, ret_code, leave_fun); - -// //sorting bnd_states and according norm_const in descending order based on magnitude of imaginary part -// for (UINT i = 0; i < K; ++i){ -// for (UINT j = i + 1; j < K; ++j){ -// if (CIMAG(bnd_states[i]) < CIMAG(bnd_states[j])){ -// tmp = bnd_states[i]; -// bnd_states[i] = bnd_states[j]; -// bnd_states[j] = tmp; -// tmp = norm_consts[i]; -// norm_consts[i] = norm_consts[j]; -// norm_consts[j] = tmp; -// } -// } -// } - - - -// // Apply normconst -// //TODO - -// // Classic Crum transformation - -// UINT N_rest = K; - -// while(N_rest > 0){ - -// // Scale trajectories -// //TODO - -// ret_code = add_two_solitons_E(q, x_grid, ); -// CHECK_RETCODE(ret_code, leave_fun); -// } - -// leave_fun: -// free(bound_states); -// free(norm_consts); - -// return ret_code; -// } - - INT add_one_soliton_E( - COMPLEX const * const bound_state, - REAL const * const theta_E1, - REAL const * const theta_E2, + COMPLEX const bound_state, + COMPLEX const * const theta_E1, + COMPLEX const * const theta_E2, REAL const * const x_grid, const UINT D, COMPLEX * const q) @@ -98,7 +38,7 @@ INT add_one_soliton_E( CHECK_NOMEM(M_min1_11, ret_code, leave_fun); // Transformation into a real number. - COMPLEX k1 = bound_state[0]*I; + COMPLEX k1 = bound_state*I; // printf("Marker1"); @@ -265,36 +205,70 @@ INT fnft_kdvv_inverse( x_grid[n]= T[0] + n*eps_t; } - // Initialize theta_E - REAL * theta_E1 = NULL; - theta_E1 = malloc(D * sizeof(REAL)); - CHECK_NOMEM(theta_E1, ret_code, leave_fun); + // Initialize new storage for bound states and normconsts for later sorting + COMPLEX tmp; - REAL * theta_E2 = NULL; - theta_E2 = malloc(D * sizeof(REAL)); - CHECK_NOMEM(theta_E2, ret_code, leave_fun); + COMPLEX * bound_states_sorted = NULL; + bound_states_sorted = malloc(K * sizeof(COMPLEX)); + CHECK_NOMEM(bound_states_sorted, ret_code, leave_fun); - //TODO: only temporary solution! - for (UINT i=0; i Date: Mon, 2 Feb 2026 09:44:16 +0100 Subject: [PATCH 05/82] Refactoring add_one_solition. Beginning implementing add_two solitions. Refactoring of fnft_kdvv_inverse() for this as preparation. --- src/fnft_kdvv_inverse.c | 308 ++++++++++++++++++++++++++++++++++------ 1 file changed, 263 insertions(+), 45 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 21a2701..6f84e0d 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -24,7 +24,7 @@ INT add_one_soliton_E( - COMPLEX const bound_state, + COMPLEX const * const bound_state, COMPLEX const * const theta_E1, COMPLEX const * const theta_E2, REAL const * const x_grid, @@ -38,9 +38,138 @@ INT add_one_soliton_E( CHECK_NOMEM(M_min1_11, ret_code, leave_fun); // Transformation into a real number. - COMPLEX k1 = bound_state*I; + COMPLEX k1 = bound_state[0]*I; - // printf("Marker1"); + fnft_printf_ptr_t printf_ptr = fnft_errwarn_getprintf(); + printf_ptr("Marker! %d \n", ret_code); + + // Determine where x>0 + UINT i_pos; + + for (UINT i = 0; i < D; i++){ + if (x_grid[i] > 0){ + i_pos = i; + break; + } + } + + COMPLEX * w_inv = NULL; + w_inv = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(w_inv, ret_code, leave_fun); + + COMPLEX * prefactor = NULL; + prefactor = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(prefactor, ret_code, leave_fun); + + + // Calculation for positive x + + for (UINT i = i_pos; i0 UINT i_pos; @@ -74,30 +203,77 @@ INT add_one_soliton_E( } // divide theta_E - COMPLEX * th1_pos = NULL; - COMPLEX * th2_pos = NULL; - COMPLEX * th1_neg = NULL; - COMPLEX * th2_neg = NULL; - - th1_pos = malloc(n_pos * sizeof(COMPLEX)); - CHECK_NOMEM(th1_pos, ret_code, leave_fun); - th2_pos = malloc(n_pos * sizeof(COMPLEX)); - CHECK_NOMEM(th2_pos, ret_code, leave_fun); - th1_neg = malloc(n_neg * sizeof(COMPLEX)); - CHECK_NOMEM(th1_neg, ret_code, leave_fun); - th2_neg = malloc(n_neg * sizeof(COMPLEX)); - CHECK_NOMEM(th2_neg, ret_code, leave_fun); + + // theta_E for 1.eigenvalue + COMPLEX * th11_pos = NULL; + COMPLEX * th12_pos = NULL; + COMPLEX * th11_neg = NULL; + COMPLEX * th12_neg = NULL; + + th11_pos = malloc(n_pos * sizeof(COMPLEX)); + CHECK_NOMEM(th11_pos, ret_code, leave_fun); + th12_pos = malloc(n_pos * sizeof(COMPLEX)); + CHECK_NOMEM(th12_pos, ret_code, leave_fun); + th11_neg = malloc(n_neg * sizeof(COMPLEX)); + CHECK_NOMEM(th11_neg, ret_code, leave_fun); + th12_neg = malloc(n_neg * sizeof(COMPLEX)); + CHECK_NOMEM(th12_neg, ret_code, leave_fun); + + // theta_E for 2.eigenvalue + COMPLEX * th21_pos = NULL; + COMPLEX * th22_pos = NULL; + COMPLEX * th21_neg = NULL; + COMPLEX * th22_neg = NULL; + + th21_pos = malloc(n_pos * sizeof(COMPLEX)); + CHECK_NOMEM(th21_pos, ret_code, leave_fun); + th22_pos = malloc(n_pos * sizeof(COMPLEX)); + CHECK_NOMEM(th22_pos, ret_code, leave_fun); + th22_neg = malloc(n_neg * sizeof(COMPLEX)); + CHECK_NOMEM(th22_neg, ret_code, leave_fun); + th21_neg = malloc(n_neg * sizeof(COMPLEX)); + CHECK_NOMEM(th21_neg, ret_code, leave_fun); + + for (UINT i=0; i 0){ + + // If the number of eigenvalues left is odd, then only one solition + // should be added + if (N_rest%2==1){ - // Assign values to theta_E dependend on next eigenvalue to add - for (UINT j=0; j Date: Mon, 2 Feb 2026 13:10:48 +0100 Subject: [PATCH 06/82] Refactoring of add_one_solition() and fnft_kdvv_inverse() finished so far. Ready for implementing add_two_solitons(). Test with adding two solitions only with add_one_solition() works, the result is quite similar to matlab where these two solitons also are added only with add_one_soliton. However there are some differs. Maybe due to numerical characteristics? --- src/fnft_kdvv_inverse.c | 198 ++++++++++++++++++++++------------------ 1 file changed, 110 insertions(+), 88 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 6f84e0d..52e5478 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -20,28 +20,39 @@ #define FNFT_ENABLE_SHORT_NAMES #include "fnft_kdvv_inverse.h" -#include // for printf - -INT add_one_soliton_E( - COMPLEX const * const bound_state, - COMPLEX const * const theta_E1, - COMPLEX const * const theta_E2, +INT add_one_soliton( + COMPLEX const * const bound_state_to_add, + COMPLEX const * const bound_states_left, + UINT const N_bound_states_left, + COMPLEX * const theta_E1, + COMPLEX * const theta_E2, REAL const * const x_grid, const UINT D, COMPLEX * const q) { + // Initialize variables and arrays INT ret_code = SUCCESS; COMPLEX * M_min1_11 = NULL; M_min1_11 = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(M_min1_11, ret_code, leave_fun); - // Transformation into a real number. - COMPLEX k1 = bound_state[0]*I; + COMPLEX * w_inv = NULL; + w_inv = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(w_inv, ret_code, leave_fun); + + COMPLEX * prefactor = NULL; + prefactor = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(prefactor, ret_code, leave_fun); - fnft_printf_ptr_t printf_ptr = fnft_errwarn_getprintf(); - printf_ptr("Marker! %d \n", ret_code); + // Defining theta vectors out of theta_E vectors, belonging to the bound states to add + // These vectors are used in the Crum-Transformation step + COMPLEX * th1 = &theta_E1[0]; + COMPLEX * th2 = &theta_E2[0]; + + // Transformation into a real number + COMPLEX k1 = bound_state_to_add[0]; // Determine where x>0 UINT i_pos; @@ -53,29 +64,18 @@ INT add_one_soliton_E( } } - COMPLEX * w_inv = NULL; - w_inv = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(w_inv, ret_code, leave_fun); - - COMPLEX * prefactor = NULL; - prefactor = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(prefactor, ret_code, leave_fun); - - + // -- Crum Transform -- // Calculation for positive x - for (UINT i = i_pos; i 0){ - // If the number of eigenvalues left is odd, then only one solition - // should be added - if (N_rest%2==1){ - - for (UINT j=0; j Date: Mon, 2 Feb 2026 16:47:59 +0100 Subject: [PATCH 07/82] Implementing add_two_solitons. Works, but differences to matlab, similar to problem with add_one_solition. --- src/fnft_kdvv_inverse.c | 349 +++++++++++++++++++++++----------------- 1 file changed, 203 insertions(+), 146 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 52e5478..076a532 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -46,7 +46,7 @@ INT add_one_soliton( prefactor = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(prefactor, ret_code, leave_fun); - // Defining theta vectors out of theta_E vectors, belonging to the bound states to add + // Defining theta vectors out of theta_E vectors, belonging to the bound state to add // These vectors are used in the Crum-Transformation step COMPLEX * th1 = &theta_E1[0]; COMPLEX * th2 = &theta_E2[0]; @@ -91,10 +91,8 @@ INT add_one_soliton( UINT const is_Jost_to_update = N_jZ != 0; // theta vectors according to the bound states left, which corresponds to the Jost solutions - COMPLEX * th1_jZ = &theta_E1[D*1]; - COMPLEX * th2_jZ = &theta_E2[D*1]; - - // TODO: sind jZ komplex oder reell gemacht worden? + COMPLEX * th1_jZ = &theta_E1[D]; + COMPLEX * th2_jZ = &theta_E2[D]; // Create C_E Matrix (for Updating Jost Solution) UINT n_C_E = 2*2*D*N_jZ; @@ -125,24 +123,24 @@ INT add_one_soliton( // Calculation for positive x for (UINT i=i_pos; i0 UINT i_pos; @@ -197,159 +221,181 @@ INT add_two_solitons( } } - UINT n_pos=D-i_pos; - UINT n_neg=i_pos; + // Defining theta vectors out of theta_E vectors, belonging to the bound states to add + // These vectors are used in the Crum-Transformation step + // 1. eigenvalue + COMPLEX * th11 = &theta_E1[0]; + COMPLEX * th12 = &theta_E2[0]; + + //2. eigenvalue + COMPLEX * th21 = &theta_E1[D]; + COMPLEX * th22 = &theta_E2[D]; - // Divide x_grid into positive and negative - REAL * x_pos = NULL; - REAL * x_neg = NULL; + // Intermediate variables + COMPLEX * t1 = NULL; + COMPLEX * t2 = NULL; + + t1 = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(t1, ret_code, leave_fun); + t2 = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(t2, ret_code, leave_fun); + + // -- Crum-Transformation -- + // for positive x + for (UINT i=i_pos; i 0){ // If the number of eigenvalues left is odd, then only one solition should be added - // if (N_rest%2==1){ - if (N_rest){ + if (N_rest%2==1){ // Crum-Transformation step ret_code = add_one_soliton(&bound_states_sorted[step_idx], &bound_states_sorted[step_idx+1], @@ -479,7 +524,19 @@ INT fnft_kdvv_inverse( // If the number of eigenvalues left is even, then two solitions // should be added else { - + // Crum-Transformation step + ret_code = add_two_solitons(&bound_states_sorted[step_idx], + &bound_states_sorted[step_idx+2], + N_rest-2, + &theta_E1[step_idx*D], + &theta_E2[step_idx*D], + x_grid, + D, + q); + CHECK_RETCODE(ret_code, leave_fun); + + step_idx += 2; + N_rest -= 2; } // Filter out errors that result in negative values From 203819dcb72606ce29b4d3c8709e13adfb50eaf4 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 2 Feb 2026 18:11:42 +0100 Subject: [PATCH 08/82] small refactorings --- src/fnft_kdvv_inverse.c | 72 ++++++++++++++++++----------------------- 1 file changed, 31 insertions(+), 41 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 076a532..ea3345e 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -67,16 +67,16 @@ INT add_one_soliton( // -- Crum Transform -- // Calculation for positive x for (UINT i = i_pos; i 0){ + // Select indices for this step and for the remaining eigenvalues to add // If the number of eigenvalues left is odd, then only one solition should be added - if (N_rest%2==1){ - // Crum-Transformation step + if (N_rest%2==1){ N_step = 1; } + else { N_step = 2; } + + // Scale trajectories (magnitudes of th1 and th2 symmetric around 0) + for (UINT i=0; i Date: Mon, 2 Feb 2026 18:28:27 +0100 Subject: [PATCH 09/82] Remove mod_mex_fnft_kdvv_inverse_example_1.m from project --- .../mod_mex_fnft_kdvv_inverse_example_1.m | 33 ------------------- 1 file changed, 33 deletions(-) delete mode 100644 examples/mod_mex_fnft_kdvv_inverse_example_1.m diff --git a/examples/mod_mex_fnft_kdvv_inverse_example_1.m b/examples/mod_mex_fnft_kdvv_inverse_example_1.m deleted file mode 100644 index cf42565..0000000 --- a/examples/mod_mex_fnft_kdvv_inverse_example_1.m +++ /dev/null @@ -1,33 +0,0 @@ -% This file is part of FNFT. -% -% FNFT is free software; you can redistribute it and/or -% modify it under the terms of the version 2 of the GNU General -% Public License as published by the Free Software Foundation. -% -% FNFT is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program. If not, see . -% -% Contributors: -% Sander Wahls (KIT) 2026. - -close all; - -contspec = []; -XI = [0 1]; -% bound_states = 1i*[1, 2, 3, 4] /2; -bound_states = 1i*sqrt([1] /2); -% normconsts = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); -normconsts = complex([1].*[10]); -D = 1001; -T = [-20 20]; -q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T); - -t = linspace(T(1), T(2), D); -plot(t, q) -xlabel('t') -ylabel('q(t)') From 8a20f8a19ba92960b74cf78d922e6516d4be2c32 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 9 Feb 2026 10:27:43 +0100 Subject: [PATCH 10/82] Fixing scaling problem in fnft_kdvv_inverse --- src/fnft_kdvv_inverse.c | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index ea3345e..5c240d8 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -27,7 +27,7 @@ INT add_one_soliton( UINT const N_bound_states_left, COMPLEX * const theta_E1, COMPLEX * const theta_E2, - REAL const * const x_grid, + COMPLEX const * const x_grid, const UINT D, COMPLEX * const q) { @@ -58,7 +58,7 @@ INT add_one_soliton( UINT i_pos; for (UINT i = 0; i < D; i++){ - if (x_grid[i] > 0){ + if (CREAL(x_grid[i]) > 0){ i_pos = i; break; } @@ -79,12 +79,12 @@ INT add_one_soliton( M_min1_11[i] = prefactor[i] * CEXP(-2*k1*x_grid[i]); } - + // Update output for (UINT i=0; i 0){ + if (CREAL(x_grid[i]) > 0){ i_pos = i; break; } @@ -436,11 +435,11 @@ INT fnft_kdvv_inverse( } // Initialize spatial grid with D points between T[0] and T[1] - REAL * x_grid = NULL; - x_grid = malloc(D * sizeof(REAL)); + COMPLEX * x_grid = NULL; + x_grid = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(x_grid, ret_code, leave_fun); - const REAL eps_t = (T[1] - T[0])/(D - 1); + const COMPLEX eps_t = (T[1] - T[0])/(D - 1); for (UINT n=0; n Date: Mon, 9 Feb 2026 12:15:56 +0100 Subject: [PATCH 11/82] Fixing error at initializing C_E in add_one_solition --- src/fnft_kdvv_inverse.c | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 5c240d8..5776c12 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -48,8 +48,8 @@ INT add_one_soliton( // Defining theta vectors out of theta_E vectors, belonging to the bound state to add // These vectors are used in the Crum-Transformation step - COMPLEX * th1 = &theta_E1[0]; - COMPLEX * th2 = &theta_E2[0]; + COMPLEX * const th1 = &theta_E1[0]; + COMPLEX * const th2 = &theta_E2[0]; // Transformation into a real number COMPLEX k1 = bound_state_to_add[0]; @@ -111,8 +111,8 @@ INT add_one_soliton( for (UINT i=0; i Date: Mon, 9 Feb 2026 15:40:58 +0100 Subject: [PATCH 12/82] fixing errors in add_two_solitions --- src/fnft_kdvv_inverse.c | 68 +++++++++++++++++++++++++---------------- 1 file changed, 42 insertions(+), 26 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 5776c12..6b693ca 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -201,8 +201,8 @@ INT add_two_solitons( CHECK_NOMEM(pm0_jZ, ret_code, leave_fun); for (UINT i=0; i Date: Tue, 10 Feb 2026 17:41:51 +0100 Subject: [PATCH 13/82] Refactoring, defining auxiliary functions --- src/fnft_kdvv_inverse.c | 497 ++++++++++++++++++++++++++-------------- 1 file changed, 328 insertions(+), 169 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 6b693ca..74e4202 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -21,50 +21,36 @@ #include "fnft_kdvv_inverse.h" -INT add_one_soliton( - COMPLEX const * const bound_state_to_add, - COMPLEX const * const bound_states_left, - UINT const N_bound_states_left, - COMPLEX * const theta_E1, - COMPLEX * const theta_E2, - COMPLEX const * const x_grid, - const UINT D, - COMPLEX * const q) +static INT find_first_positive_value( + COMPLEX const * const array, + const UINT length, + UINT * const i_pos_ptr) { - // Initialize variables and arrays INT ret_code = SUCCESS; - - COMPLEX * M_min1_11 = NULL; - M_min1_11 = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(M_min1_11, ret_code, leave_fun); - - COMPLEX * w_inv = NULL; - w_inv = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(w_inv, ret_code, leave_fun); - - COMPLEX * prefactor = NULL; - prefactor = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(prefactor, ret_code, leave_fun); - - // Defining theta vectors out of theta_E vectors, belonging to the bound state to add - // These vectors are used in the Crum-Transformation step - COMPLEX * const th1 = &theta_E1[0]; - COMPLEX * const th2 = &theta_E2[0]; - - // Transformation into a real number - COMPLEX k1 = bound_state_to_add[0]; - - // Determine where x>0 - UINT i_pos; - for (UINT i = 0; i < D; i++){ - if (CREAL(x_grid[i]) > 0){ - i_pos = i; + for (UINT i = 0; i < length; i++){ + if (CREAL(array[i]) > 0){ + *i_pos_ptr = i; break; } } - // -- Crum Transform -- + return ret_code; +} + +static INT one_solition_crum_transformation( + COMPLEX const k1, + COMPLEX const * const th1, + COMPLEX const * const th2, + const UINT D, + const UINT i_pos, + COMPLEX const * const x_grid, + COMPLEX * const w_inv, + COMPLEX * const prefactor, + COMPLEX * const M_min1_11) +{ + INT ret_code = SUCCESS; + // Calculation for positive x for (UINT i = i_pos; i0 - UINT i_pos; - - for (UINT i = 0; i < D; i++){ - if (CREAL(x_grid[i]) > 0){ - i_pos = i; - break; - } - } + // Check, if Jost has to be updated. Only when there are bound states left, that are not added yet + UINT const is_Jost_to_update = N_jZ != 0; // Defining theta vectors out of theta_E vectors, belonging to the bound states to add // These vectors are used in the Crum-Transformation step // thetas according to 1. eigenvalue - COMPLEX * th11 = &theta_E1[0]; - COMPLEX * th12 = &theta_E2[0]; + COMPLEX const * const th11 = &theta_E1[0]; + COMPLEX const * const th12 = &theta_E2[0]; // thetas according to 2. eigenvalue - COMPLEX * th21 = &theta_E1[D]; - COMPLEX * th22 = &theta_E2[D]; - - // Intermediate variables - COMPLEX * t1 = NULL; - COMPLEX * t2 = NULL; - - t1 = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(t1, ret_code, leave_fun); - t2 = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(t2, ret_code, leave_fun); - - // -- Crum-Transformation -- - // for positive x - for (UINT i=i_pos; i0 + UINT i_pos; + ret_code = find_first_positive_value(x_grid, D, &i_pos); + CHECK_RETCODE(ret_code, leave_fun); + + // -- Crum Transform -- + COMPLEX bound_state_added_here = bound_states[0]; + ret_code = one_solition_crum_transformation( bound_state_added_here, + &theta_E1[0], // thetas belonging to the bound state to add + &theta_E2[0], // thetas belonging to the bound state to add + D, + i_pos, + x_grid, + w_inv, + prefactor, + M_min1_11); + CHECK_RETCODE(ret_code, leave_fun); + + // Update output + for (UINT i=0; i0 + UINT i_pos; + ret_code = find_first_positive_value(x_grid, D, &i_pos); + CHECK_RETCODE(ret_code, leave_fun); + + // -- Crum-Transformation -- + ret_code = two_solitions_crum_transformation( k1, + k2, + N_pm0_jZ, + pm0_jZ, + theta_E1, + theta_E2, + t1, + t2, + D, + i_pos, + x_grid, + w_inv, + prefactor, + dq, + s ); + CHECK_RETCODE(ret_code, leave_fun); + + + // Update output + for (UINT i=0; i Date: Mon, 16 Feb 2026 14:25:20 +0100 Subject: [PATCH 14/82] Adding first test for inverse kdvv fnft --- .gitignore | 1 + include/fnft_kdvv_inverse.h | 2 + .../private/fnft__kdvv_inverse_testcases.h | 27 ++++ src/fnft_kdvv_inverse.c | 9 -- .../fnft_kdvv_inverse_test_1.c | 123 ++++++++++++++++++ 5 files changed, 153 insertions(+), 9 deletions(-) create mode 100644 include/private/fnft__kdvv_inverse_testcases.h create mode 100644 test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c diff --git a/.gitignore b/.gitignore index e6907a0..d6c8d60 100644 --- a/.gitignore +++ b/.gitignore @@ -35,6 +35,7 @@ !/test/fnft__poly/ !/test/fnft_version_test/ !/test/fnft_version_test.c/ +!/test/fnft_kdvv_inverse/ # Ignore MacOS hidden files in subdirectories /**/*.DS_Store diff --git a/include/fnft_kdvv_inverse.h b/include/fnft_kdvv_inverse.h index 1977d42..d558ee2 100644 --- a/include/fnft_kdvv_inverse.h +++ b/include/fnft_kdvv_inverse.h @@ -28,6 +28,8 @@ #define FNFT_KDVV_INVERSE_H #include "fnft__errwarn.h" +#include "fnft__misc.h" +#include "fnft_numtypes.h" /** * @brief Fast inverse nonlinear Fourier transform for the diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h new file mode 100644 index 0000000..d32700c --- /dev/null +++ b/include/private/fnft__kdvv_inverse_testcases.h @@ -0,0 +1,27 @@ +/* +* This file is part of FNFT. +* +* FNFT is free software; you can redistribute it and/or +* modify it under the terms of the version 2 of the GNU General +* Public License as published by the Free Software Foundation. +* +* FNFT is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Contributors: +* Sander Wahls (TU Delft) 2018. +*/ + +#ifndef FNFT__KDVV_INVERSE_TESTCASES_H +#define FNFT__KDVV_INVERSE_TESTCASES_H + +#include "fnft_kdvv.h" +#include "fnft_kdvv_inverse.h" +#include "fnft_numtypes.h" + +#endif diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 74e4202..829cf8d 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -584,8 +584,6 @@ INT add_two_solitons( - - INT fnft_kdvv_inverse( const UINT M, COMPLEX * const contspec, @@ -605,11 +603,6 @@ INT fnft_kdvv_inverse( */ - /* Tests - number of normconsts = number bound_states - */ - - INT ret_code = SUCCESS; // Initialize q @@ -687,7 +680,6 @@ INT fnft_kdvv_inverse( // Add solitions while (N_rest > 0){ - // Select indices for this step and for the remaining eigenvalues to add // If the number of eigenvalues left is odd, then only one solition should be added if (N_rest%2==1){ N_step = 1; } @@ -732,7 +724,6 @@ INT fnft_kdvv_inverse( for (UINT i=0; i. +* +* Contributors: +*/ + +#define FNFT_ENABLE_SHORT_NAMES + +#include + +#include "fnft__kdvv_inverse_testcases.h" +#include "fnft__errwarn.h" +#include "fnft__misc.h" + +INT main() +{ + INT ret_code = SUCCESS; + + // General parameters + UINT M = 0; + REAL XI[2] = {-2.0, 2.0}; + UINT D = 1001; + REAL T[2] = {-10.0, 10.0}; + + COMPLEX * q = NULL; + q = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(q, ret_code, leave_fun); + + // Initialize variables for a specific test case without continuous spectrum + COMPLEX * contspec_i = NULL; + contspec_i = malloc(M * sizeof(COMPLEX)); + const UINT K_i = 5; + COMPLEX bound_states_i[5] = {I*SQRT(1.0/2.0), I*SQRT(2.0/2.0), I*SQRT(3.0/2.0), I*SQRT(4.0/2.0), I*SQRT(5.0/2.0)}; + COMPLEX normconsts_i[5] = {1*10, -1*0.1, 1*1, -1*0.00001, 1*10000000}; + + // const UINT K_i = 1; + // COMPLEX bound_states_i[1] = {SQRT(1.0/2)*I}; + // COMPLEX normconsts_i[1] = {1*10}; + + ret_code = fnft_kdvv_inverse(M, contspec_i, XI, K_i, bound_states_i, normconsts_i, D, q, T, NULL); + CHECK_RETCODE(ret_code, leave_fun); + + // printf("function inverse: %d \n", ret_code); + + // // variables for the result of the forward fnft_kdvv + + COMPLEX * contspec_r = NULL; + contspec_r = malloc(M * sizeof(COMPLEX)); + CHECK_NOMEM(contspec_r, ret_code, leave_fun); + + UINT K_r = 5; + + COMPLEX * bound_states_r = NULL; + bound_states_r = malloc(K_r * sizeof(COMPLEX)); + CHECK_NOMEM(bound_states_r, ret_code, leave_fun); + + COMPLEX * normconsts_r = NULL; + normconsts_r = malloc(K_r * sizeof(COMPLEX)); + CHECK_NOMEM(normconsts_r, ret_code, leave_fun); + + fnft_kdvv_opts_t opts = fnft_kdvv_default_opts(); + + ret_code = fnft_kdvv(D, q, T, M, contspec_r, XI, &K_r, bound_states_r, normconsts_r, &opts); + CHECK_RETCODE(ret_code, leave_fun); + + + if (ret_code != FNFT_SUCCESS) { + printf("An error occured!\n"); + return EXIT_FAILURE; + } + + + /** Step 4: Print the results **/ + + printf("Number of samples:\n D = %u\n", (unsigned int)D); + + FNFT_REAL eps_xi = (XI[1] - XI[0]) / (M - 1); + printf("Continuous spectrum:\n"); + for (FNFT_UINT i=0; i Date: Mon, 16 Feb 2026 15:25:05 +0100 Subject: [PATCH 15/82] Adding conditions in fnft_kdvv_inverse_test_1 --- .../fnft_kdvv_inverse_test_1.c | 41 ++++++++++++++----- 1 file changed, 30 insertions(+), 11 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index 76dddb4..538d6aa 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -51,16 +51,13 @@ INT main() ret_code = fnft_kdvv_inverse(M, contspec_i, XI, K_i, bound_states_i, normconsts_i, D, q, T, NULL); CHECK_RETCODE(ret_code, leave_fun); - - // printf("function inverse: %d \n", ret_code); - // // variables for the result of the forward fnft_kdvv - + // Prepare forward fnft_kdvv COMPLEX * contspec_r = NULL; contspec_r = malloc(M * sizeof(COMPLEX)); CHECK_NOMEM(contspec_r, ret_code, leave_fun); - UINT K_r = 5; + UINT K_r = K_i; COMPLEX * bound_states_r = NULL; bound_states_r = malloc(K_r * sizeof(COMPLEX)); @@ -75,15 +72,37 @@ INT main() ret_code = fnft_kdvv(D, q, T, M, contspec_r, XI, &K_r, bound_states_r, normconsts_r, &opts); CHECK_RETCODE(ret_code, leave_fun); - - if (ret_code != FNFT_SUCCESS) { - printf("An error occured!\n"); - return EXIT_FAILURE; + // Check result + // bound_states_r is in ascending order -> start from the end + for (UINT i=0; i 0) && is_bsr_pure_imaginary; + UINT is_bsr_in_tolerance = CABS(bsr - bsi)/CABS(bsr) < 1e-3; + + UINT is_ncr_real = CABS(CIMAG(ncr)) < 1e-8; + UINT is_ncr_in_tolerance = CABS(ncr - nci)/CABS(ncr) < 1e-2; + + if (is_bsr_pure_imaginary && + is_bsr_positive_imaginary && + is_bsr_in_tolerance && + is_ncr_real && + is_ncr_in_tolerance) + { + ret_code = SUCCESS; + } + else { + ret_code = FNFT_EC_TEST_FAILED; + break; + } } - /** Step 4: Print the results **/ - + // Print the results printf("Number of samples:\n D = %u\n", (unsigned int)D); FNFT_REAL eps_xi = (XI[1] - XI[0]) / (M - 1); From b1e9e2239c46bd6ca412000c6c43305dc435233c Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 17 Feb 2026 16:45:18 +0100 Subject: [PATCH 16/82] Change in fnft__misc.c from printf to fnft_errwarn_getprinf() for output in matlab --- src/private/fnft__misc.c | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/private/fnft__misc.c b/src/private/fnft__misc.c index 1dc90b6..eb53af8 100644 --- a/src/private/fnft__misc.c +++ b/src/private/fnft__misc.c @@ -17,6 +17,7 @@ * Sander Wahls (TU Delft) 2017. * Peter J Prins (TU Delft) 2018-2020. * Shrinivas Chimmalgi (TU Delft) 2019-2020. +* Fabian Fischer (Hiwi KIT) 2026 */ #define FNFT_ENABLE_SHORT_NAMES @@ -28,27 +29,31 @@ void misc_print_buf(const UINT len, COMPLEX const * const buf, char const * const varname) { + fnft_printf_ptr_t printf_ptr = fnft_errwarn_getprintf(); + UINT i; - printf("%s = [", varname); + printf_ptr("%s = [", varname); for (i = 0; i < len; i++) { - printf("%1.12e+%1.12ej", CREAL(buf[i]), CIMAG(buf[i])); + printf_ptr("%1.12e+%1.12ej", CREAL(buf[i]), CIMAG(buf[i])); if (i != len-1) - printf(", "); + printf_ptr(", "); } - printf("];\n"); + printf_ptr("];\n"); } void misc_print_buf_real(const UINT len, REAL const * const buf, char const * const varname) { + fnft_printf_ptr_t printf_ptr = fnft_errwarn_getprintf(); + UINT i; - printf("%s = [", varname); + printf_ptr("%s = [", varname); for (i = 0; i < len; i++) { - printf("%1.12e", buf[i]); + printf_ptr("%1.12e", buf[i]); if (i != len-1) - printf(", "); + printf_ptr(", "); } - printf("];\n"); + printf_ptr("];\n"); } REAL misc_rel_err(const UINT len, COMPLEX const * const vec_numer, From 0fcc53c37af2149e36d16c9bfcced10ef7d7d11f Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 17 Feb 2026 17:10:22 +0100 Subject: [PATCH 17/82] Rename in fnft_kdvv_inverse.c x_grid to t_grid, t1 to temp_t1, t2 to temp_t2 --- src/fnft_kdvv_inverse.c | 162 ++++++++++++++++++++-------------------- 1 file changed, 83 insertions(+), 79 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 829cf8d..e446197 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -17,6 +17,8 @@ * Sander Wahls (KIT) 2026 */ +// TODO comment matlab + #define FNFT_ENABLE_SHORT_NAMES #include "fnft_kdvv_inverse.h" @@ -35,6 +37,8 @@ static INT find_first_positive_value( } } + // TODO wenn kein positive value?! + return ret_code; } @@ -44,7 +48,7 @@ static INT one_solition_crum_transformation( COMPLEX const * const th2, const UINT D, const UINT i_pos, - COMPLEX const * const x_grid, + COMPLEX const * const t_grid, COMPLEX * const w_inv, COMPLEX * const prefactor, COMPLEX * const M_min1_11) @@ -53,16 +57,16 @@ static INT one_solition_crum_transformation( // Calculation for positive x for (UINT i = i_pos; i0 UINT i_pos; - ret_code = find_first_positive_value(x_grid, D, &i_pos); + ret_code = find_first_positive_value(t_grid, D, &i_pos); CHECK_RETCODE(ret_code, leave_fun); // -- Crum Transform -- @@ -421,7 +425,7 @@ static INT add_one_soliton( &theta_E2[0], // thetas belonging to the bound state to add D, i_pos, - x_grid, + t_grid, w_inv, prefactor, M_min1_11); @@ -443,7 +447,7 @@ static INT add_one_soliton( theta_E2, D, i_pos, - x_grid, + t_grid, w_inv, prefactor, M_min1_11); @@ -464,7 +468,7 @@ INT add_two_solitons( UINT const N_bound_states, COMPLEX * const theta_E1, COMPLEX * const theta_E2, - COMPLEX const * const x_grid, + COMPLEX const * const t_grid, const UINT D, COMPLEX * const q) { @@ -509,13 +513,13 @@ INT add_two_solitons( CHECK_NOMEM(s, ret_code, leave_fun); // Intermediate variables - COMPLEX * t1 = NULL; - COMPLEX * t2 = NULL; + COMPLEX * temp_t1 = NULL; + COMPLEX * temp_t2 = NULL; - t1 = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(t1, ret_code, leave_fun); - t2 = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(t2, ret_code, leave_fun); + temp_t1 = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(temp_t1, ret_code, leave_fun); + temp_t2 = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(temp_t2, ret_code, leave_fun); // Change Sign of the bound states to add COMPLEX k1 = -bound_states[0]; @@ -523,7 +527,7 @@ INT add_two_solitons( // Determine where x>0 UINT i_pos; - ret_code = find_first_positive_value(x_grid, D, &i_pos); + ret_code = find_first_positive_value(t_grid, D, &i_pos); CHECK_RETCODE(ret_code, leave_fun); // -- Crum-Transformation -- @@ -533,11 +537,11 @@ INT add_two_solitons( pm0_jZ, theta_E1, theta_E2, - t1, - t2, + temp_t1, + temp_t2, D, i_pos, - x_grid, + t_grid, w_inv, prefactor, dq, @@ -559,11 +563,11 @@ INT add_two_solitons( pm0_jZ, theta_E1, theta_E2, - t1, - t2, + temp_t1, + temp_t2, D, i_pos, - x_grid, + t_grid, w_inv, prefactor, s ); @@ -576,8 +580,8 @@ INT add_two_solitons( free(dq); free(pm0_jZ); free(s); - free(t1); - free(t2); + free(temp_t1); + free(temp_t2); return ret_code; } @@ -611,14 +615,14 @@ INT fnft_kdvv_inverse( } // Initialize spatial grid with D points between T[0] and T[1] - COMPLEX * x_grid = NULL; - x_grid = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(x_grid, ret_code, leave_fun); + COMPLEX * t_grid = NULL; + t_grid = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(t_grid, ret_code, leave_fun); const COMPLEX eps_t = (T[1] - T[0])/(D - 1); for (UINT n=0; n Date: Tue, 17 Feb 2026 17:13:11 +0100 Subject: [PATCH 18/82] Adapt Credentials fnft_kdvv_inverse.c --- src/fnft_kdvv_inverse.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index e446197..4ffed4c 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -15,9 +15,12 @@ * * Contributors: * Sander Wahls (KIT) 2026 + * Fabian Fischer (Hiwi KIT) 2026 + * + * Following Algorithms are based on Crum Transformation Matlab project + * of P.J. Prins 2020 */ -// TODO comment matlab #define FNFT_ENABLE_SHORT_NAMES From 47a10a5fa93dbb057d55a86afd448d01cd56f425 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 17 Feb 2026 17:39:24 +0100 Subject: [PATCH 19/82] Add in find_first_positive_value condition if positive value was found --- src/fnft_kdvv_inverse.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 4ffed4c..f7fe04e 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -32,15 +32,22 @@ static INT find_first_positive_value( UINT * const i_pos_ptr) { INT ret_code = SUCCESS; + UINT is_positive_value_found = 0; for (UINT i = 0; i < length; i++){ if (CREAL(array[i]) > 0){ *i_pos_ptr = i; + is_positive_value_found = 1; break; } } - // TODO wenn kein positive value?! + if (is_positive_value_found == 0){ + ret_code = FNFT_EC_INVALID_ARGUMENT; + } + else { + ret_code = SUCCESS; + } return ret_code; } From 25e70c841f520df379e14949341ffbe720a2f5c2 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 23 Feb 2026 09:49:41 +0100 Subject: [PATCH 20/82] Add in fnft_kdvv_inverse_test_1.c function run_test and swap major procedure from main into it. Same style then fnft_kdvv_test_accounting_99sech.c --- .../fnft_kdvv_inverse_test_1.c | 94 ++++++++++++++++--- 1 file changed, 83 insertions(+), 11 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index 538d6aa..c9ed514 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -24,15 +24,20 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" -INT main() + +static INT run_test(const UINT D, + REAL err_bnd_bound_states, + REAL err_bnd_spurious_bound_states, + REAL err_bnd_normconst, + REAL err_bnd_contspec) + { - INT ret_code = SUCCESS; +INT ret_code = SUCCESS; // General parameters UINT M = 0; REAL XI[2] = {-2.0, 2.0}; - UINT D = 1001; - REAL T[2] = {-10.0, 10.0}; + REAL T[2] = {-20.0, 20.0}; COMPLEX * q = NULL; q = malloc(D * sizeof(COMPLEX)); @@ -57,7 +62,7 @@ INT main() contspec_r = malloc(M * sizeof(COMPLEX)); CHECK_NOMEM(contspec_r, ret_code, leave_fun); - UINT K_r = K_i; + UINT K_r = D; COMPLEX * bound_states_r = NULL; bound_states_r = malloc(K_r * sizeof(COMPLEX)); @@ -72,8 +77,44 @@ INT main() ret_code = fnft_kdvv(D, q, T, M, contspec_r, XI, &K_r, bound_states_r, normconsts_r, &opts); CHECK_RETCODE(ret_code, leave_fun); + // TODO: sorting Eigenvalues from bound_states_r? + // Check result // bound_states_r is in ascending order -> start from the end +// for (UINT i=0; i 0) && is_bsr_pure_imaginary; +// UINT is_bsr_in_tolerance = CABS(bsr - bsi)/CABS(bsr) < err_bnd_bound_states; + +// UINT are_spurious_bound_states_small = err_bnd_spurious_bound_states; + +// UINT is_ncr_real = CABS(CIMAG(ncr)) < 1e-9; +// UINT is_ncr_in_tolerance = CABS(ncr - nci)/CABS(ncr) < err_bnd_normconst; + +// // misc_hausdorff_dist + +// UINT is_contspec_small = err_bnd_contspec; + + +// if (is_bsr_pure_imaginary && +// is_bsr_positive_imaginary && +// is_bsr_in_tolerance && +// is_ncr_real && +// is_ncr_in_tolerance) +// { +// ret_code = SUCCESS; +// } +// else { +// ret_code = FNFT_EC_TEST_FAILED; +// break; +// } +// } + for (UINT i=0; i Date: Mon, 23 Feb 2026 11:51:12 +0100 Subject: [PATCH 21/82] Add misc_hausdorff_dist_normed based on misc_hausdorff_dist --- include/private/fnft__misc.h | 20 ++++++++++++++++++++ src/private/fnft__misc.c | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+) diff --git a/include/private/fnft__misc.h b/include/private/fnft__misc.h index 5dc9b2a..7a3ae9f 100644 --- a/include/private/fnft__misc.h +++ b/include/private/fnft__misc.h @@ -98,6 +98,25 @@ FNFT_REAL fnft__misc_hausdorff_dist(const FNFT_UINT lenA, FNFT_COMPLEX const * const vecA, const FNFT_UINT lenB, FNFT_COMPLEX const * const vecB); + +/** + * @brief Hausdorff distance between two vectors, normed by the considered element. + * + * @ingroup misc + * This function computes the Hausdorff distance between two vectors vecA and vecB. + * By calculating the maximum distance, the distance is divided by the absolute of the + * according element. + * @param[in] lenA Length of vector vecA. + * @param[in] vecA Complex vector of length lenA. + * @param[in] lenB length of vector vecB. + * @param[in] vecB Complex vector of length lenB. + * @return Returns the real valued Hausdorff distance between the vectors vecA and vecB. + */ +FNFT_REAL fnft__misc_hausdorff_dist_normed(const FNFT_UINT lenA, + FNFT_COMPLEX const * const vecA, const FNFT_UINT lenB, + FNFT_COMPLEX const * const vecB); + + /** * @brief Hyperbolic secant. * @@ -453,6 +472,7 @@ static inline FNFT_INT fnft__misc_normalize_vector(const FNFT_UINT len, FNFT_COM #define misc_rel_err(...) fnft__misc_rel_err(__VA_ARGS__) #define misc_rel_err_real(...) fnft__misc_rel_err_real(__VA_ARGS__) #define misc_hausdorff_dist(...) fnft__misc_hausdorff_dist(__VA_ARGS__) +#define misc_hausdorff_dist_normed(...) fnft__misc_hausdorff_dist_normed(__VA_ARGS__) #define misc_sech(...) fnft__misc_sech(__VA_ARGS__) #define misc_l2norm2(...) fnft__misc_l2norm2(__VA_ARGS__) #define misc_filter(...) fnft__misc_filter(__VA_ARGS__) diff --git a/src/private/fnft__misc.c b/src/private/fnft__misc.c index eb53af8..e478a1c 100644 --- a/src/private/fnft__misc.c +++ b/src/private/fnft__misc.c @@ -112,6 +112,38 @@ REAL misc_hausdorff_dist(const UINT lenA, return max_dist; } +REAL misc_hausdorff_dist_normed(const UINT lenA, + COMPLEX const * const vecA, const UINT lenB, + COMPLEX const * const vecB) +{ + UINT i, j; + double tmp, dist, max_dist = -1.0; + + for (i=0; i max_dist) + max_dist = dist; + } + + for (j=0; j max_dist) + max_dist = dist; + } + + return max_dist; +} + COMPLEX misc_sech(COMPLEX Z) { return 2.0 / (CEXP(Z) + CEXP(-Z)); From 679ee6d184c6705a34cc53c1e14be334fc454db1 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 23 Feb 2026 11:52:50 +0100 Subject: [PATCH 22/82] Add test conditions in fnft_kdvv_inverse_test_1 and test for quadratic convergence now --- .../fnft_kdvv_inverse_test_1.c | 195 +++++++++++------- 1 file changed, 117 insertions(+), 78 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index c9ed514..b03a79d 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -24,6 +24,39 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" +static INT print_test_results( COMPLEX * bound_states_r, + COMPLEX * normconsts_r, + COMPLEX * contspec_r, + REAL * XI, + const UINT M, + const UINT D, + const UINT K_r) +{ + printf("Number of samples:\n D = %u\n", (unsigned int)D); + + FNFT_REAL eps_xi = (XI[1] - XI[0]) / (M - 1); + printf("Continuous spectrum:\n"); + for (FNFT_UINT i=0; i start from the end -// for (UINT i=0; i 0) && is_bsr_pure_imaginary; -// UINT is_bsr_in_tolerance = CABS(bsr - bsi)/CABS(bsr) < err_bnd_bound_states; - -// UINT are_spurious_bound_states_small = err_bnd_spurious_bound_states; - -// UINT is_ncr_real = CABS(CIMAG(ncr)) < 1e-9; -// UINT is_ncr_in_tolerance = CABS(ncr - nci)/CABS(ncr) < err_bnd_normconst; - -// // misc_hausdorff_dist - -// UINT is_contspec_small = err_bnd_contspec; - - -// if (is_bsr_pure_imaginary && -// is_bsr_positive_imaginary && -// is_bsr_in_tolerance && -// is_ncr_real && -// is_ncr_in_tolerance) -// { -// ret_code = SUCCESS; -// } -// else { -// ret_code = FNFT_EC_TEST_FAILED; -// break; -// } -// } + // -- Check results -- + // Simple general tests for (UINT i=0; i 0) && is_bsr_pure_imaginary; - UINT is_bsr_in_tolerance = CABS(bsr - bsi)/CABS(bsr) < 1e-3; UINT is_ncr_real = CABS(CIMAG(ncr)) < 1e-8; - UINT is_ncr_in_tolerance = CABS(ncr - nci)/CABS(ncr) < 1e-2; if (is_bsr_pure_imaginary && is_bsr_positive_imaginary && - is_bsr_in_tolerance && - is_ncr_real && - is_ncr_in_tolerance) + is_ncr_real) { ret_code = SUCCESS; } @@ -141,6 +135,53 @@ INT ret_code = SUCCESS; break; } } + CHECK_RETCODE(ret_code, leave_fun); + + // Check matching of the computed bound states + // bound_states_r is sorted in ascending order, K_r is the number of found bound states by fnft_kdvv + // take only the last K_i bound states, because they are the bigger ones and more likely the ones + // which corresponds to the initial given bound states + COMPLEX * const candidate_bound_states_r_ptr = &bound_states_r[K_r-K_i]; + REAL hausdorff_dist_bound_states = misc_hausdorff_dist_normed( K_i, bound_states_i, K_i, + candidate_bound_states_r_ptr); + UINT is_bsr_in_tolerance = hausdorff_dist_bound_states < err_bnd_bound_states; + + COMPLEX * const candidate_normconsts_r_ptr = &normconsts_r[K_r-K_i]; + REAL hausdorff_dist_normconsts = misc_hausdorff_dist_normed(K_i, normconsts_i, K_i, + candidate_normconsts_r_ptr); + UINT is_ncr_in_tolerance = hausdorff_dist_normconsts < err_bnd_normconst; + + UINT is_contspec_small = 1; + for (UINT i=0; i err_bnd_contspec){ + is_contspec_small = 0; + } + } + + UINT are_spurious_bound_states_small = 1; + for (UINT i=0; i<(K_r-K_i); i++){ + if (CABS(bound_states_r[i]) > err_bnd_spurious_bound_states){ + are_spurious_bound_states_small = 0; + } + } + + if (is_bsr_in_tolerance && + is_ncr_in_tolerance && + is_contspec_small && + are_spurious_bound_states_small) + { + ret_code = SUCCESS; + } + else { + ret_code = FNFT_EC_TEST_FAILED; + } + +#ifdef DEBUG + printf("Number of bound_states:\n K_r = %u\n", (unsigned int)K_r); + printf("Hausdorff dist bound states:\n dist = %f\n", hausdorff_dist_bound_states); + printf("Hausdorff dist normconsts:\n dist = %f\n", hausdorff_dist_normconsts); + print_test_results(bound_states_r, normconsts_r, contspec_r, XI, M, D, K_r); +#endif leave_fun: @@ -157,46 +198,44 @@ INT ret_code = SUCCESS; } -static INT print_test_results( COMPLEX * bound_states_r, - COMPLEX * normconsts_r, - COMPLEX * contspec_r, - REAL * XI, - const UINT M, - const UINT D, - const UINT K_r) +INT main() { - printf("Number of samples:\n D = %u\n", (unsigned int)D); + INT ret_code = SUCCESS; + + UINT D = 256; + REAL err_bnd_bound_states = 1.4e-3; + REAL err_bnd_spurious_bound_states = 1e-2; + REAL err_bnd_normconst = 4.5e-2; + REAL err_bnd_contspec = 1e-1; - FNFT_REAL eps_xi = (XI[1] - XI[0]) / (M - 1); - printf("Continuous spectrum:\n"); - for (FNFT_UINT i=0; i Date: Mon, 2 Mar 2026 10:16:49 +0100 Subject: [PATCH 23/82] Add struct fnft_kdvv_params and implement it into functions in fnft_kdvv_inverse_test_1 --- .../fnft_kdvv_inverse_test_1.c | 100 +++++++++++------- 1 file changed, 60 insertions(+), 40 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index b03a79d..ef9d1aa 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -24,6 +24,19 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" + +typedef struct { + UINT D; + REAL T[2]; + UINT K; + COMPLEX * bound_states; + COMPLEX * normconsts; + UINT M; + REAL XI[2]; + REAL * contspec; +} fnft_kdvv_params; + + static INT print_test_results( COMPLEX * bound_states_r, COMPLEX * normconsts_r, COMPLEX * contspec_r, @@ -58,44 +71,36 @@ static INT print_test_results( COMPLEX * bound_states_r, -static INT run_test(const UINT D, +static INT run_test(const fnft_kdvv_params params_i, REAL err_bnd_bound_states, REAL err_bnd_spurious_bound_states, REAL err_bnd_normconst, REAL err_bnd_contspec) { -INT ret_code = SUCCESS; - - // General parameters - UINT M = 10; - REAL XI[2] = {-2.0, 2.0}; - REAL T[2] = {-20.0, 20.0}; + INT ret_code = SUCCESS; COMPLEX * q = NULL; - q = malloc(D * sizeof(COMPLEX)); + q = malloc(params_i.D * sizeof(COMPLEX)); CHECK_NOMEM(q, ret_code, leave_fun); - // Initialize variables for a specific test case without continuous spectrum - COMPLEX * contspec_i = NULL; - contspec_i = malloc(M * sizeof(COMPLEX)); - const UINT K_i = 5; - COMPLEX bound_states_i[5] = {I*SQRT(1.0/2.0), I*SQRT(2.0/2.0), I*SQRT(3.0/2.0), I*SQRT(4.0/2.0), I*SQRT(5.0/2.0)}; - COMPLEX normconsts_i[5] = {1*10, -1*0.1, 1*1, -1*1e-5, 1*1e7}; - - // const UINT K_i = 1; - // COMPLEX bound_states_i[1] = {SQRT(1.0/2)*I}; - // COMPLEX normconsts_i[1] = {1*10}; - - ret_code = fnft_kdvv_inverse(M, contspec_i, XI, K_i, bound_states_i, normconsts_i, D, q, T, NULL); + ret_code = fnft_kdvv_inverse( params_i.M, + params_i.contspec, + params_i.XI, + params_i.K, + params_i.bound_states, + params_i.normconsts, + params_i.D, + q, + params_i.T, NULL); CHECK_RETCODE(ret_code, leave_fun); // Prepare forward fnft_kdvv COMPLEX * contspec_r = NULL; - contspec_r = malloc(M * sizeof(COMPLEX)); + contspec_r = malloc(params_i.M * sizeof(COMPLEX)); CHECK_NOMEM(contspec_r, ret_code, leave_fun); - UINT K_r = D; + UINT K_r = params_i.D; COMPLEX * bound_states_r = NULL; bound_states_r = malloc(K_r * sizeof(COMPLEX)); @@ -107,15 +112,15 @@ INT ret_code = SUCCESS; fnft_kdvv_opts_t opts = fnft_kdvv_default_opts(); - ret_code = fnft_kdvv(D, q, T, M, contspec_r, XI, &K_r, bound_states_r, normconsts_r, &opts); + ret_code = fnft_kdvv(params_i.D, q, params_i.T, params_i.M, contspec_r, params_i.XI, &K_r, bound_states_r, normconsts_r, &opts); CHECK_RETCODE(ret_code, leave_fun); // -- Check results -- // Simple general tests - for (UINT i=0; i err_bnd_contspec){ is_contspec_small = 0; } } UINT are_spurious_bound_states_small = 1; - for (UINT i=0; i<(K_r-K_i); i++){ + for (UINT i=0; i<(K_r-params_i.K); i++){ if (CABS(bound_states_r[i]) > err_bnd_spurious_bound_states){ are_spurious_bound_states_small = 0; } @@ -186,7 +191,6 @@ INT ret_code = SUCCESS; leave_fun: free(q); - free(contspec_i); free(contspec_r); free(bound_states_r); free(normconsts_r); @@ -201,39 +205,55 @@ INT ret_code = SUCCESS; INT main() { INT ret_code = SUCCESS; + + COMPLEX bound_states_i[5] = {I*SQRT(1.0/2.0), I*SQRT(2.0/2.0), I*SQRT(3.0/2.0), I*SQRT(4.0/2.0), I*SQRT(5.0/2.0)}; + COMPLEX normconsts_i[5] = {1*10, -1*0.1, 1*1, -1*1e-5, 1*1e7}; + + COMPLEX * contspec_i = NULL; + + fnft_kdvv_params kdvv_parameters = { + .D = 256, + .T = {-20.0, 20.0}, + .K = 5, + .bound_states = bound_states_i, + .normconsts = normconsts_i, + .M = 10, + .XI = {-2.0, 2.0}, + .contspec = contspec_i + }; + - UINT D = 256; REAL err_bnd_bound_states = 1.4e-3; REAL err_bnd_spurious_bound_states = 1e-2; REAL err_bnd_normconst = 4.5e-2; REAL err_bnd_contspec = 1e-1; - ret_code = run_test(D, err_bnd_bound_states, err_bnd_spurious_bound_states, + ret_code = run_test(kdvv_parameters, err_bnd_bound_states, err_bnd_spurious_bound_states, err_bnd_normconst, err_bnd_contspec); CHECK_RETCODE(ret_code, leave_fun); // Check quadratic convergence - D *= 2; + kdvv_parameters.D *= 2; err_bnd_bound_states /= 4; err_bnd_spurious_bound_states /= 4; err_bnd_normconst = 1.1e-2; err_bnd_contspec = 5e-3; - ret_code = run_test(D, err_bnd_bound_states, err_bnd_spurious_bound_states, + ret_code = run_test(kdvv_parameters, err_bnd_bound_states, err_bnd_spurious_bound_states, err_bnd_normconst, err_bnd_contspec); CHECK_RETCODE(ret_code, leave_fun); - D *= 2; + kdvv_parameters.D *= 2; err_bnd_bound_states /= 4; err_bnd_spurious_bound_states /= 4; err_bnd_normconst /= 4; err_bnd_contspec /= 4; - ret_code = run_test(D, err_bnd_bound_states, err_bnd_spurious_bound_states, + ret_code = run_test(kdvv_parameters, err_bnd_bound_states, err_bnd_spurious_bound_states, err_bnd_normconst, err_bnd_contspec); CHECK_RETCODE(ret_code, leave_fun); From e1b5e415619df4de103d2646e54be10b3fe033f1 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 2 Mar 2026 11:33:55 +0100 Subject: [PATCH 24/82] Refactoring: Move test prcedure fnft__kdvv_inverse_testcases_get_spectrum_of_inverse, which calculates the spectrum of the inverse of some bound states to ffnft_kdvv_inverse_testcases.c, so that it can be used for several tests --- .../private/fnft__kdvv_inverse_testcases.h | 46 ++++- src/private/fnft__kdvv_inverse_testcases.c | 189 +++++++++++++++++ .../fnft_kdvv_inverse_test_1.c | 192 +----------------- 3 files changed, 243 insertions(+), 184 deletions(-) create mode 100644 src/private/fnft__kdvv_inverse_testcases.c diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h index d32700c..0592a86 100644 --- a/include/private/fnft__kdvv_inverse_testcases.h +++ b/include/private/fnft__kdvv_inverse_testcases.h @@ -20,8 +20,50 @@ #ifndef FNFT__KDVV_INVERSE_TESTCASES_H #define FNFT__KDVV_INVERSE_TESTCASES_H -#include "fnft_kdvv.h" #include "fnft_kdvv_inverse.h" -#include "fnft_numtypes.h" +#include "fnft_kdvv.h" + +/** + * @struct fnft_kdvv_params + */ +typedef struct { + UINT D; + REAL T[2]; + UINT K; + COMPLEX * bound_states; + COMPLEX * normconsts; + UINT M; + REAL XI[2]; + COMPLEX * contspec; +} fnft_kdvv_params; + + +/** + * @brief Routine to run tests for \link fnft_kdvv_inverse \endlink. + * + * This routine is used by the tests for \link fnft_kdvv_inverse \endlink. + * + * @param[in] params_i \link fnft_kdvv_params \endlink + * @param[in] err_bnd_bound_states + * @param[in] err_bnd_spurious_bound_states + * @param[in] err_bnd_normconst + * @param[in] err_bnd_contspec + * @return If all errors stay below bounds the routine + * \link FNFT_SUCCESS \endlink. Otherwise, it returns an error code + * (normally, \link FNFT_EC_TEST_FAILED \endlink). + * + * @ingroup kdv + */ +FNFT_INT fnft__kdvv_inverse_testcases_get_spectrum_of_inverse( + const fnft_kdvv_params params_i, + const FNFT_REAL err_bnd_bound_states, + const FNFT_REAL err_bnd_spurious_bound_states, + const FNFT_REAL err_bnd_normconst, + const FNFT_REAL err_bnd_contspec); + + +#ifdef FNFT_ENABLE_SHORT_NAMES +#define kdvv_testcases_get_spectrum_of_inverse(...) fnft__kdvv_inverse_testcases_get_spectrum_of_inverse(__VA_ARGS__) +#endif #endif diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c new file mode 100644 index 0000000..906c0de --- /dev/null +++ b/src/private/fnft__kdvv_inverse_testcases.c @@ -0,0 +1,189 @@ +/* +* This file is part of FNFT. +* +* FNFT is free software; you can redistribute it and/or +* modify it under the terms of the version 2 of the GNU General +* Public License as published by the Free Software Foundation. +* +* FNFT is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Contributors: +* Fabian Fischer (Hiwi KIT) 2026 +*/ + +#define FNFT_ENABLE_SHORT_NAMES + +#include + +#include "fnft__kdvv_inverse_testcases.h" +#include "fnft__errwarn.h" +#include "fnft__misc.h" + +static INT print_test_results( COMPLEX * bound_states_r, + COMPLEX * normconsts_r, + COMPLEX * contspec_r, + REAL * XI, + const UINT M, + const UINT D, + const UINT K_r) +{ + printf("Number of samples:\n D = %u\n", (unsigned int)D); + + FNFT_REAL eps_xi = (XI[1] - XI[0]) / (M - 1); + printf("Continuous spectrum:\n"); + for (FNFT_UINT i=0; i 0) && is_bsr_pure_imaginary; + + UINT is_ncr_real = CABS(CIMAG(ncr)) < 1e-8; + + if (is_bsr_pure_imaginary && + is_bsr_positive_imaginary && + is_ncr_real) + { + ret_code = SUCCESS; + } + else { + ret_code = FNFT_EC_TEST_FAILED; + break; + } + } + CHECK_RETCODE(ret_code, leave_fun); + + // Check matching of the computed bound states + // bound_states_r is sorted in ascending order, K_r is the number of found bound states by fnft_kdvv + // take only the last params_i.K bound states, because they are the bigger ones and more likely the ones + // which corresponds to the initial given bound states + COMPLEX * const candidate_bound_states_r_ptr = &bound_states_r[K_r-params_i.K]; + REAL hausdorff_dist_bound_states = misc_hausdorff_dist_normed( params_i.K, params_i.bound_states, params_i.K, + candidate_bound_states_r_ptr); + UINT is_bsr_in_tolerance = hausdorff_dist_bound_states < err_bnd_bound_states; + + COMPLEX * const candidate_normconsts_r_ptr = &normconsts_r[K_r-params_i.K]; + REAL hausdorff_dist_normconsts = misc_hausdorff_dist_normed(params_i.K, params_i.normconsts, params_i.K, + candidate_normconsts_r_ptr); + UINT is_ncr_in_tolerance = hausdorff_dist_normconsts < err_bnd_normconst; + + UINT is_contspec_small = 1; + for (UINT i=0; i err_bnd_contspec){ + is_contspec_small = 0; + } + } + + UINT are_spurious_bound_states_small = 1; + for (UINT i=0; i<(K_r-params_i.K); i++){ + if (CABS(bound_states_r[i]) > err_bnd_spurious_bound_states){ + are_spurious_bound_states_small = 0; + } + } + + if (is_bsr_in_tolerance && + is_ncr_in_tolerance && + is_contspec_small && + are_spurious_bound_states_small) + { + ret_code = SUCCESS; + } + else { + ret_code = FNFT_EC_TEST_FAILED; + } + +#ifdef DEBUG + printf("Number of bound_states:\n K_r = %u\n", (unsigned int)K_r); + printf("Hausdorff dist bound states:\n dist = %f\n", hausdorff_dist_bound_states); + printf("Hausdorff dist normconsts:\n dist = %f\n", hausdorff_dist_normconsts); + print_test_results(bound_states_r, normconsts_r, contspec_r, XI, M, D, K_r); +#endif + + +leave_fun: + free(q); + free(contspec_r); + free(bound_states_r); + free(normconsts_r); + + if (ret_code != SUCCESS) + return EXIT_FAILURE; + else + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index ef9d1aa..9d39545 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -14,6 +14,7 @@ * along with this program. If not, see . * * Contributors: +* Fabian Fischer (Hiwi KIT) 2026 */ #define FNFT_ENABLE_SHORT_NAMES @@ -25,182 +26,6 @@ #include "fnft__misc.h" -typedef struct { - UINT D; - REAL T[2]; - UINT K; - COMPLEX * bound_states; - COMPLEX * normconsts; - UINT M; - REAL XI[2]; - REAL * contspec; -} fnft_kdvv_params; - - -static INT print_test_results( COMPLEX * bound_states_r, - COMPLEX * normconsts_r, - COMPLEX * contspec_r, - REAL * XI, - const UINT M, - const UINT D, - const UINT K_r) -{ - printf("Number of samples:\n D = %u\n", (unsigned int)D); - - FNFT_REAL eps_xi = (XI[1] - XI[0]) / (M - 1); - printf("Continuous spectrum:\n"); - for (FNFT_UINT i=0; i 0) && is_bsr_pure_imaginary; - - UINT is_ncr_real = CABS(CIMAG(ncr)) < 1e-8; - - if (is_bsr_pure_imaginary && - is_bsr_positive_imaginary && - is_ncr_real) - { - ret_code = SUCCESS; - } - else { - ret_code = FNFT_EC_TEST_FAILED; - break; - } - } - CHECK_RETCODE(ret_code, leave_fun); - - // Check matching of the computed bound states - // bound_states_r is sorted in ascending order, K_r is the number of found bound states by fnft_kdvv - // take only the last params_i.K bound states, because they are the bigger ones and more likely the ones - // which corresponds to the initial given bound states - COMPLEX * const candidate_bound_states_r_ptr = &bound_states_r[K_r-params_i.K]; - REAL hausdorff_dist_bound_states = misc_hausdorff_dist_normed( params_i.K, params_i.bound_states, params_i.K, - candidate_bound_states_r_ptr); - UINT is_bsr_in_tolerance = hausdorff_dist_bound_states < err_bnd_bound_states; - - COMPLEX * const candidate_normconsts_r_ptr = &normconsts_r[K_r-params_i.K]; - REAL hausdorff_dist_normconsts = misc_hausdorff_dist_normed(params_i.K, params_i.normconsts, params_i.K, - candidate_normconsts_r_ptr); - UINT is_ncr_in_tolerance = hausdorff_dist_normconsts < err_bnd_normconst; - - UINT is_contspec_small = 1; - for (UINT i=0; i err_bnd_contspec){ - is_contspec_small = 0; - } - } - - UINT are_spurious_bound_states_small = 1; - for (UINT i=0; i<(K_r-params_i.K); i++){ - if (CABS(bound_states_r[i]) > err_bnd_spurious_bound_states){ - are_spurious_bound_states_small = 0; - } - } - - if (is_bsr_in_tolerance && - is_ncr_in_tolerance && - is_contspec_small && - are_spurious_bound_states_small) - { - ret_code = SUCCESS; - } - else { - ret_code = FNFT_EC_TEST_FAILED; - } - -#ifdef DEBUG - printf("Number of bound_states:\n K_r = %u\n", (unsigned int)K_r); - printf("Hausdorff dist bound states:\n dist = %f\n", hausdorff_dist_bound_states); - printf("Hausdorff dist normconsts:\n dist = %f\n", hausdorff_dist_normconsts); - print_test_results(bound_states_r, normconsts_r, contspec_r, XI, M, D, K_r); -#endif - - -leave_fun: - free(q); - free(contspec_r); - free(bound_states_r); - free(normconsts_r); - - if (ret_code != SUCCESS) - return EXIT_FAILURE; - else - return EXIT_SUCCESS; -} - INT main() { @@ -229,8 +54,9 @@ INT main() REAL err_bnd_contspec = 1e-1; - ret_code = run_test(kdvv_parameters, err_bnd_bound_states, err_bnd_spurious_bound_states, - err_bnd_normconst, err_bnd_contspec); + ret_code = kdvv_testcases_get_spectrum_of_inverse( kdvv_parameters, err_bnd_bound_states, + err_bnd_spurious_bound_states, + err_bnd_normconst, err_bnd_contspec); CHECK_RETCODE(ret_code, leave_fun); @@ -242,8 +68,9 @@ INT main() err_bnd_normconst = 1.1e-2; err_bnd_contspec = 5e-3; - ret_code = run_test(kdvv_parameters, err_bnd_bound_states, err_bnd_spurious_bound_states, - err_bnd_normconst, err_bnd_contspec); + ret_code = kdvv_testcases_get_spectrum_of_inverse( kdvv_parameters, err_bnd_bound_states, + err_bnd_spurious_bound_states, + err_bnd_normconst, err_bnd_contspec); CHECK_RETCODE(ret_code, leave_fun); @@ -253,8 +80,9 @@ INT main() err_bnd_normconst /= 4; err_bnd_contspec /= 4; - ret_code = run_test(kdvv_parameters, err_bnd_bound_states, err_bnd_spurious_bound_states, - err_bnd_normconst, err_bnd_contspec); + ret_code = kdvv_testcases_get_spectrum_of_inverse( kdvv_parameters, err_bnd_bound_states, + err_bnd_spurious_bound_states, + err_bnd_normconst, err_bnd_contspec); CHECK_RETCODE(ret_code, leave_fun); From 136c7d6a13fe3d725a1f829e5ebb72247e800d61 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 2 Mar 2026 15:50:16 +0100 Subject: [PATCH 25/82] Add fnft_kdvv_inverse_example.c --- examples/fnft_kdvv_inverse_example.c | 79 ++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) create mode 100644 examples/fnft_kdvv_inverse_example.c diff --git a/examples/fnft_kdvv_inverse_example.c b/examples/fnft_kdvv_inverse_example.c new file mode 100644 index 0000000..7e0a91d --- /dev/null +++ b/examples/fnft_kdvv_inverse_example.c @@ -0,0 +1,79 @@ +/* +* This file is part of FNFT. +* +* FNFT is free software; you can redistribute it and/or +* modify it under the terms of the version 2 of the GNU General +* Public License as published by the Free Software Foundation. +* +* FNFT is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Contributors: +* Fabian Fischer (Hiwi KIT) 2026 +*/ + +#define FNFT_ENABLE_SHORT_NAMES + +#include + +#include "fnft_kdvv_inverse.h" +#include "fnft__errwarn.h" +#include "fnft__misc.h" + + +#define K_I 25 + + +INT main() +{ + INT ret_code = SUCCESS; + + + // General parameters + UINT D = 256; + UINT M = 10; + REAL XI[2] = {-2.0, 2.0}; + REAL T[2] = {-20.0, 20.0}; + + COMPLEX * q = NULL; + q = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(q, ret_code, leave_fun); + + // Initialize variables for a specific test case without continuous spectrum + COMPLEX * contspec_i = NULL; + contspec_i = malloc(M * sizeof(COMPLEX)); + const UINT K_i = 25; + + + + COMPLEX bound_states_i[K_I] = { 50.0, 49.0, 48.0, 47.0, 46.0, + 40.0, 39.0, 38.0, 37.0, 36.0, + 30.0, 29.0, 28.0, 27.0, 26.0, + 20.0, 19.0, 18.0, 17.0, 16.0, + 10.0, 9.0, 8.0, 7.0, 6.0}; + + COMPLEX normconsts_i[K_I] = { 1e20, -1e-7, 1e5, -1e3, 1e1, + -1e0, 1e2, -1e4, 1e-6, -1e8, + 1e2, -1e4, 1e6, -1e8, 1e-10, + -1e7, 1e-6, -1e5, 1e-9, -1e-11, + 1e-3, -1e3, 1e2, -1e2, 1e1}; + + + for (UINT i = 0; i Date: Mon, 2 Mar 2026 15:52:46 +0100 Subject: [PATCH 26/82] fnft__misc.c: Fix Indexing Error in misc_hausdorff_dist_normed --- src/private/fnft__misc.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/private/fnft__misc.c b/src/private/fnft__misc.c index e478a1c..d2b31ad 100644 --- a/src/private/fnft__misc.c +++ b/src/private/fnft__misc.c @@ -133,7 +133,7 @@ REAL misc_hausdorff_dist_normed(const UINT lenA, for (j=0; j Date: Mon, 2 Mar 2026 15:53:46 +0100 Subject: [PATCH 27/82] fnft__kdvv_inverse_testcases: Improve Debugging Messages --- src/private/fnft__kdvv_inverse_testcases.c | 26 +++++++++++++++------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c index 906c0de..768a3f8 100644 --- a/src/private/fnft__kdvv_inverse_testcases.c +++ b/src/private/fnft__kdvv_inverse_testcases.c @@ -19,6 +19,8 @@ #define FNFT_ENABLE_SHORT_NAMES +#define DEBUG + #include #include "fnft__kdvv_inverse_testcases.h" @@ -67,6 +69,11 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, { INT ret_code = SUCCESS; + #ifdef DEBUG + // print initial parameters + print_test_results(params_i.bound_states, params_i.normconsts, params_i.contspec, params_i.XI, params_i.M, params_i.D, params_i.K); + #endif + COMPLEX * q = NULL; q = malloc(params_i.D * sizeof(COMPLEX)); CHECK_NOMEM(q, ret_code, leave_fun); @@ -102,6 +109,10 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, ret_code = fnft_kdvv(params_i.D, q, params_i.T, params_i.M, contspec_r, params_i.XI, &K_r, bound_states_r, normconsts_r, &opts); CHECK_RETCODE(ret_code, leave_fun); + #ifdef DEBUG + print_test_results(bound_states_r, normconsts_r, contspec_r, params_i.XI, params_i.M, params_i.D, K_r); + #endif + // -- Check results -- // Simple general tests @@ -141,7 +152,13 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, COMPLEX * const candidate_normconsts_r_ptr = &normconsts_r[K_r-params_i.K]; REAL hausdorff_dist_normconsts = misc_hausdorff_dist_normed(params_i.K, params_i.normconsts, params_i.K, candidate_normconsts_r_ptr); - UINT is_ncr_in_tolerance = hausdorff_dist_normconsts < err_bnd_normconst; + UINT is_ncr_in_tolerance = hausdorff_dist_normconsts < err_bnd_normconst; + + #ifdef DEBUG + printf("Number of bound_states:\n K_r = %u\n", (unsigned int)K_r); + printf("Hausdorff dist bound states:\n dist = %f\n", hausdorff_dist_bound_states); + printf("Hausdorff dist normconsts:\n dist = %f\n", hausdorff_dist_normconsts); + #endif UINT is_contspec_small = 1; for (UINT i=0; i Date: Mon, 2 Mar 2026 15:55:19 +0100 Subject: [PATCH 28/82] fnft_kdvv_inverse_test_1.c: Now allocate for contspec of kdvv_parameters an empty array --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index 9d39545..7e9a797 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -34,8 +34,6 @@ INT main() COMPLEX bound_states_i[5] = {I*SQRT(1.0/2.0), I*SQRT(2.0/2.0), I*SQRT(3.0/2.0), I*SQRT(4.0/2.0), I*SQRT(5.0/2.0)}; COMPLEX normconsts_i[5] = {1*10, -1*0.1, 1*1, -1*1e-5, 1*1e7}; - COMPLEX * contspec_i = NULL; - fnft_kdvv_params kdvv_parameters = { .D = 256, .T = {-20.0, 20.0}, @@ -44,7 +42,7 @@ INT main() .normconsts = normconsts_i, .M = 10, .XI = {-2.0, 2.0}, - .contspec = contspec_i + .contspec = malloc(10*sizeof(COMPLEX)) }; From 99089a8751c0768132ebce7118b38df1045dc401 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 2 Mar 2026 15:56:26 +0100 Subject: [PATCH 29/82] Add fnft_kdvv_inverse_test_2.c (Test with 25 Eigenvalues) --- .../fnft_kdvv_inverse_test_2.c | 109 ++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c new file mode 100644 index 0000000..25c518f --- /dev/null +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c @@ -0,0 +1,109 @@ +/* +* This file is part of FNFT. +* +* FNFT is free software; you can redistribute it and/or +* modify it under the terms of the version 2 of the GNU General +* Public License as published by the Free Software Foundation. +* +* FNFT is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Contributors: +* Fabian Fischer (Hiwi KIT) 2026 +*/ + +#define FNFT_ENABLE_SHORT_NAMES + +#include + +#include "fnft__kdvv_inverse_testcases.h" +#include "fnft__errwarn.h" +#include "fnft__misc.h" + +#define K_I 25 + +INT main() +{ + INT ret_code = SUCCESS; + + COMPLEX bound_states_i[K_I] = { 50.0, 49.0, 48.0, 47.0, 46.0, + 40.0, 39.0, 38.0, 37.0, 36.0, + 30.0, 29.0, 28.0, 27.0, 26.0, + 20.0, 19.0, 18.0, 17.0, 16.0, + 10.0, 9.0, 8.0, 7.0, 6.0}; + + COMPLEX normconsts_i[K_I] = { 1e20, -1e-7, 1e5, -1e3, 1e1, + -1e0, 1e2, -1e4, 1e-6, -1e8, + 1e2, -1e4, 1e6, -1e8, 1e-10, + -1e7, 1e-6, -1e5, 1e-9, -1e-11, + 1e-3, -1e3, 1e2, -1e2, 1e1}; + + for (UINT i = 0; i Date: Mon, 16 Mar 2026 09:38:27 +0100 Subject: [PATCH 30/82] Change error bounds ind fnft_kdvv_inverse_test_1.c --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index 7e9a797..1fdb940 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -63,8 +63,8 @@ INT main() kdvv_parameters.D *= 2; err_bnd_bound_states /= 4; err_bnd_spurious_bound_states /= 4; - err_bnd_normconst = 1.1e-2; - err_bnd_contspec = 5e-3; + err_bnd_normconst /= 4; + err_bnd_contspec /= 4; ret_code = kdvv_testcases_get_spectrum_of_inverse( kdvv_parameters, err_bnd_bound_states, err_bnd_spurious_bound_states, From 2f0db168f14a3769ec11020e80c0b003af332c3f Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 16 Mar 2026 10:19:26 +0100 Subject: [PATCH 31/82] Changing number of Eigenvalues in fnft_kdvv_inverse_test_2 and error bounds --- .../fnft_kdvv_inverse_test_2.c | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c index 25c518f..d9f82f7 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c @@ -25,23 +25,21 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" -#define K_I 25 +#define K_I 19 INT main() { INT ret_code = SUCCESS; - COMPLEX bound_states_i[K_I] = { 50.0, 49.0, 48.0, 47.0, 46.0, - 40.0, 39.0, 38.0, 37.0, 36.0, + COMPLEX bound_states_i[K_I] = { 40.0, 39.0, 38.0, 37.0, 36.0, 30.0, 29.0, 28.0, 27.0, 26.0, 20.0, 19.0, 18.0, 17.0, 16.0, - 10.0, 9.0, 8.0, 7.0, 6.0}; + 10.0, 9.0, 8.0, 7.0}; COMPLEX normconsts_i[K_I] = { 1e20, -1e-7, 1e5, -1e3, 1e1, -1e0, 1e2, -1e4, 1e-6, -1e8, 1e2, -1e4, 1e6, -1e8, 1e-10, - -1e7, 1e-6, -1e5, 1e-9, -1e-11, - 1e-3, -1e3, 1e2, -1e2, 1e1}; + -1e7, 1e-6, -1e5, 1e-9}; for (UINT i = 0; i Date: Mon, 16 Mar 2026 10:22:07 +0100 Subject: [PATCH 32/82] Add short test description --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c | 2 +- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index 1fdb940..7f8478b 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -25,7 +25,7 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" - +/* This is a simple test case */ INT main() { diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c index d9f82f7..edbd4ad 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c @@ -25,6 +25,8 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" +/* This is a testcase with many eigenvalues */ + #define K_I 19 INT main() From bc0dca071b2245ae9152be09b1d535124fe2f09b Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 16 Mar 2026 10:49:06 +0100 Subject: [PATCH 33/82] Rename print_test_results into kdvv_print_spectrum and make it a public function. --- .../private/fnft__kdvv_inverse_testcases.h | 9 ++++++++ src/private/fnft__kdvv_inverse_testcases.c | 21 +++++++------------ 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h index 0592a86..a718663 100644 --- a/include/private/fnft__kdvv_inverse_testcases.h +++ b/include/private/fnft__kdvv_inverse_testcases.h @@ -61,9 +61,18 @@ FNFT_INT fnft__kdvv_inverse_testcases_get_spectrum_of_inverse( const FNFT_REAL err_bnd_normconst, const FNFT_REAL err_bnd_contspec); +FNFT_INT fnft__kdvv_print_spectrum( FNFT_COMPLEX * bound_states_r, + FNFT_COMPLEX * normconsts_r, + FNFT_COMPLEX * contspec_r, + FNFT_REAL * XI, + const UINT M, + const UINT D, + const UINT K_r); + #ifdef FNFT_ENABLE_SHORT_NAMES #define kdvv_testcases_get_spectrum_of_inverse(...) fnft__kdvv_inverse_testcases_get_spectrum_of_inverse(__VA_ARGS__) +#define kdvv_print_spectrum(...) fnft__kdvv_print_spectrum(__VA_ARGS__) #endif #endif diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c index 768a3f8..acd6557 100644 --- a/src/private/fnft__kdvv_inverse_testcases.c +++ b/src/private/fnft__kdvv_inverse_testcases.c @@ -27,13 +27,13 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" -static INT print_test_results( COMPLEX * bound_states_r, - COMPLEX * normconsts_r, - COMPLEX * contspec_r, - REAL * XI, - const UINT M, - const UINT D, - const UINT K_r) +INT kdvv_print_spectrum(COMPLEX * bound_states_r, + COMPLEX * normconsts_r, + COMPLEX * contspec_r, + REAL * XI, + const UINT M, + const UINT D, + const UINT K_r) { printf("Number of samples:\n D = %u\n", (unsigned int)D); @@ -69,11 +69,6 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, { INT ret_code = SUCCESS; - #ifdef DEBUG - // print initial parameters - print_test_results(params_i.bound_states, params_i.normconsts, params_i.contspec, params_i.XI, params_i.M, params_i.D, params_i.K); - #endif - COMPLEX * q = NULL; q = malloc(params_i.D * sizeof(COMPLEX)); CHECK_NOMEM(q, ret_code, leave_fun); @@ -110,7 +105,7 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, CHECK_RETCODE(ret_code, leave_fun); #ifdef DEBUG - print_test_results(bound_states_r, normconsts_r, contspec_r, params_i.XI, params_i.M, params_i.D, K_r); + kdvv_print_spectrum(bound_states_r, normconsts_r, contspec_r, params_i.XI, params_i.M, params_i.D, K_r); #endif // -- Check results -- From 198b7e104bb3fdd74b3129d7db176c43033b1f66 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 16 Mar 2026 10:50:17 +0100 Subject: [PATCH 34/82] Add fnft_kdvv_inverse_test_3.c for test with assymetric window --- .../fnft_kdvv_inverse_test_3.c | 125 ++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c new file mode 100644 index 0000000..6360ce3 --- /dev/null +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c @@ -0,0 +1,125 @@ +/* +* This file is part of FNFT. +* +* FNFT is free software; you can redistribute it and/or +* modify it under the terms of the version 2 of the GNU General +* Public License as published by the Free Software Foundation. +* +* FNFT is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Contributors: +* Fabian Fischer (Hiwi KIT) 2026 +*/ + +#define FNFT_ENABLE_SHORT_NAMES + +#include + +#include "fnft__kdvv_inverse_testcases.h" +#include "fnft__errwarn.h" +#include "fnft__misc.h" + +/* This is a testcase with an assymetric window */ + +#define K_I 8 +#define DEBUG + +INT main() +{ + INT ret_code = SUCCESS; + + COMPLEX bound_states_i[K_I] = { 40.0, 35.0, 28.0, 23.0, 19.0, + 16.0, 12.0, 4.0}; + + COMPLEX normconsts_i[K_I] = { 1e-3, -1e-17, 1e-15, -1e-20, 1e-1, + -1e-5, 1e-12, -1e-14}; + + for (UINT i = 0; i Date: Mon, 16 Mar 2026 11:00:02 +0100 Subject: [PATCH 35/82] Changing error bounds in fnft_kdvv_inverse_test_3.c --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c index 6360ce3..63b2695 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c @@ -57,20 +57,22 @@ INT main() }; #ifdef DEBUG - // print initial parameters + printf("Initial parameters: \n"); kdvv_print_spectrum(kdvv_parameters.bound_states, kdvv_parameters.normconsts, kdvv_parameters.contspec, kdvv_parameters.XI, kdvv_parameters.M, kdvv_parameters.D, kdvv_parameters.K); #endif + + #ifdef DEBUG printf("-- 1. Iteration: -- \n"); #endif - REAL err_bnd_bound_states = 2e-3; - REAL err_bnd_spurious_bound_states = 0.1; - REAL err_bnd_normconst = 1.6e-1; + REAL err_bnd_bound_states = 3.8e-3; + REAL err_bnd_spurious_bound_states = 0.12; + REAL err_bnd_normconst = 2e-1; REAL err_bnd_contspec = 2.5e-1; ret_code = kdvv_testcases_get_spectrum_of_inverse( kdvv_parameters, err_bnd_bound_states, @@ -80,6 +82,7 @@ INT main() CHECK_RETCODE(ret_code, leave_fun); + #ifdef DEBUG printf("-- 2. Iteration: -- \n"); #endif @@ -97,7 +100,8 @@ INT main() CHECK_RETCODE(ret_code, leave_fun); - + + #ifdef DEBUG printf("-- 3. Iteration: -- \n"); #endif From f1e47c1194cf2128eb50b32cd236792f344749c6 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 16 Mar 2026 11:14:38 +0100 Subject: [PATCH 36/82] Changing parameters in fnft_kdvv_inverse_example.c --- examples/fnft_kdvv_inverse_example.c | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/examples/fnft_kdvv_inverse_example.c b/examples/fnft_kdvv_inverse_example.c index 7e0a91d..432b14c 100644 --- a/examples/fnft_kdvv_inverse_example.c +++ b/examples/fnft_kdvv_inverse_example.c @@ -26,7 +26,7 @@ #include "fnft__misc.h" -#define K_I 25 +#define K_I 19 INT main() @@ -38,7 +38,7 @@ INT main() UINT D = 256; UINT M = 10; REAL XI[2] = {-2.0, 2.0}; - REAL T[2] = {-20.0, 20.0}; + REAL T[2] = {-18.0, 12.0}; COMPLEX * q = NULL; q = malloc(D * sizeof(COMPLEX)); @@ -47,28 +47,24 @@ INT main() // Initialize variables for a specific test case without continuous spectrum COMPLEX * contspec_i = NULL; contspec_i = malloc(M * sizeof(COMPLEX)); - const UINT K_i = 25; - - COMPLEX bound_states_i[K_I] = { 50.0, 49.0, 48.0, 47.0, 46.0, - 40.0, 39.0, 38.0, 37.0, 36.0, + COMPLEX bound_states_i[K_I] = { 40.0, 39.0, 38.0, 37.0, 36.0, 30.0, 29.0, 28.0, 27.0, 26.0, 20.0, 19.0, 18.0, 17.0, 16.0, - 10.0, 9.0, 8.0, 7.0, 6.0}; + 10.0, 9.0, 8.0, 7.0}; COMPLEX normconsts_i[K_I] = { 1e20, -1e-7, 1e5, -1e3, 1e1, -1e0, 1e2, -1e4, 1e-6, -1e8, 1e2, -1e4, 1e6, -1e8, 1e-10, - -1e7, 1e-6, -1e5, 1e-9, -1e-11, - 1e-3, -1e3, 1e2, -1e2, 1e1}; + -1e7, 1e-6, -1e5, 1e-9}; for (UINT i = 0; i Date: Mon, 16 Mar 2026 11:19:46 +0100 Subject: [PATCH 37/82] Change Parameters in mex_fnft_kdvv_inverse_example_1.m --- examples/mex_fnft_kdvv_inverse_example_1.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index 09e5cf0..7d8f821 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -20,10 +20,10 @@ contspec = []; XI = [0 1]; -bound_states = 1i*[1, 2, 3, 4] /2; +bound_states = 1i*sqrt([1, 2, 3, 4] /2); normconsts = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); D = 1001; -T = [-20 20]; +T = [-10 10]; q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T) t = linspace(T(1), T(2), D); From d2d96578592f025181faa56502d27ee36793e920 Mon Sep 17 00:00:00 2001 From: FF Date: Mon, 16 Mar 2026 11:20:28 +0100 Subject: [PATCH 38/82] Refactoring fnft_kdvv_inverse_test_3.c --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c index 63b2695..4c508cd 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c @@ -57,7 +57,7 @@ INT main() }; #ifdef DEBUG - printf("Initial parameters: \n"); + printf("\n Initial parameters: \n"); kdvv_print_spectrum(kdvv_parameters.bound_states, kdvv_parameters.normconsts, kdvv_parameters.contspec, kdvv_parameters.XI, kdvv_parameters.M, kdvv_parameters.D, @@ -65,9 +65,9 @@ INT main() #endif - + #ifdef DEBUG - printf("-- 1. Iteration: -- \n"); + printf("\n -- 1. Iteration: -- \n"); #endif REAL err_bnd_bound_states = 3.8e-3; @@ -84,7 +84,7 @@ INT main() #ifdef DEBUG - printf("-- 2. Iteration: -- \n"); + printf("\n -- 2. Iteration: -- \n"); #endif // Check quadratic convergence @@ -103,7 +103,7 @@ INT main() #ifdef DEBUG - printf("-- 3. Iteration: -- \n"); + printf("\n -- 3. Iteration: -- \n"); #endif kdvv_parameters.D *= 2; From 1d2afecb86d14013594f0d01898c6d7cfc9028b8 Mon Sep 17 00:00:00 2001 From: FF Date: Fri, 27 Mar 2026 14:55:15 +0100 Subject: [PATCH 39/82] Changing variable names in fnft_kdvv_inverse_example.c --- examples/fnft_kdvv_inverse_example.c | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/fnft_kdvv_inverse_example.c b/examples/fnft_kdvv_inverse_example.c index 432b14c..55141be 100644 --- a/examples/fnft_kdvv_inverse_example.c +++ b/examples/fnft_kdvv_inverse_example.c @@ -26,7 +26,7 @@ #include "fnft__misc.h" -#define K_I 19 +#define K 19 INT main() @@ -45,26 +45,26 @@ INT main() CHECK_NOMEM(q, ret_code, leave_fun); // Initialize variables for a specific test case without continuous spectrum - COMPLEX * contspec_i = NULL; - contspec_i = malloc(M * sizeof(COMPLEX)); + COMPLEX * contspec = NULL; + contspec = malloc(M * sizeof(COMPLEX)); - COMPLEX bound_states_i[K_I] = { 40.0, 39.0, 38.0, 37.0, 36.0, + COMPLEX bound_states[K] = { 40.0, 39.0, 38.0, 37.0, 36.0, 30.0, 29.0, 28.0, 27.0, 26.0, 20.0, 19.0, 18.0, 17.0, 16.0, 10.0, 9.0, 8.0, 7.0}; - COMPLEX normconsts_i[K_I] = { 1e20, -1e-7, 1e5, -1e3, 1e1, - -1e0, 1e2, -1e4, 1e-6, -1e8, - 1e2, -1e4, 1e6, -1e8, 1e-10, - -1e7, 1e-6, -1e5, 1e-9}; + COMPLEX normconsts[K] = { 1e20, -1e-7, 1e5, -1e3, 1e1, + -1e0, 1e2, -1e4, 1e-6, -1e8, + 1e2, -1e4, 1e6, -1e8, 1e-10, + -1e7, 1e-6, -1e5, 1e-9}; - for (UINT i = 0; i Date: Tue, 31 Mar 2026 18:13:36 +0200 Subject: [PATCH 40/82] Add an initial version of an analytic test case --- .../fnft_kdvv_inverse_analytic_test.c | 143 ++++++++++++++++++ 1 file changed, 143 insertions(+) create mode 100644 test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c new file mode 100644 index 0000000..47b6a72 --- /dev/null +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -0,0 +1,143 @@ +/* +* This file is part of FNFT. +* +* FNFT is free software; you can redistribute it and/or +* modify it under the terms of the version 2 of the GNU General +* Public License as published by the Free Software Foundation. +* +* FNFT is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Contributors: +* Fabian Fischer (Hiwi KIT) 2026 +*/ + +#define FNFT_ENABLE_SHORT_NAMES + +#include + +#include "fnft__kdvv_inverse_testcases.h" +#include "fnft__errwarn.h" +#include "fnft__misc.h" + +/* This is a testcase with analytic reference */ + + +/* This testcase is based on the example discussed on p. 74-78 in "Solitons - an introduction" by +* P. G. Drazin, R. S. Johnson, Cambridge University Press (1989) +* +* - the eigenvalues are k_1 = 1 and k_2 = 2 +* - the eigenfunctions are: +* > psi_1 ~ sqrt(6)*exp(-t) as t->+infinity and psi_1 ~ sqrt(6)*exp(t) as t->-infinity +* ==> normconst_1 = sqrt(6) + +* the values for the normconsts are +* +* Note: The usage of x and t is inverted in the example towards the usage in this library +* - the resulting signal has to be multiplied by -1 +*/ + + +#define K 2 +#define DEBUG + +REAL analytic_signal(REAL t, REAL x){ + return 12.0*(3.0 + 4.0 * COSH(2.0*t-8.0*x) + COSH(4.0*t-64.0*x))/POW((3.0*COSH(t-28.0*x) + COSH(3.0*t-36.0*x)), 2); +} + + +INT main() +{ + INT ret_code = SUCCESS; + + REAL const x = 0.1; + + UINT D = 256; + UINT M = 10; + REAL XI[2] = {-2.0, 2.0}; + REAL T[2] = {-10.0, 10.0}; + + COMPLEX * contspec = NULL; + contspec = malloc(10*sizeof(COMPLEX)); + CHECK_NOMEM(contspec, ret_code, leave_fun); + + COMPLEX * q = NULL; + q = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(q, ret_code, leave_fun); + + COMPLEX bound_states[K] = { I*2.0, I*1.0 }; + // COMPLEX normconsts[K] = { 2.0*SQRT(3.0)*CEXP(32.0*x), -SQRT(6.0)*CEXP(4.0*x) }; + COMPLEX normconsts[K] = { CEXP(64.0*x), -CEXP(8.0*x) }; + // COMPLEX normconsts[K] = { 1, -1}; + + ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states, normconsts, D, q, T, NULL); + CHECK_RETCODE(ret_code, leave_fun); + + COMPLEX * t_grid = NULL; + t_grid = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(t_grid, ret_code, leave_fun); + + const COMPLEX eps_t = (T[1] - T[0])/(D - 1); + + for (UINT n=0; n max_quadratic_error) {max_quadratic_error = error;} + } + + #ifdef DEBUG + printf("resulting max quadratic error: %e \n", max_quadratic_error); + printf("resulting quadratic error sum: %e \n", quadratic_error_sum); + #endif + + UINT is_quadratic_error_sum_in_tolerance = quadratic_error_sum < 1e-27; + UINT is_max_quadratic_error_in_tolerance = max_quadratic_error < 1e-28; + + + if (is_max_quadratic_error_in_tolerance && + is_quadratic_error_sum_in_tolerance) + { + ret_code = SUCCESS; + } + else { + ret_code = FNFT_EC_TEST_FAILED; + } + +leave_fun: + free(contspec); + free(q); + free(analytic_q_a); + + if (ret_code != SUCCESS) + return EXIT_FAILURE; + else + return EXIT_SUCCESS; +} From 7c72bf64f1c60ded1cdd7088ff7c178b38f06bcc Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 10:16:31 +0200 Subject: [PATCH 41/82] Minor refactoring in analytic test case --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index 47b6a72..aee9409 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -40,6 +40,8 @@ * * Note: The usage of x and t is inverted in the example towards the usage in this library * - the resulting signal has to be multiplied by -1 +why COMPLEX normconsts[K] = { CEXP(64.0*x), -CEXP(8.0*x) } instead of +COMPLEX normconsts[K] = { 2.0*SQRT(3.0)*CEXP(32.0*x), -SQRT(6.0)*CEXP(4.0*x) };???? */ @@ -71,10 +73,9 @@ INT main() CHECK_NOMEM(q, ret_code, leave_fun); COMPLEX bound_states[K] = { I*2.0, I*1.0 }; - // COMPLEX normconsts[K] = { 2.0*SQRT(3.0)*CEXP(32.0*x), -SQRT(6.0)*CEXP(4.0*x) }; - COMPLEX normconsts[K] = { CEXP(64.0*x), -CEXP(8.0*x) }; - // COMPLEX normconsts[K] = { 1, -1}; - + COMPLEX normconsts[K] = { CEXP(8.0*I*CPOW(bound_states[0], 3)*x), + -CEXP(8.0*I*CPOW(bound_states[1], 3)*x) }; + ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states, normconsts, D, q, T, NULL); CHECK_RETCODE(ret_code, leave_fun); From fbb45e5d53dc8f3fffb3b3cd0534404fc8439599 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 11:28:20 +0200 Subject: [PATCH 42/82] Add comments and plots to mex_fnft_kdvv_inverse_example_1.m --- examples/mex_fnft_kdvv_inverse_example_1.m | 51 +++++++++++++++++++--- 1 file changed, 44 insertions(+), 7 deletions(-) diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index 7d8f821..a269b04 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -14,19 +14,56 @@ % % Contributors: % Sander Wahls (KIT) 2026. +% Fabian Fischer (Hiwi KIT) 2026. + +% This example should illustrate how to use the inverse kdvv transform clear all; close all; +% defining the continuous spectrum: +% out of function, just for seek of completeness (state 04/2026) contspec = []; -XI = [0 1]; -bound_states = 1i*sqrt([1, 2, 3, 4] /2); -normconsts = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); +XI = [1e-6 1]; + +% desired height of solitions: +desired_solitions = [1, 2, 3, 4]; + +% resulting bound states out of desired heights of solitions: +bound_states = 1i*sqrt(desired_solitions ./2); + +% desired norming constants +% defines how much the solitions are shifted towards each other +% signs have to be alternating regards to the order of the bound state +% sign of the normconst for the biggest eigenvalue has to be positive +norming_constants = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); + +% Number of samples of the output of the inverse kdvv D = 1001; + +% Area, for which the output has to be computed +% Can be asymmetric T = [-10 10]; -q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T) +% calls function of c-library FNFT +q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, norming_constants, D, T); + +% --- Plot the results --- t = linspace(T(1), T(2), D); -plot(t, q) -xlabel('t') -ylabel('q(t)') + +% the output of the inverse kdvv is defined as a complex number (state 04/2026), +% however only the imaginary part should be zero and can be neglected +q = double(real(q(:)')); + +figure; +plot(t, q); +title('output of inverse kdvv: time-domain'); +xlabel('t'); +ylabel('q(t)'); +legend('q(t)'); + +figure; +stem(imag(bound_states),real(norming_constants), 'x'); +title('Bound states and norming constants'); +xlabel('bound states'); +ylabel('norming constants'); \ No newline at end of file From 403da4ce00ce8017ca796d847a97b6ab58af9d34 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 11:29:26 +0200 Subject: [PATCH 43/82] Add mex_fnft_kdvv_inverse_example_2.m as illustration for the relation between inverse and forward kdvv --- examples/mex_fnft_kdvv_inverse_example_2.m | 85 ++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 examples/mex_fnft_kdvv_inverse_example_2.m diff --git a/examples/mex_fnft_kdvv_inverse_example_2.m b/examples/mex_fnft_kdvv_inverse_example_2.m new file mode 100644 index 0000000..57b14ec --- /dev/null +++ b/examples/mex_fnft_kdvv_inverse_example_2.m @@ -0,0 +1,85 @@ +% This file is part of FNFT. +% +% FNFT is free software; you can redistribute it and/or +% modify it under the terms of the version 2 of the GNU General +% Public License as published by the Free Software Foundation. +% +% FNFT is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% +% Contributors: +% Sander Wahls (KIT) 2026. +% Fabian Fischer (Hiwi KIT) 2026. + +% This example should illustrate, how to use the inverse kdvv in combination +% with the forward kdvv transform. For basic usage of the inverse kdvv, please +% please have a look at mex_fnft_kdvv_inverse_example_1.m (state 04/2026). + +clear all; +close all; + +contspec_initial = []; +XI = [1e-6 10]; + +bound_states_initial = 1i*sqrt([1, 2, 3, 4] ./2); +norming_constants_initial = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); + +D = 1001; +T = [-10 10]; + +% inverse kdvv transform +q = mex_fnft_kdvv_inverse(contspec_initial, XI, bound_states_initial, norming_constants_initial, D, T); + +q = double(real(q(:)')); + +% compute the nonlinear Fourier transform of the output of the inverse kdvv +[contspec_computed, bound_states_computed, norming_constants_computed] = mex_fnft_kdvv(q, T, XI); + + +% --- Plot the results --- +t = linspace(T(1), T(2), D); + +figure; +plot(t, q); +title('output of inverse kdvv: time-domain'); +xlabel('t'); +ylabel('q(t)'); +legend('q(t)'); + +% plotting the initial and the computed discrete spectrum +% the difference between the initial and the computed spectrum should be very small + +figure; +stem(imag(bound_states_initial),real(norming_constants_initial), 'x', 'color', 'r'); +hold on; +stem(imag(bound_states_computed),real(norming_constants_computed), 'x', 'color', 'g'); +hold off; +title('initial and computed bound states and norming constants'); +xlabel('bound states'); +ylabel('norming constants'); +legend('initial', 'computed'); + +% plotting of continuous spectrum +% values are small because initial continuous spectrum was assumed as zero +% and the computed continuous spectrum is only the result of numerical errors + +ep_xi = (XI(2) - XI(1)) / (D - 1); +xi = XI(1):ep_xi:XI(2); + +figure; +hmag=subplot(2,1,1); +semilogy(xi, abs(contspec_computed)); +title('computed continuous spectrum'); +xlabel('\xi'); +ylabel('|r(\xi)|'); +hang=subplot(2,1,2); +plot(xi, angle(contspec_computed)); +xlabel('\xi'); +ylabel('\angle r(\xi)'); +linkaxes([hmag,hang],'x'); + From cfcac8586e6baf61ed729429e0ed3961e01ec71c Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 12:14:55 +0200 Subject: [PATCH 44/82] Add mex_fnft_kdvv_inverse_example_3.m as a fast way to try different eigenvalues and normconsts --- examples/mex_fnft_kdvv_inverse_example_3.m | 79 ++++++++++++++++++++++ 1 file changed, 79 insertions(+) create mode 100644 examples/mex_fnft_kdvv_inverse_example_3.m diff --git a/examples/mex_fnft_kdvv_inverse_example_3.m b/examples/mex_fnft_kdvv_inverse_example_3.m new file mode 100644 index 0000000..9a2497c --- /dev/null +++ b/examples/mex_fnft_kdvv_inverse_example_3.m @@ -0,0 +1,79 @@ +% This file is part of FNFT. +% +% FNFT is free software; you can redistribute it and/or +% modify it under the terms of the version 2 of the GNU General +% Public License as published by the Free Software Foundation. +% +% FNFT is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% +% Contributors: +% Fabian Fischer (Hiwi KIT) 2026. + +% This example shows a fast way to try different values (bound states, norming constants) +% for the inverse kdvv transform, because this script sorts and considers the alternating +% signs of the norming constants automatically. + +clear all; +close all; + +% === modify here ==================================================================================== +% a solition in the output of inverse kdvv is defined by the elements of both arrays at the same index +desired_solitions = [9, 2, 4, 3, 7]; +desired_shifting = [10000, 0.1, 1, 0.00001, 2]; +% ==================================================================================================== + + +% Check if both arrays contains the same number of elements +if length(desired_shifting) ~= length(desired_solitions) + error('Number of desired solitions is not equal to number of values given for shifting!'); +end + +% --- determining sorted bound states and norming constants out of desired values --- + +number_values = length(desired_solitions); + +tmp_bound_states = 1i*sqrt(desired_solitions ./2); + +% sorting +[bound_states, indices] = sort(tmp_bound_states, "ascend"); +tmp_norming_constants = desired_shifting(indices); + +% alternating sequence 1 and -1 as elements and with 1 as last element +tmp_seq = (-1).^( 0:(number_values-1) ); +tmp_seq = flip(tmp_seq); + +norming_constants = complex(tmp_seq.*tmp_norming_constants); + +% inverse kdvv +contspec = []; +XI = [1e-6 1]; +D = 1001; +T = [-20 20]; +q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, norming_constants, D, T); + + +% --- Plot the results --- +t = linspace(T(1), T(2), D); + +% the output of the inverse kdvv is defined as a complex number (state 04/2026), +% however only the imaginary part should be zero and can be neglected +q = double(real(q(:)')); + +figure; +plot(t, q); +title('output of inverse kdvv: time-domain'); +xlabel('t'); +ylabel('q(t)'); +legend('q(t)'); + +figure; +stem(imag(bound_states),real(norming_constants), 'x'); +title('Bound states and norming constants'); +xlabel('bound states'); +ylabel('norming constants'); \ No newline at end of file From e574198d540fc830ba8554d2be409543eac9e6cc Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 12:18:02 +0200 Subject: [PATCH 45/82] Add mex_fnft_kdvv_inverse_example_4.m as a very simple example --- examples/mex_fnft_kdvv_inverse_example_4.m | 37 ++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 examples/mex_fnft_kdvv_inverse_example_4.m diff --git a/examples/mex_fnft_kdvv_inverse_example_4.m b/examples/mex_fnft_kdvv_inverse_example_4.m new file mode 100644 index 0000000..c036dcb --- /dev/null +++ b/examples/mex_fnft_kdvv_inverse_example_4.m @@ -0,0 +1,37 @@ +%% +% This file is part of FNFT. +% +% FNFT is free software; you can redistribute it and/or +% modify it under the terms of the version 2 of the GNU General +% Public License as published by the Free Software Foundation. +% +% FNFT is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% +% Contributors: +% Sander Wahls (KIT) 2026. + +% This is a simple example of how to use the inverse kdvv. + +close all; + +contspec = []; +XI = [0 1]; + +bound_states = 1i*sqrt( [1, 2, 3, 4, 5] /2); + +normconsts = complex([1, -1, 1, -1, 1].*[10, 0.1, 1, 0.00001, 10000000]); + +D = 1001; +T = [-10 10]; +q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T); + +t = linspace(T(1), T(2), D); +plot(t, q) +xlabel('t') +ylabel('q(t)') From d737ffc893ce533174415c2f450cd0d611c36844 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 12:19:55 +0200 Subject: [PATCH 46/82] Minor refactoring in mex_fnft_kdvv_inverse_example_1.m --- examples/mex_fnft_kdvv_inverse_example_1.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index a269b04..ca19ab5 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -21,11 +21,6 @@ clear all; close all; -% defining the continuous spectrum: -% out of function, just for seek of completeness (state 04/2026) -contspec = []; -XI = [1e-6 1]; - % desired height of solitions: desired_solitions = [1, 2, 3, 4]; @@ -45,6 +40,11 @@ % Can be asymmetric T = [-10 10]; +% defining the continuous spectrum: +% out of function, just for seek of completeness (state 04/2026) +contspec = []; +XI = [1e-6 1]; + % calls function of c-library FNFT q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, norming_constants, D, T); From 3853ca127cb14cebc6681464ca67db15757b8447 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 14:11:08 +0200 Subject: [PATCH 47/82] Refactoring examples for inverse kdvv because of consistency in the order of provided bound states to add --- examples/mex_fnft_kdvv_inverse_example_1.m | 7 +++++-- examples/mex_fnft_kdvv_inverse_example_2.m | 4 ++-- examples/mex_fnft_kdvv_inverse_example_3.m | 5 ++--- examples/mex_fnft_kdvv_inverse_example_4.m | 2 +- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index ca19ab5..c8620fd 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -22,16 +22,17 @@ close all; % desired height of solitions: -desired_solitions = [1, 2, 3, 4]; +desired_solitions = [4, 3, 2, 1]; % resulting bound states out of desired heights of solitions: +% - bound states have to be true and positive imaginary bound_states = 1i*sqrt(desired_solitions ./2); % desired norming constants % defines how much the solitions are shifted towards each other % signs have to be alternating regards to the order of the bound state % sign of the normconst for the biggest eigenvalue has to be positive -norming_constants = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); +norming_constants = complex([1, -1, 1, -1].*[10, 0.1, 1, 0.00001]); % Number of samples of the output of the inverse kdvv D = 1001; @@ -47,6 +48,8 @@ % calls function of c-library FNFT q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, norming_constants, D, T); +% the result consists of solitions with desired height (if the solitions are far +% away enough to each other) % --- Plot the results --- t = linspace(T(1), T(2), D); diff --git a/examples/mex_fnft_kdvv_inverse_example_2.m b/examples/mex_fnft_kdvv_inverse_example_2.m index 57b14ec..bacc36a 100644 --- a/examples/mex_fnft_kdvv_inverse_example_2.m +++ b/examples/mex_fnft_kdvv_inverse_example_2.m @@ -26,8 +26,8 @@ contspec_initial = []; XI = [1e-6 10]; -bound_states_initial = 1i*sqrt([1, 2, 3, 4] ./2); -norming_constants_initial = complex([-1, 1, -1, 1].*[10, 0.1, 1, 0.00001]); +bound_states_initial = 1i*sqrt([4, 3, 2, 1] ./2); +norming_constants_initial = complex([1, -1, 1, -1].*[10, 0.1, 1, 0.00001]); D = 1001; T = [-10 10]; diff --git a/examples/mex_fnft_kdvv_inverse_example_3.m b/examples/mex_fnft_kdvv_inverse_example_3.m index 9a2497c..a71f2eb 100644 --- a/examples/mex_fnft_kdvv_inverse_example_3.m +++ b/examples/mex_fnft_kdvv_inverse_example_3.m @@ -41,12 +41,11 @@ tmp_bound_states = 1i*sqrt(desired_solitions ./2); % sorting -[bound_states, indices] = sort(tmp_bound_states, "ascend"); +[bound_states, indices] = sort(tmp_bound_states, "descend"); tmp_norming_constants = desired_shifting(indices); -% alternating sequence 1 and -1 as elements and with 1 as last element +% alternating sequence 1 and -1 as elements and with 1 as first element tmp_seq = (-1).^( 0:(number_values-1) ); -tmp_seq = flip(tmp_seq); norming_constants = complex(tmp_seq.*tmp_norming_constants); diff --git a/examples/mex_fnft_kdvv_inverse_example_4.m b/examples/mex_fnft_kdvv_inverse_example_4.m index c036dcb..30f74d5 100644 --- a/examples/mex_fnft_kdvv_inverse_example_4.m +++ b/examples/mex_fnft_kdvv_inverse_example_4.m @@ -23,7 +23,7 @@ contspec = []; XI = [0 1]; -bound_states = 1i*sqrt( [1, 2, 3, 4, 5] /2); +bound_states = 1i*sqrt( [5, 4, 3, 2, 1] /2); normconsts = complex([1, -1, 1, -1, 1].*[10, 0.1, 1, 0.00001, 10000000]); From cb5a2edc648bfd44c8e13009a666bff4d3d6d482 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 14:12:33 +0200 Subject: [PATCH 48/82] Improving comment in mex_fnft_kdvv_inverse_example_1.m --- examples/mex_fnft_kdvv_inverse_example_1.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index c8620fd..7889e2b 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -30,8 +30,8 @@ % desired norming constants % defines how much the solitions are shifted towards each other -% signs have to be alternating regards to the order of the bound state -% sign of the normconst for the biggest eigenvalue has to be positive +% - signs have to alternate regards to the order of the bound states +% - sign of the normconst for the biggest eigenvalue has to be positive norming_constants = complex([1, -1, 1, -1].*[10, 0.1, 1, 0.00001]); % Number of samples of the output of the inverse kdvv From bb2d5abad417f2ade75a21d48c0c0a07beb78597 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 14:13:27 +0200 Subject: [PATCH 49/82] Add fnft_kdvv_inverse_example_1.c as a simple example for the usage of inverse kdvv in C --- examples/fnft_kdvv_inverse_example_1.c | 73 ++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 examples/fnft_kdvv_inverse_example_1.c diff --git a/examples/fnft_kdvv_inverse_example_1.c b/examples/fnft_kdvv_inverse_example_1.c new file mode 100644 index 0000000..5f303c8 --- /dev/null +++ b/examples/fnft_kdvv_inverse_example_1.c @@ -0,0 +1,73 @@ +/* +* This file is part of FNFT. +* +* FNFT is free software; you can redistribute it and/or +* modify it under the terms of the version 2 of the GNU General +* Public License as published by the Free Software Foundation. +* +* FNFT is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* Contributors: +* Fabian Fischer (Hiwi KIT) 2026 +*/ + +// This is a simple example for the usage of the inverse kdvv transform in C + +#define FNFT_ENABLE_SHORT_NAMES + +#include + +#include "fnft_kdvv_inverse.h" +#include "fnft__errwarn.h" +#include "fnft__misc.h" + +#define K 5 // number of bound states to add + +INT main() +{ + INT ret_code = SUCCESS; + + // five desired bound states and norming constants + // - bound states needs to be true imaginary positive + // - every norming constants belong to the bound state at the same index + // - the signs of the norming constants have to alternate regarding the order + // of bound states. The norming constant of the biggest bound state has to be positive. + COMPLEX bound_states[K] = { I*5.0, I*4.0, I*3.0, I*2.0, I*1.0 }; + COMPLEX norming_constants[K] = { 1e5, -1e-2, 1e0, -1e3, 1e1 }; + + // General parameters + UINT D = 256; // Number of samples of the computed output + REAL T[2] = {-10.0, 10.0}; // area for which the output should be computed + + // allocation of memory for the computed output + COMPLEX * q = NULL; + q = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(q, ret_code, leave_fun); + + // Define continuous spectrum + // continuous spectrum is out of function, but needs to be defined/allocated for using + // the inverse kdvv (state 04/2026) + UINT M = 10; + REAL XI[2] = {-2.0, 2.0}; + + COMPLEX * contspec = NULL; + contspec = malloc(M * sizeof(COMPLEX)); + CHECK_NOMEM(contspec, ret_code, leave_fun); + + // call of the inverse kdvv + ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states, norming_constants, D, q, T, NULL); + CHECK_RETCODE(ret_code, leave_fun); + + // prints the results in the console + misc_print_buf(D, q, "output"); + +leave_fun: + free(q); + free(contspec); +} From d4391a3448ce1790acc7bbd2dc550c45bddbb02d Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 14:17:30 +0200 Subject: [PATCH 50/82] Change fnft_kdvv_inverse_example.c to fnft_kdvv_inverse_example_2.c. This is an example for the usage of the inverse_kdvv with many eigenvalues. Adding explanatory comments. --- ...xample.c => fnft_kdvv_inverse_example_2.c} | 60 +++++++++++-------- 1 file changed, 36 insertions(+), 24 deletions(-) rename examples/{fnft_kdvv_inverse_example.c => fnft_kdvv_inverse_example_2.c} (56%) diff --git a/examples/fnft_kdvv_inverse_example.c b/examples/fnft_kdvv_inverse_example_2.c similarity index 56% rename from examples/fnft_kdvv_inverse_example.c rename to examples/fnft_kdvv_inverse_example_2.c index 55141be..b0179c4 100644 --- a/examples/fnft_kdvv_inverse_example.c +++ b/examples/fnft_kdvv_inverse_example_2.c @@ -17,6 +17,9 @@ * Fabian Fischer (Hiwi KIT) 2026 */ +// This is an example for the usage of the inverse kdvv transform with many +// bound states + #define FNFT_ENABLE_SHORT_NAMES #include @@ -25,51 +28,60 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" - -#define K 19 - +#define K 19 // Number of bound states to add INT main() { INT ret_code = SUCCESS; + COMPLEX desired_solitions_height[K] = { 40.0, 39.0, 38.0, 37.0, 36.0, + 30.0, 29.0, 28.0, 27.0, 26.0, + 20.0, 19.0, 18.0, 17.0, 16.0, + 10.0, 9.0, 8.0, 7.0}; + + COMPLEX normconsts[K] = { 1e20, -1e-7, 1e5, -1e3, 1e1, + -1e0, 1e2, -1e4, 1e-6, -1e8, + 1e2, -1e4, 1e6, -1e8, 1e-10, + -1e7, 1e-6, -1e5, 1e-9}; + + // resulting bound states out of desired solitions height + COMPLEX * bound_states = NULL; + bound_states = malloc(K * sizeof(COMPLEX)); + CHECK_NOMEM(bound_states, ret_code, leave_fun); + + for (UINT i = 0; i Date: Tue, 7 Apr 2026 14:23:32 +0200 Subject: [PATCH 51/82] Adapting bound states and windows for correct scaling --- examples/mex_fnft_kdvv_inverse_example_1.m | 2 +- examples/mex_fnft_kdvv_inverse_example_2.m | 4 ++-- examples/mex_fnft_kdvv_inverse_example_3.m | 2 +- examples/mex_fnft_kdvv_inverse_example_4.m | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index 7889e2b..da1120a 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -32,7 +32,7 @@ % defines how much the solitions are shifted towards each other % - signs have to alternate regards to the order of the bound states % - sign of the normconst for the biggest eigenvalue has to be positive -norming_constants = complex([1, -1, 1, -1].*[10, 0.1, 1, 0.00001]); +norming_constants = complex([1, -1, 1, -1].*[0.00001, 0.1, 1, 10]); % Number of samples of the output of the inverse kdvv D = 1001; diff --git a/examples/mex_fnft_kdvv_inverse_example_2.m b/examples/mex_fnft_kdvv_inverse_example_2.m index bacc36a..fcc8811 100644 --- a/examples/mex_fnft_kdvv_inverse_example_2.m +++ b/examples/mex_fnft_kdvv_inverse_example_2.m @@ -27,10 +27,10 @@ XI = [1e-6 10]; bound_states_initial = 1i*sqrt([4, 3, 2, 1] ./2); -norming_constants_initial = complex([1, -1, 1, -1].*[10, 0.1, 1, 0.00001]); +norming_constants_initial = complex([1, -1, 1, -1].*[10000, 0.1, 1, 0.001]); D = 1001; -T = [-10 10]; +T = [-15 10]; % inverse kdvv transform q = mex_fnft_kdvv_inverse(contspec_initial, XI, bound_states_initial, norming_constants_initial, D, T); diff --git a/examples/mex_fnft_kdvv_inverse_example_3.m b/examples/mex_fnft_kdvv_inverse_example_3.m index a71f2eb..9334d3a 100644 --- a/examples/mex_fnft_kdvv_inverse_example_3.m +++ b/examples/mex_fnft_kdvv_inverse_example_3.m @@ -53,7 +53,7 @@ contspec = []; XI = [1e-6 1]; D = 1001; -T = [-20 20]; +T = [-12 8]; q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, norming_constants, D, T); diff --git a/examples/mex_fnft_kdvv_inverse_example_4.m b/examples/mex_fnft_kdvv_inverse_example_4.m index 30f74d5..0ca29a4 100644 --- a/examples/mex_fnft_kdvv_inverse_example_4.m +++ b/examples/mex_fnft_kdvv_inverse_example_4.m @@ -25,7 +25,7 @@ bound_states = 1i*sqrt( [5, 4, 3, 2, 1] /2); -normconsts = complex([1, -1, 1, -1, 1].*[10, 0.1, 1, 0.00001, 10000000]); +normconsts = complex([1, -1, 1, -1, 1].*[10000000, 0.1, 1, 0.001, 10]); D = 1001; T = [-10 10]; From f5013f145b0671526a656a938c927300ddd05074 Mon Sep 17 00:00:00 2001 From: FF Date: Tue, 7 Apr 2026 15:19:29 +0200 Subject: [PATCH 52/82] Adding checks for invalid inputs in fnft_kdvv_inverse.c --- src/fnft_kdvv_inverse.c | 54 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 48 insertions(+), 6 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index f7fe04e..fe4aefe 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -473,7 +473,7 @@ static INT add_one_soliton( } -INT add_two_solitons( +static INT add_two_solitons( COMPLEX const * const bound_states, UINT const N_bound_states, COMPLEX * const theta_E1, @@ -619,6 +619,30 @@ INT fnft_kdvv_inverse( INT ret_code = SUCCESS; + + // checks for valid input --------------------------------------------------------------------------------- + + for (UINT i=0; i T[0] + if (T[1] <= T[0]) {E_INVALID_ARGUMENT(T);} + CHECK_RETCODE(ret_code, leave_fun); + + // ---------------------------------------------------------------------------------------------------------- + + + // Initialize q for (UINT n=0; n 0){ + ret_code = E_INVALID_ARGUMENT(normconsts_or_residues) + } + } + CHECK_RETCODE(ret_code, leave_fun); + + // ---------------------------------------------------------------------------------------------------------- + + // shift sign of norm_consts if K even if (K%2 == 0){ for (UINT i=0; i Date: Fri, 10 Apr 2026 10:19:10 +0200 Subject: [PATCH 53/82] Fix crash of matlab after invalid arguments because of wrong deallocating of memory: declaration of all pointers (which are later used with malloc) needs to be before first usage of goto leave_fun (where all free() statements are). --- src/fnft_kdvv_inverse.c | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index fe4aefe..26e8f16 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -164,6 +164,9 @@ static INT one_solition_update_Jost( COMPLEX const * const M_min1_11) { INT ret_code = SUCCESS; + + COMPLEX * C_E = NULL; + COMPLEX * prefactor_C_E = NULL; // Check, if Jost has to be updated. Only when there are bound states left, that are not added yet UINT const is_Jost_to_update = N_jZ != 0; @@ -178,7 +181,6 @@ static INT one_solition_update_Jost( // Create C_E Matrix (for Updating Jost Solution) UINT n_C_E = 2*2*D*N_jZ; - COMPLEX * C_E = NULL; C_E = malloc(n_C_E*sizeof(COMPLEX)); CHECK_NOMEM(C_E, ret_code, leave_fun); @@ -199,7 +201,6 @@ static INT one_solition_update_Jost( } // Compute a new prefactor for C_E - COMPLEX * prefactor_C_E = NULL; prefactor_C_E = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(prefactor_C_E, ret_code, leave_fun); @@ -266,6 +267,7 @@ static INT two_solitions_update_Jost( COMPLEX const * const s) { INT ret_code = SUCCESS; + COMPLEX * C_E = NULL; // Check, if Jost has to be updated. Only when there are bound states left, that are not added yet UINT const is_Jost_to_update = N_jZ != 0; @@ -286,7 +288,6 @@ static INT two_solitions_update_Jost( // Create C_E Matrix (for Updating Jost Solution) UINT n_C_E = 2*2*D*N_jZ; - COMPLEX * C_E = NULL; C_E = malloc(n_C_E*sizeof(COMPLEX)); CHECK_NOMEM(C_E, ret_code, leave_fun); @@ -416,6 +417,7 @@ static INT add_one_soliton( COMPLEX * M_min1_11 = NULL; COMPLEX * w_inv = NULL; COMPLEX * prefactor = NULL; + M_min1_11 = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(M_min1_11, ret_code, leave_fun); w_inv = malloc(D * sizeof(COMPLEX)); @@ -487,6 +489,11 @@ static INT add_two_solitons( COMPLEX * w_inv = NULL; COMPLEX * prefactor = NULL; COMPLEX * dq = NULL; + COMPLEX * pm0_jZ = NULL; + COMPLEX * s = NULL; + COMPLEX * temp_t1 = NULL; + COMPLEX * temp_t2 = NULL; + w_inv = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(w_inv, ret_code, leave_fun); prefactor = malloc(D * sizeof(COMPLEX)); @@ -506,7 +513,6 @@ static INT add_two_solitons( // Initialize pm0_jZ UINT N_pm0_jZ = 2 * N_bound_states_left + 1; - COMPLEX * pm0_jZ = NULL; pm0_jZ = malloc((2*N_bound_states_left+1) * sizeof(COMPLEX)); CHECK_NOMEM(pm0_jZ, ret_code, leave_fun); @@ -518,14 +524,10 @@ static INT add_two_solitons( pm0_jZ[2*N_bound_states_left] = 0; // dimension of s is dependend on number of bound states left - COMPLEX * s = NULL; s = malloc(D * N_pm0_jZ * sizeof(COMPLEX)); CHECK_NOMEM(s, ret_code, leave_fun); // Intermediate variables - COMPLEX * temp_t1 = NULL; - COMPLEX * temp_t2 = NULL; - temp_t1 = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(temp_t1, ret_code, leave_fun); temp_t2 = malloc(D * sizeof(COMPLEX)); @@ -618,8 +620,13 @@ INT fnft_kdvv_inverse( */ INT ret_code = SUCCESS; - - + + COMPLEX * t_grid = NULL; + COMPLEX * bound_states_sorted = NULL; + COMPLEX * normconsts_sorted = NULL; + COMPLEX * theta_E1 = NULL; + COMPLEX * theta_E2 = NULL; + // checks for valid input --------------------------------------------------------------------------------- for (UINT i=0; i Date: Fri, 10 Apr 2026 11:51:49 +0200 Subject: [PATCH 54/82] Changing in tests and examples fo fnft_kdvv_inversey so that pointers are declared at the beginning. Idea analog to last commit. --- examples/fnft_kdvv_inverse_example_1.c | 5 +++-- examples/fnft_kdvv_inverse_example_2.c | 7 ++++--- .../fnft_kdvv_inverse_analytic_test.c | 12 +++++++----- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c | 6 ++++-- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c | 6 ++++-- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c | 6 +++++- 6 files changed, 27 insertions(+), 15 deletions(-) diff --git a/examples/fnft_kdvv_inverse_example_1.c b/examples/fnft_kdvv_inverse_example_1.c index 5f303c8..4285e12 100644 --- a/examples/fnft_kdvv_inverse_example_1.c +++ b/examples/fnft_kdvv_inverse_example_1.c @@ -32,6 +32,9 @@ INT main() { INT ret_code = SUCCESS; + + COMPLEX * q = NULL; + COMPLEX * contspec = NULL; // five desired bound states and norming constants // - bound states needs to be true imaginary positive @@ -46,7 +49,6 @@ INT main() REAL T[2] = {-10.0, 10.0}; // area for which the output should be computed // allocation of memory for the computed output - COMPLEX * q = NULL; q = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(q, ret_code, leave_fun); @@ -56,7 +58,6 @@ INT main() UINT M = 10; REAL XI[2] = {-2.0, 2.0}; - COMPLEX * contspec = NULL; contspec = malloc(M * sizeof(COMPLEX)); CHECK_NOMEM(contspec, ret_code, leave_fun); diff --git a/examples/fnft_kdvv_inverse_example_2.c b/examples/fnft_kdvv_inverse_example_2.c index b0179c4..04de867 100644 --- a/examples/fnft_kdvv_inverse_example_2.c +++ b/examples/fnft_kdvv_inverse_example_2.c @@ -34,6 +34,10 @@ INT main() { INT ret_code = SUCCESS; + COMPLEX * bound_states = NULL; + COMPLEX * q = NULL; + COMPLEX * contspec = NULL; + COMPLEX desired_solitions_height[K] = { 40.0, 39.0, 38.0, 37.0, 36.0, 30.0, 29.0, 28.0, 27.0, 26.0, 20.0, 19.0, 18.0, 17.0, 16.0, @@ -45,7 +49,6 @@ INT main() -1e7, 1e-6, -1e5, 1e-9}; // resulting bound states out of desired solitions height - COMPLEX * bound_states = NULL; bound_states = malloc(K * sizeof(COMPLEX)); CHECK_NOMEM(bound_states, ret_code, leave_fun); @@ -58,7 +61,6 @@ INT main() REAL T[2] = {-18.0, 12.0}; // allocating memory for the output of the inverse kdvv - COMPLEX * q = NULL; q = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(q, ret_code, leave_fun); @@ -69,7 +71,6 @@ INT main() UINT M = 10; REAL XI[2] = {-2.0, 2.0}; - COMPLEX * contspec = NULL; contspec = malloc(M * sizeof(COMPLEX)); CHECK_NOMEM(contspec, ret_code, leave_fun); diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index aee9409..382a9a3 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -57,18 +57,21 @@ INT main() { INT ret_code = SUCCESS; - REAL const x = 0.1; + COMPLEX * contspec = NULL; + COMPLEX * q = NULL; + COMPLEX * t_grid = NULL; + COMPLEX * analytic_q_a = NULL; + + REAL const x = 0.5; UINT D = 256; UINT M = 10; REAL XI[2] = {-2.0, 2.0}; REAL T[2] = {-10.0, 10.0}; - COMPLEX * contspec = NULL; contspec = malloc(10*sizeof(COMPLEX)); CHECK_NOMEM(contspec, ret_code, leave_fun); - COMPLEX * q = NULL; q = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(q, ret_code, leave_fun); @@ -79,7 +82,6 @@ INT main() ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states, normconsts, D, q, T, NULL); CHECK_RETCODE(ret_code, leave_fun); - COMPLEX * t_grid = NULL; t_grid = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(t_grid, ret_code, leave_fun); @@ -89,7 +91,6 @@ INT main() t_grid[n]= T[0] + n*eps_t; } - COMPLEX * analytic_q_a = NULL; analytic_q_a = malloc(D*sizeof(COMPLEX)); for (UINT i=0; i < D; i++) { @@ -135,6 +136,7 @@ INT main() leave_fun: free(contspec); free(q); + free(t_grid); free(analytic_q_a); if (ret_code != SUCCESS) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index 7f8478b..57410bb 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -42,9 +42,11 @@ INT main() .normconsts = normconsts_i, .M = 10, .XI = {-2.0, 2.0}, - .contspec = malloc(10*sizeof(COMPLEX)) + .contspec = NULL, }; + kdvv_parameters.contspec = malloc(kdvv_parameters.M*sizeof(COMPLEX)); + CHECK_NOMEM(kdvv_parameters.contspec, ret_code, leave_fun); REAL err_bnd_bound_states = 1.4e-3; REAL err_bnd_spurious_bound_states = 1e-2; @@ -86,7 +88,7 @@ INT main() leave_fun: - + free(kdvv_parameters.contspec); if (ret_code != SUCCESS) return EXIT_FAILURE; diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c index edbd4ad..ef22dfe 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c @@ -56,9 +56,11 @@ INT main() .normconsts = normconsts_i, .M = 10, .XI = {-2.0, 2.0}, - .contspec = malloc(10*sizeof(COMPLEX)) + .contspec = NULL, }; + kdvv_parameters.contspec = malloc(kdvv_parameters.M*sizeof(COMPLEX)); + CHECK_NOMEM(kdvv_parameters.contspec, ret_code, leave_fun); REAL err_bnd_bound_states = 2e-3; REAL err_bnd_spurious_bound_states = 0.1; @@ -100,7 +102,7 @@ INT main() leave_fun: - + free(kdvv_parameters.contspec); if (ret_code != SUCCESS) return EXIT_FAILURE; diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c index 4c508cd..64bc384 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c @@ -53,9 +53,12 @@ INT main() .normconsts = normconsts_i, .M = 10, .XI = {-2.0, 2.0}, - .contspec = malloc(10*sizeof(COMPLEX)) + .contspec = NULL, }; + kdvv_parameters.contspec = malloc(kdvv_parameters.M*sizeof(COMPLEX)); + CHECK_NOMEM(kdvv_parameters.contspec, ret_code, leave_fun); + #ifdef DEBUG printf("\n Initial parameters: \n"); kdvv_print_spectrum(kdvv_parameters.bound_states, kdvv_parameters.normconsts, @@ -120,6 +123,7 @@ INT main() leave_fun: + free(kdvv_parameters.contspec); if (ret_code != SUCCESS) From 35f85ce8064f109fb63a827611faef3250e16d29 Mon Sep 17 00:00:00 2001 From: FF Date: Thu, 16 Apr 2026 16:02:58 +0200 Subject: [PATCH 55/82] Remove sorting of the handed over bound states and norming constants in fnft_kdvv_inverse. As a conseqeunce, adapt checks vor valid inputs and add E_INVALID_ARGUMENT_MSG as error code, which can handle descriptive error messages. --- include/private/fnft__errwarn.h | 14 ++++ src/fnft_kdvv_inverse.c | 125 ++++++++++++++------------------ 2 files changed, 67 insertions(+), 72 deletions(-) diff --git a/include/private/fnft__errwarn.h b/include/private/fnft__errwarn.h index 1bf1e80..9c5643e 100644 --- a/include/private/fnft__errwarn.h +++ b/include/private/fnft__errwarn.h @@ -47,6 +47,12 @@ */ #define FNFT__E_INVALID_ARGUMENT(name) FNFT__ERRMSG(FNFT_EC_INVALID_ARGUMENT, FNFT__E_INVALID_ARGUMENT_(name)) +/** + * Macro for raising an invalid argument error with detailed description. + * @ingroup private_errwarn + */ +#define FNFT__E_INVALID_ARGUMENT_MSG(name,msg) FNFT__ERRMSG(FNFT_EC_INVALID_ARGUMENT, FNFT__E_INVALID_ARGUMENT_MSG_(name, msg)) + /** * Macro for raising an subroutine error, where ec is the error code returned * by the subroutine. The error code returned by this macro is -abs(ec). Since @@ -111,6 +117,7 @@ #define WARN(msg) FNFT__WARN(msg) #define E_NOMEM FNFT__E_NOMEM #define E_INVALID_ARGUMENT(name) FNFT__E_INVALID_ARGUMENT(name) +#define E_INVALID_ARGUMENT_MSG(name,msg) FNFT__E_INVALID_ARGUMENT_MSG(name, msg) #define E_SUBROUTINE(ec) FNFT__E_SUBROUTINE(ec) #define E_DIV_BY_ZERO FNFT__E_DIV_BY_ZERO #define E_TEST_FAILED FNFT__E_TEST_FAILED @@ -139,6 +146,13 @@ */ #define FNFT__E_INVALID_ARGUMENT_(name) "Invalid argument "#name"." +/** + * Auxiliary macro that is used to stringify the input to + * \link FNFT__E_INVALID_ARGUMENT_MSG \endlink. Do not call directly. + * @ingroup private_errwarn + */ +#define FNFT__E_INVALID_ARGUMENT_MSG_(name,msg) "Invalid argument "#name". "#msg + /** * Auxiliary macro used to stringify the input to * \link FNFT__E_NOT_YET_IMPLEMENTED \endlink. Do not call directly. diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 26e8f16..c457b5e 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -417,7 +417,7 @@ static INT add_one_soliton( COMPLEX * M_min1_11 = NULL; COMPLEX * w_inv = NULL; COMPLEX * prefactor = NULL; - + M_min1_11 = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(M_min1_11, ret_code, leave_fun); w_inv = malloc(D * sizeof(COMPLEX)); @@ -619,36 +619,51 @@ INT fnft_kdvv_inverse( */ - INT ret_code = SUCCESS; - - COMPLEX * t_grid = NULL; - COMPLEX * bound_states_sorted = NULL; - COMPLEX * normconsts_sorted = NULL; - COMPLEX * theta_E1 = NULL; - COMPLEX * theta_E2 = NULL; + + // checks for valid input --------------------------------------------------------------------------------- + + if (T[1] <= T[0]) { + return E_INVALID_ARGUMENT_MSG(T,Position of the first sample is bigger than of the last sample (i.e. T[1] <= T[0])!); + } - // checks for valid input --------------------------------------------------------------------------------- - for (UINT i=0; i 0){ + return E_INVALID_ARGUMENT_MSG(normconsts_or_residues,The signs of the norming constants does not alternate!); + } + } + + // first checking if correct order of bound states exist, than checking if the first norming constant is positive! + if (CREAL(normconsts_or_residues[0]) < 0){ + return E_INVALID_ARGUMENT_MSG(normconsts_or_residues,The norming constant of the biggest bound state is negative!); } - - // check if T[1] > T[0] - if (T[1] <= T[0]) {E_INVALID_ARGUMENT(T);} - CHECK_RETCODE(ret_code, leave_fun); // ---------------------------------------------------------------------------------------------------------- - - + + INT ret_code = SUCCESS; + + COMPLEX * t_grid = NULL; + COMPLEX * bound_states_i = NULL; + COMPLEX * normconsts_i = NULL; + COMPLEX * theta_E1 = NULL; + COMPLEX * theta_E2 = NULL; // Initialize q for (UINT n=0; n 0){ - ret_code = E_INVALID_ARGUMENT(normconsts_or_residues) - } - } - CHECK_RETCODE(ret_code, leave_fun); - - // ---------------------------------------------------------------------------------------------------------- - // shift sign of norm_consts if K even if (K%2 == 0){ for (UINT i=0; i Date: Thu, 16 Apr 2026 17:08:16 +0200 Subject: [PATCH 56/82] Improve test description in fnft_kdvv_inverse_analytic_test --- .../fnft_kdvv_inverse_analytic_test.c | 35 ++++++++++--------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index 382a9a3..a4bae37 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -25,23 +25,24 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" -/* This is a testcase with analytic reference */ - - -/* This testcase is based on the example discussed on p. 74-78 in "Solitons - an introduction" by -* P. G. Drazin, R. S. Johnson, Cambridge University Press (1989) -* -* - the eigenvalues are k_1 = 1 and k_2 = 2 -* - the eigenfunctions are: -* > psi_1 ~ sqrt(6)*exp(-t) as t->+infinity and psi_1 ~ sqrt(6)*exp(t) as t->-infinity -* ==> normconst_1 = sqrt(6) - -* the values for the normconsts are -* -* Note: The usage of x and t is inverted in the example towards the usage in this library -* - the resulting signal has to be multiplied by -1 -why COMPLEX normconsts[K] = { CEXP(64.0*x), -CEXP(8.0*x) } instead of -COMPLEX normconsts[K] = { 2.0*SQRT(3.0)*CEXP(32.0*x), -SQRT(6.0)*CEXP(4.0*x) };???? + +/* This is a testcase with analytic reference +* +* This test case is based on the example discussed on p. 74-78 in [1]. Similar to this example +* is the example discussed in section 5.2 (p. 13) in [2]. Note that we use different notation and +* the definition of the norming constants in [2]. Thus we have to do the following adaptions of the +* formulas given in [1]: +* +* - substitute x with t and t with x +* - the analytic signal needs to be multiplied by -1 +* - norming constants: +* b(2i, x) = exp(8i*(2i)^3*x) --> b(2i, 0) = 1 +* b(1i, x) = -exp(8i*(1i)^3*x) --> b(1i, 0) = -1 +* +* [1] P. G. Drazin, R. S. Johnson (1989). Solitons - an introduction. Cambridge University Press +* [2] Prins, P. J., & Wahls, S. (2021). An accurate O(N 2) floating point algorithm for the Crum +* transform of the KdV equation. Communications in Nonlinear Science and Numerical Simulation, +* 102, Article 105782. https://doi.org/10.1016/j.cnsns.2021.105782 */ From f42051341af3ad0958de059fa36a06d40c48b423 Mon Sep 17 00:00:00 2001 From: FF Date: Fri, 17 Apr 2026 14:47:33 +0200 Subject: [PATCH 57/82] Define error tolerances as constants and adapt window position in analytic test case --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index a4bae37..1c40376 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -48,6 +48,8 @@ #define K 2 #define DEBUG +#define QUADRATIC_ERROR_SUM_TOLERANCE 1e-27 +#define MAX_QUADRATIC_ERROR 1e-28 REAL analytic_signal(REAL t, REAL x){ return 12.0*(3.0 + 4.0 * COSH(2.0*t-8.0*x) + COSH(4.0*t-64.0*x))/POW((3.0*COSH(t-28.0*x) + COSH(3.0*t-36.0*x)), 2); @@ -68,7 +70,7 @@ INT main() UINT D = 256; UINT M = 10; REAL XI[2] = {-2.0, 2.0}; - REAL T[2] = {-10.0, 10.0}; + REAL T[2] = {-8.0, 12.0}; contspec = malloc(10*sizeof(COMPLEX)); CHECK_NOMEM(contspec, ret_code, leave_fun); @@ -121,8 +123,8 @@ INT main() printf("resulting quadratic error sum: %e \n", quadratic_error_sum); #endif - UINT is_quadratic_error_sum_in_tolerance = quadratic_error_sum < 1e-27; - UINT is_max_quadratic_error_in_tolerance = max_quadratic_error < 1e-28; + UINT is_quadratic_error_sum_in_tolerance = quadratic_error_sum < QUADRATIC_ERROR_SUM_TOLERANCE; + UINT is_max_quadratic_error_in_tolerance = max_quadratic_error < MAX_QUADRATIC_ERROR; if (is_max_quadratic_error_in_tolerance && From 4db0a8b18756dffe124329be5420011901834161 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Mon, 20 Apr 2026 19:41:18 +0200 Subject: [PATCH 58/82] Fix typos in comments --- include/fnft_kdvv.h | 2 +- .../fnft_kdvv_inverse_analytic_test.c | 34 +++++++++---------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/include/fnft_kdvv.h b/include/fnft_kdvv.h index 542ac97..cedb0f2 100644 --- a/include/fnft_kdvv.h +++ b/include/fnft_kdvv.h @@ -270,7 +270,7 @@ fnft_kdvv_opts_t fnft_kdvv_default_opts(); * - fnft_kdv_discretization_TES4(_VANILLA) * * The accuray of the computed quantities for a given signal depends primarily on the number of samples \f$ D\f$ and the numerical method. When the exact spectrum is - * is know, the accuracy can be quantified by defining a suitable error. The error usually decreases with increasing \f$ D\f$ assuming everthing else remains the same. + * known, the accuracy can be quantified by defining a suitable error. The error usually decreases with increasing \f$ D\f$ assuming everthing else remains the same. * The rate at which the error decreases with increase in \f$ D\f$ is known as the order of the method. The orders of the various discretizations can be found at \link fnft_kdv_discretization_t \endlink. * The orders of the discretizations which use exponential splitting schemes should be the same as their base methods but can deviate when accuracy of the splitting scheme is low. * diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index 1c40376..3383d14 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -27,23 +27,23 @@ /* This is a testcase with analytic reference -* -* This test case is based on the example discussed on p. 74-78 in [1]. Similar to this example -* is the example discussed in section 5.2 (p. 13) in [2]. Note that we use different notation and -* the definition of the norming constants in [2]. Thus we have to do the following adaptions of the -* formulas given in [1]: -* -* - substitute x with t and t with x -* - the analytic signal needs to be multiplied by -1 -* - norming constants: -* b(2i, x) = exp(8i*(2i)^3*x) --> b(2i, 0) = 1 -* b(1i, x) = -exp(8i*(1i)^3*x) --> b(1i, 0) = -1 -* -* [1] P. G. Drazin, R. S. Johnson (1989). Solitons - an introduction. Cambridge University Press -* [2] Prins, P. J., & Wahls, S. (2021). An accurate O(N 2) floating point algorithm for the Crum -* transform of the KdV equation. Communications in Nonlinear Science and Numerical Simulation, -* 102, Article 105782. https://doi.org/10.1016/j.cnsns.2021.105782 -*/ + * + * This test case is based on the example discussed on p. 74-78 in [1]. Similar to this example + * is the example discussed in section 5.2 (p. 13) in [2]. Note that we use different notation and + * the definition of the norming constants in [2]. Thus we have to do the following adaptions of the + * formulas given in [1]: + * + * - substitute x with t and t with x + * - the analytic signal needs to be multiplied by -1 + * - norming constants: + * b(2i, x) = exp(8i*(2i)^3*x) --> b(2i, 0) = 1 + * b(1i, x) = -exp(8i*(1i)^3*x) --> b(1i, 0) = -1 + * + * [1] P. G. Drazin, R. S. Johnson (1989). Solitons - an introduction. Cambridge University Press + * [2] Prins, P. J., & Wahls, S. (2021). An accurate O(N^2) floating point algorithm for the Crum + * transform of the KdV equation. Communications in Nonlinear Science and Numerical Simulation, + * 102, Article 105782. https://doi.org/10.1016/j.cnsns.2021.105782 + */ #define K 2 From a8ba21f3d18d69211bef4a484945a5d4d8d06eec Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Mon, 20 Apr 2026 19:44:11 +0200 Subject: [PATCH 59/82] Clarification in comments of examples to inverse nft --- examples/fnft_kdvv_inverse_example_1.c | 3 ++- examples/mex_fnft_kdvv_inverse_example_1.m | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/fnft_kdvv_inverse_example_1.c b/examples/fnft_kdvv_inverse_example_1.c index 4285e12..52e1d26 100644 --- a/examples/fnft_kdvv_inverse_example_1.c +++ b/examples/fnft_kdvv_inverse_example_1.c @@ -37,7 +37,8 @@ INT main() COMPLEX * contspec = NULL; // five desired bound states and norming constants - // - bound states needs to be true imaginary positive + // - bound states have to be positive, purely imaginary numbers + // - bound states have to be in descending order // - every norming constants belong to the bound state at the same index // - the signs of the norming constants have to alternate regarding the order // of bound states. The norming constant of the biggest bound state has to be positive. diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m index da1120a..6bb45f9 100644 --- a/examples/mex_fnft_kdvv_inverse_example_1.m +++ b/examples/mex_fnft_kdvv_inverse_example_1.m @@ -25,7 +25,8 @@ desired_solitions = [4, 3, 2, 1]; % resulting bound states out of desired heights of solitions: -% - bound states have to be true and positive imaginary +% - bound states have to be positive, purely imaginary numbers +% - the bound states have to be in descendent order bound_states = 1i*sqrt(desired_solitions ./2); % desired norming constants From 40f123cd5be4a9c3db0c1d15e50b3838d75c9a7b Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Mon, 20 Apr 2026 19:45:21 +0200 Subject: [PATCH 60/82] Improve comments for doxygen in fnft_kdvv_inverse.h --- include/fnft_kdvv_inverse.h | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/include/fnft_kdvv_inverse.h b/include/fnft_kdvv_inverse.h index d558ee2..d50bf55 100644 --- a/include/fnft_kdvv_inverse.h +++ b/include/fnft_kdvv_inverse.h @@ -15,6 +15,7 @@ * * Contributors: * Sander Wahls (KIT) 2026 + * Fabian Fischer (Hiwi KIT) 2026 */ /** @@ -34,6 +35,17 @@ /** * @brief Fast inverse nonlinear Fourier transform for the * Korteweg-de Vries equation with vanishing boundary conditions. + * + * This routine computes the inverse nonlinear Fourier transform for the + * Korteweg-de Vries equation + * \f[ q_x + 6qq_{t} + q_{ttt}=0, \quad q=q(x,t), \f] + * of Gardner et al. ( + * Phys. Rev. Lett., 1967) + * + * The Fast inverse nonlinear Fourier transform for the Korteweg-de Vries equation + * uses the crum transformation. The main references are: + * - Prins and Wahls, "An accurate O(N^2) floating point algorithm for the Crum transform of the KdV equation," Communications in Nonlinear Science and Numerical Simulation 102, Article 105782, 2021. + * - The matlab project of P. Prins to the Crum Transformation * * @param[in] M Number of samples of the continuous spectrum. * @param[in,out] contspec Array of length M, contains samples @@ -41,29 +53,40 @@ * and \f$n=0,1,\dots,M-1\f$, of the to-be-inverted continuous spectrum in * ascending order (i.e., * \f$ \hat{q}(\xi_0), \hat{q}(\xi_1), \dots, \hat{q}(\xi_{M-1}) \f$). - * Please note that the routine currently might overwrite this array. + * Note: contspec functionality currently out of function (state 04/2026)! + * It will be implicitly assumed to \f$ 0 \f$ for all \f$ \xi \f$. * @param[in] XI Array of length 2, contains the position of the first and the last * sample of the continuous spectrum. + * Note: contspec functionality currently out of function (state 04/2026)! It is + * equal which values are chosen at the moment. * @param[in] K Number of discrete spectrum points. - * @param[in] bound_states Complex array of length K. Complex roots of - * \f$ a(\xi) \f$ in the upper half of the complex-plane. + * @param[in] bound_states Complex array of length K. Bound states have to be positive, + * purely imaginary numbers (lie on the upper half of the imaginary axis). The bound + * states have to be in descending order. To add a solition with height \f$ h_i \f$ + * the bound state have to be \f$ \gamma_i = \sqrt{ h_i/2} \f$. * @param[in] normconsts_or_residues Complex array of length K. Values of * either the norming constants \f$ b(\xi) \f$ or the residues * \f$ \frac{b(\xi)}{\partial{a(\xi)}/\partial{\xi}}\f$ at the values bound_states. + * The signs of the norming constants have to alternate regards to the order of the + * bound states. The sign of the normconst for the biggest eigenvalue has to be positive + * Note: currently this array is always interpreted as norming constants (state 04/2026)! + * Residues functionality is not implemented yet! * @param[in] D Number of samples of the to be generated signal q. Should be a * positive power of two. * @param[out] q Array of length D. Is filled with samples * \f$ q(t_n) \f$, where \f$ t_n = T[0] + n(T[1]-T[0])/(D-1) \f$ * and \f$n=0,1,\dots,D-1\f$, of the to-be-generated signal in ascending order * (i.e., \f$ q(t_0), q(t_1), \dots, q(t_{D-1}) \f$). + * Has to be preallocated by the user. * @param[in] T Array of length 2, contains the position in time of the first and - * of the last sample of q. It should be T[0] Date: Wed, 22 Apr 2026 17:54:50 +0200 Subject: [PATCH 61/82] Minor Refactoring in fnft_kdvv_inverse.c --- src/fnft_kdvv_inverse.c | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index c457b5e..a8ba2bd 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -611,15 +611,7 @@ INT fnft_kdvv_inverse( COMPLEX * const q, REAL const * const T, void *opts_ptr) -{ - - /* - D number of spatial points - out output of transformation - - */ - - +{ // checks for valid input --------------------------------------------------------------------------------- if (T[1] <= T[0]) { @@ -667,7 +659,7 @@ INT fnft_kdvv_inverse( // Initialize q for (UINT n=0; n Date: Wed, 22 Apr 2026 18:54:28 +0200 Subject: [PATCH 62/82] Refactoring of fnft_kdvv_inverse_testcases.c due to warning messages --- src/private/fnft__kdvv_inverse_testcases.c | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c index acd6557..cef222a 100644 --- a/src/private/fnft__kdvv_inverse_testcases.c +++ b/src/private/fnft__kdvv_inverse_testcases.c @@ -19,18 +19,16 @@ #define FNFT_ENABLE_SHORT_NAMES -#define DEBUG - #include #include "fnft__kdvv_inverse_testcases.h" #include "fnft__errwarn.h" #include "fnft__misc.h" -INT kdvv_print_spectrum(COMPLEX * bound_states_r, - COMPLEX * normconsts_r, - COMPLEX * contspec_r, - REAL * XI, +kdvv_print_spectrum(COMPLEX const * const bound_states_r, + COMPLEX const * const normconsts_r, + COMPLEX const * const contspec_r, + REAL const * const XI, const UINT M, const UINT D, const UINT K_r) @@ -112,8 +110,6 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, // Simple general tests for (UINT i=0; i Date: Wed, 22 Apr 2026 19:11:14 +0200 Subject: [PATCH 63/82] Adapting fnft__kdvv_inverse_testcases.h due to last commit --- include/private/fnft__kdvv_inverse_testcases.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h index a718663..edf46e6 100644 --- a/include/private/fnft__kdvv_inverse_testcases.h +++ b/include/private/fnft__kdvv_inverse_testcases.h @@ -61,10 +61,10 @@ FNFT_INT fnft__kdvv_inverse_testcases_get_spectrum_of_inverse( const FNFT_REAL err_bnd_normconst, const FNFT_REAL err_bnd_contspec); -FNFT_INT fnft__kdvv_print_spectrum( FNFT_COMPLEX * bound_states_r, - FNFT_COMPLEX * normconsts_r, - FNFT_COMPLEX * contspec_r, - FNFT_REAL * XI, +FNFT_INT fnft__kdvv_print_spectrum( FNFT_COMPLEX const * const bound_states_r, + FNFT_COMPLEX const * const normconsts_r, + FNFT_COMPLEX const * const contspec_r, + FNFT_REAL const * const XI, const UINT M, const UINT D, const UINT K_r); From 35cf355dc999d48ae81e433cab6ad2c218c2ff37 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Wed, 22 Apr 2026 19:12:28 +0200 Subject: [PATCH 64/82] Rename normconst_or_residues to norming_constants in context of inverse kdvv --- include/fnft_kdvv_inverse.h | 4 ++-- matlab/mex_fnft_kdvv_inverse.c | 18 +++++++++--------- matlab/mex_fnft_kdvv_inverse.m | 2 +- src/fnft_kdvv_inverse.c | 18 ++++++++---------- 4 files changed, 20 insertions(+), 22 deletions(-) diff --git a/include/fnft_kdvv_inverse.h b/include/fnft_kdvv_inverse.h index d50bf55..5429edc 100644 --- a/include/fnft_kdvv_inverse.h +++ b/include/fnft_kdvv_inverse.h @@ -64,7 +64,7 @@ * purely imaginary numbers (lie on the upper half of the imaginary axis). The bound * states have to be in descending order. To add a solition with height \f$ h_i \f$ * the bound state have to be \f$ \gamma_i = \sqrt{ h_i/2} \f$. - * @param[in] normconsts_or_residues Complex array of length K. Values of + * @param[in] norming_constants Complex array of length K. Values of * either the norming constants \f$ b(\xi) \f$ or the residues * \f$ \frac{b(\xi)}{\partial{a(\xi)}/\partial{\xi}}\f$ at the values bound_states. * The signs of the norming constants have to alternate regards to the order of the @@ -98,7 +98,7 @@ FNFT_INT fnft_kdvv_inverse( FNFT_REAL const * const XI, FNFT_UINT const K, FNFT_COMPLEX const * const bound_states, - FNFT_COMPLEX const * const normconsts_or_residues, + FNFT_COMPLEX const * const norming_constants, const FNFT_UINT D, FNFT_COMPLEX * const q, FNFT_REAL const * const T, diff --git a/matlab/mex_fnft_kdvv_inverse.c b/matlab/mex_fnft_kdvv_inverse.c index a30221f..c52187e 100644 --- a/matlab/mex_fnft_kdvv_inverse.c +++ b/matlab/mex_fnft_kdvv_inverse.c @@ -33,7 +33,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) FNFT_UINT M, K; FNFT_COMPLEX * contspec = NULL; FNFT_COMPLEX * bound_states = NULL; - FNFT_COMPLEX * normconsts_or_residues = NULL; + FNFT_COMPLEX * norming_constants = NULL; FNFT_REAL * XI; FNFT_UINT i; double *re, *im; @@ -45,7 +45,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) return; /* Check types and dimensions of the first seven inputs: contspec, XI, - bound_states, normconsts_or_residues, D, T */ + bound_states, norming_constants, D, T */ if ( nrhs < 6 ) mexErrMsgTxt("At least seven inputs expected."); @@ -56,7 +56,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if ( !mxIsEmpty(prhs[2]) && (!mxIsDouble(prhs[2]) || !mxIsComplex(prhs[2]) || mxGetM(prhs[2]) != 1) ) mexErrMsgTxt("Third input bound_states should be a complex row vector (double precision) or []. Try passing complex(double(bound_states(:)'))."); if ( !mxIsEmpty(prhs[3]) && (!mxIsDouble(prhs[3]) || !mxIsComplex(prhs[3]) || mxGetM(prhs[3]) != 1) ) - mexErrMsgTxt("Fourth input normconsts_or_residues should be a complex row vector (double precision) or []. Try passing complex(double(normconsts_or_residues(:)'))."); + mexErrMsgTxt("Fourth input norming_constants should be a complex row vector (double precision) or []. Try passing complex(double(norming_constants(:)'))."); if ( mxIsComplex(prhs[4]) || !mxIsDouble(prhs[4]) || mxGetNumberOfElements(prhs[4]) != 1 ) mexErrMsgTxt("Fifth input D should be a real scalar (double precision)."); if ( mxIsComplex(prhs[5]) || !mxIsDouble(prhs[5]) || mxGetM(prhs[5]) != 1 || mxGetN(prhs[5]) != 2 ) @@ -71,7 +71,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* Check values of first four inputs */ if ( K != mxGetNumberOfElements(prhs[3]) ) - mexErrMsgTxt("bound_states and normconsts_or_residues should have the same lengths."); + mexErrMsgTxt("bound_states and norming_constants should have the same lengths."); if ( T[0] >= T[1] ) mexErrMsgTxt("T(1) >= T(2)."); if ( XI[0] >= XI[1] ) @@ -118,10 +118,10 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) contspec = mxMalloc(M * sizeof(FNFT_COMPLEX)); if (K>0) { bound_states = mxMalloc(K * sizeof(FNFT_COMPLEX)); - normconsts_or_residues = mxMalloc(K * sizeof(FNFT_COMPLEX)); + norming_constants = mxMalloc(K * sizeof(FNFT_COMPLEX)); } if ( q == NULL || (M>0 && contspec == NULL) || (K>0 && bound_states == NULL) - || (K>0 && normconsts_or_residues == NULL) ) { + || (K>0 && norming_constants == NULL) ) { snprintf(msg, sizeof msg, "Out of memory."); goto on_error; } @@ -139,12 +139,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) re = mxGetPr(prhs[3]); im = mxGetPi(prhs[3]); for (i=0; i 0){ - return E_INVALID_ARGUMENT_MSG(normconsts_or_residues,The signs of the norming constants does not alternate!); + if (CREAL(norming_constants[i]*norming_constants[i+1]) > 0){ + return E_INVALID_ARGUMENT_MSG(norming_constants,The signs of the norming constants does not alternate!); } } // first checking if correct order of bound states exist, than checking if the first norming constant is positive! - if (CREAL(normconsts_or_residues[0]) < 0){ - return E_INVALID_ARGUMENT_MSG(normconsts_or_residues,The norming constant of the biggest bound state is negative!); + if (CREAL(norming_constants[0]) < 0){ + return E_INVALID_ARGUMENT_MSG(norming_constants,The norming constant of the biggest bound state is negative!); } // ---------------------------------------------------------------------------------------------------------- @@ -681,7 +679,7 @@ INT fnft_kdvv_inverse( for (UINT i=0; i Date: Wed, 22 Apr 2026 19:35:41 +0200 Subject: [PATCH 65/82] Fix error with fnft_kdvv_print_spectrum --- include/private/fnft__kdvv_inverse_testcases.h | 14 +++++++------- src/private/fnft__kdvv_inverse_testcases.c | 14 +++++++------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h index edf46e6..0922351 100644 --- a/include/private/fnft__kdvv_inverse_testcases.h +++ b/include/private/fnft__kdvv_inverse_testcases.h @@ -61,13 +61,13 @@ FNFT_INT fnft__kdvv_inverse_testcases_get_spectrum_of_inverse( const FNFT_REAL err_bnd_normconst, const FNFT_REAL err_bnd_contspec); -FNFT_INT fnft__kdvv_print_spectrum( FNFT_COMPLEX const * const bound_states_r, - FNFT_COMPLEX const * const normconsts_r, - FNFT_COMPLEX const * const contspec_r, - FNFT_REAL const * const XI, - const UINT M, - const UINT D, - const UINT K_r); +void fnft__kdvv_print_spectrum( FNFT_COMPLEX const * const bound_states_r, + FNFT_COMPLEX const * const normconsts_r, + FNFT_COMPLEX const * const contspec_r, + FNFT_REAL const * const XI, + const UINT M, + const UINT D, + const UINT K_r); #ifdef FNFT_ENABLE_SHORT_NAMES diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c index cef222a..599cf5c 100644 --- a/src/private/fnft__kdvv_inverse_testcases.c +++ b/src/private/fnft__kdvv_inverse_testcases.c @@ -25,13 +25,13 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" -kdvv_print_spectrum(COMPLEX const * const bound_states_r, - COMPLEX const * const normconsts_r, - COMPLEX const * const contspec_r, - REAL const * const XI, - const UINT M, - const UINT D, - const UINT K_r) +void kdvv_print_spectrum( COMPLEX const * const bound_states_r, + COMPLEX const * const normconsts_r, + COMPLEX const * const contspec_r, + REAL const * const XI, + const UINT M, + const UINT D, + const UINT K_r) { printf("Number of samples:\n D = %u\n", (unsigned int)D); From 348acc46c468e2b142829f1916c207d9a31d3ca1 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Wed, 22 Apr 2026 19:52:09 +0200 Subject: [PATCH 66/82] Formatting fnft_kdvv_inverse.c --- src/fnft_kdvv_inverse.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 435201a..195d2fd 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -599,16 +599,16 @@ static INT add_two_solitons( INT fnft_kdvv_inverse( - const UINT M, + UINT const M, COMPLEX * const contspec, REAL const * const XI, UINT const K, COMPLEX const * const bound_states, COMPLEX const * const norming_constants, - const UINT D, + UINT const D, COMPLEX * const q, REAL const * const T, - void *opts_ptr) + void * opts_ptr) { // checks for valid input --------------------------------------------------------------------------------- From a7d191d58f4746d1595ce0913e92633cfa5be8d0 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 30 Apr 2026 12:30:57 +0200 Subject: [PATCH 67/82] Fixing error with free of allocated storage: the pointers of all allocated storage needs to be declared before the first potential jump to leave_fun --- src/private/fnft__kdvv_inverse_testcases.c | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c index 599cf5c..3e7ea72 100644 --- a/src/private/fnft__kdvv_inverse_testcases.c +++ b/src/private/fnft__kdvv_inverse_testcases.c @@ -66,8 +66,13 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, const REAL err_bnd_contspec) { INT ret_code = SUCCESS; - + COMPLEX * q = NULL; + COMPLEX * contspec_r = NULL; + COMPLEX * bound_states_r = NULL; + COMPLEX * normconsts_r = NULL; + + q = malloc(params_i.D * sizeof(COMPLEX)); CHECK_NOMEM(q, ret_code, leave_fun); @@ -83,17 +88,14 @@ INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i, CHECK_RETCODE(ret_code, leave_fun); // Prepare forward fnft_kdvv - COMPLEX * contspec_r = NULL; contspec_r = malloc(params_i.M * sizeof(COMPLEX)); CHECK_NOMEM(contspec_r, ret_code, leave_fun); UINT K_r = params_i.D; - COMPLEX * bound_states_r = NULL; bound_states_r = malloc(K_r * sizeof(COMPLEX)); CHECK_NOMEM(bound_states_r, ret_code, leave_fun); - COMPLEX * normconsts_r = NULL; normconsts_r = malloc(K_r * sizeof(COMPLEX)); CHECK_NOMEM(normconsts_r, ret_code, leave_fun); From edecdd89d716f0fddb90e9b1e13cd6b55237f419 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 30 Apr 2026 12:35:43 +0200 Subject: [PATCH 68/82] Fixing false order of eigenvalues. in fnft_kdvv_inverse_test_1.c. fnft_kdvv_inverse requires eigenvalues in descending order due to a recently change. --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c index 57410bb..bc831b1 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c @@ -31,8 +31,8 @@ INT main() { INT ret_code = SUCCESS; - COMPLEX bound_states_i[5] = {I*SQRT(1.0/2.0), I*SQRT(2.0/2.0), I*SQRT(3.0/2.0), I*SQRT(4.0/2.0), I*SQRT(5.0/2.0)}; - COMPLEX normconsts_i[5] = {1*10, -1*0.1, 1*1, -1*1e-5, 1*1e7}; + COMPLEX bound_states_i[5] = {I*SQRT(5.0/2.0), I*SQRT(4.0/2.0), I*SQRT(3.0/2.0), I*SQRT(2.0/2.0), I*SQRT(1.0/2.0)}; + COMPLEX normconsts_i[5] = {1*1e7, -1*1e-5, 1*1, -1*0.1, 1*10}; fnft_kdvv_params kdvv_parameters = { .D = 256, From 7ec4205da04d0c047067260a349e46ca5d7baf5b Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 30 Apr 2026 13:10:32 +0200 Subject: [PATCH 69/82] Change in fnft_kdvv_inverse.c numerical constants from 0 to 0.0 --- src/fnft_kdvv_inverse.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 195d2fd..cd4ccb3 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -751,7 +751,7 @@ INT fnft_kdvv_inverse( // Filter out errors that result in negative values for (UINT i=0; i Date: Mon, 4 May 2026 20:08:24 +0200 Subject: [PATCH 70/82] Adapt fnft_kdvv_inverse_analytic_test so that output is compared with analytic signal for two x-values. --- .../fnft_kdvv_inverse_analytic_test.c | 120 +++++++++++------- 1 file changed, 76 insertions(+), 44 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index 3383d14..2f7243e 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -51,7 +51,7 @@ #define QUADRATIC_ERROR_SUM_TOLERANCE 1e-27 #define MAX_QUADRATIC_ERROR 1e-28 -REAL analytic_signal(REAL t, REAL x){ +static REAL analytic_signal(REAL t, REAL x){ return 12.0*(3.0 + 4.0 * COSH(2.0*t-8.0*x) + COSH(4.0*t-64.0*x))/POW((3.0*COSH(t-28.0*x) + COSH(3.0*t-36.0*x)), 2); } @@ -60,71 +60,101 @@ INT main() { INT ret_code = SUCCESS; - COMPLEX * contspec = NULL; - COMPLEX * q = NULL; COMPLEX * t_grid = NULL; - COMPLEX * analytic_q_a = NULL; - - REAL const x = 0.5; + COMPLEX * q_1 = NULL; + COMPLEX * q_2 = NULL; + COMPLEX * analytic_q_1 = NULL; + COMPLEX * analytic_q_2 = NULL; UINT D = 256; - UINT M = 10; - REAL XI[2] = {-2.0, 2.0}; - REAL T[2] = {-8.0, 12.0}; - - contspec = malloc(10*sizeof(COMPLEX)); - CHECK_NOMEM(contspec, ret_code, leave_fun); - - q = malloc(D * sizeof(COMPLEX)); - CHECK_NOMEM(q, ret_code, leave_fun); - - COMPLEX bound_states[K] = { I*2.0, I*1.0 }; - COMPLEX normconsts[K] = { CEXP(8.0*I*CPOW(bound_states[0], 3)*x), - -CEXP(8.0*I*CPOW(bound_states[1], 3)*x) }; - - ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states, normconsts, D, q, T, NULL); - CHECK_RETCODE(ret_code, leave_fun); + REAL const T[2] = {-8.0, 12.0}; t_grid = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(t_grid, ret_code, leave_fun); - - const COMPLEX eps_t = (T[1] - T[0])/(D - 1); - + + COMPLEX const eps_t = (T[1] - T[0])/(D - 1); + for (UINT n=0; n max_quadratic_error_1) {max_quadratic_error_1 = error_1;} } #ifdef DEBUG - misc_print_buf(D, t_grid, "t_grid"); - misc_print_buf(D, q, "output_inverse"); - misc_print_buf(D, analytic_q_a, "analytic_q"); + // misc_print_buf(K, bound_states_1, "bs1"); + // misc_print_buf(D, q_1, "q_1"); + printf("resulting max quadratic error: %e \n", max_quadratic_error_1); + printf("resulting quadratic error sum: %e \n", quadratic_error_sum_1); #endif - REAL analytic_q = 0.0; - REAL error = 0.0; - REAL max_quadratic_error = 0.0; - REAL quadratic_error_sum = 0.0; + + // Value #2 for x + q_2 = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(q_2, ret_code, leave_fun); + + analytic_q_2 = malloc(D * sizeof(COMPLEX)); + CHECK_NOMEM(analytic_q_2, ret_code, leave_fun); + + ret_code = fnft_kdvv_inverse(0, NULL, NULL, K, bound_states_1, normconsts_2, D, q_2, T, NULL); + CHECK_RETCODE(ret_code, leave_fun); + for (UINT i=0; i < D; i++) { + analytic_q_2[i] = analytic_signal(t_grid[i], x_2); + } + REAL error_2 = 0.0; + REAL max_quadratic_error_2 = 0.0; + REAL quadratic_error_sum_2 = 0.0; + for (UINT i=0; i < D; i++) { - error = POW((analytic_q_a[i] - q[i]), 2); - quadratic_error_sum += error; - if (error > max_quadratic_error) {max_quadratic_error = error;} + error_2 = CABS(CPOW((analytic_q_2[i] - q_2[i]), 2)); + quadratic_error_sum_2 += error_2; + if (error_2 > max_quadratic_error_2) {max_quadratic_error_2 = error_2;} } #ifdef DEBUG - printf("resulting max quadratic error: %e \n", max_quadratic_error); - printf("resulting quadratic error sum: %e \n", quadratic_error_sum); + // misc_print_buf(D, q_2, "q_2"); + printf("resulting max quadratic error: %e \n", max_quadratic_error_2); + printf("resulting quadratic error sum: %e \n", quadratic_error_sum_2); #endif - UINT is_quadratic_error_sum_in_tolerance = quadratic_error_sum < QUADRATIC_ERROR_SUM_TOLERANCE; - UINT is_max_quadratic_error_in_tolerance = max_quadratic_error < MAX_QUADRATIC_ERROR; + UINT is_quadratic_error_sum_in_tolerance = (quadratic_error_sum_1 < QUADRATIC_ERROR_SUM_TOLERANCE) && + (quadratic_error_sum_2 < QUADRATIC_ERROR_SUM_TOLERANCE); + UINT is_max_quadratic_error_in_tolerance = (max_quadratic_error_1 < MAX_QUADRATIC_ERROR) && + (max_quadratic_error_2 < MAX_QUADRATIC_ERROR); if (is_max_quadratic_error_in_tolerance && @@ -135,12 +165,14 @@ INT main() else { ret_code = FNFT_EC_TEST_FAILED; } + leave_fun: - free(contspec); - free(q); free(t_grid); - free(analytic_q_a); + free(q_1); + free(q_2); + free(analytic_q_1); + free(analytic_q_2); if (ret_code != SUCCESS) return EXIT_FAILURE; From 80882c228e54618c519dfd4d90de17967c91e225 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Mon, 4 May 2026 20:17:15 +0200 Subject: [PATCH 71/82] Add explanation to fnft_kdvv_inverse_analytic_test due to own experience of rounded intermediate results --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index 2f7243e..a27c2bd 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -77,6 +77,11 @@ INT main() for (UINT n=0; n Date: Mon, 4 May 2026 20:28:27 +0200 Subject: [PATCH 72/82] Adapt doxygen comment for fnft_kdvv_inverse.h --- include/fnft_kdvv_inverse.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/fnft_kdvv_inverse.h b/include/fnft_kdvv_inverse.h index 5429edc..0963681 100644 --- a/include/fnft_kdvv_inverse.h +++ b/include/fnft_kdvv_inverse.h @@ -71,8 +71,7 @@ * bound states. The sign of the normconst for the biggest eigenvalue has to be positive * Note: currently this array is always interpreted as norming constants (state 04/2026)! * Residues functionality is not implemented yet! - * @param[in] D Number of samples of the to be generated signal q. Should be a - * positive power of two. + * @param[in] D Number of samples of the to be generated signal q. * @param[out] q Array of length D. Is filled with samples * \f$ q(t_n) \f$, where \f$ t_n = T[0] + n(T[1]-T[0])/(D-1) \f$ * and \f$n=0,1,\dots,D-1\f$, of the to-be-generated signal in ascending order From 9671b046e7169c425e8dba67092c0574911b7cf7 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 4 Jun 2026 17:53:04 +0200 Subject: [PATCH 73/82] Adapt matlab help descriptions in mex-file for inverse kdvv --- matlab/mex_fnft_kdvv_inverse.m | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/matlab/mex_fnft_kdvv_inverse.m b/matlab/mex_fnft_kdvv_inverse.m index 2651830..29245bc 100644 --- a/matlab/mex_fnft_kdvv_inverse.m +++ b/matlab/mex_fnft_kdvv_inverse.m @@ -2,7 +2,7 @@ % Korteweg-de Vries equation with vanishing boundaries. % % q = MEX_FNFT_KDVVV_INVERSE(contspec, XI, bound_states, ... -% normconsts_or_residuals, D, T); +% norming_constants, D, T); % % DESCRIPTION % Provides an interface to the C routine fnft_kdvv_inverse. @@ -14,19 +14,20 @@ % b-scattering coefficient on an equidistant grid. % Pass [] if the continuous spectrum is zero % (i.e., a multi-soliton is desired) -% XI Real 1x2 vector, contains the location of the first and -% the last sample in contspec +% Note: contspec functionality currently out of function +% (state 04/2026)! +% It will be implicitly assumed to \f$ 0 \f$ for all \f$ \xi \f$. +% XI Real 1x2 vector, contains the position of the first and the last +% sample of the continuous spectrum. +% Note: contspec functionality currently out of function +% (state 04/2026)! +% It does not matter which values are chosen at the moment. % bound_states Complex row vector, contains the desired bound states. -% Pass [] if the discrete spectrum is empty. % norming_constants Complex row vector, same length as bound_states. -% Contains the corresponding norming constants (default) or, -% if the corresponding option is passed, residues. Pass [] -% if the discrete spectrum is empty. -% D Real scalar, number of time domain samples; must be a -% positive power of two +% Contains the corresponding norming constants. +% D Real scalar, number of time domain samples. % T Real 1x2 vector, contains the location of the first and % the last sample in q -% kappa +1.0 or -1.0 % % OUTPUTS % q Complex row vector of length D @@ -46,4 +47,5 @@ % along with this program. If not, see . % % Contributors: -% Sander Wahls (KIT) 2026. \ No newline at end of file +% Sander Wahls (KIT) 2026. +% Fabian Fischer (Hiwi KIT) 2026. \ No newline at end of file From 03a8e051be8309ee0d081cd5796089b36f9ed7a2 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 4 Jun 2026 17:55:55 +0200 Subject: [PATCH 74/82] Clarification in doxygen comment in fnft_kdvv_inverse.h --- include/fnft_kdvv_inverse.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/fnft_kdvv_inverse.h b/include/fnft_kdvv_inverse.h index 0963681..b566788 100644 --- a/include/fnft_kdvv_inverse.h +++ b/include/fnft_kdvv_inverse.h @@ -57,8 +57,8 @@ * It will be implicitly assumed to \f$ 0 \f$ for all \f$ \xi \f$. * @param[in] XI Array of length 2, contains the position of the first and the last * sample of the continuous spectrum. - * Note: contspec functionality currently out of function (state 04/2026)! It is - * equal which values are chosen at the moment. + * Note: contspec functionality currently out of function (state 04/2026)! It does not + * matter which values are chosen at the moment. * @param[in] K Number of discrete spectrum points. * @param[in] bound_states Complex array of length K. Bound states have to be positive, * purely imaginary numbers (lie on the upper half of the imaginary axis). The bound From b9fb84b3954d762fc79f01aaf0fbfa4f0c0e994e Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 4 Jun 2026 18:08:35 +0200 Subject: [PATCH 75/82] Adding check for valid input (check for duplicity of bound states) in fnft_kdvv_inverse. Therefore the constant THRESHOLD_EQUALITY_BOUND_STATES is introduced. --- src/fnft_kdvv_inverse.c | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index cd4ccb3..0352690 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -24,6 +24,8 @@ #define FNFT_ENABLE_SHORT_NAMES +#define THRESHOLD_EQUALITY_BOUND_STATES 1e-6 + #include "fnft_kdvv_inverse.h" static INT find_first_positive_value( @@ -638,6 +640,12 @@ INT fnft_kdvv_inverse( if (CREAL(norming_constants[i]*norming_constants[i+1]) > 0){ return E_INVALID_ARGUMENT_MSG(norming_constants,The signs of the norming constants does not alternate!); } + + for (UINT j=i+1; j Date: Thu, 4 Jun 2026 18:13:35 +0200 Subject: [PATCH 76/82] Adaption on comment to the observed wrong behavior in fnft_kdvv_inverse_analytic_test --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index a27c2bd..d67f73f 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -80,8 +80,9 @@ INT main() // Calculating norming constants for both x-values. // Attention: both norming constants arrays need to be calculated here! If there are further - // commands between the calculation of both norming constants arrays, the processor calculates the - // 2. norming constants with rounded intermediate results. + // commands between the calculation of both norming constants arrays and the compiler optimize + // drastically, there may be errors. Is is speculated, that in that case the processor + // calculates the 2. norming constants with rounded intermediate results. COMPLEX const bound_states_1[K] = { I*2.0, I*1.0 }; From ed269cfa3f3b2dfdcd336299d3dd9e3206df1d1d Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 4 Jun 2026 18:17:53 +0200 Subject: [PATCH 77/82] Changing Value of constant THRESHOLD_EQUALITY_BOUND_STATES --- src/fnft_kdvv_inverse.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 0352690..01a03b5 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -24,7 +24,7 @@ #define FNFT_ENABLE_SHORT_NAMES -#define THRESHOLD_EQUALITY_BOUND_STATES 1e-6 +#define THRESHOLD_EQUALITY_BOUND_STATES 1e-15 #include "fnft_kdvv_inverse.h" From 00654495426111b1128da054a2335ebd97315d98 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Mon, 15 Jun 2026 18:52:26 +0200 Subject: [PATCH 78/82] Fixing compiler errors due to wrong usage of CABS() --- src/fnft_kdvv_inverse.c | 2 +- src/private/fnft__kdvv_inverse_testcases.c | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c index 01a03b5..7c18024 100644 --- a/src/fnft_kdvv_inverse.c +++ b/src/fnft_kdvv_inverse.c @@ -642,7 +642,7 @@ INT fnft_kdvv_inverse( } for (UINT j=i+1; j 0) && is_bsr_pure_imaginary; - UINT is_ncr_real = CABS(CIMAG(ncr)) < 1e-8; + UINT is_ncr_real = FABS(CIMAG(ncr)) < 1e-8; if (is_bsr_pure_imaginary && is_bsr_positive_imaginary && From 237cf1cf13909223751913a059192808073f4198 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Mon, 15 Jun 2026 18:54:10 +0200 Subject: [PATCH 79/82] Refactoring fnft_kdvv_inverse_analytic_test, adapting error bounds due to failed test at linux computer --- .../fnft_kdvv_inverse_analytic_test.c | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index d67f73f..aabc319 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -48,8 +48,8 @@ #define K 2 #define DEBUG -#define QUADRATIC_ERROR_SUM_TOLERANCE 1e-27 -#define MAX_QUADRATIC_ERROR 1e-28 +#define QUADRATIC_ERROR_SUM_TOLERANCE 1e-26 +#define MAX_QUADRATIC_ERROR 1e-27 static REAL analytic_signal(REAL t, REAL x){ return 12.0*(3.0 + 4.0 * COSH(2.0*t-8.0*x) + COSH(4.0*t-64.0*x))/POW((3.0*COSH(t-28.0*x) + COSH(3.0*t-36.0*x)), 2); @@ -84,22 +84,22 @@ INT main() // drastically, there may be errors. Is is speculated, that in that case the processor // calculates the 2. norming constants with rounded intermediate results. - COMPLEX const bound_states_1[K] = { I*2.0, I*1.0 }; + COMPLEX const bound_states[K] = { I*2.0, I*1.0 }; REAL const x_1 = -0.1; - COMPLEX normconsts_1[K] = { CREAL(CEXP(8.0*I*CPOW(bound_states_1[0], 3)*x_1)), - -CREAL(CEXP(8.0*I*CPOW(bound_states_1[1], 3)*x_1)) }; + COMPLEX normconsts_1[K] = { CREAL(CEXP(8.0*I*CPOW(bound_states[0], 3)*x_1)), + -CREAL(CEXP(8.0*I*CPOW(bound_states[1], 3)*x_1)) }; REAL const x_2 = 0.5; - COMPLEX normconsts_2[K] = { CREAL(CEXP(8.0*I*CPOW(bound_states_1[0], 3)*x_2)), - -CREAL(CEXP(8.0*I*CPOW(bound_states_1[1], 3)*x_2)) }; + COMPLEX normconsts_2[K] = { CREAL(CEXP(8.0*I*CPOW(bound_states[0], 3)*x_2)), + -CREAL(CEXP(8.0*I*CPOW(bound_states[1], 3)*x_2)) }; // Value #1 for x q_1 = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(q_1, ret_code, leave_fun); - ret_code = fnft_kdvv_inverse(0, NULL, NULL, K, bound_states_1, normconsts_1, D, q_1, T, NULL); + ret_code = fnft_kdvv_inverse(0, NULL, NULL, K, bound_states, normconsts_1, D, q_1, T, NULL); CHECK_RETCODE(ret_code, leave_fun); analytic_q_1 = malloc(D * sizeof(COMPLEX)); @@ -120,7 +120,7 @@ INT main() } #ifdef DEBUG - // misc_print_buf(K, bound_states_1, "bs1"); + // misc_print_buf(K, bound_states, "bs1"); // misc_print_buf(D, q_1, "q_1"); printf("resulting max quadratic error: %e \n", max_quadratic_error_1); printf("resulting quadratic error sum: %e \n", quadratic_error_sum_1); @@ -134,7 +134,7 @@ INT main() analytic_q_2 = malloc(D * sizeof(COMPLEX)); CHECK_NOMEM(analytic_q_2, ret_code, leave_fun); - ret_code = fnft_kdvv_inverse(0, NULL, NULL, K, bound_states_1, normconsts_2, D, q_2, T, NULL); + ret_code = fnft_kdvv_inverse(0, NULL, NULL, K, bound_states, normconsts_2, D, q_2, T, NULL); CHECK_RETCODE(ret_code, leave_fun); for (UINT i=0; i < D; i++) { From 644f6b5616345916137af81c9c7cd1d13d47c646 Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 18 Jun 2026 17:27:14 +0200 Subject: [PATCH 80/82] Changing names of arguments --- .../private/fnft__kdvv_inverse_testcases.h | 10 ++++----- src/private/fnft__kdvv_inverse_testcases.c | 22 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h index 0922351..7046719 100644 --- a/include/private/fnft__kdvv_inverse_testcases.h +++ b/include/private/fnft__kdvv_inverse_testcases.h @@ -60,14 +60,14 @@ FNFT_INT fnft__kdvv_inverse_testcases_get_spectrum_of_inverse( const FNFT_REAL err_bnd_spurious_bound_states, const FNFT_REAL err_bnd_normconst, const FNFT_REAL err_bnd_contspec); - -void fnft__kdvv_print_spectrum( FNFT_COMPLEX const * const bound_states_r, - FNFT_COMPLEX const * const normconsts_r, - FNFT_COMPLEX const * const contspec_r, + +void fnft__kdvv_print_spectrum( FNFT_COMPLEX const * const bound_states, + FNFT_COMPLEX const * const normconsts, + FNFT_COMPLEX const * const contspec, FNFT_REAL const * const XI, const UINT M, const UINT D, - const UINT K_r); + const UINT K); #ifdef FNFT_ENABLE_SHORT_NAMES diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c index 3e7ea72..df4095a 100644 --- a/src/private/fnft__kdvv_inverse_testcases.c +++ b/src/private/fnft__kdvv_inverse_testcases.c @@ -25,13 +25,13 @@ #include "fnft__errwarn.h" #include "fnft__misc.h" -void kdvv_print_spectrum( COMPLEX const * const bound_states_r, - COMPLEX const * const normconsts_r, - COMPLEX const * const contspec_r, +void kdvv_print_spectrum( COMPLEX const * const bound_states, + COMPLEX const * const normconsts, + COMPLEX const * const contspec, REAL const * const XI, const UINT M, const UINT D, - const UINT K_r) + const UINT K) { printf("Number of samples:\n D = %u\n", (unsigned int)D); @@ -41,18 +41,18 @@ void kdvv_print_spectrum( COMPLEX const * const bound_states_r, FNFT_REAL xi = XI[0] + i*eps_xi; printf(" continuous_spectrum(xi=%f) \t= %g + %gI\n", (double)xi, - (double)FNFT_CREAL(contspec_r[i]), - (double)FNFT_CIMAG(contspec_r[i]) + (double)FNFT_CREAL(contspec[i]), + (double)FNFT_CIMAG(contspec[i]) ); } printf("Discrete spectrum:\n"); - for (FNFT_UINT i=0; i Date: Thu, 18 Jun 2026 17:27:46 +0200 Subject: [PATCH 81/82] Adding doxygen comments --- .../private/fnft__kdvv_inverse_testcases.h | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h index 7046719..db131c8 100644 --- a/include/private/fnft__kdvv_inverse_testcases.h +++ b/include/private/fnft__kdvv_inverse_testcases.h @@ -60,7 +60,26 @@ FNFT_INT fnft__kdvv_inverse_testcases_get_spectrum_of_inverse( const FNFT_REAL err_bnd_spurious_bound_states, const FNFT_REAL err_bnd_normconst, const FNFT_REAL err_bnd_contspec); - + + +/** + * @brief Routine to print the continuous spectrum and discrete spectrum, consisting of + * bound states and norming constants, as a result of \link fnft_kdvv \endlink. + * + * This routine is used by the tests for \link fnft_kdvv_inverse \endlink. + * + * @param[in] bound_states + * @param[in] normconsts + * @param[in] contspec + * @param[in] XI Array of length 2, contains the position of the first and the last + * sample of the continuous spectrum. + * @param[in] M Number of points at which the continuous spectrum is computed. + * @param[in] D Number of samples of the potential. + * @param[in] K Number of bound states (same than number of norming constants) + * @return void + * + * @ingroup kdv + */ void fnft__kdvv_print_spectrum( FNFT_COMPLEX const * const bound_states, FNFT_COMPLEX const * const normconsts, FNFT_COMPLEX const * const contspec, From 360541470939241e7026ae8f659726c1c4c3966f Mon Sep 17 00:00:00 2001 From: ffischer7 Date: Thu, 18 Jun 2026 17:28:21 +0200 Subject: [PATCH 82/82] Remove DEBUG constant --- test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c | 1 - test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c | 1 - 2 files changed, 2 deletions(-) diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c index d67f73f..b1bd7f7 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c @@ -47,7 +47,6 @@ #define K 2 -#define DEBUG #define QUADRATIC_ERROR_SUM_TOLERANCE 1e-27 #define MAX_QUADRATIC_ERROR 1e-28 diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c index 64bc384..b825ff0 100644 --- a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c +++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c @@ -28,7 +28,6 @@ /* This is a testcase with an assymetric window */ #define K_I 8 -#define DEBUG INT main() {