Merge pull request #35 from FastNFT/release-0.2

Release 0.2.1

Author: wahls

CHANGELOG.md

@@ -1,6 +1,12 @@
 # Changelog
 
-## [0.2.0] -- 2018-09
+## [0.2.1] -- 2018-09-28
+
+### Fixed
+
+- Computation of norming constants from residues in fnft_nsev_inverse was incorrect when continuous spectrum was present.
+
+## [0.2.0] -- 2018-09-21
 
 ### Added
 

CMakeLists.txt

@@ -21,7 +21,7 @@ project(fnft C)
 
 set(FNFT_VERSION_MAJOR 0)
 set(FNFT_VERSION_MINOR 2)
-set(FNFT_VERSION_PATCH 0)
+set(FNFT_VERSION_PATCH 1)
 set(FNFT_VERSION_SUFFIX "") # should not be longer than FNFT_SUFFIX_MAXLEN
 set(FNFT_VERSION ${FNFT_VERSION_MAJOR}.${FNFT_VERSION_MINOR}.${FNFT_VERSION_PATCH}${FNFT_VERSION_SUFFIX})
 

README.md

@@ -9,10 +9,12 @@ FNFT is a software library for the fast numerical computation of (inverse) nonli
 ### Forward Transforms
 
 * Nonlinear Schroedinger equation
+
     * Vanishing boundary conditions
       * Reflection coefficient and/or scattering coefficients (a and b)
       * Bound states (eigenvalues)
       * Norming constants and/or residues
+
     * Periodic boundary conditions
       * Main spectrum
       * Auxiliary spectrum
@@ -23,6 +25,7 @@ FNFT is a software library for the fast numerical computation of (inverse) nonli
 ### Inverse Transforms
 
 * Nonlinear Schroedinger equation
+
     * Vanishing boundary conditions
       * Inversion of reflection coefficients, b-scattering coefficients or the inverse Fourier transform of the b-coefficient
       * Bound states (eigenvalues) can be added with arbitrary norming constants/residuals

src/fnft_nsev_inverse.c

@@ -29,6 +29,8 @@
 #include "fnft__fft_wrapper.h"
 #include "fnft__poly_specfact.h"
 #include "fnft__misc.h"
+#include "fnft__nse_scatter.h"
+
 
 static fnft_nsev_inverse_opts_t default_opts = {
     .discretization = nse_discretization_2SPLIT2A,
@@ -704,6 +706,7 @@ static INT add_discrete_spectrum(
     COMPLEX * bnd_states_conj = NULL;
     COMPLEX * bnd_states_diff = NULL;
     COMPLEX tmp;
+    COMPLEX * acoeff_cs = NULL;
     COMPLEX * phi = NULL;
     COMPLEX * psi = NULL;
     UINT i,j;
@@ -716,9 +719,10 @@ static INT add_discrete_spectrum(
     norm_consts = malloc(K * sizeof(COMPLEX));
     bnd_states_conj = malloc(K * sizeof(COMPLEX));
     bnd_states_diff = malloc(K * sizeof(COMPLEX));
+    acoeff_cs = malloc(3*K * sizeof(COMPLEX));
     if (t == NULL || bnd_states == NULL ||
             bnd_states_conj == NULL || bnd_states_diff == NULL ||
-            norm_consts == NULL) {
+            norm_consts == NULL || acoeff_cs == NULL) {
         ret_code = E_NOMEM;
         goto leave_fun;
     }
@@ -768,16 +772,29 @@ static INT add_discrete_spectrum(
     //CDT works with norming constants so residues need to be converted to
     // norming constants
     if (opts_ptr->discspec_type == fnft_nsev_inverse_dstype_RESIDUES){
+        if (contspec_flag != 0){
+            //When continuous spectrum exists, the non-solitonic component
+            // of the potential contributes to the residues and needs to be
+            // removed.
+            UINT trunc_index;
+            trunc_index = D;
+            ret_code = nse_scatter_bound_states(D, q, T, &trunc_index, K,
+                    bnd_states, acoeff_cs, acoeff_cs+K, acoeff_cs+2*K, nse_discretization_BO);
+            CHECK_RETCODE(ret_code, leave_fun);
+        }
+        else {
+            for (i = 0; i < K; i++)
+                acoeff_cs[i] = 1;
+        }
         for (i = 0; i < K; i++){
-            tmp = 1;
+            tmp = acoeff_cs[i];
             for (j = 0; j < K; j++){
                 if (j != i){
                     tmp = tmp*(bnd_states[i]-bnd_states[j])/(bnd_states[i]-CONJ(bnd_states[j]));
                 }
             }
             norm_consts[i]=(norm_consts[i]/bnd_states_diff[i])*tmp;
-        }
-
+        }        
     }
 
     // In absence of contspec some tricks can be used to
@@ -896,6 +913,7 @@ static INT add_discrete_spectrum(
         free(bnd_states_conj);
         free(bnd_states_diff);
         free(norm_consts);
+        free(acoeff_cs);
         return ret_code;
 }
 

test/fnft_nsev_inverse/fnft_nsev_inverse_test_against_forward/fnft_nsev_inverse_test_against_forward_w_discrete.inc

@@ -30,7 +30,7 @@
 #include <stdio.h>
 
 static INT fnft_nsev_inverse_test(const UINT D, const UINT M,
-   const REAL error_bound, fnft_nsev_inverse_opts_t * opts_ptr)
+   const REAL error_bound, fnft_nsev_inverse_opts_t * opts_inv_ptr)
 {
     if (M < D)
         return E_INVALID_ARGUMENT(M);
@@ -70,9 +70,13 @@ static INT fnft_nsev_inverse_test(const UINT D, const UINT M,
     //Setting up forward transform options
     fnft_nsev_opts_t opts_fwd = fnft_nsev_default_opts();
     opts_fwd.discretization = nse_discretization_2SPLIT4B;
+    if (opts_inv_ptr->discspec_type == fnft_nsev_inverse_dstype_RESIDUES)
+       opts_fwd.discspec_type = nsev_dstype_RESIDUES;
+    else if (opts_inv_ptr->discspec_type == fnft_nsev_inverse_dstype_NORMING_CONSTANTS)
+       opts_fwd.discspec_type = nsev_dstype_NORMING_CONSTANTS;
 
     //Obtain XI grid
-    ret_code = fnft_nsev_inverse_XI(D, T, M, XI, opts_ptr->discretization);
+    ret_code = fnft_nsev_inverse_XI(D, T, M, XI, opts_inv_ptr->discretization);
     CHECK_RETCODE(ret_code, leave_fun);
 
     //Perform full forward NFT
@@ -82,7 +86,7 @@ static INT fnft_nsev_inverse_test(const UINT D, const UINT M,
 
     //Perform full inverse NFT
     ret_code = fnft_nsev_inverse(M, contspec, XI, *K_ptr, bound_states, normconsts_or_residues, D, q, T,
-            kappa, opts_ptr);
+            kappa, opts_inv_ptr);
     CHECK_RETCODE(ret_code, leave_fun);
 
     REAL error = misc_rel_err(D, q, q_exact);

test/fnft_nsev_inverse/fnft_nsev_inverse_test_against_forward/fnft_nsev_inverse_test_against_forward_w_discrete_2split2A.c

@@ -27,21 +27,34 @@ int main()
     REAL error_bound;
     INT ret_code = SUCCESS;
 
-    fnft_nsev_inverse_opts_t opts = fnft_nsev_inverse_default_opts();
-    opts.discretization = nse_discretization_2SPLIT2A;
-    opts.discspec_type = fnft_nsev_inverse_dstype_NORMING_CONSTANTS;
+    fnft_nsev_inverse_opts_t opts_inv = fnft_nsev_inverse_default_opts();
+    opts_inv.discretization = nse_discretization_2SPLIT2A;
+    opts_inv.discspec_type = fnft_nsev_inverse_dstype_NORMING_CONSTANTS;      
 
+    D = 512;
+    M = 2*D;
+    error_bound = 0.014;
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+    D = D*2;
+    M = 2*D;
+    error_bound = error_bound/4;
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    opts_inv.discspec_type = fnft_nsev_inverse_dstype_RESIDUES;      
 
     D = 512;
     M = 2*D;
     error_bound = 0.014;
-    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts);
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
     CHECK_RETCODE(ret_code, leave_fun);
 
     D = D*2;
     M = 2*D;
     error_bound = error_bound/4;
-    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts);
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
     CHECK_RETCODE(ret_code, leave_fun);
 
 

test/fnft_nsev_inverse/fnft_nsev_inverse_test_against_forward/fnft_nsev_inverse_test_against_forward_w_discrete_2split2_modal.c

@@ -27,21 +27,34 @@ int main()
     REAL error_bound;
     INT ret_code = SUCCESS;
 
-    fnft_nsev_inverse_opts_t opts = fnft_nsev_inverse_default_opts();
-    opts.discretization = nse_discretization_2SPLIT2_MODAL;
-    opts.discspec_type = fnft_nsev_inverse_dstype_NORMING_CONSTANTS;
+    fnft_nsev_inverse_opts_t opts_inv = fnft_nsev_inverse_default_opts();
+    opts_inv.discretization = nse_discretization_2SPLIT2_MODAL;
+    opts_inv.discspec_type = fnft_nsev_inverse_dstype_NORMING_CONSTANTS;
 
+    D = 512;
+    M = 2*D;
+    error_bound = 0.014;
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+    D = D*2;
+    M = 2*D;
+    error_bound = error_bound/4;
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    opts_inv.discspec_type = fnft_nsev_inverse_dstype_RESIDUES;
 
     D = 512;
     M = 2*D;
     error_bound = 0.014;
-    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts);
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
     CHECK_RETCODE(ret_code, leave_fun);
 
     D = D*2;
     M = 2*D;
     error_bound = error_bound/4;
-    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts);
+    ret_code = fnft_nsev_inverse_test(D, M, error_bound, &opts_inv);
     CHECK_RETCODE(ret_code, leave_fun);