jfm dev

README.md

@@ -13,45 +13,65 @@ Compute various exponential sums involving arbitrary point distributions in opti
 
 ### Dependencies
 
-The basic libraries need a C++ compiler, GNU make, FFTW, and optionally OpenMP (the makefile can be adjusted for single-threaded).
-The fortran wrappers need a fortran compiler.
-See settings in the `makefile`.
+For the basic libraries
+
+- C++ compiler
+- GNU make
+- FFTW
+- Optional OpenMP (the makefile can be adjusted for single-threaded)
+
+For the fotran wrappers
+
+- Fortran compiler (see settings in the makefile)
 
 On a fedora linux system, the dependencies can be installed as follows:
 
-`sudo yum install git fftw3 fftw3-devel libgomp`
+```bash
+sudo yum install git fftw3 fftw3-devel libgomp
+```
+
+On Ubuntu:
+
+```bash
+sudo apt-get install git libfftw3-dev
+```
 
 ### Installation
 
-1. Download using `git`, `svn`, or as a zip (see green button above).
-1. `cp makefile.dist makefile`
-1. edit `makefile` for your system
-1. `make`. This will compile the library then run a set of multi-threaded and single-threaded speed tests  
+- Clone using git (or checkout using svn, or download as a zip -- see green button above)
+- Copy makefile.dist to makefile and edit for your system
+- Compile using:
+
+```bash
+make
+```
+This will compile the library then run a set of multi-threaded and single-threaded speed tests  
 
 Other useful make modes include:
 
-  `make test1d` : small accuracy test for components in 1D. Analogously for 2D, 3D  
-  `make spreadtestnd` : benchmark the spreader routines, all dimensions  
-  `make examples/testutils` : test various low-level utilities  
-  `make nufft` : compile the library without testing  
-  `make fortran` : compile and test the fortran interfaces  
+```bash
+make test1d # small accuracy test for components in 1D. Analogously for 2D, 3D  
+make spreadtestnd # benchmark the spreader routines, all dimensions  
+make examples/testutils # test various low-level utilities  
+make nufft # compile the library without testing  
+make fortran # compile and test the fortran interfaces  
+```
 
 ### Contents of this package
 
-  `src` : main library source and headers.  
-  `examples` : test codes (drivers) which verify libaries are working correctly, perform speed tests, and show how to call them. In this directory are the useful scripts:
-  - `nuffttestnd.sh` : benchmark and display accuracy for all types and dimensions (3x3 = 9 in total) of NUFFT at fixed requested tolerance  
-  - `checkallaccs.sh dim` (where `dim` is 1, 2, or 3) : sweep over all tolerances checking the spreader and NUFFT at a single dimension  
-
-  `examples/results` : accuracy and timing outputs.  
-  `contrib` : 3rd-party code.  
-  `fortran` : wrappers and drivers for Fortran.
-  `matlab` : wrappers and examples for MATLAB. (Not yet working)  
-  `devel` : various obsolete or in-development codes (experts only)  
-  `doc` : the manual (not yet there)  
-  `README.md`  
-  `LICENSE`  
-  `makefile.dist` : GNU makefile (user should first copy to `makefile`)  
+- `src` : main library source and headers.  
+- `examples` : test codes (drivers) which verify libaries are working correctly, perform speed tests, and show how to call them. 
+- `examples/nuffttestnd.sh` : benchmark and display accuracy for all types and dimensions (3x3 = 9 in total) of NUFFT at fixed requested tolerance  
+- `examples/checkallaccs.sh [dim]` : sweep over all tolerances checking the spreader and NUFFT at a single dimension;  [dim] is 1, 2, or 3
+- `examples/results` : accuracy and timing outputs.  
+- `contrib` : 3rd-party code.  
+- `fortran` : wrappers and drivers for Fortran.
+- `matlab` : wrappers and examples for MATLAB. (Not yet working)  
+- `devel` : various obsolete or in-development codes (experts only)  
+- `doc` : the manual (not yet there)  
+- `README.md`
+- `LICENSE`
+- `makefile.dist` : GNU makefile (user should first copy to `makefile`)  
 
 ### Notes
 

contrib/testi0.cpp

@@ -6,6 +6,8 @@ int main(int argc, char* argv[])
 // speed and accuracy tests for I_0 modified Bessel of 1st kind in besseli.cpp.
 // Barnett 2/2/17
 {
+  (void)argc; //tell compiler this variable is unused
+  (void)argv; //tell compiler this variable is unused
   double R=100.0;    // upper lim
   int n=1e7;
    CNTime timer; timer.start();

devel/qtpro/finufft.pro

@@ -0,0 +1,48 @@
+QT += core network
+QT -= gui
+
+CONFIG += c++11
+
+CONFIG -= app_bundle #Please apple, don't make a bundle
+
+DESTDIR = bin
+OBJECTS_DIR = build
+MOC_DIR=build
+TARGET = finufft_qtpro
+TEMPLATE = app
+
+DEFINES += NEED_EXTERN_C
+
+HEADERS += \
+    ../../src/cnufftspread.h \
+    ../../src/common.h \
+    ../../src/dirft.h \
+    ../../src/finufft.h \
+    ../../src/twopispread.h \
+    ../../src/utils.h \
+    ../../contrib/besseli.h \
+    ../../contrib/legendre_rule_fast.h
+
+SOURCES += \
+    ../../src/cnufftspread.cpp \
+    ../../src/common.cpp \
+    ../../src/dirft1d.cpp \
+    ../../src/dirft2d.cpp \
+    ../../src/dirft3d.cpp \
+    ../../src/finufft1d.cpp \
+    ../../src/finufft2d.cpp \
+    ../../src/finufft3d.cpp \
+    ../../src/twopispread.cpp \
+    ../../src/utils.cpp \
+    ../../contrib/besseli.cpp \
+    ../../contrib/testi0.cpp \
+    ../../contrib/legendre_rule_fast.c
+
+#FFTW
+LIBS += -lfftw3 -lfftw3_threads
+
+#OPENMP
+!macx {
+  QMAKE_LFLAGS += -fopenmp
+  QMAKE_CXXFLAGS += -fopenmp
+}

src/cnufftspread.cpp

@@ -333,7 +333,7 @@ bool set_thread_index_box(long *i1th,long *i2th,long *i3th,long N1,long N2,long
   i2th[0] = 0; i2th[1] = N2-1;
   i3th[0] = 0; i3th[1] = N3-1;
   // set Ntop the lowest (in grid order) nontrivial dimension, ie Nd for d = #dims...
-  BIGINT ith[2], Ntop = N1, dims=1;
+  BIGINT ith[2], Ntop = N1;
   if (ndims==2) Ntop = N2;
   if (ndims==3) Ntop = N3;
   if (N2==1 || !opts.checkerboard) {  // slice only along one dim
@@ -408,6 +408,8 @@ std::vector<long> compute_sort_indices(long M,double *kx, double *ky, double *kz
    * Magland, Dec 2016; Barnett tweaked so doesn't examine ky in 1d, or kz in 1d or 2d.
    */
 {
+  (void)kx; //tell compiler this is an unused variable
+  (void)N1; //tell compiler this is an unused variable
   bool isky=(N2>1), iskz=(N3>1);           // are ky,kz available? cannot access if not!
   std::vector<long> counts(N2*N3);
   for (long j=0; j<N2*N3; j++)

src/common.cpp

@@ -1,5 +1,12 @@
 #include "common.h"
+
+#ifdef NEED_EXTERN_C
+extern "C" {
+    #include "../contrib/legendre_rule_fast.h"
+}
+#else
 #include "../contrib/legendre_rule_fast.h"
+#endif
 #include <fftw3.h>
 #include <math.h>
 #include <stdio.h>

src/finufft1d.cpp

@@ -45,6 +45,7 @@ int finufft1d1(BIGINT nj,double* xj,dcomplex* cj,int iflag,double eps,BIGINT ms,
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1; set_nf_type12(ms,opts,spopts,&nf1);
   cout << scientific << setprecision(15);  // for debug
 
@@ -132,6 +133,7 @@ int finufft1d2(BIGINT nj,double* xj,dcomplex* cj,int iflag,double eps,BIGINT ms,
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1; set_nf_type12(ms,opts,spopts,&nf1);
   cout << scientific << setprecision(15);  // for debug
 
@@ -221,6 +223,7 @@ int finufft1d3(BIGINT nj,double* xj,dcomplex* cj,int iflag, double eps, BIGINT n
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1;
   double X1,C1,S1,D1,h1,gam1;
   cout << scientific << setprecision(15);  // for debug

src/finufft2d.cpp

@@ -50,6 +50,7 @@ int finufft2d1(BIGINT nj,double* xj,double *yj,dcomplex* cj,int iflag,
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1; set_nf_type12(ms,opts,spopts,&nf1);
   BIGINT nf2; set_nf_type12(mt,opts,spopts,&nf2);
   cout << scientific << setprecision(15);  // for debug
@@ -134,6 +135,7 @@ int finufft2d2(BIGINT nj,double* xj,double *yj,dcomplex* cj,int iflag,double eps
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1; set_nf_type12(ms,opts,spopts,&nf1);
   BIGINT nf2; set_nf_type12(mt,opts,spopts,&nf2);
   cout << scientific << setprecision(15);  // for debug
@@ -225,6 +227,7 @@ int finufft2d3(BIGINT nj,double* xj,double* yj,dcomplex* cj,int iflag, double ep
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1,nf2;
   double X1,C1,S1,D1,h1,gam1,X2,C2,S2,D2,h2,gam2;
   cout << scientific << setprecision(15);  // for debug

src/finufft3d.cpp

@@ -52,6 +52,7 @@ int finufft3d1(BIGINT nj,double* xj,double *yj,double *zj,dcomplex* cj,int iflag
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1; set_nf_type12(ms,opts,spopts,&nf1);
   BIGINT nf2; set_nf_type12(mt,opts,spopts,&nf2);
   BIGINT nf3; set_nf_type12(mu,opts,spopts,&nf3);
@@ -145,6 +146,7 @@ int finufft3d2(BIGINT nj,double* xj,double *yj,double *zj,dcomplex* cj,
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1; set_nf_type12(ms,opts,spopts,&nf1);
   BIGINT nf2; set_nf_type12(mt,opts,spopts,&nf2);
   BIGINT nf3; set_nf_type12(mu,opts,spopts,&nf3);
@@ -246,6 +248,7 @@ int finufft3d3(BIGINT nj,double* xj,double* yj,double *zj, dcomplex* cj,
 {
   spread_opts spopts;
   int ier_set = setup_kernel(spopts,eps,opts.R);
+  if (ier_set) exit(ier_set);
   BIGINT nf1,nf2,nf3;
   double X1,C1,S1,D1,h1,gam1,X2,C2,S2,D2,h2,gam2,X3,C3,S3,D3,h3,gam3;
   cout << scientific << setprecision(15);  // for debug

src/twopispread.cpp

@@ -32,7 +32,7 @@
 int twopispread1d(BIGINT nf1,dcomplex *fw,BIGINT nj,double* xj,dcomplex* cj,
 		  spread_opts opts)
 {
-  double *dummy;   // note this should never be read from!
+  double *dummy=0;   // note this should never be read from!
   double *xjscal = (double*)malloc(sizeof(double)*nj);
   double s = nf1/(2*M_PI);
   for (BIGINT i=0;i<nj;++i)
@@ -46,7 +46,7 @@ int twopispread1d(BIGINT nf1,dcomplex *fw,BIGINT nj,double* xj,dcomplex* cj,
 int twopispread2d(long nf1,long nf2, dcomplex *fw,BIGINT nj,double* xj,
 		  double *yj,dcomplex* cj,spread_opts opts)
 {
-  double *dummy;
+  double *dummy=0;   // note this should never be read from!
   double *xjscal = (double*)malloc(sizeof(double)*nj);
   double *yjscal = (double*)malloc(sizeof(double)*nj);
   double s1 = nf1/(2*M_PI);