[nfft] 02/04: Imported Upstream version 3.3.0~beta1
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Fri Apr 17 15:36:35 UTC 2015
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to branch debian/experimental
in repository nfft.
commit 8b764041f63182ff5a5f6cc97b9fd7c2dba70d86
Author: Ghislain Antony Vaillant <ghisvail at gmail.com>
Date: Fri Apr 17 16:17:37 2015 +0100
Imported Upstream version 3.3.0~beta1
---
Makefile.am | 3 +
README | 37 +++--
.../fastgauss/{fastgauss.c => fastgauss.c.in} | 178 ++++++++++-----------
applications/fastsum/fastsum.c | 2 +-
applications/mri/mri2d/construct_data_2d.c | 1 -
applications/mri/mri2d/construct_data_inh_2d1d.c | 7 +-
applications/mri/mri2d/construct_data_inh_3d.c | 7 +-
applications/mri/mri2d/reconstruct_data_2d.c | 13 +-
applications/mri/mri2d/reconstruct_data_gridding.c | 1 -
applications/mri/mri2d/reconstruct_data_inh_2d1d.c | 19 ++-
applications/mri/mri2d/reconstruct_data_inh_3d.c | 19 ++-
.../mri/mri2d/reconstruct_data_inh_nnfft.c | 19 ++-
applications/mri/mri3d/construct_data_2d1d.c | 1 -
applications/mri/mri3d/construct_data_3d.c | 1 -
applications/mri/mri3d/reconstruct_data_2d1d.c | 1 -
applications/mri/mri3d/reconstruct_data_3d.c | 1 -
applications/mri/mri3d/reconstruct_data_gridding.c | 1 -
...{linogram_fft_test.c => linogram_fft_test.c.in} | 158 +++++++++---------
.../{mpolar_fft_test.c => mpolar_fft_test.c.in} | 162 +++++++++----------
.../{polar_fft_test.c => polar_fft_test.c.in} | 86 +++++-----
.../radon/{inverse_radon.c => inverse_radon.c.in} | 78 +++++----
applications/radon/{radon.c => radon.c.in} | 74 +++++----
configure.ac | 28 +++-
examples/fpt/simple_test.c | 9 +-
examples/nfct/{simple_test.c => simple_test.c.in} | 19 ++-
examples/nfft/{simple_test.c => simple_test.c.in} | 70 ++++----
...ple_test_threads.c => simple_test_threads.c.in} | 35 ++--
examples/nfsft/simple_test.c | 14 +-
examples/nfsft/simple_test_threads.c | 16 +-
examples/nfsoft/simple_test.c | 46 +++---
examples/nfst/{simple_test.c => simple_test.c.in} | 19 ++-
examples/nnfft/accuracy.c | 13 +-
examples/nnfft/simple_test.c | 22 +--
examples/nsfft/nsfft_test.c | 38 ++---
examples/nsfft/simple_test.c | 8 +-
examples/solver/{glacier.c => glacier.c.in} | 60 ++++---
.../solver/{simple_test.c => simple_test.c.in} | 33 ++--
include/Makefile.am | 2 +-
include/infft.h | 32 ----
include/nfft3.h | 129 +++++++++++----
include/nfft3mp.h | 104 ++++++++++++
kernel/nfct/nfct.c | 2 +-
kernel/nfft/nfft.c | 2 +-
kernel/nfsft/nfsft.c | 2 +-
kernel/nfsoft/nfsoft.c | 2 +-
kernel/nfst/nfst.c | 2 +-
kernel/nnfft/nnfft.c | 4 +-
kernel/util/time.c | 21 +++
matlab/nfft/Makefile.am | 3 +-
.../FG_PSI.m => nfft/NFFT_OMP_BLOCKWISE_ADJOINT.m} | 14 +-
matlab/nfft/nfft.m | 6 +-
matlab/nfft/simple_test.m | 4 +-
matlab/nfft/test_nfft1d.m | 2 +-
matlab/nfft/test_nfft2d.m | 2 +-
matlab/nfft/test_nfft3d.m | 2 +-
matlab/nnfft/FFTW_ESTIMATE.m | 25 ---
matlab/nnfft/FFTW_MEASURE.m | 25 ---
matlab/nnfft/FFT_OUT_OF_PLACE.m | 25 ---
matlab/nnfft/Makefile.am | 13 +-
matlab/nnfft/PRE_FG_PSI.m | 28 ----
matlab/nnfft/nnfft.m | 72 ++++-----
matlab/nnfft/out.txt | 32 ----
matlab/nnfft/test_nfft1d.m | 70 --------
matlab/nnfft/test_nnfft1d.m | 19 ++-
matlab/nnfft/test_nnfft2d.m | 55 +++++--
matlab/nnfft/test_nnfft2d_N215.m | 86 ++++++++++
nfft3.pc.in | 4 +-
support/Portfile | 59 +++++++
68 files changed, 1140 insertions(+), 1007 deletions(-)
diff --git a/Makefile.am b/Makefile.am
index aea6dae..6060d70 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -63,8 +63,11 @@ endif
EXTRA_DIST = bootstrap.sh doxygen.dox nfft3.pc.in doc
+if HAVE_NON_DOUBLE_PRECISION
nfft3 at PREC_SUFFIX@.pc: nfft3.pc
cp -f nfft3.pc nfft3 at PREC_SUFFIX@.pc
+endif
+
pkgconfigdir = $(libdir)/pkgconfig
pkgconfig_DATA = nfft3 at PREC_SUFFIX@.pc
diff --git a/README b/README
index 7bc5619..5ac8aa0 100644
--- a/README
+++ b/README
@@ -2,7 +2,7 @@ NFFT3 - Nonequispaced FFT, generalisations, inversion, and applications
Overview
--------
-NFFT3 is a software library written in C for computing nonequispaced fast Fourier
+NFFT3 is a software library written in C for computing nonequispaced fast Fourier
and related transformations. In detail, NFFT3 implements
1) The nonequispaced fast Fourier transform (NFFT)
@@ -14,7 +14,7 @@ and related transformations. In detail, NFFT3 implements
- to the hyperbolic cross (NSFFT)
- to real-valued data, i.e. (co)sine transforms, (NFCT, NFST)
- to the rotation group (NFSOFT)
- 3) Generalised inverses based on iterative methods, e.g. CGNR, CGNE
+ 3) Generalised inverses based on iterative methods, e.g. CGNR, CGNE
4) Applications in
- medical imaging
(i) magnetic resonance imaging
@@ -25,17 +25,24 @@ and related transformations. In detail, NFFT3 implements
(iii) zonal kernels
- polar FFT, discrete Radon transform, ridgelet transform
-For an introduction, please read the "NFFT 3.0 - Tutorial" first. It is
-available in the subdirectory doc/tutorial/tutorial.pdf. Detailed API
-documentation can be found in HTML format in /doc/api/html/index.html
-or in LaTeX format for self compilation in /doc/api/latex/refman.tex.
+Detailed API documentation can be found in HTML format in
+/doc/api/html/index.html.
+
For installation instructions, you can also refer to the file INSTALL
-in this directory.
+in this directory. The installation of NFFT3 follows the typical steps
+ ./configure
+ make
+ make install
+Optionally, some test programs can be built and run with
+ make check
+If the file configure is missing, please run
+ ./bootstrap.sh
+first.
-The most current general paper, and the one that we recommend if you wish
-to cite NFFT, is: The paper by Keiner, J., Kunis, S., and Potts, D.
-''Using NFFT 3 - a software library for various nonequispaced fast Fourier transforms''
-ACM Trans. Math. Software,36, Article 19, 1-30, 2009
+The most current general paper, and the one that we recommend if you wish
+to cite NFFT, is: The paper by Keiner, J., Kunis, S., and Potts, D.
+''Using NFFT 3 - a software library for various nonequispaced fast Fourier transforms''
+ACM Trans. Math. Software,36, Article 19, 1-30, 2009
Directory structure
-------------------
3rdparty (dir) Third-party source code
@@ -46,7 +53,7 @@ bootstrap.sh Bootstrap shell script that call Autoconf and friends
ChangeLog A short version history
config.guess Used by configure script
config.sub Used by configure script
-configure Configure script
+configure Configure script (created by calling ./bootstrap.sh)
configure.in Autoconf configure script template
CONVENTIONS Internal coding conventions
COPYING Information about redistributing NFFT3
@@ -61,11 +68,12 @@ ltmain.sh Used by configure script
Makefile.am Automake Makefile template
Makefile.in Makefile template generated from Makefile.am,
processed by configure script
+matlab (dir) Matlab MEX interfaces for nfft, nfsft, nfsoft, nfft
missing Used by configure script
NEWS New and noteworthy
README This file
+tests (dir) CUnit tests
TODO Work to be done
-util (dir) Source code for auxilliary routines
Feedback
--------
@@ -98,7 +106,7 @@ NFFT3 is free software. You can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version. If not stated otherwise, this applies to all files contained in this
-package and its sub-directories.
+package and its sub-directories.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -108,4 +116,3 @@ 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, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-
diff --git a/applications/fastgauss/fastgauss.c b/applications/fastgauss/fastgauss.c.in
similarity index 70%
rename from applications/fastgauss/fastgauss.c
rename to applications/fastgauss/fastgauss.c.in
index 8f1063a..a77366c 100644
--- a/applications/fastgauss/fastgauss.c
+++ b/applications/fastgauss/fastgauss.c.in
@@ -24,17 +24,14 @@
* \{
*/
-#include "config.h"
-
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
-#ifdef HAVE_COMPLEX_H
- #include <complex.h>
-#endif
+#include <complex.h>
+
+#define @NFFT_PRECISION_MACRO@
-#include "nfft3.h"
-#include "infft.h"
+#include "nfft3mp.h"
/**
* If this flag is set, the whole matrix is precomputed and stored for the
@@ -73,23 +70,23 @@ typedef struct
int N; /**< number of source nodes */
int M; /**< number of target nodes */
- C *alpha; /**< source coefficients */
- C *f; /**< target evaluations */
+ NFFT_C *alpha; /**< source coefficients */
+ NFFT_C *f; /**< target evaluations */
unsigned flags; /**< flags for precomputation and
approximation type */
- C sigma; /**< parameter of the Gaussian */
+ NFFT_C sigma; /**< parameter of the Gaussian */
- R *x; /**< source nodes in \f$[-1/4,1/4]\f$*/
- R *y; /**< target nodes in \f$[-1/4,1/4]\f$*/
+ NFFT_R *x; /**< source nodes in \f$[-1/4,1/4]\f$*/
+ NFFT_R *y; /**< target nodes in \f$[-1/4,1/4]\f$*/
- C *pre_cexp; /**< precomputed values for dgt */
+ NFFT_C *pre_cexp; /**< precomputed values for dgt */
int n; /**< expansion degree */
- R p; /**< period, at least 1 */
+ NFFT_R p; /**< period, at least 1 */
- C *b; /**< expansion coefficients */
+ NFFT_C *b; /**< expansion coefficients */
NFFT(plan) *nplan1; /**< source nfft plan */
NFFT(plan) *nplan2; /**< target nfft plan */
@@ -105,10 +102,10 @@ typedef struct
*/
static void dgt_trafo(fgt_plan *ths)
{
- INT j, k, l;
+ NFFT_INT j, k, l;
for (j = 0; j < ths->M; j++)
- ths->f[j] = K(0.0);
+ ths->f[j] = NFFT_K(0.0);
if (ths->flags & DGT_PRE_CEXP)
for (j = 0, l = 0; j < ths->M; j++)
@@ -118,7 +115,7 @@ static void dgt_trafo(fgt_plan *ths)
for (j = 0; j < ths->M; j++)
for (k = 0; k < ths->N; k++)
ths->f[j] += ths->alpha[k]
- * CEXP(
+ * NFFT_M(cexp)(
-ths->sigma * (ths->y[j] - ths->x[k])
* (ths->y[j] - ths->x[k]));
}
@@ -132,7 +129,7 @@ static void dgt_trafo(fgt_plan *ths)
*/
static void fgt_trafo(fgt_plan *ths)
{
- INT l;
+ NFFT_INT l;
if (ths->flags & FGT_NDFT)
{
@@ -168,7 +165,7 @@ static void fgt_trafo(fgt_plan *ths)
*
* \author Stefan Kunis
*/
-static void fgt_init_guru(fgt_plan *ths, int N, int M, C sigma, int n, R p, int m,
+static void fgt_init_guru(fgt_plan *ths, int N, int M, NFFT_C sigma, int n, NFFT_R p, int m,
unsigned flags)
{
int j, n_fftw;
@@ -179,15 +176,15 @@ static void fgt_init_guru(fgt_plan *ths, int N, int M, C sigma, int n, R p, int
ths->sigma = sigma;
ths->flags = flags;
- ths->x = (R*) NFFT(malloc)((size_t)(ths->N) * sizeof(R));
- ths->y = (R*) NFFT(malloc)((size_t)(ths->M) * sizeof(R));
- ths->alpha = (C*) NFFT(malloc)((size_t)(ths->N) * sizeof(C));
- ths->f = (C*) NFFT(malloc)((size_t)(ths->M) * sizeof(C));
+ ths->x = (NFFT_R*) NFFT(malloc)((size_t)(ths->N) * sizeof(NFFT_R));
+ ths->y = (NFFT_R*) NFFT(malloc)((size_t)(ths->M) * sizeof(NFFT_R));
+ ths->alpha = (NFFT_C*) NFFT(malloc)((size_t)(ths->N) * sizeof(NFFT_C));
+ ths->f = (NFFT_C*) NFFT(malloc)((size_t)(ths->M) * sizeof(NFFT_C));
ths->n = n;
ths->p = p;
- ths->b = (C*) NFFT(malloc)((size_t)(ths->n) * sizeof(C));
+ ths->b = (NFFT_C*) NFFT(malloc)((size_t)(ths->n) * sizeof(NFFT_C));
ths->nplan1 = (NFFT(plan)*) NFFT(malloc)(sizeof(NFFT(plan)));
ths->nplan2 = (NFFT(plan)*) NFFT(malloc)(sizeof(NFFT(plan)));
@@ -211,9 +208,9 @@ static void fgt_init_guru(fgt_plan *ths, int N, int M, C sigma, int n, R p, int
FFTW_MEASURE);
for (j = 0; j < ths->n; j++)
- ths->b[j] = CEXP(
- -ths->p * ths->p * ths->sigma * ((R)(j) - (R)(ths->n) / K(2.0)) * ((R)(j) - (R)(ths->n) / K(2.0))
- / ((R) (ths->n * ths->n))) / ((R)(ths->n));
+ ths->b[j] = NFFT_M(cexp)(
+ -ths->p * ths->p * ths->sigma * ((NFFT_R)(j) - (NFFT_R)(ths->n) / NFFT_K(2.0)) * ((NFFT_R)(j) - (NFFT_R)(ths->n) / NFFT_K(2.0))
+ / ((NFFT_R) (ths->n * ths->n))) / ((NFFT_R)(ths->n));
NFFT(fftshift_complex_int)(ths->b, 1, &ths->n);
FFTW(execute)(fplan);
@@ -224,9 +221,9 @@ static void fgt_init_guru(fgt_plan *ths, int N, int M, C sigma, int n, R p, int
else
{
for (j = 0; j < ths->n; j++)
- ths->b[j] = K(1.0) / ths->p * CSQRT(KPI / ths->sigma)
- * CEXP(
- -KPI * KPI * ((R)(j) - (R)(ths->n) / K(2.0)) * ((R)(j) - (R)(ths->n) / K(2.0))
+ ths->b[j] = NFFT_K(1.0) / ths->p * NFFT_M(csqrt)(NFFT_KPI / ths->sigma)
+ * NFFT_M(cexp)(
+ -NFFT_KPI * NFFT_KPI * ((NFFT_R)(j) - (NFFT_R)(ths->n) / NFFT_K(2.0)) * ((NFFT_R)(j) - (NFFT_R)(ths->n) / NFFT_K(2.0))
/ (ths->p * ths->p * ths->sigma));
}
}
@@ -242,19 +239,19 @@ static void fgt_init_guru(fgt_plan *ths, int N, int M, C sigma, int n, R p, int
*
* \author Stefan Kunis
*/
-static void fgt_init(fgt_plan *ths, int N, int M, C sigma, R eps)
+static void fgt_init(fgt_plan *ths, int N, int M, NFFT_C sigma, NFFT_R eps)
{
- R p;
+ NFFT_R p;
int n;
- p = K(0.5) + SQRT(-LOG(eps) / CREAL(sigma));
+ p = NFFT_K(0.5) + NFFT_M(sqrt)(-NFFT_M(log)(eps) / NFFT_M(creal)(sigma));
- if (p < K(1.0))
- p = K(1.0);
+ if (p < NFFT_K(1.0))
+ p = NFFT_K(1.0);
- n = (int)(2 * (LRINT(CEIL(p * CABS(sigma) / KPI * SQRT(-LOG(eps) / CREAL(sigma))))));
+ n = (int)(2 * (NFFT_M(lrint)(NFFT_M(ceil)(p * NFFT_M(cabs)(sigma) / NFFT_KPI * NFFT_M(sqrt)(-NFFT_M(log)(eps) / NFFT_M(creal)(sigma))))));
- fgt_init_guru(ths, N, M, sigma, n, p, 7, N * M <= ((INT) (1U << 20)) ? DGT_PRE_CEXP : 0);
+ fgt_init_guru(ths, N, M, sigma, n, p, 7, N * M <= ((NFFT_INT) (1U << 20)) ? DGT_PRE_CEXP : 0);
}
/**
@@ -269,11 +266,11 @@ static void fgt_init_node_dependent(fgt_plan *ths)
if (ths->flags & DGT_PRE_CEXP)
{
- ths->pre_cexp = (C*) NFFT(malloc)((size_t)(ths->M * ths->N) * sizeof(C));
+ ths->pre_cexp = (NFFT_C*) NFFT(malloc)((size_t)(ths->M * ths->N) * sizeof(NFFT_C));
for (j = 0, l = 0; j < ths->M; j++)
for (k = 0; k < ths->N; k++, l++)
- ths->pre_cexp[l] = CEXP(
+ ths->pre_cexp[l] = NFFT_M(cexp)(
-ths->sigma * (ths->y[j] - ths->x[k]) * (ths->y[j] - ths->x[k]));
}
@@ -319,17 +316,17 @@ static void fgt_finalize(fgt_plan *ths)
*/
static void fgt_test_init_rand(fgt_plan *ths)
{
- INT j, k;
+ NFFT_INT j, k;
for (k = 0; k < ths->N; k++)
- ths->x[k] = NFFT(drand48)() / K(2.0) - K(1.0) / K(4.0);
+ ths->x[k] = NFFT(drand48)() / NFFT_K(2.0) - NFFT_K(1.0) / NFFT_K(4.0);
for (j = 0; j < ths->M; j++)
- ths->y[j] = NFFT(drand48)() / K(2.0) - K(1.0) / K(4.0);
+ ths->y[j] = NFFT(drand48)() / NFFT_K(2.0) - NFFT_K(1.0) / NFFT_K(4.0);
for (k = 0; k < ths->N; k++)
- ths->alpha[k] = NFFT(drand48)() - K(1.0) / K(2.0)
- + II * (NFFT(drand48)() - K(1.0) / K(2.0));
+ ths->alpha[k] = NFFT(drand48)() - NFFT_K(1.0) / NFFT_K(2.0)
+ + _Complex_I * (NFFT(drand48)() - NFFT_K(1.0) / NFFT_K(2.0));
}
/**
@@ -340,27 +337,28 @@ static void fgt_test_init_rand(fgt_plan *ths)
*
* \author Stefan Kunis
*/
-static R fgt_test_measure_time(fgt_plan *ths, unsigned dgt)
+static NFFT_R fgt_test_measure_time(fgt_plan *ths, unsigned dgt)
{
int r;
- ticks t0, t1;
- R t_out;
- R tau = K(0.01);
+ NFFT_R t0, t1, time_diff;
+ NFFT_R t_out;
+ NFFT_R tau = NFFT_K(0.01);
- t_out = K(0.0);
+ t_out = NFFT_K(0.0);
r = 0;
while (t_out < tau)
{
r++;
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
if (dgt)
dgt_trafo(ths);
else
fgt_trafo(ths);
- t1 = getticks();
- t_out += NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
+ time_diff = t1 - t0;
+ t_out += time_diff;
}
- t_out /= (R)(r);
+ t_out /= (NFFT_R)(r);
return t_out;
}
@@ -374,26 +372,26 @@ static R fgt_test_measure_time(fgt_plan *ths, unsigned dgt)
*
* \author Stefan Kunis
*/
-static void fgt_test_simple(int N, int M, C sigma, R eps)
+static void fgt_test_simple(int N, int M, NFFT_C sigma, NFFT_R eps)
{
fgt_plan my_plan;
- C *swap_dgt;
+ NFFT_C *swap_dgt;
fgt_init(&my_plan, N, M, sigma, eps);
- swap_dgt = (C*) NFFT(malloc)((size_t)(my_plan.M) * sizeof(C));
+ swap_dgt = (NFFT_C*) NFFT(malloc)((size_t)(my_plan.M) * sizeof(NFFT_C));
fgt_test_init_rand(&my_plan);
fgt_init_node_dependent(&my_plan);
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
dgt_trafo(&my_plan);
NFFT(vpr_complex)(my_plan.f, my_plan.M, "discrete gauss transform");
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
fgt_trafo(&my_plan);
NFFT(vpr_complex)(my_plan.f, my_plan.M, "fast gauss transform");
- printf("\n relative error: %1.3" __FES__ "\n",
+ printf("\n relative error: %1.3" NFFT__FES__ "\n",
NFFT(error_l_infty_1_complex)(swap_dgt, my_plan.f, my_plan.M,
my_plan.alpha, my_plan.N));
@@ -413,23 +411,23 @@ static void fgt_test_simple(int N, int M, C sigma, R eps)
static void fgt_test_andersson(void)
{
fgt_plan my_plan;
- C *swap_dgt;
+ NFFT_C *swap_dgt;
int N;
- C sigma = K(4.0) * (K(138.0) + II * K(100.0));
+ NFFT_C sigma = NFFT_K(4.0) * (NFFT_K(138.0) + _Complex_I * NFFT_K(100.0));
int n = 128;
int N_dgt_pre_exp = (int) (1U << 11);
int N_dgt = (int) (1U << 19);
- printf("n=%d, sigma=%1.3" __FES__ "+i%1.3" __FES__ "\n", n, CREAL(sigma), CIMAG(sigma));
+ printf("n=%d, sigma=%1.3" NFFT__FES__ "+i%1.3" NFFT__FES__ "\n", n, NFFT_M(creal)(sigma), NFFT_M(cimag)(sigma));
- for (N = ((INT) (1U << 6)); N < ((INT) (1U << 22)); N = N << 1)
+ for (N = ((NFFT_INT) (1U << 6)); N < ((NFFT_INT) (1U << 22)); N = N << 1)
{
printf("$%d$\t & ", N);
fgt_init_guru(&my_plan, N, N, sigma, n, 1, 7, N < N_dgt_pre_exp ? DGT_PRE_CEXP : 0);
- swap_dgt = (C*) NFFT(malloc)((size_t)(my_plan.M) * sizeof(C));
+ swap_dgt = (NFFT_C*) NFFT(malloc)((size_t)(my_plan.M) * sizeof(NFFT_C));
fgt_test_init_rand(&my_plan);
@@ -437,30 +435,30 @@ static void fgt_test_andersson(void)
if (N < N_dgt)
{
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
if (N < N_dgt_pre_exp)
my_plan.flags ^= DGT_PRE_CEXP;
- printf("$%1.1" __FES__ "$\t & ", fgt_test_measure_time(&my_plan, 1));
+ printf("$%1.1" NFFT__FES__ "$\t & ", fgt_test_measure_time(&my_plan, 1));
if (N < N_dgt_pre_exp)
my_plan.flags ^= DGT_PRE_CEXP;
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
}
else
printf("\t\t & ");
if (N < N_dgt_pre_exp)
- printf("$%1.1" __FES__ "$\t & ", fgt_test_measure_time(&my_plan, 1));
+ printf("$%1.1" NFFT__FES__ "$\t & ", fgt_test_measure_time(&my_plan, 1));
else
printf("\t\t & ");
my_plan.flags ^= FGT_NDFT;
- printf("$%1.1" __FES__ "$\t & ", fgt_test_measure_time(&my_plan, 0));
+ printf("$%1.1" NFFT__FES__ "$\t & ", fgt_test_measure_time(&my_plan, 0));
my_plan.flags ^= FGT_NDFT;
- printf("$%1.1" __FES__ "$\t & ", fgt_test_measure_time(&my_plan, 0));
+ printf("$%1.1" NFFT__FES__ "$\t & ", fgt_test_measure_time(&my_plan, 0));
- printf("$%1.1" __FES__ "$\t \\\\ \n",
+ printf("$%1.1" NFFT__FES__ "$\t \\\\ \n",
NFFT(error_l_infty_1_complex)(swap_dgt, my_plan.f, my_plan.M, my_plan.alpha, my_plan.N));
fflush(stdout);
@@ -479,19 +477,19 @@ static void fgt_test_andersson(void)
static void fgt_test_error(void)
{
fgt_plan my_plan;
- C *swap_dgt;
+ NFFT_C *swap_dgt;
int n, mi;
- C sigma = K(4.0) * (K(138.0) + II * K(100.0));
+ NFFT_C sigma = NFFT_K(4.0) * (NFFT_K(138.0) + _Complex_I * NFFT_K(100.0));
int N = 1000;
int M = 1000;
int m[2] = { 7, 3 };
- printf("N=%d;\tM=%d;\nsigma=%1.3" __FES__ "+i*%1.3" __FES__ ";\n", N, M, CREAL(sigma),
- CIMAG(sigma));
+ printf("N=%d;\tM=%d;\nsigma=%1.3" NFFT__FES__ "+i*%1.3" NFFT__FES__ ";\n", N, M, NFFT_M(creal)(sigma),
+ NFFT_M(cimag)(sigma));
printf("error=[\n");
- swap_dgt = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
+ swap_dgt = (NFFT_C*) NFFT(malloc)((size_t)(M) * sizeof(NFFT_C));
for (n = 8; n <= 128; n += 4)
{
@@ -502,13 +500,13 @@ static void fgt_test_error(void)
fgt_test_init_rand(&my_plan);
fgt_init_node_dependent(&my_plan);
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
dgt_trafo(&my_plan);
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
fgt_trafo(&my_plan);
- printf("%1.3" __FES__ "\t",
+ printf("%1.3" NFFT__FES__ "\t",
NFFT(error_l_infty_1_complex)(swap_dgt, my_plan.f, my_plan.M, my_plan.alpha, my_plan.N));
fflush(stdout);
@@ -531,19 +529,19 @@ static void fgt_test_error(void)
static void fgt_test_error_p(void)
{
fgt_plan my_plan;
- C *swap_dgt;
+ NFFT_C *swap_dgt;
int n, pi;
- C sigma = K(20.0) + K(40.0) * II;
+ NFFT_C sigma = NFFT_K(20.0) + NFFT_K(40.0) * _Complex_I;
int N = 1000;
int M = 1000;
- R p[3] = {K(1.0), K(1.5), K(2.0)};
+ NFFT_R p[3] = {NFFT_K(1.0), NFFT_K(1.5), NFFT_K(2.0)};
- printf("N=%d;\tM=%d;\nsigma=%1.3" __FES__ "+i*%1.3" __FES__ ";\n", N, M, CREAL(sigma),
- CIMAG(sigma));
+ printf("N=%d;\tM=%d;\nsigma=%1.3" NFFT__FES__ "+i*%1.3" NFFT__FES__ ";\n", N, M, NFFT_M(creal)(sigma),
+ NFFT_M(cimag)(sigma));
printf("error=[\n");
- swap_dgt = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
+ swap_dgt = (NFFT_C*) NFFT(malloc)((size_t)(M) * sizeof(NFFT_C));
for (n = 8; n <= 128; n += 4)
{
@@ -554,13 +552,13 @@ static void fgt_test_error_p(void)
fgt_test_init_rand(&my_plan);
fgt_init_node_dependent(&my_plan);
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
dgt_trafo(&my_plan);
- CSWAP(swap_dgt, my_plan.f);
+ NFFT_CSWAP(swap_dgt, my_plan.f);
fgt_trafo(&my_plan);
- printf("%1.3" __FES__ "\t",
+ printf("%1.3" NFFT__FES__ "\t",
NFFT(error_l_infty_1_complex)(swap_dgt, my_plan.f, my_plan.M, my_plan.alpha, my_plan.N));
fflush(stdout);
@@ -594,7 +592,7 @@ int main(int argc, char **argv)
}
if (atoi(argv[1]) == 0)
- fgt_test_simple(10, 10, K(5.0) + K(3.0) * II, K(0.001));
+ fgt_test_simple(10, 10, NFFT_K(5.0) + NFFT_K(3.0) * _Complex_I, NFFT_K(0.001));
if (atoi(argv[1]) == 1)
fgt_test_andersson();
diff --git a/applications/fastsum/fastsum.c b/applications/fastsum/fastsum.c
index 5aaf49e..ea21862 100644
--- a/applications/fastsum/fastsum.c
+++ b/applications/fastsum/fastsum.c
@@ -36,7 +36,7 @@
#include "fastsum.h"
#include "infft.h"
-/** Required for test if (ths->k == one_over_x) */
+// Required for test if (ths->k == one_over_x)
#include "kernels.h"
/**
diff --git a/applications/mri/mri2d/construct_data_2d.c b/applications/mri/mri2d/construct_data_2d.c
index 0abc733..ab1546f 100644
--- a/applications/mri/mri2d/construct_data_2d.c
+++ b/applications/mri/mri2d/construct_data_2d.c
@@ -27,7 +27,6 @@
#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri2d_construct_data_2d construct_data_2d
diff --git a/applications/mri/mri2d/construct_data_inh_2d1d.c b/applications/mri/mri2d/construct_data_inh_2d1d.c
index 6126818..deccda2 100644
--- a/applications/mri/mri2d/construct_data_inh_2d1d.c
+++ b/applications/mri/mri2d/construct_data_inh_2d1d.c
@@ -27,7 +27,10 @@
#endif
#include "nfft3.h"
-#include "infft.h"
+
+#ifndef MAX
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#endif
/**
* \defgroup applications_mri2d_construct_data_inh_2d1d construct_data__inh_2d1d
@@ -123,7 +126,7 @@ static void construct(char * file, int N, int M)
for(j=0;j<N*N;j++)
{
fscanf(fi,"%le ",&real);
- my_plan.f_hat[j] = real*cexp(2.0*_Complex_I*KPI*Ts*my_plan.w[j]*W);
+ my_plan.f_hat[j] = real*cexp(2.0*_Complex_I*M_PI*Ts*my_plan.w[j]*W);
}
if(my_plan.plan.flags & PRE_PSI)
diff --git a/applications/mri/mri2d/construct_data_inh_3d.c b/applications/mri/mri2d/construct_data_inh_3d.c
index 2025eb0..b0dc66a 100644
--- a/applications/mri/mri2d/construct_data_inh_3d.c
+++ b/applications/mri/mri2d/construct_data_inh_3d.c
@@ -27,7 +27,10 @@
#endif
#include "nfft3.h"
-#include "infft.h"
+
+#ifndef MAX
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#endif
/**
* \defgroup applications_mri2d_construct_data_inh_3d construct_data_inh_3d
@@ -122,7 +125,7 @@ static void construct(char * file, int N, int M)
for(j=0;j<N*N;j++)
{
fscanf(fi,"%le ",&real);
- my_plan.f_hat[j] = real*cexp(2.0*_Complex_I*KPI*Ts*my_plan.w[j]*W);
+ my_plan.f_hat[j] = real*cexp(2.0*_Complex_I*M_PI*Ts*my_plan.w[j]*W);
}
if(my_plan.plan.flags & PRE_PSI)
diff --git a/applications/mri/mri2d/reconstruct_data_2d.c b/applications/mri/mri2d/reconstruct_data_2d.c
index 44f9557..d60ef2b 100644
--- a/applications/mri/mri2d/reconstruct_data_2d.c
+++ b/applications/mri/mri2d/reconstruct_data_2d.c
@@ -17,16 +17,11 @@
*/
/* $Id$ */
-#include "config.h"
-
#include <math.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri2d_reconstruct_data_2d reconstruct_data_2d
@@ -40,7 +35,7 @@
static void reconstruct(char* filename,int N,int M,int iteration, int weight)
{
int j,k,l; /* some variables */
- ticks t0, t1;
+ double t0, t1;
double real,imag,t; /* to read the real and imag part of a complex number */
nfft_plan my_plan; /* plan for the two dimensional nfft */
solver_plan_complex my_iplan; /* plan for the two dimensional infft */
@@ -124,7 +119,7 @@ static void reconstruct(char* filename,int N,int M,int iteration, int weight)
for(k=0;k<my_plan.N_total;k++)
my_iplan.f_hat_iter[k]=0.0;
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
/* inverse trafo */
solver_before_loop_complex(&my_iplan);
@@ -139,8 +134,8 @@ static void reconstruct(char* filename,int N,int M,int iteration, int weight)
}
- t1 = getticks();
- t=nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t=t1-t0;
fout_real=fopen("output_real.dat","w");
fout_imag=fopen("output_imag.dat","w");
diff --git a/applications/mri/mri2d/reconstruct_data_gridding.c b/applications/mri/mri2d/reconstruct_data_gridding.c
index 0f74548..671909b 100644
--- a/applications/mri/mri2d/reconstruct_data_gridding.c
+++ b/applications/mri/mri2d/reconstruct_data_gridding.c
@@ -26,7 +26,6 @@
#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri2d_construct_data_gridding construct_data_gridding
diff --git a/applications/mri/mri2d/reconstruct_data_inh_2d1d.c b/applications/mri/mri2d/reconstruct_data_inh_2d1d.c
index c273667..d91cf2e 100644
--- a/applications/mri/mri2d/reconstruct_data_inh_2d1d.c
+++ b/applications/mri/mri2d/reconstruct_data_inh_2d1d.c
@@ -17,17 +17,16 @@
*/
/* $Id$ */
-#include "config.h"
-
#include <stdlib.h>
#include <math.h>
#include <limits.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
+
+#ifndef MAX
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#endif
/**
* \defgroup applications_mri2d_reconstruct_data_inh_2d1d reconstruct_data__inh_2d1d
@@ -39,7 +38,7 @@ static void reconstruct(char* filename,int N,int M,int iteration , int weight)
{
int j,k,l;
double time,min_time,max_time,min_inh,max_inh;
- ticks t0, t1;
+ double t0, t1;
double t,real,imag;
double w,epsilon=0.0000003; /* epsilon is a the break criterium for
the iteration */;
@@ -177,7 +176,7 @@ static void reconstruct(char* filename,int N,int M,int iteration , int weight)
my_iplan.f_hat_iter[j]=0.0;
}
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
/* inverse trafo */
solver_before_loop_complex(&my_iplan);
@@ -191,15 +190,15 @@ static void reconstruct(char* filename,int N,int M,int iteration , int weight)
solver_loop_one_step_complex(&my_iplan);
}
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = t1-t0;
fout_real=fopen("output_real.dat","w");
fout_imag=fopen("output_imag.dat","w");
for (j=0;j<N*N;j++) {
/* Verschiebung wieder herausrechnen */
- my_iplan.f_hat_iter[j]*=cexp(-2.0*_Complex_I*KPI*Ts*my_plan.w[j]*W);
+ my_iplan.f_hat_iter[j]*=cexp(-2.0*_Complex_I*M_PI*Ts*my_plan.w[j]*W);
fprintf(fout_real,"%le ",creal(my_iplan.f_hat_iter[j]));
fprintf(fout_imag,"%le ",cimag(my_iplan.f_hat_iter[j]));
diff --git a/applications/mri/mri2d/reconstruct_data_inh_3d.c b/applications/mri/mri2d/reconstruct_data_inh_3d.c
index decaad6..dee0f35 100644
--- a/applications/mri/mri2d/reconstruct_data_inh_3d.c
+++ b/applications/mri/mri2d/reconstruct_data_inh_3d.c
@@ -17,17 +17,16 @@
*/
/* $Id$ */
-#include "config.h"
-
#include <stdlib.h>
#include <math.h>
#include <limits.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
+
+#ifndef MAX
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#endif
/**
* \defgroup applications_mri2d_reconstruct_data_inh_3d reconstruct_data_inh_3d
@@ -38,7 +37,7 @@
static void reconstruct(char* filename,int N,int M,int iteration , int weight)
{
int j,k,l;
- ticks t0, t1;
+ double t0, t1;
double time,min_time,max_time,min_inh,max_inh;
double t,real,imag;
double w,epsilon=0.0000003; /* epsilon is a the break criterium for
@@ -171,7 +170,7 @@ static void reconstruct(char* filename,int N,int M,int iteration , int weight)
my_iplan.f_hat_iter[j]=0.0;
}
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
/* inverse trafo */
solver_before_loop_complex(&my_iplan);
@@ -185,15 +184,15 @@ static void reconstruct(char* filename,int N,int M,int iteration , int weight)
solver_loop_one_step_complex(&my_iplan);
}
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = t1-t0;
fout_real=fopen("output_real.dat","w");
fout_imag=fopen("output_imag.dat","w");
for (j=0;j<N*N;j++) {
/* Verschiebung wieder herausrechnen */
- my_iplan.f_hat_iter[j]*=cexp(-2.0*_Complex_I*KPI*Ts*my_plan.w[j]*W);
+ my_iplan.f_hat_iter[j]*=cexp(-2.0*_Complex_I*M_PI*Ts*my_plan.w[j]*W);
fprintf(fout_real,"%le ",creal(my_iplan.f_hat_iter[j]));
fprintf(fout_imag,"%le ",cimag(my_iplan.f_hat_iter[j]));
diff --git a/applications/mri/mri2d/reconstruct_data_inh_nnfft.c b/applications/mri/mri2d/reconstruct_data_inh_nnfft.c
index e16e5ca..846866e 100644
--- a/applications/mri/mri2d/reconstruct_data_inh_nnfft.c
+++ b/applications/mri/mri2d/reconstruct_data_inh_nnfft.c
@@ -17,17 +17,16 @@
*/
/* $Id$ */
-#include "config.h"
-
#include <stdlib.h>
#include <math.h>
#include <limits.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
+
+#ifndef MAX
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#endif
/**
* \defgroup applications_mri2d_construct_data_inh_nnfft construct_data_inh_nnfft
@@ -49,7 +48,7 @@ static void reconstruct(char* filename,int N,int M,int iteration, int weight)
FILE* fout_real; /* output file */
FILE* fout_imag; /* output file */
int my_N[3],my_n[3]; /* to init the nfft */
- ticks t0, t1;
+ double t0, t1;
double t,epsilon=0.0000003; /* epsilon is a the break criterium for
the iteration */
unsigned infft_flags = CGNR | PRECOMPUTE_DAMP; /* flags for the infft*/
@@ -184,7 +183,7 @@ static void reconstruct(char* filename,int N,int M,int iteration, int weight)
my_iplan.f_hat_iter[k]=0.0;
}
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
/* inverse trafo */
solver_before_loop_complex(&my_iplan);
@@ -198,15 +197,15 @@ static void reconstruct(char* filename,int N,int M,int iteration, int weight)
solver_loop_one_step_complex(&my_iplan);
}
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = t1-t0;
fout_real=fopen("output_real.dat","w");
fout_imag=fopen("output_imag.dat","w");
for(k=0;k<my_plan.N_total;k++) {
- my_iplan.f_hat_iter[k]*=cexp(2.0*_Complex_I*KPI*Ts*w[k]);
+ my_iplan.f_hat_iter[k]*=cexp(2.0*_Complex_I*M_PI*Ts*w[k]);
fprintf(fout_real,"%le ", creal(my_iplan.f_hat_iter[k]));
fprintf(fout_imag,"%le ", cimag(my_iplan.f_hat_iter[k]));
diff --git a/applications/mri/mri3d/construct_data_2d1d.c b/applications/mri/mri3d/construct_data_2d1d.c
index 95562ff..0e981f4 100644
--- a/applications/mri/mri3d/construct_data_2d1d.c
+++ b/applications/mri/mri3d/construct_data_2d1d.c
@@ -26,7 +26,6 @@
#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri3d_construct_data_1d2d construct_data_1d2d
diff --git a/applications/mri/mri3d/construct_data_3d.c b/applications/mri/mri3d/construct_data_3d.c
index 69d4624..a8ecb4e 100644
--- a/applications/mri/mri3d/construct_data_3d.c
+++ b/applications/mri/mri3d/construct_data_3d.c
@@ -26,7 +26,6 @@
#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri3d_construct_data_3d construct_data_3d
diff --git a/applications/mri/mri3d/reconstruct_data_2d1d.c b/applications/mri/mri3d/reconstruct_data_2d1d.c
index 2b77f8f..e0e98ff 100644
--- a/applications/mri/mri3d/reconstruct_data_2d1d.c
+++ b/applications/mri/mri3d/reconstruct_data_2d1d.c
@@ -26,7 +26,6 @@
#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri3d_reconstruct_data_1d2d reconstruct_data_1d2d
diff --git a/applications/mri/mri3d/reconstruct_data_3d.c b/applications/mri/mri3d/reconstruct_data_3d.c
index 0904498..b705cad 100644
--- a/applications/mri/mri3d/reconstruct_data_3d.c
+++ b/applications/mri/mri3d/reconstruct_data_3d.c
@@ -26,7 +26,6 @@
#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri3d_reconstruct_data_3d reconstruct_data_3d
diff --git a/applications/mri/mri3d/reconstruct_data_gridding.c b/applications/mri/mri3d/reconstruct_data_gridding.c
index d6c5e8a..f5021d2 100644
--- a/applications/mri/mri3d/reconstruct_data_gridding.c
+++ b/applications/mri/mri3d/reconstruct_data_gridding.c
@@ -26,7 +26,6 @@
#endif
#include "nfft3.h"
-#include "infft.h"
/**
* \defgroup applications_mri3d_reconstruct_data_gridding reconstruct_data_gridding
diff --git a/applications/polarFFT/linogram_fft_test.c b/applications/polarFFT/linogram_fft_test.c.in
similarity index 70%
rename from applications/polarFFT/linogram_fft_test.c
rename to applications/polarFFT/linogram_fft_test.c.in
index 309189e..72fad97 100644
--- a/applications/polarFFT/linogram_fft_test.c
+++ b/applications/polarFFT/linogram_fft_test.c.in
@@ -26,16 +26,14 @@
* \author Markus Fenn
* \date 2006
*/
-#include "config.h"
#include <math.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
/**
* \defgroup applications_polarFFT_linogramm linogram_fft_test
@@ -43,15 +41,15 @@
* \{
*/
-R GLOBAL_elapsed_time;
+NFFT_R GLOBAL_elapsed_time;
/** Generates the points x with weights w
* for the linogram grid with T slopes and R offsets.
*/
-static int linogram_grid(int T, int rr, R *x, R *w)
+static int linogram_grid(int T, int rr, NFFT_R *x, NFFT_R *w)
{
int t, r;
- R W = (R) T * (((R) rr / K(2.0)) * ((R) rr / K(2.0)) + K(1.0) / K(4.0));
+ NFFT_R W = (NFFT_R) T * (((NFFT_R) rr / NFFT_K(2.0)) * ((NFFT_R) rr / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
for (t = -T / 2; t < T / 2; t++)
{
@@ -59,20 +57,20 @@ static int linogram_grid(int T, int rr, R *x, R *w)
{
if (t < 0)
{
- x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = (R) (r) / (R)(rr);
- x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = K(4.0) * ((R)(t) + (R)(T) / K(4.0)) / (R)(T)
- * (R)(r) / (R)(rr);
+ x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = (NFFT_R) (r) / (NFFT_R)(rr);
+ x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = NFFT_K(4.0) * ((NFFT_R)(t) + (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
+ * (NFFT_R)(r) / (NFFT_R)(rr);
}
else
{
- x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = -K(4.0) * ((R)(t) - (R)(T) / K(4.0)) / (R)(T)
- * (R)(r) / (R)(rr);
- x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = (R) r / (R)(rr);
+ x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = -NFFT_K(4.0) * ((NFFT_R)(t) - (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
+ * (NFFT_R)(r) / (NFFT_R)(rr);
+ x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = (NFFT_R) r / (NFFT_R)(rr);
}
if (r == 0)
- w[(t + T / 2) * rr + (r + rr / 2)] = K(1.0) / K(4.0) / W;
+ w[(t + T / 2) * rr + (r + rr / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
else
- w[(t + T / 2) * rr + (r + rr / 2)] = FABS((R) r) / W;
+ w[(t + T / 2) * rr + (r + rr / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
}
}
@@ -80,13 +78,13 @@ static int linogram_grid(int T, int rr, R *x, R *w)
}
/** discrete pseudo-polar FFT */
-static int linogram_dft(C *f_hat, int NN, C *f, int T, int rr, int m)
+static int linogram_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
{
- ticks t0, t1;
+ double t0, t1;
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int N[2], n[2];
int M = T * rr; /**< number of knots */
@@ -96,11 +94,11 @@ static int linogram_dft(C *f_hat, int NN, C *f, int T, int rr, int m)
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -123,10 +121,10 @@ static int linogram_dft(C *f_hat, int NN, C *f, int T, int rr, int m)
my_nfft_plan.f_hat[k] = f_hat[k];
/** NFFT-2D */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
NFFT(trafo_direct)(&my_nfft_plan);
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
/** copy result */
for (j = 0; j < my_nfft_plan.M_total; j++)
@@ -141,13 +139,13 @@ static int linogram_dft(C *f_hat, int NN, C *f, int T, int rr, int m)
}
/** NFFT-based pseudo-polar FFT */
-static int linogram_fft(C *f_hat, int NN, C *f, int T, int rr, int m)
+static int linogram_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
{
- ticks t0, t1;
+ double t0, t1;
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int N[2], n[2];
int M = T * rr; /**< number of knots */
@@ -157,11 +155,11 @@ static int linogram_fft(C *f_hat, int NN, C *f, int T, int rr, int m)
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -194,10 +192,10 @@ static int linogram_fft(C *f_hat, int NN, C *f, int T, int rr, int m)
my_nfft_plan.f_hat[k] = f_hat[k];
/** NFFT-2D */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
NFFT(trafo)(&my_nfft_plan);
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
/** copy result */
for (j = 0; j < my_nfft_plan.M_total; j++)
@@ -212,15 +210,15 @@ static int linogram_fft(C *f_hat, int NN, C *f, int T, int rr, int m)
}
/** NFFT-based inverse pseudo-polar FFT */
-static int inverse_linogram_fft(C *f, int T, int rr, C *f_hat, int NN,
+static int inverse_linogram_fft(NFFT_C *f, int T, int rr, NFFT_C *f_hat, int NN,
int max_i, int m)
{
- ticks t0, t1;
+ double t0, t1;
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
SOLVER(plan_complex) my_infft_plan; /**< plan for the inverse nfft */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int l; /**< index for iterations */
int N[2], n[2];
@@ -231,11 +229,11 @@ static int inverse_linogram_fft(C *f, int T, int rr, C *f_hat, int NN,
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -275,17 +273,17 @@ static int inverse_linogram_fft(C *f, int T, int rr, C *f_hat, int NN,
for (k = 0; k < my_nfft_plan.N[1]; k++)
{
my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
- SQRT(
- POW((R)(j - my_nfft_plan.N[0] / 2), K(2.0))
- + POW((R)(k - my_nfft_plan.N[1] / 2), K(2.0)))
- > (R)(my_nfft_plan.N[0] / 2) ? 0 : 1);
+ NFFT_M(sqrt)(
+ NFFT_M(pow)((NFFT_R)(j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
+ + NFFT_M(pow)((NFFT_R)(k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
+ > (NFFT_R)(my_nfft_plan.N[0] / 2) ? 0 : 1);
}
/** initialise some guess f_hat_0 */
for (k = 0; k < my_nfft_plan.N_total; k++)
- my_infft_plan.f_hat_iter[k] = K(0.0) + II * K(0.0);
+ my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
/** solve the system */
SOLVER(before_loop_complex)(&my_infft_plan);
@@ -303,8 +301,8 @@ static int inverse_linogram_fft(C *f, int T, int rr, C *f_hat, int NN,
}
}
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
/** copy result */
for (k = 0; k < my_nfft_plan.N_total; k++)
@@ -322,30 +320,30 @@ static int inverse_linogram_fft(C *f, int T, int rr, C *f_hat, int NN,
/** Comparison of the FFTW, linogram FFT, and inverse linogram FFT */
static int comparison_fft(FILE *fp, int N, int T, int rr)
{
- ticks t0, t1;
+ double t0, t1;
FFTW(plan) my_fftw_plan;
- C *f_hat, *f;
+ NFFT_C *f_hat, *f;
int m, k;
- R t_fft, t_dft_linogram;
+ NFFT_R t_fft, t_dft_linogram;
- f_hat = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
- f = (C *) NFFT(malloc)(sizeof(C) * (size_t)((T * rr / 4) * 5));
+ f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
+ f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)((T * rr / 4) * 5));
my_fftw_plan = FFTW(plan_dft_2d)(N, N, f_hat, f, FFTW_BACKWARD, FFTW_MEASURE);
for (k = 0; k < N * N; k++)
- f_hat[k] = NFFT(drand48)() + II * NFFT(drand48)();
+ f_hat[k] = NFFT(drand48)() + _Complex_I * NFFT(drand48)();
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
for (m = 0; m < 65536 / N; m++)
{
FFTW(execute)(my_fftw_plan);
/* touch */
- f_hat[2] = K(2.0) * f_hat[0];
+ f_hat[2] = NFFT_K(2.0) * f_hat[0];
}
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
- t_fft = (R)(N) * GLOBAL_elapsed_time / K(65536.0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
+ t_fft = (NFFT_R)(N) * GLOBAL_elapsed_time / NFFT_K(65536.0);
if (N < 256)
{
@@ -356,25 +354,25 @@ static int comparison_fft(FILE *fp, int N, int T, int rr)
for (m = 3; m <= 9; m += 3)
{
if ((m == 3) && (N < 256))
- fprintf(fp, "%d\t&\t&\t%1.1" __FES__ "&\t%1.1" __FES__ "&\t%d\t", N, t_fft, t_dft_linogram,
+ fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t%1.1" NFFT__FES__ "&\t%d\t", N, t_fft, t_dft_linogram,
m);
else if (m == 3)
- fprintf(fp, "%d\t&\t&\t%1.1" __FES__ "&\t &\t%d\t", N, t_fft, m);
+ fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t &\t%d\t", N, t_fft, m);
else
fprintf(fp, " \t&\t&\t &\t &\t%d\t", m);
- printf("N=%d\tt_fft=%1.1" __FES__ "\tt_dft_linogram=%1.1" __FES__ "\tm=%d\t", N, t_fft,
+ printf("N=%d\tt_fft=%1.1" NFFT__FES__ "\tt_dft_linogram=%1.1" NFFT__FES__ "\tm=%d\t", N, t_fft,
t_dft_linogram, m);
linogram_fft(f_hat, N, f, T, rr, m);
- fprintf(fp, "%1.1" __FES__ "&\t", GLOBAL_elapsed_time);
- printf("t_linogram=%1.1" __FES__ "\t", GLOBAL_elapsed_time);
+ fprintf(fp, "%1.1" NFFT__FES__ "&\t", GLOBAL_elapsed_time);
+ printf("t_linogram=%1.1" NFFT__FES__ "\t", GLOBAL_elapsed_time);
inverse_linogram_fft(f, T, rr, f_hat, N, m + 3, m);
if (m == 9)
- fprintf(fp, "%1.1" __FES__ "\\\\\\hline\n", GLOBAL_elapsed_time);
+ fprintf(fp, "%1.1" NFFT__FES__ "\\\\\\hline\n", GLOBAL_elapsed_time);
else
- fprintf(fp, "%1.1" __FES__ "\\\\\n", GLOBAL_elapsed_time);
- printf("t_ilinogram=%1.1" __FES__ "\n", GLOBAL_elapsed_time);
+ fprintf(fp, "%1.1" NFFT__FES__ "\\\\\n", GLOBAL_elapsed_time);
+ printf("t_ilinogram=%1.1" NFFT__FES__ "\n", GLOBAL_elapsed_time);
}
fflush(fp);
@@ -391,12 +389,12 @@ int main(int argc, char **argv)
int N; /**< linogram FFT size NxN */
int T, rr; /**< number of directions/offsets */
int M; /**< number of knots of linogram grid */
- R *x, *w; /**< knots and associated weights */
- C *f_hat, *f, *f_direct, *f_tilde;
+ NFFT_R *x, *w; /**< knots and associated weights */
+ NFFT_C *f_hat, *f, *f_direct, *f_tilde;
int k;
int max_i; /**< number of iterations */
int m;
- R temp1, temp2, E_max = K(0.0);
+ NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
FILE *fp1, *fp2;
char filename[30];
int logN;
@@ -427,13 +425,13 @@ int main(int argc, char **argv)
rr = atoi(argv[3]);
printf("N=%d, linogram grid with T=%d, R=%d => ", N, T, rr);
- x = (R *) NFFT(malloc)((size_t)(5 * T * rr / 2) * (sizeof(R)));
- w = (R *) NFFT(malloc)((size_t)(5 * T * rr / 4) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * rr / 2) * (sizeof(NFFT_R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * rr / 4) * (sizeof(NFFT_R)));
- f_hat = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
- f = (C *) NFFT(malloc)(sizeof(C) * (size_t)(5 * T * rr / 4)); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
- f_direct = (C *) NFFT(malloc)(sizeof(C) * (size_t)(5 * T * rr / 4));
- f_tilde = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
+ f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
+ f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(5 * T * rr / 4)); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
+ f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(5 * T * rr / 4));
+ f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
/** generate knots of linogram grid */
M = linogram_grid(T, rr, x, w);
@@ -446,9 +444,9 @@ int main(int argc, char **argv)
return EXIT_FAILURE;
for (k = 0; k < N * N; k++)
{
- fscanf(fp1, __FR__ " ", &temp1);
- fscanf(fp2, __FR__ " ", &temp2);
- f_hat[k] = temp1 + II * temp2;
+ fscanf(fp1, NFFT__FR__ " ", &temp1);
+ fscanf(fp2, NFFT__FR__ " ", &temp2);
+ f_hat[k] = temp1 + _Complex_I * temp2;
}
fclose(fp1);
fclose(fp2);
@@ -469,8 +467,8 @@ int main(int argc, char **argv)
E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
//E_max=NFFT(error_l_2_complex)(f_direct,f,M);
- printf("m=%2d: E_max = %" __FES__ "\n", m, E_max);
- fprintf(fp1, "%" __FES__ "\n", E_max);
+ printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max);
+ fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
}
fclose(fp1);
@@ -487,8 +485,8 @@ int main(int argc, char **argv)
/** compute maximum error */
E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
- printf("%3d iterations: E_max = %" __FES__ "\n", max_i, E_max);
- fprintf(fp1, "%" __FES__ "\n", E_max);
+ printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max);
+ fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
}
fclose(fp1);
}
diff --git a/applications/polarFFT/mpolar_fft_test.c b/applications/polarFFT/mpolar_fft_test.c.in
similarity index 71%
rename from applications/polarFFT/mpolar_fft_test.c
rename to applications/polarFFT/mpolar_fft_test.c.in
index 1fba37a..827ac22 100644
--- a/applications/polarFFT/mpolar_fft_test.c
+++ b/applications/polarFFT/mpolar_fft_test.c.in
@@ -27,16 +27,14 @@
* \author Markus Fenn
* \date 2006
*/
-#include "config.h"
#include <math.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
/**
* \defgroup applications_polarFFT_mpolar mpolar_fft_test
@@ -44,7 +42,7 @@
* \{
*/
-R GLOBAL_elapsed_time;
+NFFT_R GLOBAL_elapsed_time;
/** Generates the points \f$x_{t,j}\f$ with weights \f$w_{t,j}\f$
* for the modified polar grid with \f$T\f$ angles and \f$R\f$ offsets.
@@ -59,32 +57,32 @@ R GLOBAL_elapsed_time;
* The number of nodes for the modified polar grid can be estimated as
* \f$M \approx \frac{4}{\pi}\log(1+\sqrt{2}) T R\f$.
*/
-static int mpolar_grid(int T, int S, R *x, R *w)
+static int mpolar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
{
int t, r;
- R W;
- int R2 = 2 * (int)(LRINT(CEIL(SQRT(K(2.0)) * (R)(S) / K(2.0))));
- R xx, yy;
+ NFFT_R W;
+ int R2 = 2 * (int)(NFFT_M(lrint)(NFFT_M(ceil)(NFFT_M(sqrt)(NFFT_K(2.0)) * (NFFT_R)(S) / NFFT_K(2.0))));
+ NFFT_R xx, yy;
int M = 0;
for (t = -T / 2; t < T / 2; t++)
{
for (r = -R2 / 2; r < R2 / 2; r++)
{
- xx = (R) (r) / (R)(S) * COS(KPI * (R)(t) / (R)(T));
- yy = (R) (r) / (R)(S) * SIN(KPI * (R)(t) / (R)(T));
+ xx = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
+ yy = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
- if (((-K(0.5) - K(1.0) / (R) S) <= xx) & (xx <= (K(0.5) + K(1.0) / (R) S))
- & ((-K(0.5) - K(1.0) / (R) S) <= yy)
- & (yy <= (K(0.5) + K(1.0) / (R) S)))
+ if (((-NFFT_K(0.5) - NFFT_K(1.0) / (NFFT_R) S) <= xx) & (xx <= (NFFT_K(0.5) + NFFT_K(1.0) / (NFFT_R) S))
+ & ((-NFFT_K(0.5) - NFFT_K(1.0) / (NFFT_R) S) <= yy)
+ & (yy <= (NFFT_K(0.5) + NFFT_K(1.0) / (NFFT_R) S)))
{
x[2 * M + 0] = xx;
x[2 * M + 1] = yy;
if (r == 0)
- w[M] = K(1.0) / K(4.0);
+ w[M] = NFFT_K(1.0) / NFFT_K(4.0);
else
- w[M] = FABS((R) r);
+ w[M] = NFFT_M(fabs)((NFFT_R) r);
M++; /** count the knots */
}
@@ -92,7 +90,7 @@ static int mpolar_grid(int T, int S, R *x, R *w)
}
/** normalize the weights */
- W = K(0.0);
+ W = NFFT_K(0.0);
for (t = 0; t < M; t++)
W += w[t];
@@ -103,13 +101,13 @@ static int mpolar_grid(int T, int S, R *x, R *w)
}
/** discrete mpolar FFT */
-static int mpolar_dft(C *f_hat, int NN, C *f, int T, int S, int m)
+static int mpolar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
{
- ticks t0, t1;
+ double t0, t1;
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int N[2], n[2];
int M; /**< number of knots */
@@ -119,11 +117,11 @@ static int mpolar_dft(C *f_hat, int NN, C *f, int T, int S, int m)
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(5 * (T / 2) * S) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * S) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(5 * (T * S) / 4) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T * S) / 4) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -145,13 +143,13 @@ static int mpolar_dft(C *f_hat, int NN, C *f, int T, int S, int m)
for (k = 0; k < my_nfft_plan.N_total; k++)
my_nfft_plan.f_hat[k] = f_hat[k];
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
/** NDFT-2D */
NFFT(trafo_direct)(&my_nfft_plan);
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
/** copy result */
for (j = 0; j < my_nfft_plan.M_total; j++)
@@ -166,13 +164,13 @@ static int mpolar_dft(C *f_hat, int NN, C *f, int T, int S, int m)
}
/** NFFT-based mpolar FFT */
-static int mpolar_fft(C *f_hat, int NN, C *f, int T, int S, int m)
+static int mpolar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, int m)
{
- ticks t0, t1;
+ double t0, t1;
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int N[2], n[2];
int M; /**< number of knots */
@@ -182,11 +180,11 @@ static int mpolar_fft(C *f_hat, int NN, C *f, int T, int S, int m)
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -219,13 +217,13 @@ static int mpolar_fft(C *f_hat, int NN, C *f, int T, int S, int m)
for (k = 0; k < my_nfft_plan.N_total; k++)
my_nfft_plan.f_hat[k] = f_hat[k];
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
/** NFFT-2D */
NFFT(trafo)(&my_nfft_plan);
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
/** copy result */
for (j = 0; j < my_nfft_plan.M_total; j++)
@@ -240,15 +238,15 @@ static int mpolar_fft(C *f_hat, int NN, C *f, int T, int S, int m)
}
/** inverse NFFT-based mpolar FFT */
-static int inverse_mpolar_fft(C *f, int T, int S, C *f_hat, int NN, int max_i,
+static int inverse_mpolar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, int NN, int max_i,
int m)
{
- ticks t0, t1;
+ double t0, t1;
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
SOLVER(plan_complex) my_infft_plan; /**< plan for the inverse nfft */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int l; /**< index for iterations */
int N[2], n[2];
@@ -259,11 +257,11 @@ static int inverse_mpolar_fft(C *f, int T, int S, C *f_hat, int NN, int max_i,
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -303,17 +301,17 @@ static int inverse_mpolar_fft(C *f, int T, int S, C *f_hat, int NN, int max_i,
for (k = 0; k < my_nfft_plan.N[1]; k++)
{
my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
- SQRT(
- POW((R)(j - my_nfft_plan.N[0] / 2), K(2.0))
- + POW((R)(k - my_nfft_plan.N[1] / 2), K(2.0)))
- > (R)(my_nfft_plan.N[0] / 2) ? K(0.0) : K(1.0));
+ NFFT_M(sqrt)(
+ NFFT_M(pow)((NFFT_R)(j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
+ + NFFT_M(pow)((NFFT_R)(k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
+ > (NFFT_R)(my_nfft_plan.N[0] / 2) ? NFFT_K(0.0) : NFFT_K(1.0));
}
/** initialise some guess f_hat_0 */
for (k = 0; k < my_nfft_plan.N_total; k++)
- my_infft_plan.f_hat_iter[k] = K(0.0) + II * K(0.0);
+ my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
/** solve the system */
SOLVER(before_loop_complex)(&my_infft_plan);
@@ -332,8 +330,8 @@ static int inverse_mpolar_fft(C *f, int T, int S, C *f_hat, int NN, int max_i,
}
}
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
/** copy result */
for (k = 0; k < my_nfft_plan.N_total; k++)
@@ -351,30 +349,30 @@ static int inverse_mpolar_fft(C *f, int T, int S, C *f_hat, int NN, int max_i,
/** Comparison of the FFTW, mpolar FFT, and inverse mpolar FFT */
static int comparison_fft(FILE *fp, int N, int T, int S)
{
- ticks t0, t1;
+ double t0, t1;
FFTW(plan) my_fftw_plan;
- C *f_hat, *f;
+ NFFT_C *f_hat, *f;
int m, k;
- R t_fft, t_dft_mpolar;
+ NFFT_R t_fft, t_dft_mpolar;
- f_hat = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
- f = (C *) NFFT(malloc)(sizeof(C) * (size_t)((T * S / 4) * 5));
+ f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
+ f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)((T * S / 4) * 5));
my_fftw_plan = FFTW(plan_dft_2d)(N, N, f_hat, f, FFTW_BACKWARD, FFTW_MEASURE);
for (k = 0; k < N * N; k++)
- f_hat[k] = NFFT(drand48)() + II * NFFT(drand48)();
+ f_hat[k] = NFFT(drand48)() + _Complex_I * NFFT(drand48)();
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
for (m = 0; m < 65536 / N; m++)
{
FFTW(execute)(my_fftw_plan);
/* touch */
- f_hat[2] = K(2.0) * f_hat[0];
+ f_hat[2] = NFFT_K(2.0) * f_hat[0];
}
- t1 = getticks();
- GLOBAL_elapsed_time = NFFT(elapsed_seconds)(t1, t0);
- t_fft = (R)(N) * GLOBAL_elapsed_time / K(65536.0);
+ t1 = NFFT(clock_gettime_seconds)();
+ GLOBAL_elapsed_time = (t1 - t0);
+ t_fft = (NFFT_R)(N) * GLOBAL_elapsed_time / NFFT_K(65536.0);
if (N < 256)
{
@@ -385,24 +383,24 @@ static int comparison_fft(FILE *fp, int N, int T, int S)
for (m = 3; m <= 9; m += 3)
{
if ((m == 3) && (N < 256))
- fprintf(fp, "%d\t&\t&\t%1.1" __FES__ "&\t%1.1" __FES__ "&\t%d\t", N, t_fft, t_dft_mpolar, m);
+ fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t%1.1" NFFT__FES__ "&\t%d\t", N, t_fft, t_dft_mpolar, m);
else if (m == 3)
- fprintf(fp, "%d\t&\t&\t%1.1" __FES__ "&\t &\t%d\t", N, t_fft, m);
+ fprintf(fp, "%d\t&\t&\t%1.1" NFFT__FES__ "&\t &\t%d\t", N, t_fft, m);
else
fprintf(fp, " \t&\t&\t &\t &\t%d\t", m);
- printf("N=%d\tt_fft=%1.1" __FES__ "\tt_dft_mpolar=%1.1" __FES__ "\tm=%d\t", N, t_fft,
+ printf("N=%d\tt_fft=%1.1" NFFT__FES__ "\tt_dft_mpolar=%1.1" NFFT__FES__ "\tm=%d\t", N, t_fft,
t_dft_mpolar, m);
mpolar_fft(f_hat, N, f, T, S, m);
- fprintf(fp, "%1.1" __FES__ "&\t", GLOBAL_elapsed_time);
- printf("t_mpolar=%1.1" __FES__ "\t", GLOBAL_elapsed_time);
+ fprintf(fp, "%1.1" NFFT__FES__ "&\t", GLOBAL_elapsed_time);
+ printf("t_mpolar=%1.1" NFFT__FES__ "\t", GLOBAL_elapsed_time);
inverse_mpolar_fft(f, T, S, f_hat, N, 2 * m, m);
if (m == 9)
- fprintf(fp, "%1.1" __FES__ "\\\\\\hline\n", GLOBAL_elapsed_time);
+ fprintf(fp, "%1.1" NFFT__FES__ "\\\\\\hline\n", GLOBAL_elapsed_time);
else
- fprintf(fp, "%1.1" __FES__ "\\\\\n", GLOBAL_elapsed_time);
- printf("t_impolar=%1.1" __FES__ "\n", GLOBAL_elapsed_time);
+ fprintf(fp, "%1.1" NFFT__FES__ "\\\\\n", GLOBAL_elapsed_time);
+ printf("t_impolar=%1.1" NFFT__FES__ "\n", GLOBAL_elapsed_time);
}
fflush(fp);
@@ -419,12 +417,12 @@ int main(int argc, char **argv)
int N; /**< mpolar FFT size NxN */
int T, S; /**< number of directions/offsets */
int M; /**< number of knots of mpolar grid */
- R *x, *w; /**< knots and associated weights */
- C *f_hat, *f, *f_direct, *f_tilde;
+ NFFT_R *x, *w; /**< knots and associated weights */
+ NFFT_C *f_hat, *f, *f_direct, *f_tilde;
int k;
int max_i; /**< number of iterations */
int m;
- R temp1, temp2, E_max = K(0.0);
+ NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
FILE *fp1, *fp2;
char filename[30];
int logN;
@@ -455,13 +453,13 @@ int main(int argc, char **argv)
S = atoi(argv[3]);
printf("N=%d, modified polar grid with T=%d, R=%d => ", N, T, S);
- x = (R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(R)));
- w = (R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 2) * (sizeof(NFFT_R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * S / 4) * (sizeof(NFFT_R)));
- f_hat = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
- f = (C *) NFFT(malloc)(sizeof(C) * (size_t)(1.25 * T * S)); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
- f_direct = (C *) NFFT(malloc)(sizeof(C) * (size_t)(1.25 * T * S));
- f_tilde = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
+ f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
+ f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(1.25 * T * S)); /* 4/pi*log(1+sqrt(2)) = 1.122... < 1.25 */
+ f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(1.25 * T * S));
+ f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
/** generate knots of mpolar grid */
M = mpolar_grid(T, S, x, w);
@@ -474,9 +472,9 @@ int main(int argc, char **argv)
return (-1);
for (k = 0; k < N * N; k++)
{
- fscanf(fp1, __FR__ " ", &temp1);
- fscanf(fp2, __FR__ " ", &temp2);
- f_hat[k] = temp1 + II * temp2;
+ fscanf(fp1, NFFT__FR__ " ", &temp1);
+ fscanf(fp2, NFFT__FR__ " ", &temp2);
+ f_hat[k] = temp1 + _Complex_I * temp2;
}
fclose(fp1);
fclose(fp2);
@@ -495,8 +493,8 @@ int main(int argc, char **argv)
/** compute error of fast mpolar FFT */
E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
- printf("m=%2d: E_max = %" __FES__ "\n", m, E_max);
- fprintf(fp1, "%" __FES__ "\n", E_max);
+ printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max);
+ fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
}
fclose(fp1);
@@ -513,8 +511,8 @@ int main(int argc, char **argv)
/** compute maximum relativ error */
E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
- printf("%3d iterations: E_max = %" __FES__ "\n", max_i, E_max);
- fprintf(fp1, "%" __FES__ "\n", E_max);
+ printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max);
+ fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
}
fclose(fp1);
}
diff --git a/applications/polarFFT/polar_fft_test.c b/applications/polarFFT/polar_fft_test.c.in
similarity index 80%
rename from applications/polarFFT/polar_fft_test.c
rename to applications/polarFFT/polar_fft_test.c.in
index 6ed51dc..3e2c9a4 100644
--- a/applications/polarFFT/polar_fft_test.c
+++ b/applications/polarFFT/polar_fft_test.c.in
@@ -26,16 +26,14 @@
* \author Markus Fenn
* \date 2006
*/
-#include "config.h"
#include <math.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
/**
* \defgroup applications_polarFFT_polar polar_fft_test
@@ -74,21 +72,21 @@
* Thus, the sum of all weights is \f$\frac{\pi}{4}(1+\frac{1}{R^2})\f$ and
* we divide by this value for normalization.
*/
-static int polar_grid(int T, int S, R *x, R *w)
+static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
{
int t, r;
- R W = (R) T * (((R) S / K(2.0)) * ((R) S / K(2.0)) + K(1.0) / K(4.0));
+ NFFT_R W = (NFFT_R) T * (((NFFT_R) S / NFFT_K(2.0)) * ((NFFT_R) S / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
for (t = -T / 2; t < T / 2; t++)
{
for (r = -S / 2; r < S / 2; r++)
{
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (R) (r) / (R)(S) * COS(KPI * (R)(t) / (R)(T));
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (R) (r) / (R)(S) * SIN(KPI * (R)(t) / (R)(T));
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) (r) / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
if (r == 0)
- w[(t + T / 2) * S + (r + S / 2)] = K(1.0) / K(4.0) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
else
- w[(t + T / 2) * S + (r + S / 2)] = FABS((R) r) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
}
}
@@ -96,13 +94,13 @@ static int polar_grid(int T, int S, R *x, R *w)
}
/** discrete polar FFT */
-static int polar_dft(C *f_hat, int NN, C *f, int T, int S,
+static int polar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S,
int m)
{
int j, k; /**< index for nodes and frequencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int N[2], n[2];
int M = T * S; /**< number of knots */
@@ -112,11 +110,11 @@ static int polar_dft(C *f_hat, int NN, C *f, int T, int S,
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * S) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -154,13 +152,13 @@ static int polar_dft(C *f_hat, int NN, C *f, int T, int S,
}
/** NFFT-based polar FFT */
-static int polar_fft(C *f_hat, int NN, C *f, int T, int S,
+static int polar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S,
int m)
{
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int N[2], n[2];
int M = T * S; /**< number of knots */
@@ -170,11 +168,11 @@ static int polar_fft(C *f_hat, int NN, C *f, int T, int S,
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * S) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -222,14 +220,14 @@ static int polar_fft(C *f_hat, int NN, C *f, int T, int S,
}
/** inverse NFFT-based polar FFT */
-static int inverse_polar_fft(C *f, int T, int S, C *f_hat,
+static int inverse_polar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat,
int NN, int max_i, int m)
{
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
SOLVER(plan_complex) my_infft_plan; /**< plan for the inverse nfft */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int l; /**< index for iterations */
int N[2], n[2];
@@ -240,11 +238,11 @@ static int inverse_polar_fft(C *f, int T, int S, C *f_hat,
N[1] = NN;
n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */
- x = (R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * S) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -284,15 +282,15 @@ static int inverse_polar_fft(C *f, int T, int S, C *f_hat,
for (k = 0; k < my_nfft_plan.N[1]; k++)
{
my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
- SQRT(
- POW((R) (j - my_nfft_plan.N[0] / 2), K(2.0))
- + POW((R) (k - my_nfft_plan.N[1] / 2), K(2.0)))
- > ((R) (my_nfft_plan.N[0] / 2)) ? 0 : 1);
+ NFFT_M(sqrt)(
+ NFFT_M(pow)((NFFT_R) (j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
+ + NFFT_M(pow)((NFFT_R) (k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
+ > ((NFFT_R) (my_nfft_plan.N[0] / 2)) ? 0 : 1);
}
/** initialise some guess f_hat_0 */
for (k = 0; k < my_nfft_plan.N_total; k++)
- my_infft_plan.f_hat_iter[k] = K(0.0) + II * K(0.0);
+ my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
/** solve the system */
SOLVER(before_loop_complex)(&my_infft_plan);
@@ -330,12 +328,12 @@ int main(int argc, char **argv)
int N; /**< mpolar FFT size NxN */
int T, S; /**< number of directions/offsets */
int M; /**< number of knots of mpolar grid */
- R *x, *w; /**< knots and associated weights */
- C *f_hat, *f, *f_direct, *f_tilde;
+ NFFT_R *x, *w; /**< knots and associated weights */
+ NFFT_C *f_hat, *f, *f_direct, *f_tilde;
int k;
int max_i; /**< number of iterations */
int m = 1;
- R temp1, temp2, E_max = K(0.0);
+ NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
FILE *fp1, *fp2;
char filename[30];
@@ -354,13 +352,13 @@ int main(int argc, char **argv)
S = atoi(argv[3]);
printf("N=%d, polar grid with T=%d, R=%d => ", N, T, S);
- x = (R *) NFFT(malloc)((size_t)(2 * 5 * (T / 2) * (S / 2)) * (sizeof(R)));
- w = (R *) NFFT(malloc)((size_t)(5 * (T / 2) * (S / 2)) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * 5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R)));
- f_hat = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
- f = (C *) NFFT(malloc)(sizeof(C) * (size_t)(T * S));
- f_direct = (C *) NFFT(malloc)(sizeof(C) * (size_t)(T * S));
- f_tilde = (C *) NFFT(malloc)(sizeof(C) * (size_t)(N * N));
+ f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
+ f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S));
+ f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S));
+ f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N));
/** generate knots of mpolar grid */
M = polar_grid(T, S, x, w);
@@ -373,9 +371,9 @@ int main(int argc, char **argv)
return (-1);
for (k = 0; k < N * N; k++)
{
- fscanf(fp1, __FR__ " ", &temp1);
- fscanf(fp2, __FR__ " ", &temp2);
- f_hat[k] = temp1 + II * temp2;
+ fscanf(fp1, NFFT__FR__ " ", &temp1);
+ fscanf(fp2, NFFT__FR__ " ", &temp2);
+ f_hat[k] = temp1 + _Complex_I * temp2;
}
fclose(fp1);
fclose(fp2);
@@ -394,8 +392,8 @@ int main(int argc, char **argv)
/** compute error of fast polar FFT */
E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
- printf("m=%2d: E_max = %" __FES__ "\n", m, E_max);
- fprintf(fp1, "%" __FES__ "\n", E_max);
+ printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max);
+ fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
}
fclose(fp1);
@@ -419,8 +417,8 @@ int main(int argc, char **argv)
}
*/
E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
- printf("%3d iterations: E_max = %" __FES__ "\n", max_i, E_max);
- fprintf(fp1, "%" __FES__ "\n", E_max);
+ printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max);
+ fprintf(fp1, "%" NFFT__FES__ "\n", E_max);
}
fclose(fp1);
}
diff --git a/applications/radon/inverse_radon.c b/applications/radon/inverse_radon.c.in
similarity index 72%
rename from applications/radon/inverse_radon.c
rename to applications/radon/inverse_radon.c.in
index b0ec803..80e4fc0 100644
--- a/applications/radon/inverse_radon.c
+++ b/applications/radon/inverse_radon.c.in
@@ -36,41 +36,39 @@
* \author Markus Fenn
* \date 2005
*/
-#include "config.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
/** define weights of kernel function for discrete Radon transform */
/*#define KERNEL(r) 1.0 */
-#define KERNEL(r) (K(1.0)-FABS((R)(r))/((R)S/2))
+#define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))
/** generates the points x with weights w
* for the polar grid with T angles and R offsets
*/
-static int polar_grid(int T, int S, R *x, R *w)
+static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
{
int t, r;
- R W = (R) T * (((R) S / K(2.0)) * ((R) S / K(2.0)) + K(1.0) / K(4.0));
+ NFFT_R W = (NFFT_R) T * (((NFFT_R) S / NFFT_K(2.0)) * ((NFFT_R) S / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
for (t = -T / 2; t < T / 2; t++)
{
for (r = -S / 2; r < S / 2; r++)
{
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (R) r / (R)(S) * COS(KPI * (R)(t) / (R)(T));
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (R) r / (R)(S) * SIN(KPI * (R)(t) / (R)(T));
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
if (r == 0)
- w[(t + T / 2) * S + (r + S / 2)] = K(1.0) / K(4.0) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
else
- w[(t + T / 2) * S + (r + S / 2)] = FABS((R) r) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
}
}
@@ -80,10 +78,10 @@ static int polar_grid(int T, int S, R *x, R *w)
/** generates the points x with weights w
* for the linogram grid with T slopes and R offsets
*/
-static int linogram_grid(int T, int S, R *x, R *w)
+static int linogram_grid(int T, int S, NFFT_R *x, NFFT_R *w)
{
int t, r;
- R W = (R) T * (((R) S / K(2.0)) * ((R) S / K(2.0)) + K(1.0) / K(4.0));
+ NFFT_R W = (NFFT_R) T * (((NFFT_R) S / NFFT_K(2.0)) * ((NFFT_R) S / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
for (t = -T / 2; t < T / 2; t++)
{
@@ -91,20 +89,20 @@ static int linogram_grid(int T, int S, R *x, R *w)
{
if (t < 0)
{
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (R) r / (R)(S);
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = K(4.0) * ((R)(t) + (R)(T) / K(4.0)) / (R)(T) * (R)(r)
- / (R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = NFFT_K(4.0) * ((NFFT_R)(t) + (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T) * (NFFT_R)(r)
+ / (NFFT_R)(S);
}
else
{
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = -K(4.0) * ((R)(t) - (R)(T) / K(4.0)) / (R)(T)
- * (R)(r) / (R)(S);
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (R) r / (R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = -NFFT_K(4.0) * ((NFFT_R)(t) - (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
+ * (NFFT_R)(r) / (NFFT_R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S);
}
if (r == 0)
- w[(t + T / 2) * S + (r + S / 2)] = K(1.0) / K(4.0) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
else
- w[(t + T / 2) * S + (r + S / 2)] = FABS((R) r) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
}
}
@@ -115,18 +113,18 @@ static int linogram_grid(int T, int S, R *x, R *w)
* on the grid given by gridfcn() with T angles and R offsets
* by a NFFT-based CG-type algorithm
*/
-static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, R *Rf, int NN, R *f,
+static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, NFFT_R *Rf, int NN, NFFT_R *f,
int max_i)
{
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
SOLVER(plan_complex) my_infft_plan; /**< plan for the inverse nfft */
- C *fft; /**< variable for the fftw-1Ds */
+ NFFT_C *fft; /**< variable for the fftw-1Ds */
FFTW(plan) my_fftw_plan; /**< plan for the fftw-1Ds */
int t, r; /**< index for directions and offsets */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int l; /**< index for iterations */
int N[2], n[2];
@@ -137,14 +135,14 @@ static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, R *Rf, int NN, R
N[1] = NN;
n[1] = 2 * N[1];
- fft = (C *) NFFT(malloc)((size_t)(S) * sizeof(C));
+ fft = (NFFT_C *) NFFT(malloc)((size_t)(S) * sizeof(NFFT_C));
my_fftw_plan = FFTW(plan_dft_1d)(S, fft, fft, FFTW_FORWARD, FFTW_MEASURE);
- x = (R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * S) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -167,7 +165,7 @@ static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, R *Rf, int NN, R
if (j % S)
my_infft_plan.w[j] = w[j];
else
- my_infft_plan.w[j] = K(0.0);
+ my_infft_plan.w[j] = NFFT_K(0.0);
}
/** precompute psi, the entries of the matrix B */
@@ -184,26 +182,26 @@ static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, R *Rf, int NN, R
for (t = 0; t < T; t++)
{
/* for(r=0; r<R/2; r++)
- fft[r] = cexp(I*KPI*r)*Rf[t*R+(r+R/2)];
+ fft[r] = cexp(I*NFFT_KPI*r)*Rf[t*R+(r+R/2)];
for(r=0; r<R/2; r++)
- fft[r+R/2] = cexp(I*KPI*r)*Rf[t*R+r];
+ fft[r+R/2] = cexp(I*NFFT_KPI*r)*Rf[t*R+r];
*/
for (r = 0; r < S; r++)
- fft[r] = Rf[t * S + r] + II * K(0.0);
+ fft[r] = Rf[t * S + r] + _Complex_I * NFFT_K(0.0);
NFFT(fftshift_complex_int)(fft, 1, &S);
FFTW(execute)(my_fftw_plan);
NFFT(fftshift_complex_int)(fft, 1, &S);
- my_infft_plan.y[t * S] = K(0.0);
+ my_infft_plan.y[t * S] = NFFT_K(0.0);
for (r = -S / 2 + 1; r < S / 2; r++)
my_infft_plan.y[t * S + (r + S / 2)] = fft[r + S / 2] / KERNEL(r);
}
/** initialise some guess f_hat_0 */
for (k = 0; k < my_nfft_plan.N_total; k++)
- my_infft_plan.f_hat_iter[k] = K(0.0) + II * K(0.0);
+ my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
/** solve the system */
SOLVER(before_loop_complex)(&my_infft_plan);
@@ -226,7 +224,7 @@ static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, R *Rf, int NN, R
/** copy result */
for (k = 0; k < my_nfft_plan.N_total; k++)
- f[k] = CREAL(my_infft_plan.f_hat_iter[k]);
+ f[k] = NFFT_M(creal)(my_infft_plan.f_hat_iter[k]);
/** finalise the plans and free the variables */
FFTW(destroy_plan)(my_fftw_plan);
@@ -246,7 +244,7 @@ int main(int argc, char **argv)
int T, S; /**< number of directions/offsets */
FILE *fp;
int N; /**< image size */
- R *Rf, *iRf;
+ NFFT_R *Rf, *iRf;
int max_i; /**< number of iterations */
if (argc != 6)
@@ -272,14 +270,14 @@ int main(int argc, char **argv)
/*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
max_i = atoi(argv[5]);
- Rf = (R *) NFFT(malloc)((size_t)(T * S) * (sizeof(R)));
- iRf = (R *) NFFT(malloc)((size_t)(N * N) * (sizeof(R)));
+ Rf = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
+ iRf = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (sizeof(NFFT_R)));
/** load data */
fp = fopen("sinogram_data.bin", "rb");
if (fp == NULL)
return EXIT_FAILURE;
- fread(Rf, sizeof(R), (size_t)(T * S), fp);
+ fread(Rf, sizeof(NFFT_R), (size_t)(T * S), fp);
fclose(fp);
/** inverse Radon transform */
@@ -289,7 +287,7 @@ int main(int argc, char **argv)
fp = fopen("output_data.bin", "wb+");
if (fp == NULL)
return EXIT_FAILURE;
- fwrite(iRf, sizeof(R), (size_t)(N * N), fp);
+ fwrite(iRf, sizeof(NFFT_R), (size_t)(N * N), fp);
fclose(fp);
/** free the variables */
diff --git a/applications/radon/radon.c b/applications/radon/radon.c.in
similarity index 69%
rename from applications/radon/radon.c
rename to applications/radon/radon.c.in
index 30a8640..1956428 100644
--- a/applications/radon/radon.c
+++ b/applications/radon/radon.c.in
@@ -37,41 +37,39 @@
* \author Markus Fenn
* \date 2005
*/
-#include "config.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
/** define weights of kernel function for discrete Radon transform */
/*#define KERNEL(r) 1.0 */
-#define KERNEL(r) (K(1.0)-FABS((R)(r))/((R)S/2))
+#define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))
/** generates the points x with weights w
* for the polar grid with T angles and R offsets
*/
-static int polar_grid(int T, int S, R *x, R *w)
+static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
{
int t, r;
- R W = (R) T * (((R) S / K(2.0)) * ((R) S / K(2.0)) + K(1.0) / K(4.0));
+ NFFT_R W = (NFFT_R) T * (((NFFT_R) S / NFFT_K(2.0)) * ((NFFT_R) S / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
for (t = -T / 2; t < T / 2; t++)
{
for (r = -S / 2; r < S / 2; r++)
{
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (R) r / (R)(S) * COS(KPI * (R)(t) / (R)(T));
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (R) r / (R)(S) * SIN(KPI * (R)(t) / (R)(T));
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
if (r == 0)
- w[(t + T / 2) * S + (r + S / 2)] = K(1.0) / K(4.0) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
else
- w[(t + T / 2) * S + (r + S / 2)] = FABS((R) r) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
}
}
@@ -81,10 +79,10 @@ static int polar_grid(int T, int S, R *x, R *w)
/** generates the points x with weights w
* for the linogram grid with T slopes and R offsets
*/
-static int linogram_grid(int T, int S, R *x, R *w)
+static int linogram_grid(int T, int S, NFFT_R *x, NFFT_R *w)
{
int t, r;
- R W = (R) T * (((R) S / K(2.0)) * ((R) S / K(2.0)) + K(1.0) / K(4.0));
+ NFFT_R W = (NFFT_R) T * (((NFFT_R) S / NFFT_K(2.0)) * ((NFFT_R) S / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
for (t = -T / 2; t < T / 2; t++)
{
@@ -92,20 +90,20 @@ static int linogram_grid(int T, int S, R *x, R *w)
{
if (t < 0)
{
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (R) r / (R)(S);
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = K(4.0) * ((R)(t) + (R)(T) / K(4.0)) / (R)(T) * (R)(r)
- / (R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = NFFT_K(4.0) * ((NFFT_R)(t) + (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T) * (NFFT_R)(r)
+ / (NFFT_R)(S);
}
else
{
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = -K(4.0) * ((R)(t) - (R)(T) / K(4.0)) / (R)(T)
- * (R)(r) / (R)(S);
- x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (R) r / (R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = -NFFT_K(4.0) * ((NFFT_R)(t) - (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
+ * (NFFT_R)(r) / (NFFT_R)(S);
+ x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S);
}
if (r == 0)
- w[(t + T / 2) * S + (r + S / 2)] = K(1.0) / K(4.0) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
else
- w[(t + T / 2) * S + (r + S / 2)] = FABS((R) r) / W;
+ w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
}
}
@@ -115,16 +113,16 @@ static int linogram_grid(int T, int S, R *x, R *w)
/** computes the NFFT-based discrete Radon transform of f
* on the grid given by gridfcn() with T angles and R offsets
*/
-static int Radon_trafo(int (*gridfcn)(), int T, int S, R *f, int NN, R *Rf)
+static int Radon_trafo(int (*gridfcn)(), int T, int S, NFFT_R *f, int NN, NFFT_R *Rf)
{
int j, k; /**< index for nodes and freqencies */
NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */
- C *fft; /**< variable for the fftw-1Ds */
+ NFFT_C *fft; /**< variable for the fftw-1Ds */
FFTW(plan) my_fftw_plan; /**< plan for the fftw-1Ds */
int t, r; /**< index for directions and offsets */
- R *x, *w; /**< knots and associated weights */
+ NFFT_R *x, *w; /**< knots and associated weights */
int N[2], n[2];
int M = T * S;
@@ -134,14 +132,14 @@ static int Radon_trafo(int (*gridfcn)(), int T, int S, R *f, int NN, R *Rf)
N[1] = NN;
n[1] = 2 * N[1];
- fft = (C *) NFFT(malloc)((size_t)(S) * sizeof(C));
+ fft = (NFFT_C *) NFFT(malloc)((size_t)(S) * sizeof(NFFT_C));
my_fftw_plan = FFTW(plan_dft_1d)(S, fft, fft, FFTW_BACKWARD, FFTW_MEASURE);
- x = (R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(R)));
+ x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
if (x == NULL)
return EXIT_FAILURE;
- w = (R *) NFFT(malloc)((size_t)(T * S) * (sizeof(R)));
+ w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
if (w == NULL)
return EXIT_FAILURE;
@@ -171,7 +169,7 @@ static int Radon_trafo(int (*gridfcn)(), int T, int S, R *f, int NN, R *Rf)
/** init Fourier coefficients from given image */
for (k = 0; k < my_nfft_plan.N_total; k++)
- my_nfft_plan.f_hat[k] = f[k] + II * K(0.0);
+ my_nfft_plan.f_hat[k] = f[k] + _Complex_I * NFFT_K(0.0);
/** NFFT-2D */
NFFT(trafo)(&my_nfft_plan);
@@ -179,7 +177,7 @@ static int Radon_trafo(int (*gridfcn)(), int T, int S, R *f, int NN, R *Rf)
/** FFTW-1Ds */
for (t = 0; t < T; t++)
{
- fft[0] = K(0.0);
+ fft[0] = NFFT_K(0.0);
for (r = -S / 2 + 1; r < S / 2; r++)
fft[r + S / 2] = KERNEL(r) * my_nfft_plan.f[t * S + (r + S / 2)];
@@ -188,12 +186,12 @@ static int Radon_trafo(int (*gridfcn)(), int T, int S, R *f, int NN, R *Rf)
NFFT(fftshift_complex_int)(fft, 1, &S);
for (r = 0; r < S; r++)
- Rf[t * S + r] = CREAL(fft[r]) / (R)(S);
+ Rf[t * S + r] = NFFT_M(creal)(fft[r]) / (NFFT_R)(S);
/* for(r=0; r<R/2; r++)
- Rf[t*R+(r+R/2)] = creal(cexp(-I*KPI*r)*fft[r]);
+ Rf[t*R+(r+R/2)] = creal(cexp(-I*NFFT_KPI*r)*fft[r]);
for(r=0; r<R/2; r++)
- Rf[t*R+r] = creal(cexp(-I*KPI*r)*fft[r+R/2]);
+ Rf[t*R+r] = creal(cexp(-I*NFFT_KPI*r)*fft[r+R/2]);
*/
}
@@ -214,7 +212,7 @@ int main(int argc, char **argv)
int T, S; /**< number of directions/offsets */
FILE *fp;
int N; /**< image size */
- R *f, *Rf;
+ NFFT_R *f, *Rf;
if (argc != 5)
{
@@ -237,14 +235,14 @@ int main(int argc, char **argv)
S = atoi(argv[4]);
/*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
- f = (R *) NFFT(malloc)((size_t)(N * N) * (sizeof(R)));
- Rf = (R *) NFFT(malloc)((size_t)(T * S) * (sizeof(R)));
+ f = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (sizeof(NFFT_R)));
+ Rf = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
/** load data */
fp = fopen("input_data.bin", "rb");
if (fp == NULL)
return EXIT_FAILURE;
- fread(f, sizeof(R), (size_t)(N * N), fp);
+ fread(f, sizeof(NFFT_R), (size_t)(N * N), fp);
fclose(fp);
/** Radon transform */
@@ -254,7 +252,7 @@ int main(int argc, char **argv)
fp = fopen("sinogram_data.bin", "wb+");
if (fp == NULL)
return EXIT_FAILURE;
- fwrite(Rf, sizeof(R), (size_t)(T * S), fp);
+ fwrite(Rf, sizeof(NFFT_R), (size_t)(T * S), fp);
fclose(fp);
/** free the variables */
diff --git a/configure.ac b/configure.ac
index 0eb5183..e42b02e 100644
--- a/configure.ac
+++ b/configure.ac
@@ -124,6 +124,16 @@ case "$PRECISION" in
esac
AC_SUBST(PREC_SUFFIX)
+# Library suffix
+case "$PRECISION" in
+ s) NFFT_PRECISION_MACRO=NFFT_PRECISION_SINGLE;;
+ d) NFFT_PRECISION_MACRO=NFFT_PRECISION_DOUBLE;;
+ l) NFFT_PRECISION_MACRO=NFFT_PRECISION_LONG_DOUBLE;;
+esac
+AC_SUBST(NFFT_PRECISION_MACRO)
+
+AM_CONDITIONAL(HAVE_NON_DOUBLE_PRECISION, test "x$PRECISION" != "xd" )
+
# enable or disable parts of NFFT
# whether we need the fpt module
@@ -161,8 +171,10 @@ AX_NFFT_MODULE([fpt],[FPT],[fast polynomial transform],["no"],[],[],
[test "x$ok" = "xyes" -o "x$need_fpt" = "xyes"])
# multithreaded code
+#AC_ARG_ENABLE(openmp, [AC_HELP_STRING([--enable-openmp],
+# [enable OpenMP multithreaded code])], [enable_threads=$enableval; AC_DEFINE(ENABLE_OPENMP, 1, ["Define to enable OpenMP code."])], enable_threads=no)
AC_ARG_ENABLE(openmp, [AC_HELP_STRING([--enable-openmp],
- [enable OpenMP multithreaded code])], [enable_threads=$enableval; AC_DEFINE(ENABLE_OPENMP, 1, ["Define to enable OpenMP code."])], enable_threads=no)
+ [enable OpenMP multithreaded code])], enable_threads=$enableval, enable_threads=no)
AM_CONDITIONAL(HAVE_THREADS, test "x$enable_threads" = "xyes" )
# debug mode
@@ -548,15 +560,22 @@ AC_CONFIG_FILES(Makefile \
examples/fpt/Makefile \
examples/mri/Makefile \
examples/nfct/Makefile \
+ examples/nfct/simple_test.c \
examples/nfft/Makefile \
+ examples/nfft/simple_test.c \
+ examples/nfft/simple_test_threads.c \
examples/nfsft/Makefile \
examples/nfsoft/Makefile \
examples/nfst/Makefile \
+ examples/nfst/simple_test.c \
examples/nnfft/Makefile \
examples/nsfft/Makefile \
examples/solver/Makefile \
+ examples/solver/glacier.c \
+ examples/solver/simple_test.c \
applications/Makefile \
applications/fastgauss/Makefile \
+ applications/fastgauss/fastgauss.c \
applications/fastsum/Makefile \
applications/fastsumS2/Makefile \
applications/quadratureS2/Makefile \
@@ -564,8 +583,13 @@ AC_CONFIG_FILES(Makefile \
applications/mri/mri2d/Makefile \
applications/mri/mri3d/Makefile \
applications/polarFFT/Makefile \
+ applications/polarFFT/linogram_fft_test.c \
+ applications/polarFFT/mpolar_fft_test.c \
+ applications/polarFFT/polar_fft_test.c \
applications/radon/Makefile \
- applications/iterS2/Makefile
+ applications/radon/inverse_radon.c \
+ applications/radon/radon.c \
+ applications/iterS2/Makefile \
matlab/Makefile \
matlab/nfsft/Makefile \
matlab/nfft/Makefile \
diff --git a/examples/fpt/simple_test.c b/examples/fpt/simple_test.c
index f335100..c0be8d1 100644
--- a/examples/fpt/simple_test.c
+++ b/examples/fpt/simple_test.c
@@ -18,23 +18,18 @@
/* $Id$ */
-#include "config.h"
-
/* standard headers */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
-/* It is important to include complex.h before nfft3.h. */
-#ifdef HAVE_COMPLEX_H
+/* It is important to include complex.h before nfft3.h and fftw3.h. */
#include <complex.h>
-#endif
#include <fftw3.h>
/* NFFT3 header */
#include "nfft3.h"
-#include "infft.h"
int main(void)
{
@@ -113,7 +108,7 @@ int main(void)
printf("\n2) Random Fourier coefficients a_k, k=0,1,...,N:\n");
for (k = 0; k <= N; k++)
{
- a[k] = 2.0*X(drand48)() - 1.0; /* for debugging: use k+1 */
+ a[k] = 2.0*nfft_drand48() - 1.0; /* for debugging: use k+1 */
printf(" a_%-2d = %+5.3lE\n",k,creal(a[k]));
}
}
diff --git a/examples/nfct/simple_test.c b/examples/nfct/simple_test.c.in
similarity index 83%
rename from examples/nfct/simple_test.c
rename to examples/nfct/simple_test.c.in
index 3b738da..e264f65 100644
--- a/examples/nfct/simple_test.c
+++ b/examples/nfct/simple_test.c.in
@@ -23,13 +23,16 @@
#include <string.h>
#include <stdlib.h>
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
static void simple_test_nfct_1d(void)
{
NFCT(plan) p;
+ const char *error_str;
+
int N = 14;
int M = 19;
@@ -37,16 +40,24 @@ static void simple_test_nfct_1d(void)
NFCT(init_1d)(&p,N,M);
/** init pseudo random nodes */
- NFFT(vrand_real)(p.x, p.M_total, K(0.0), K(0.5));
+ NFFT(vrand_real)(p.x, p.M_total, NFFT_K(0.0), NFFT_K(0.5));
/** precompute psi, the entries of the matrix B */
if( p.flags & PRE_ONE_PSI)
NFCT(precompute_one_psi)(&p);
/** init pseudo random Fourier coefficients and show them */
- NFFT(vrand_real)(p.f_hat, p.N_total, K(0.0), K(1.0));
+ NFFT(vrand_real)(p.f_hat, p.N_total, NFFT_K(0.0), NFFT_K(1.0));
NFFT(vpr_double)(p.f_hat,p.N_total,"given Fourier coefficients, vector f_hat");
+ /** check for valid parameters before calling any trafo/adjoint method */
+ error_str = NFCT(check)(&p);
+ if (error_str != 0)
+ {
+ printf("Error in nfct module: %s\n", error_str);
+ return;
+ }
+
/** direct trafo and show the result */
NFCT(trafo_direct)(&p);
NFFT(vpr_double)(p.f,p.M_total,"ndct, vector f");
diff --git a/examples/nfft/simple_test.c b/examples/nfft/simple_test.c.in
similarity index 74%
rename from examples/nfft/simple_test.c
rename to examples/nfft/simple_test.c.in
index f5b61a0..9d6f944 100644
--- a/examples/nfft/simple_test.c
+++ b/examples/nfft/simple_test.c.in
@@ -16,19 +16,15 @@
* Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
-/* $Id$ */
-#include "config.h"
-
+/* $Id: simple_test.c.in 4298 2015-01-15 10:24:37Z tovo $ */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
-#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
static void simple_test_nfft_1d(void)
{
@@ -37,6 +33,8 @@ static void simple_test_nfft_1d(void)
int N = 14;
int M = 19;
+ const char *error_str;
+
/** init an one dimensional plan */
NFFT(init_1d)(&p, N, M);
@@ -51,6 +49,14 @@ static void simple_test_nfft_1d(void)
NFFT(vrand_unit_complex)(p.f_hat,p.N_total);
NFFT(vpr_complex)(p.f_hat, p.N_total, "given Fourier coefficients, vector f_hat");
+ /** check for valid parameters before calling any trafo/adjoint method */
+ error_str = NFFT(check)(&p);
+ if (error_str != 0)
+ {
+ printf("Error in nfft module: %s\n", error_str);
+ return;
+ }
+
/** direct trafo and show the result */
NFFT(trafo_direct)(&p);
NFFT(vpr_complex)(p.f,p.M_total,"ndft, vector f");
@@ -74,17 +80,18 @@ static void simple_test_nfft_1d(void)
static void simple_test_nfft_2d(void)
{
int K, N[2], n[2], M;
- R t;
- ticks t0, t1;
+ NFFT_R t0, t1;
NFFT(plan) p;
+ const char *error_str;
+
N[0] = 32; n[0] = 64;
N[1] = 14; n[1] = 32;
M = N[0] * N[1];
K = 16;
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
/** init a two dimensional plan */
NFFT(init_guru)(&p, 2, N, M, n, 7,
PRE_PHI_HUT| PRE_FULL_PSI| MALLOC_F_HAT| MALLOC_X| MALLOC_F |
@@ -101,42 +108,45 @@ static void simple_test_nfft_2d(void)
/** init pseudo random Fourier coefficients and show them */
NFFT(vrand_unit_complex)(p.f_hat, p.N_total);
- t1 = getticks();
- t = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
NFFT(vpr_complex)(p.f_hat,K, "given Fourier coefficients, vector f_hat (first few entries)");
- printf(" ... initialisation took %.2" __FES__ " seconds.\n",t);
+ printf(" ... initialisation took %.2" NFFT__FES__ " seconds.\n",t1-t0);
+
+ /** check for valid parameters before calling any trafo/adjoint method */
+ error_str = NFFT(check)(&p);
+ if (error_str != 0)
+ {
+ printf("Error in nfft module: %s\n", error_str);
+ return;
+ }
/** direct trafo and show the result */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
NFFT(trafo_direct)(&p);
- t1 = getticks();
- t = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
NFFT(vpr_complex)(p.f, K, "ndft, vector f (first few entries)");
- printf(" took %.2" __FES__ " seconds.\n",t);
+ printf(" took %.2" NFFT__FES__ " seconds.\n",t1-t0);
/** approx. trafo and show the result */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
NFFT(trafo)(&p);
- t1 = getticks();
- t = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
NFFT(vpr_complex)(p.f, K, "nfft, vector f (first few entries)");
- printf(" took %.2" __FES__ " seconds.\n",t);
+ printf(" took %.2" NFFT__FES__ " seconds.\n",t1-t0);
/** direct adjoint and show the result */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
NFFT(adjoint_direct)(&p);
- t1 = getticks();
- t = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
NFFT(vpr_complex)(p.f_hat, K, "adjoint ndft, vector f_hat (first few entries)");
- printf(" took %.2" __FES__ " seconds.\n",t);
+ printf(" took %.2" NFFT__FES__ " seconds.\n",t1-t0);
/** approx. adjoint and show the result */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
NFFT(adjoint)(&p);
- t1 = getticks();
- t = NFFT(elapsed_seconds)(t1, t0);
+ t1 = NFFT(clock_gettime_seconds)();
NFFT(vpr_complex)(p.f_hat, K, "adjoint nfft, vector f_hat (first few entries)");
- printf(" took %.2" __FES__ " seconds.\n",t);
+ printf(" took %.2" NFFT__FES__ " seconds.\n",t1-t0);
/** finalise the two dimensional plan */
NFFT(finalize)(&p);
diff --git a/examples/nfft/simple_test_threads.c b/examples/nfft/simple_test_threads.c.in
similarity index 71%
rename from examples/nfft/simple_test_threads.c
rename to examples/nfft/simple_test_threads.c.in
index c7ad6b6..1688226 100644
--- a/examples/nfft/simple_test_threads.c
+++ b/examples/nfft/simple_test_threads.c.in
@@ -16,31 +16,27 @@
* Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
-/* $Id: simple_test.c 3198 2009-05-27 14:16:50Z keiner $ */
-#include "config.h"
-
+/* $Id$ */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
-#include <complex.h>
-#endif
#include <omp.h>
-#include "nfft3.h"
-#include "infft.h"
-//#include <time.h>
+#include <sys/time.h>
+
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
int main(void)
{
NFFT(plan) p;
const int N = 1000000;
const int M = 1000000;
- ticks t0, t1;
- double t;
+ NFFT_R t0, t1;
- printf("nthreads = %d\n", NFFT(get_num_threads)());
+ printf("nthreads = %td\n", NFFT(get_num_threads)());
/* init */
FFTW(init_threads)();
@@ -50,24 +46,21 @@ int main(void)
NFFT(vrand_shifted_unit_double)(p.x,p.M_total);
/* precompute psi, that is, the entries of the matrix B */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
if(p.flags & PRE_ONE_PSI)
NFFT(precompute_one_psi)(&p);
- t1 = getticks();
- t = NFFT(elapsed_seconds(t1,t0));
- fprintf(stderr,"precompute elapsed time: %.3f seconds\n",t);
+ t1 = NFFT(clock_gettime_seconds)();
+ fprintf(stderr,"precompute elapsed time: %.3" NFFT__FIS__ " seconds\n",t1-t0);
/* pseudo random Fourier coefficients */
NFFT(vrand_unit_complex)(p.f_hat,p.N_total);
/* transformation */
- t0 = getticks();
+ t0 = NFFT(clock_gettime_seconds)();
NFFT(trafo)(&p);
- t1 = getticks();
- t = NFFT(elapsed_seconds)(t1,t0);
- fprintf(stderr,"compute elapsed time: %.3f seconds\n",t);
+ t1 = NFFT(clock_gettime_seconds)();
+ fprintf(stderr,"compute elapsed time: %.3" NFFT__FIS__ " seconds\n",t1-t0);
fflush(stderr);
-// NFFT(vpr_complex)(p.f,p.M_total,"ndft, vector f");
/* cleanup */
NFFT(finalize)(&p);
diff --git a/examples/nfsft/simple_test.c b/examples/nfsft/simple_test.c
index 339f442..a8edcd2 100644
--- a/examples/nfsft/simple_test.c
+++ b/examples/nfsft/simple_test.c
@@ -18,20 +18,18 @@
/* $Id$ */
-#include "config.h"
-
/* standard headers */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
/* It is important to include complex.h before nfft3.h. */
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h" /* NFFT3 header */
-#include "infft.h" /* NFFT3 internal header */
+
+#define __FES__ "E"
+#define K(x) ((double) x)
static void simple_test_nfsft(void)
{
@@ -56,8 +54,8 @@ static void simple_test_nfsft(void)
/* pseudo-random nodes */
for (j = 0; j < plan.M_total; j++)
{
- plan.x[2*j]= X(drand48)() - K(0.5);
- plan.x[2*j+1]= K(0.5) * X(drand48)();
+ plan.x[2*j]= nfft_drand48() - K(0.5);
+ plan.x[2*j+1]= K(0.5) * nfft_drand48();
}
/* precomputation (for NFFT, node-dependent) */
@@ -67,7 +65,7 @@ static void simple_test_nfsft(void)
for (k = 0; k <= plan.N; k++)
for (n = -k; n <= k; n++)
plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
- X(drand48)() - K(0.5) + _Complex_I*(X(drand48)() - K(0.5));
+ nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
/* Direct transformation, display result. */
nfsft_trafo_direct(&plan);
diff --git a/examples/nfsft/simple_test_threads.c b/examples/nfsft/simple_test_threads.c
index 7a32178..6ca10aa 100644
--- a/examples/nfsft/simple_test_threads.c
+++ b/examples/nfsft/simple_test_threads.c
@@ -18,21 +18,19 @@
/* $Id: simple_test.c 3498 2010-05-07 18:46:08Z keiner $ */
-#include "config.h"
-
/* standard headers */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
/* It is important to include complex.h before nfft3.h. */
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include <omp.h>
#include "nfft3.h" /* NFFT3 header */
-#include "infft.h" /* NFFT3 internal header */
+
+#define __FES__ "E"
+#define K(x) ((double) x)
static void simple_test_nfsft(void)
{
@@ -57,8 +55,8 @@ static void simple_test_nfsft(void)
/* pseudo-random nodes */
for (j = 0; j < plan.M_total; j++)
{
- plan.x[2*j]= X(drand48)() - K(0.5);
- plan.x[2*j+1]= K(0.5) * X(drand48)();
+ plan.x[2*j]= nfft_drand48() - K(0.5);
+ plan.x[2*j+1]= K(0.5) * nfft_drand48();
}
/* precomputation (for NFFT, node-dependent) */
@@ -68,7 +66,7 @@ static void simple_test_nfsft(void)
for (k = 0; k <= plan.N; k++)
for (n = -k; n <= k; n++)
plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
- X(drand48)() - K(0.5) + _Complex_I*(X(drand48)() - K(0.5));
+ nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
/* Direct transformation, display result. */
nfsft_trafo_direct(&plan);
@@ -120,7 +118,7 @@ static void simple_test_nfsft(void)
int main(void)
{
- printf("nthreads = %d\n", X(get_num_threads)());
+ printf("nthreads = %d\n", nfft_get_num_threads());
/* init */
fftw_init_threads();
diff --git a/examples/nfsoft/simple_test.c b/examples/nfsoft/simple_test.c
index f9a31c9..094a4c7 100644
--- a/examples/nfsoft/simple_test.c
+++ b/examples/nfsoft/simple_test.c
@@ -18,28 +18,22 @@
/* $Id$ */
-#include "config.h"
-
/* Include standard C headers. */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-/* Include NFFT 3 utilities headers. */
/* Include NFFT3 library header. */
#include "nfft3.h"
-#include "infft.h"
static void simple_test_nfsoft(int bw, int M)
{
nfsoft_plan plan_nfsoft; /**< Plan for the NFSOFT */
nfsoft_plan plan_ndsoft; /**< Plan for the NDSOFT */
- ticks t0, t1;
+ double t0, t1;
int j; /** just an index*/
int k, m; /** the two parameters controlling the accuracy of the NFSOFT*/
double d1, d2, d3; /** indeces for initializing the Euler angles*/
@@ -67,9 +61,9 @@ static void simple_test_nfsoft(int bw, int M)
/** Init random nodes (for both plans, the same). */
for (j = 0; j < plan_nfsoft.M_total; j++)
{
- d1 = ((R) rand()) / RAND_MAX - 0.5;
- d2 = 0.5 * ((R) rand()) / RAND_MAX;
- d3 = ((R) rand()) / RAND_MAX - 0.5;
+ d1 = ((double) rand()) / RAND_MAX - 0.5;
+ d2 = 0.5 * ((double) rand()) / RAND_MAX;
+ d3 = ((double) rand()) / RAND_MAX - 0.5;
plan_nfsoft.x[3* j ] = d1; /**alpha*/
plan_nfsoft.x[3* j + 1] = d2; /**beta*/
@@ -83,8 +77,8 @@ static void simple_test_nfsoft(int bw, int M)
/** init random Fourier coefficients (again the same for both plans) and display them*/
for (j = 0; j < (bw + 1) * (4* (bw +1)*(bw+1)-1)/3;j++)
{
- d1=((R)rand())/RAND_MAX - 0.5;
- d2=((R)rand())/RAND_MAX - 0.5;
+ d1=((double)rand())/RAND_MAX - 0.5;
+ d2=((double)rand())/RAND_MAX - 0.5;
plan_nfsoft.f_hat[j]=d1 + I*d2;
plan_ndsoft.f_hat[j]=d1 + I*d2;
}
@@ -102,10 +96,10 @@ static void simple_test_nfsoft(int bw, int M)
/** Compute NFSOFT and display the time needed. */
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nfsoft_trafo(&plan_nfsoft);
- t1 = getticks();
- time = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ time = t1-t0;
if (plan_nfsoft.M_total<=20)
nfft_vpr_complex(plan_nfsoft.f,plan_nfsoft.M_total,"NFSOFT, function samples");
else
@@ -113,10 +107,10 @@ static void simple_test_nfsoft(int bw, int M)
printf(" computed in %11le seconds\n",time);
/** Compute NDSOFT and display the time needed. */
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nfsoft_trafo(&plan_ndsoft);
- t1 = getticks();
-time = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ time = t1-t0;
if (plan_ndsoft.M_total<=20)
nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total,"NDSOFT, function samples");
else
@@ -124,7 +118,7 @@ time = nfft_elapsed_seconds(t1,t0);
printf(" computed in %11le seconds\n",time);
/**compute the error between the NFSOFT and NDSOFT and display it*/
- error= X(error_l_infty_complex)(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
+ error= nfft_error_l_infty_complex(plan_ndsoft.f,plan_nfsoft.f, plan_nfsoft.M_total);
printf("\n The NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
printf("\n---------------------------------------------\n");
@@ -134,10 +128,10 @@ time = nfft_elapsed_seconds(t1,t0);
nfft_vpr_complex(plan_ndsoft.f,plan_ndsoft.M_total, "function samples to start adjoint trafo");
/** Compute the adjoint NFSOFT and display the time needed.*/
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nfsoft_adjoint(&plan_nfsoft);
- t1 = getticks();
-time = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ time = t1-t0;
if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
nfft_vpr_complex(plan_nfsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
else
@@ -145,10 +139,10 @@ time = nfft_elapsed_seconds(t1,t0);
printf(" computed in %11le seconds\n",time);
/** Compute adjoint NDSOFT and display the time needed.*/
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nfsoft_adjoint(&plan_ndsoft);
- t1 = getticks();
-time = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ time = t1-t0;
if ((bw+1)*(4*(bw+1)*(bw+1)-1)/3<=20)
nfft_vpr_complex(plan_ndsoft.f_hat,(bw+1)*(4*(bw+1)*(bw+1)-1)/3,"SO(3) Fourier coefficients");
else
@@ -157,7 +151,7 @@ time = nfft_elapsed_seconds(t1,t0);
/**compute the error between the adjoint NFSOFT and NDSOFT and display it*/
- error=X(error_l_infty_complex)(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
+ error=nfft_error_l_infty_complex(plan_ndsoft.f_hat,plan_nfsoft.f_hat, (bw+1)*(4*(bw+1)*(bw+1)-1)/3);
printf("\n The adjoint NFSOFT of bandwidth=%d for %d rotations has infty-error %11le \n",bw, M,error);
printf("\n---------------------------------------------\n");
diff --git a/examples/nfst/simple_test.c b/examples/nfst/simple_test.c.in
similarity index 83%
rename from examples/nfst/simple_test.c
rename to examples/nfst/simple_test.c.in
index 286e166..23c85f4 100644
--- a/examples/nfst/simple_test.c
+++ b/examples/nfst/simple_test.c.in
@@ -23,13 +23,16 @@
#include <string.h>
#include <stdlib.h>
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
static void simple_test_nfst_1d(void)
{
NFST(plan) p;
+ const char *error_str;
+
int N = 14;
int M = 19;
@@ -37,16 +40,24 @@ static void simple_test_nfst_1d(void)
NFST(init_1d)(&p,N,M);
/** init pseudo random nodes */
- NFFT(vrand_real)(p.x, p.M_total, K(0.0), K(0.5));
+ NFFT(vrand_real)(p.x, p.M_total, NFFT_K(0.0), NFFT_K(0.5));
/** precompute psi, the entries of the matrix B */
if( p.flags & PRE_ONE_PSI)
NFST(precompute_one_psi)(&p);
/** init pseudo random Fourier coefficients and show them */
- NFFT(vrand_real)(p.f_hat, p.N_total, K(0.0), K(1.0));
+ NFFT(vrand_real)(p.f_hat, p.N_total, NFFT_K(0.0), NFFT_K(1.0));
NFFT(vpr_double)(p.f_hat,p.N_total,"given Fourier coefficients, vector f_hat");
+ /** check for valid parameters before calling any trafo/adjoint method */
+ error_str = NFST(check)(&p);
+ if (error_str != 0)
+ {
+ printf("Error in nfst module: %s\n", error_str);
+ return;
+ }
+
/** direct trafo and show the result */
NFST(trafo_direct)(&p);
NFFT(vpr_double)(p.f,p.M_total,"ndst, vector f");
diff --git a/examples/nnfft/accuracy.c b/examples/nnfft/accuracy.c
index 8fecb04..32eb91a 100644
--- a/examples/nnfft/accuracy.c
+++ b/examples/nnfft/accuracy.c
@@ -18,16 +18,15 @@
/* $Id$ */
-#include "config.h"
-
#include <math.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
+
+/** Swap two vectors. */
+#define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
+ NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
void accuracy(int d)
{
@@ -86,8 +85,8 @@ void accuracy(int d)
nnfft_trafo(&my_plan);
printf("%e, %e\n",
- X(error_l_infty_complex)(slow, my_plan.f, M_total),
- X(error_l_infty_1_complex)(slow, my_plan.f, M_total, my_plan.f_hat,
+ nfft_error_l_infty_complex(slow, my_plan.f, M_total),
+ nfft_error_l_infty_1_complex(slow, my_plan.f, M_total, my_plan.f_hat,
my_plan.N_total));
/** finalise the one dimensional plan */
diff --git a/examples/nnfft/simple_test.c b/examples/nnfft/simple_test.c
index 900910c..a56cc4c 100644
--- a/examples/nnfft/simple_test.c
+++ b/examples/nnfft/simple_test.c
@@ -17,16 +17,16 @@
*/
/* $Id$ */
-#include "config.h"
#include <stdlib.h>
#include <math.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
+
+/** Swap two vectors. */
+#define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
+ NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
void simple_test_nnfft_1d(void)
{
@@ -267,7 +267,7 @@ static void measure_time_nnfft_1d(void)
nnfft_plan my_plan; /**< plan for the nfft */
int my_N,N=4;
double t;
- ticks t0, t1;
+ double t0, t1;
for(my_N=16; my_N<=16384; my_N*=2)
{
@@ -291,16 +291,16 @@ static void measure_time_nnfft_1d(void)
for(k=0;k<my_plan.N_total;k++)
my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nnfft_trafo_direct(&my_plan);
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = t1 - t0;
printf("t_nndft=%e,\t",t);
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nnfft_trafo(&my_plan);
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = t1 - t0;
printf("t_nnfft=%e\t",t);
printf("(N=M=%d)\n",my_N);
diff --git a/examples/nsfft/nsfft_test.c b/examples/nsfft/nsfft_test.c
index 1cab619..a005eec 100644
--- a/examples/nsfft/nsfft_test.c
+++ b/examples/nsfft/nsfft_test.c
@@ -17,18 +17,18 @@
*/
/* $Id$ */
-#include "config.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
+
+/** Swap two vectors. */
+#define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
+ NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
static void accuracy_nsfft(int d, int J, int M, int m)
{
@@ -53,7 +53,7 @@ static void accuracy_nsfft(int d, int J, int M, int m)
nsfft_trafo(&p);
printf("%5d\t %+.5E\t",J,
- X(error_l_infty_1_complex)(swap_sndft_trafo, p.f, p.M_total,
+ nfft_error_l_infty_1_complex(swap_sndft_trafo, p.f, p.M_total,
p.f_hat, p.N_total));
fflush(stdout);
@@ -68,7 +68,7 @@ static void accuracy_nsfft(int d, int J, int M, int m)
nsfft_adjoint(&p);
printf("%+.5E\n",
- X(error_l_infty_1_complex)(swap_sndft_adjoint, p.f_hat,
+ nfft_error_l_infty_1_complex(swap_sndft_adjoint, p.f_hat,
p.N_total,
p.f, p.M_total));
fflush(stdout);
@@ -84,14 +84,14 @@ static void time_nsfft(int d, int J, int M, unsigned test_nsdft, unsigned test_n
{
int r, N[d], n[d];
double t, t_nsdft, t_nfft, t_nsfft;
- ticks t0, t1;
+ double t0, t1;
nsfft_plan p;
nfft_plan np;
for(r=0;r<d;r++)
{
- N[r]= X(exp2i)(J+2);
+ N[r]= nfft_exp2i(J+2);
n[r]=(3*N[r])/2;
/*n[r]=2*N[r];*/
}
@@ -108,10 +108,10 @@ static void time_nsfft(int d, int J, int M, unsigned test_nsdft, unsigned test_n
while(t_nsdft<0.1)
{
r++;
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nsfft_trafo_direct(&p);
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = (t1 - t0);
t_nsdft+=t;
}
t_nsdft/=r;
@@ -133,10 +133,10 @@ static void time_nsfft(int d, int J, int M, unsigned test_nsdft, unsigned test_n
while(t_nfft<0.1)
{
r++;
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nfft_trafo(&np);
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = (t1 - t0);
t_nfft+=t;
}
t_nfft/=r;
@@ -153,10 +153,10 @@ static void time_nsfft(int d, int J, int M, unsigned test_nsdft, unsigned test_n
while(t_nsfft<0.1)
{
r++;
- t0 = getticks();
+ t0 = nfft_clock_gettime_seconds();
nsfft_trafo(&p);
- t1 = getticks();
- t = nfft_elapsed_seconds(t1,t0);
+ t1 = nfft_clock_gettime_seconds();
+ t = (t1 - t0);
t_nsfft+=t;
}
t_nsfft/=r;
@@ -198,9 +198,9 @@ int main(int argc,char **argv)
for(J=atoi(argv[3]); J<=atoi(argv[4]); J++)
{
if(d==2)
- M=(J+4)*X(exp2i)(J+1);
+ M=(J+4)*nfft_exp2i(J+1);
else
- M=6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+X(exp2i)(3*(J/2+1));
+ M=6*nfft_exp2i(J)*(nfft_exp2i((J+1)/2+1)-1)+nfft_exp2i(3*(J/2+1));
if(d*(J+2)<=24)
time_nsfft(d, J, M, 1, 1);
diff --git a/examples/nsfft/simple_test.c b/examples/nsfft/simple_test.c
index 74658e2..3f751b1 100644
--- a/examples/nsfft/simple_test.c
+++ b/examples/nsfft/simple_test.c
@@ -17,18 +17,14 @@
*/
/* $Id$ */
-#include "config.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
#include "nfft3.h"
-#include "infft.h"
static void simple_test_nsfft(int d, int J, int M)
{
@@ -69,7 +65,7 @@ int main(int argc,char **argv)
printf("1) computing a two dimensional nsdft, nsfft and adjoints\n\n");
d=2;
J=5;
- M=(J+4)*X(exp2i)(J+1);
+ M=(J+4)*nfft_exp2i(J+1);
simple_test_nsfft(d,J,M);
getc(stdin);
@@ -77,7 +73,7 @@ int main(int argc,char **argv)
printf("2) computing a three dimensional nsdft, nsfft and adjoints\n\n");
d=3;
J=5;
- M=6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+X(exp2i)(3*(J/2+1));
+ M=6*nfft_exp2i(J)*(nfft_exp2i((J+1)/2+1)-1)+nfft_exp2i(3*(J/2+1));
simple_test_nsfft(d,J,M);
return 1;
diff --git a/examples/solver/glacier.c b/examples/solver/glacier.c.in
similarity index 74%
rename from examples/solver/glacier.c
rename to examples/solver/glacier.c.in
index 0a7de73..bfd4769 100644
--- a/examples/solver/glacier.c
+++ b/examples/solver/glacier.c.in
@@ -17,18 +17,16 @@
*/
/* $Id$ */
-#include "config.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
/**
* \defgroup examples_solver_glacier Reconstruction of a glacier from \
@@ -38,16 +36,16 @@
*/
/** Generalised Sobolev weight */
-static R my_weight(R z, R a, R b, R c)
+static NFFT_R my_weight(NFFT_R z, NFFT_R a, NFFT_R b, NFFT_R c)
{
- return POW(K(0.25) - z * z, b) / (c + POW(FABS(z), K(2.0) * a));
+ return NFFT_M(pow)(NFFT_K(0.25) - z * z, b) / (c + NFFT_M(pow)(NFFT_M(fabs)(z), NFFT_K(2.0) * a));
}
/** Reconstruction routine */
static void glacier(int N, int M)
{
int j, k, k0, k1, l, my_N[2], my_n[2];
- R tmp_y;
+ NFFT_R tmp_y;
NFFT(plan) p;
SOLVER(plan_complex) ip;
FILE* fp;
@@ -72,7 +70,7 @@ static void glacier(int N, int M)
fp = fopen("input_data.dat", "r");
for (j = 0; j < p.M_total; j++)
{
- fscanf(fp, "%" __FES__ " %" __FES__ " %" __FES__, &p.x[2 * j + 0], &p.x[2 * j + 1], &tmp_y);
+ fscanf(fp, "%" NFFT__FES__ " %" NFFT__FES__ " %" NFFT__FES__, &p.x[2 * j + 0], &p.x[2 * j + 1], &tmp_y);
ip.y[j] = tmp_y;
}
fclose(fp);
@@ -85,25 +83,25 @@ static void glacier(int N, int M)
if (ip.flags & PRECOMPUTE_DAMP)
for (k0 = 0; k0 < p.N[0]; k0++)
for (k1 = 0; k1 < p.N[1]; k1++)
- ip.w_hat[k0 * p.N[1] + k1] = my_weight(((R) ((R)(k0) - (R)(p.N[0]) / K(2.0))) / ((R)(p.N[0])),
- K(0.5), K(3.0), K(0.001))
- * my_weight((((R)(k1) - (R)(p.N[1]) / K(2.0))) / ((R)(p.N[1])), K(0.5), K(3.0), K(0.001));
+ ip.w_hat[k0 * p.N[1] + k1] = my_weight(((NFFT_R) ((NFFT_R)(k0) - (NFFT_R)(p.N[0]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[0])),
+ NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001))
+ * my_weight((((NFFT_R)(k1) - (NFFT_R)(p.N[1]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[1])), NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001));
/* init some guess */
for (k = 0; k < p.N_total; k++)
- ip.f_hat_iter[k] = K(0.0);
+ ip.f_hat_iter[k] = NFFT_K(0.0);
/* inverse trafo */
SOLVER(before_loop_complex)(&ip);
for (l = 0; l < 40; l++)
{
- fprintf(stderr, "Residual ||r||=%" __FES__ ",\n", SQRT(ip.dot_r_iter));
+ fprintf(stderr, "Residual ||r||=%" NFFT__FES__ ",\n", NFFT_M(sqrt)(ip.dot_r_iter));
SOLVER(loop_one_step_complex)(&ip);
}
for (k = 0; k < p.N_total; k++)
- printf("%" __FES__ " %" __FES__ "\n", CREAL(ip.f_hat_iter[k]), CIMAG(ip.f_hat_iter[k]));
+ printf("%" NFFT__FES__ " %" NFFT__FES__ "\n", NFFT_M(creal)(ip.f_hat_iter[k]), NFFT_M(cimag)(ip.f_hat_iter[k]));
SOLVER(finalize_complex)(&ip);
NFFT(finalize)(&p);
@@ -113,10 +111,10 @@ static void glacier(int N, int M)
static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags)
{
int j, k, k0, k1, l, my_N[2], my_n[2];
- R tmp_y, r;
+ NFFT_R tmp_y, r;
NFFT(plan) p, cp;
SOLVER(plan_complex) ip;
- C* cp_y;
+ NFFT_C* cp_y;
FILE* fp;
int M_re = M - M_cv;
@@ -135,7 +133,7 @@ static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags)
SOLVER(init_advanced_complex)(&ip, (NFFT(mv_plan_complex)*) (&p), solver_flags);
/* initialise cp for validation */
- cp_y = (C*) NFFT(malloc)((size_t)(M) * sizeof(C));
+ cp_y = (NFFT_C*) NFFT(malloc)((size_t)(M) * sizeof(NFFT_C));
NFFT(init_guru)(&cp, 2, my_N, M, my_n, 6,
PRE_PHI_HUT | PRE_FULL_PSI |
MALLOC_X | MALLOC_F |
@@ -148,7 +146,7 @@ static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags)
fp = fopen("input_data.dat", "r");
for (j = 0; j < cp.M_total; j++)
{
- fscanf(fp, "%" __FES__ " %" __FES__ " %" __FES__, &cp.x[2 * j + 0], &cp.x[2 * j + 1], &tmp_y);
+ fscanf(fp, "%" NFFT__FES__ " %" NFFT__FES__ " %" NFFT__FES__, &cp.x[2 * j + 0], &cp.x[2 * j + 1], &tmp_y);
cp_y[j] = tmp_y;
}
fclose(fp);
@@ -173,13 +171,13 @@ static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags)
if (ip.flags & PRECOMPUTE_DAMP)
for (k0 = 0; k0 < p.N[0]; k0++)
for (k1 = 0; k1 < p.N[1]; k1++)
- ip.w_hat[k0 * p.N[1] + k1] = my_weight((((R)(k0) - (R)(p.N[0]) / K(2.0))) / ((R)(p.N[0])),
- K(0.5), K(3.0), K(0.001))
- * my_weight((((R)(k1) - (R)(p.N[1]) / K(2.0))) / ((R)(p.N[1])), K(0.5), K(3.0), K(0.001));
+ ip.w_hat[k0 * p.N[1] + k1] = my_weight((((NFFT_R)(k0) - (NFFT_R)(p.N[0]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[0])),
+ NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001))
+ * my_weight((((NFFT_R)(k1) - (NFFT_R)(p.N[1]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[1])), NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001));
/* init some guess */
for (k = 0; k < p.N_total; k++)
- ip.f_hat_iter[k] = K(0.0);
+ ip.f_hat_iter[k] = NFFT_K(0.0);
/* inverse trafo */
SOLVER(before_loop_complex)(&ip);
@@ -189,19 +187,19 @@ static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags)
//fprintf(stderr,"r=%1.2e, ",sqrt(ip.dot_r_iter)/M_re);
- CSWAP(p.f_hat, ip.f_hat_iter);
+ NFFT_CSWAP(p.f_hat, ip.f_hat_iter);
NFFT(trafo)(&p);
- CSWAP(p.f_hat, ip.f_hat_iter);
+ NFFT_CSWAP(p.f_hat, ip.f_hat_iter);
NFFT(upd_axpy_complex)(p.f, -1, ip.y, M_re);
- r = SQRT(NFFT(dot_complex)(p.f, M_re) / NFFT(dot_complex)(cp_y, M));
- fprintf(stderr, "r=%1.2" __FES__ ", ", r);
- printf("$%1.1" __FES__ "$ & ", r);
+ r = NFFT_M(sqrt)(NFFT(dot_complex)(p.f, M_re) / NFFT(dot_complex)(cp_y, M));
+ fprintf(stderr, "r=%1.2" NFFT__FES__ ", ", r);
+ printf("$%1.1" NFFT__FES__ "$ & ", r);
NFFT(trafo)(&cp);
NFFT(upd_axpy_complex)(&cp.f[M_re], -1, &cp_y[M_re], M_cv);
- r = SQRT(NFFT(dot_complex)(&cp.f[M_re], M_cv) / NFFT(dot_complex)(cp_y, M));
- fprintf(stderr, "r_1=%1.2" __FES__ "\t", r);
- printf("$%1.1" __FES__ "$ & ", r);
+ r = NFFT_M(sqrt)(NFFT(dot_complex)(&cp.f[M_re], M_cv) / NFFT(dot_complex)(cp_y, M));
+ fprintf(stderr, "r_1=%1.2" NFFT__FES__ "\t", r);
+ printf("$%1.1" NFFT__FES__ "$ & ", r);
NFFT(finalize)(&cp);
SOLVER(finalize_complex)(&ip);
diff --git a/examples/solver/simple_test.c b/examples/solver/simple_test.c.in
similarity index 87%
rename from examples/solver/simple_test.c
rename to examples/solver/simple_test.c.in
index 8162090..b4010f5 100644
--- a/examples/solver/simple_test.c
+++ b/examples/solver/simple_test.c.in
@@ -17,18 +17,16 @@
*/
/* $Id$ */
-#include "config.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
-#ifdef HAVE_COMPLEX_H
#include <complex.h>
-#endif
-#include "nfft3.h"
-#include "infft.h"
+#define @NFFT_PRECISION_MACRO@
+
+#include "nfft3mp.h"
/* void simple_test_infft_1d(int N, int M, int iter) */
/* { */
@@ -92,6 +90,7 @@ static void simple_test_solver_nfft_1d(int N, int M, int iter)
int k, l; /**< index for nodes, freqencies,iter*/
NFFT(plan) p; /**< plan for the nfft */
SOLVER(plan_complex) ip; /**< plan for the inverse nfft */
+ const char *error_str;
/** initialise an one dimensional plan */
NFFT(init_1d)(&p, N, M);
@@ -112,18 +111,26 @@ static void simple_test_solver_nfft_1d(int N, int M, int iter)
/** initialise some guess f_hat_0 and solve */
for (k = 0; k < p.N_total; k++)
- ip.f_hat_iter[k] = K(0.0);
+ ip.f_hat_iter[k] = NFFT_K(0.0);
NFFT(vpr_complex)(ip.f_hat_iter, p.N_total,
"Initial guess, vector f_hat_iter");
- CSWAP(ip.f_hat_iter, p.f_hat);
+ /** check for valid parameters before calling any trafo/adjoint method */
+ error_str = NFFT(check)(&p);
+ if (error_str != 0)
+ {
+ printf("Error in nfft module: %s\n", error_str);
+ return;
+ }
+
+ NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
NFFT(trafo)(&p);
NFFT(vpr_complex)(p.f, p.M_total, "Data fit, vector f");
- CSWAP(ip.f_hat_iter, p.f_hat);
+ NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
SOLVER(before_loop_complex)(&ip);
- printf("\n Residual r=%" __FES__ "\n", ip.dot_r_iter);
+ printf("\n Residual r=%" NFFT__FES__ "\n", ip.dot_r_iter);
for (l = 0; l < iter; l++)
{
@@ -132,12 +139,12 @@ static void simple_test_solver_nfft_1d(int N, int M, int iter)
NFFT(vpr_complex)(ip.f_hat_iter, p.N_total,
"Approximate solution, vector f_hat_iter");
- CSWAP(ip.f_hat_iter, p.f_hat);
+ NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
NFFT(trafo)(&p);
NFFT(vpr_complex)(p.f, p.M_total, "Data fit, vector f");
- CSWAP(ip.f_hat_iter, p.f_hat);
+ NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
- printf("\n Residual r=%" __FES__ "\n", ip.dot_r_iter);
+ printf("\n Residual r=%" NFFT__FES__ "\n", ip.dot_r_iter);
}
SOLVER(finalize_complex)(&ip);
@@ -149,7 +156,7 @@ int main(void)
{
printf("\n Computing a one dimensional inverse nfft\n");
- simple_test_solver_nfft_1d(8, 4, 5);
+ simple_test_solver_nfft_1d(16, 8, 9);
return EXIT_SUCCESS;
}
diff --git a/include/Makefile.am b/include/Makefile.am
index d580dbf..34a90fd 100644
--- a/include/Makefile.am
+++ b/include/Makefile.am
@@ -1,3 +1,3 @@
-include_HEADERS = nfft3.h
+include_HEADERS = nfft3.h nfft3mp.h
EXTRA_DIST = infft.h cycle.h api.h
diff --git a/include/infft.h b/include/infft.h
index 5fb9d37..47791a0 100644
--- a/include/infft.h
+++ b/include/infft.h
@@ -1452,16 +1452,11 @@ R Y(float_property)(float_property);
R Y(prod_real)(R *vec, INT d);
/* int.c: */
-INT Y(exp2i)(const INT a);
INT Y(log2i)(const INT m);
-INT Y(next_power_of_2)(const INT N);
void Y(next_power_of_2_exp)(const INT N, INT *N2, INT *t);
/* error.c: */
-R Y(error_l_infty_complex)(const C *x, const C *y, const INT n);
/* not used */ R Y(error_l_infty_double)(const R *x, const R *y, const INT n);
-R Y(error_l_infty_1_complex)(const C *x, const C *y, const INT n,
- const C *z, const INT m);
/* not used */ R Y(error_l_infty_1_double)(const R *x, const R *y, const INT n, const R *z,
const INT m);
R Y(error_l_2_complex)(const C *x, const C *y, const INT n);
@@ -1474,21 +1469,8 @@ void Y(sort_node_indices_radix_lsdf)(INT n, INT *keys0, INT *keys1, INT rhigh);
/* assert.c */
void Y(assertion_failed)(const char *s, int line, const char *file);
-/* rand.c */
-R Y(drand48)(void);
-void Y(srand48)(long int seed);
-/** Inits a vector of random complex numbers in \f$[0,1]\times[0,1]{\rm i}\f$.
- */
-void Y(vrand_unit_complex)(C *x, const INT n);
-/** Inits a vector of random double numbers in \f$[-1/2,1/2]\f$.
- */
-void Y(vrand_shifted_unit_double)(R *x, const INT n);
-void Y(vrand_real)(R *x, const INT n, const R a, const R b);
-
/* vector1.c */
/** Computes the inner/dot product \f$x^H x\f$. */
-R Y(dot_complex)(C *x, INT n);
-/** Computes the inner/dot product \f$x^H x\f$. */
R Y(dot_double)(R *x, INT n);
/** Computes the weighted inner/dot product \f$x^H (w \odot x)\f$. */
R Y(dot_w_complex)(C *x, R *w, INT n);
@@ -1515,8 +1497,6 @@ void Y(cp_w_double)(R *x, R *w, R *y, INT n);
/* vector3.c */
/** Updates \f$x \leftarrow a x + y\f$. */
-void Y(upd_axpy_complex)(C *x, R a, C *y, INT n);
-/** Updates \f$x \leftarrow a x + y\f$. */
void Y(upd_axpy_double)(R *x, R a, R *y, INT n);
/** Updates \f$x \leftarrow x + a y\f$. */
void Y(upd_xpay_complex)(C *x, R a, C *y, INT n);
@@ -1534,15 +1514,6 @@ void Y(upd_xpawy_double)(R *x, R a, R *w, R *y, INT n);
void Y(upd_axpwy_complex)(C *x, R a, R *w, C *y, INT n);
/** Updates \f$x \leftarrow a x + w\odot y\f$. */
void Y(upd_axpwy_double)(R *x, R a, R *w, R *y, INT n);
-/** Swaps each half over N[d]/2. */
-void Y(fftshift_complex)(C *x, INT d, INT* N);
-void Y(fftshift_complex_int)(C *x, int d, int* N);
-
-/* print.c */
-/** Print real vector to standard output. */
-void Y(vpr_double)(R *x, const INT n, const char *text);
-/** Print complex vector to standard output. */
-void Y(vpr_complex)(C *x, const INT n, const char *text);
/* voronoi.c */
void Y(voronoi_weights_1d)(R *w, R *x, const INT M);
@@ -1562,9 +1533,6 @@ R Y(modified_sobolev)(const R mu, const INT kk);
/** Comput damping factor for modified multiquadric kernel. */
R Y(modified_multiquadric)(const R mu, const R c, const INT kk);
-/* thread.c */
-INT Y(get_num_threads)(void);
-
/* always check */
#define CK(ex) \
(void)((ex) || (Y(assertion_failed)(#ex, __LINE__, __FILE__), 0))
diff --git a/include/nfft3.h b/include/nfft3.h
index 02a4a5a..b4cba63 100644
--- a/include/nfft3.h
+++ b/include/nfft3.h
@@ -48,12 +48,14 @@ extern "C"
# define NFFT_EXTERN extern
#endif
-typedef ptrdiff_t _INT;
+/* Integral type large enough to contain a stride (what ``int'' should have been
+ * in the first place) */
+typedef ptrdiff_t NFFT_INT;
/* Members inherited by all plans. */
#define MACRO_MV_PLAN(RC) \
- _INT N_total; /**< Total number of Fourier coefficients. */\
- _INT M_total; /**< Total number of samples. */\
+ NFFT_INT N_total; /**< Total number of Fourier coefficients. */\
+ NFFT_INT M_total; /**< Total number of samples. */\
RC *f_hat; /**< Fourier coefficients. */\
RC *f; /**< Samples. */\
void (*mv_trafo)(void*); /**< Transform. */\
@@ -105,24 +107,25 @@ typedef struct \
MACRO_MV_PLAN(R) \
} X(mv_plan_double); \
\
+/** data structure for an NFFT (nonequispaced fast Fourier transform) plan with R precision */ \
typedef struct\
{\
MACRO_MV_PLAN(C)\
\
- _INT d; /**< Dimension (rank). */\
- _INT *N; /**< Multi-bandwidth. */\
+ NFFT_INT d; /**< Dimension (rank). */\
+ NFFT_INT *N; /**< Multi-bandwidth. */\
R *sigma; /**< Oversampling factor. */\
- _INT *n; /**< Length of FFTW transforms. This is equal to sigma*N. The default
+ NFFT_INT *n; /**< Length of FFTW transforms. This is equal to sigma*N. The default
is to use a power of two that satifies \f$2\le\sigma<4\f$. */\
- _INT n_total; /**< Combined total length of FFTW transforms. */\
- _INT m; /**< Cut-off parameter for window function. Default values for the
+ NFFT_INT n_total; /**< Combined total length of FFTW transforms. */\
+ NFFT_INT m; /**< Cut-off parameter for window function. Default values for the
different window functions are
- 6 (KAISER_BESSEL),
- 9 (SINC_POWER),
- 11 (B_SPLINE),
- 12 (GAUSSIAN) */\
R *b; /**< Shape parameter for window function */\
- _INT K; /**< Number of equispaced samples of window function. Used for flag
+ NFFT_INT K; /**< Number of equispaced samples of window function. Used for flag
PRE_LIN_PSI. */\
\
unsigned flags; /**< Flags for precomputation, (de)allocation, and FFTW
@@ -133,7 +136,7 @@ typedef struct\
unsigned fftw_flags; /**< Flags for the FFTW, default is
FFTW_ESTIMATE | FFTW_DESTROY_INPUT */\
\
- R *x; /**< Nodes in time/spatial domain, size is \f$dM\f$ doubles */\
+ R *x; /**< Nodes in time/spatial domain, size is \f$dM\f$ R ## s */\
\
R MEASURE_TIME_t[3]; /**< Measured time for each step if MEASURE_TIME is
set */\
@@ -146,8 +149,8 @@ typedef struct\
is \f$N_0+\hdots+N_{d-1}\f$ doubles*/\
R *psi; /**< Precomputed data for the sparse matrix \f$B\f$, size depends
on precomputation scheme */\
- _INT *psi_index_g; /**< Indices in source/target vector for \ref PRE_FULL_PSI */\
- _INT *psi_index_f; /**< Indices in source/target vector for \ref PRE_FULL_PSI */\
+ NFFT_INT *psi_index_g; /**< Indices in source/target vector for \ref PRE_FULL_PSI */\
+ NFFT_INT *psi_index_f; /**< Indices in source/target vector for \ref PRE_FULL_PSI */\
\
C *g; /**< Oversampled vector of samples, size is \ref n_total double complex */\
C *g_hat; /**< Zero-padded vector of Fourier coefficients, size is \ref n_total fftw_complex */\
@@ -156,7 +159,7 @@ typedef struct\
\
R *spline_coeffs; /**< Input for de Boor algorithm if B_SPLINE or SINC_POWER is defined */\
\
- _INT *index_x; /**< Index array for nodes x used when flag \ref NFFT_SORT_NODES is set. */\
+ NFFT_INT *index_x; /**< Index array for nodes x used when flag \ref NFFT_SORT_NODES is set. */\
} X(plan); \
\
NFFT_EXTERN void X(trafo_direct)(const X(plan) *ths);\
@@ -221,20 +224,21 @@ NFFT_DEFINE_API(NFFT_MANGLE_LONG_DOUBLE,FFTW_MANGLE_LONG_DOUBLE,long double,fftw
* C: complex data type
*/
#define NFCT_DEFINE_API(X,Y,R,C) \
+/** data structure for an NFCT (nonequispaced fast cosine transform) plan with R precision */ \
typedef struct\
{\
/* api */\
MACRO_MV_PLAN(R)\
\
- _INT d; /**< dimension, rank */\
- _INT *N; /**< cut-off-frequencies (kernel) */\
- _INT *n; /**< length of DCT-I */\
- _INT n_total; /**< Combined total length of FFTW transforms. */\
+ NFFT_INT d; /**< dimension, rank */\
+ NFFT_INT *N; /**< cut-off-frequencies (kernel) */\
+ NFFT_INT *n; /**< length of DCT-I */\
+ NFFT_INT n_total; /**< Combined total length of FFTW transforms. */\
R *sigma; /**< oversampling-factor */\
- _INT m; /**< cut-off parameter in time-domain */\
+ NFFT_INT m; /**< cut-off parameter in time-domain */\
\
R *b; /**< shape parameters */\
- _INT K; /**< Number of equispaced samples of window function. Used for flag
+ NFFT_INT K; /**< Number of equispaced samples of window function. Used for flag
PRE_LIN_PSI. */\
\
unsigned flags; /**< flags for precomputation, malloc */\
@@ -250,9 +254,9 @@ typedef struct\
\
R **c_phi_inv; /**< precomputed data, matrix D */\
R *psi; /**< precomputed data, matrix B */\
- _INT size_psi; /**< only for thin B */\
- _INT *psi_index_g; /**< only for thin B */\
- _INT *psi_index_f; /**< only for thin B */\
+ NFFT_INT size_psi; /**< only for thin B */\
+ NFFT_INT *psi_index_g; /**< only for thin B */\
+ NFFT_INT *psi_index_f; /**< only for thin B */\
\
R *g;\
R *g_hat;\
@@ -300,20 +304,21 @@ NFCT_DEFINE_API(NFCT_MANGLE_LONG_DOUBLE,FFTW_MANGLE_LONG_DOUBLE,long double,fftw
* C: complex data type
*/
#define NFST_DEFINE_API(X,Y,R,C) \
+/** data structure for an NFST (nonequispaced fast sine transform) plan with R precision */ \
typedef struct\
{\
/* api */\
MACRO_MV_PLAN(R)\
\
- _INT d; /**< dimension, rank */\
- _INT *N; /**< bandwidth */\
- _INT *n; /**< length of DST-I */\
- _INT n_total; /**< Combined total length of FFTW transforms. */\
+ NFFT_INT d; /**< dimension, rank */\
+ NFFT_INT *N; /**< bandwidth */\
+ NFFT_INT *n; /**< length of DST-I */\
+ NFFT_INT n_total; /**< Combined total length of FFTW transforms. */\
R *sigma; /**< oversampling-factor */\
- _INT m; /**< cut-off parameter in time-domain */\
+ NFFT_INT m; /**< cut-off parameter in time-domain */\
\
R *b; /**< shape parameters */\
- _INT K; /**< Number of equispaced samples of window function. Used for flag
+ NFFT_INT K; /**< Number of equispaced samples of window function. Used for flag
PRE_LIN_PSI. */\
\
unsigned flags; /**< flags for precomputation, malloc */\
@@ -329,9 +334,9 @@ typedef struct\
\
R **c_phi_inv; /**< precomputed data, matrix D */\
R *psi; /**< precomputed data, matrix B */\
- _INT size_psi; /**< only for thin B */\
- _INT *psi_index_g; /**< only for thin B */\
- _INT *psi_index_f; /**< only for thin B */\
+ NFFT_INT size_psi; /**< only for thin B */\
+ NFFT_INT *psi_index_g; /**< only for thin B */\
+ NFFT_INT *psi_index_f; /**< only for thin B */\
\
R *g;\
R *g_hat;\
@@ -382,6 +387,7 @@ NFST_DEFINE_API(NFST_MANGLE_LONG_DOUBLE,FFTW_MANGLE_LONG_DOUBLE,long double,fftw
* C: complex data type
*/
#define NNFFT_DEFINE_API(X,Y,Z,R,C) \
+/** data structure for an NNFFT (nonequispaced in time and frequency fast Fourier transform) plan with R precision */ \
typedef struct\
{\
/* api */\
@@ -449,6 +455,7 @@ NNFFT_DEFINE_API(NNFFT_MANGLE_LONG_DOUBLE,FFTW_MANGLE_LONG_DOUBLE,NFFT_MANGLE_LO
* C: complex data type
*/
#define NSFFT_DEFINE_API(X,Y,Z,R,C) \
+/** data structure for an NSFFT (nonequispaced sparse fast Fourier transform) plan with R precision */ \
typedef struct\
{\
MACRO_MV_PLAN(C)\
@@ -554,6 +561,7 @@ MRI_DEFINE_API(MRI_MANGLE_LONG_DOUBLE,NFFT_MANGLE_LONG_DOUBLE,long double,fftwl_
* C: complex data type
*/
#define NFSFT_DEFINE_API(X,Z,R,C) \
+/** data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with R precision */ \
typedef struct\
{\
MACRO_MV_PLAN(C)\
@@ -735,7 +743,7 @@ NFSOFT_DEFINE_API(NFSOFT_MANGLE_LONG_DOUBLE,NFFT_MANGLE_LONG_DOUBLE,FPT_MANGLE_L
#define NFSOFT_INDEX_TWO(m,n,l,B) ((B+1)*(B+1)+(B+1)*(B+1)*(m+B)-((m-1)*m*(2*m-1)+(B+1)*(B+2)*(2*B+3))/6)+(posN(n,m,B))+(l-MAX(ABS(m),ABS(n)))
#define NFSOFT_F_HAT_SIZE(B) (((B)+1)*(4*((B)+1)*((B)+1)-1)/3)
-/*solver */
+/* solver */
/* name mangling macros */
#define SOLVER_MANGLE_DOUBLE(name) NFFT_CONCAT(solver_, name)
@@ -750,6 +758,7 @@ NFSOFT_DEFINE_API(NFSOFT_MANGLE_LONG_DOUBLE,NFFT_MANGLE_LONG_DOUBLE,FPT_MANGLE_L
* C: complex data type
*/
#define SOLVER_DEFINE_API(X,Y,R,C)\
+/** data structure for an inverse NFFT plan with R precision */ \
typedef struct\
{\
Y(mv_plan_complex) *mv; /**< matrix vector multiplication */\
@@ -778,6 +787,7 @@ NFFT_EXTERN void X(before_loop_complex)(X(plan_complex)* ths);\
NFFT_EXTERN void X(loop_one_step_complex)(X(plan_complex) *ths);\
NFFT_EXTERN void X(finalize_complex)(X(plan_complex) *ths);\
\
+/** data structure for an inverse NFFT plan with R precision */ \
typedef struct\
{\
Y(mv_plan_double) *mv; /**< matrix vector multiplication */\
@@ -820,6 +830,61 @@ SOLVER_DEFINE_API(SOLVER_MANGLE_LONG_DOUBLE,NFFT_MANGLE_LONG_DOUBLE,long double,
#define PRECOMPUTE_WEIGHT (1U<< 5)
#define PRECOMPUTE_DAMP (1U<< 6)
+/* util */
+
+/* huge second-order macro that defines prototypes for all utility API functions.
+ * We expand this macro for each supported precision.
+ * Y: nfft name-mangling macro
+ * R: real data type
+ * C: complex data type
+ */
+#define NFFT_DEFINE_UTIL_API(Y,R,C) \
+/* rand.c */ \
+R Y(drand48)(void); \
+void Y(srand48)(long int seed); \
+\
+/** Inits a vector of random complex numbers in \f$[0,1]\times[0,1]{\rm i}\f$. \
+ */ \
+void Y(vrand_unit_complex)(C *x, const NFFT_INT n); \
+\
+/** Inits a vector of random double numbers in \f$[-1/2,1/2]\f$. \
+ */ \
+void Y(vrand_shifted_unit_double)(R *x, const NFFT_INT n); \
+\
+void Y(vrand_real)(R *x, const NFFT_INT n, const R a, const R b); \
+\
+/* print.c */ \
+/** Print real vector to standard output. */ \
+void Y(vpr_double)(R *x, const NFFT_INT n, const char *text); \
+\
+/** Print complex vector to standard output. */ \
+void Y(vpr_complex)(C *x, const NFFT_INT n, const char *text); \
+/* thread.c */ \
+NFFT_INT Y(get_num_threads)(void); \
+/* time.c */ \
+R Y(clock_gettime_seconds)(void); \
+/* error.c: */ \
+R Y(error_l_infty_complex)(const C *x, const C *y, const NFFT_INT n); \
+R Y(error_l_infty_1_complex)(const C *x, const C *y, const NFFT_INT n, \
+ const C *z, const NFFT_INT m); \
+/* int.c: */ \
+NFFT_INT Y(exp2i)(const NFFT_INT a); \
+NFFT_INT Y(next_power_of_2)(const NFFT_INT N); \
+/* vector1.c */ \
+/** Computes the inner/dot product \f$x^H x\f$. */ \
+R Y(dot_complex)(C *x, NFFT_INT n); \
+/* vector3.c */ \
+/** Updates \f$x \leftarrow a x + y\f$. */ \
+void Y(upd_axpy_complex)(C *x, R a, C *y, NFFT_INT n); \
+/** Swaps each half over N[d]/2. */ \
+void Y(fftshift_complex)(C *x, NFFT_INT d, NFFT_INT* N); \
+void Y(fftshift_complex_int)(C *x, int d, int* N);
+
+
+NFFT_DEFINE_UTIL_API(NFFT_MANGLE_FLOAT,float,fftwf_complex)
+NFFT_DEFINE_UTIL_API(NFFT_MANGLE_DOUBLE,double,fftw_complex)
+NFFT_DEFINE_UTIL_API(NFFT_MANGLE_LONG_DOUBLE,long double,fftwl_complex)
+
#ifdef __cplusplus
} /* extern "C" */
#endif /* __cplusplus */
diff --git a/include/nfft3mp.h b/include/nfft3mp.h
new file mode 100644
index 0000000..7da7754
--- /dev/null
+++ b/include/nfft3mp.h
@@ -0,0 +1,104 @@
+/*
+ * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
+ *
+ * This program is free software; you can redistribute it and/or modify it under
+ * the terms of the GNU General Public License as published by the Free Software
+ * Foundation; either version 2 of the License, or (at your option) any later
+ * version.
+ *
+ * This program 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, write to the Free Software Foundation, Inc., 51
+ * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* $Id$ */
+
+#ifndef __NFFT3MP_H__
+#define __NFFT3MP_H__
+
+#include "nfft3.h"
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+#if defined(NFFT_PRECISION_SINGLE)
+typedef float NFFT_R;
+typedef float _Complex NFFT_C;
+#define NFFT_K(x) ((NFFT_R) x)
+#define NFFT_M(name) NFFT_CONCAT(name,f)
+#define FFTW(name) NFFT_CONCAT(fftwf_,name)
+#define NFFT(name) NFFT_CONCAT(nfftf_,name)
+#define NFCT(name) NFFT_CONCAT(nfctf_,name)
+#define NFST(name) NFFT_CONCAT(nfstf_,name)
+#define NFSFT(name) NFFT_CONCAT(nfsftf_,name)
+#define SOLVER(name) NFFT_CONCAT(solverf_,name)
+#elif defined(NFFT_PRECISION_LONG_DOUBLE)
+typedef long double NFFT_R;
+typedef long double _Complex NFFT_C;
+#define NFFT_K(x) ((NFFT_R) x##L)
+#define NFFT_M(name) NFFT_CONCAT(name,l)
+#define FFTW(name) NFFT_CONCAT(fftwl_,name)
+#define NFFT(name) NFFT_CONCAT(nfftl_,name)
+#define NFCT(name) NFFT_CONCAT(nfctl_,name)
+#define NFST(name) NFFT_CONCAT(nfstl_,name)
+#define NFSFT(name) NFFT_CONCAT(nfsftl_,name)
+#define SOLVER(name) NFFT_CONCAT(solverl_,name)
+#elif defined(NFFT_PRECISION_DOUBLE)
+typedef double NFFT_R;
+typedef double _Complex NFFT_C;
+#define NFFT_K(x) ((NFFT_R) x)
+#define NFFT_M(name) name
+#define FFTW(name) NFFT_CONCAT(fftw_,name)
+#define NFFT(name) NFFT_CONCAT(nfft_,name)
+#define NFCT(name) NFFT_CONCAT(nfct_,name)
+#define NFST(name) NFFT_CONCAT(nfst_,name)
+#define NFSFT(name) NFFT_CONCAT(nfsft_,name)
+#define SOLVER(name) NFFT_CONCAT(solver_,name)
+#else
+#error Either define macro NFFT_PRECISION_SINGLE, NFFT_PRECISION_DOUBLE or NFFT_PRECISION_LONG_DOUBLE for single, double or long double precision
+#endif
+
+/* format strings */
+#if defined(NFFT_PRECISION_LONG_DOUBLE)
+# define NFFT__FGS__ "Lg"
+# define NFFT__FES__ "LE"
+# define NFFT__FE__ "% 36.32LE"
+# define NFFT__FI__ "%Lf"
+# define NFFT__FIS__ "Lf"
+# define NFFT__FR__ "%La"
+#elif defined(NFFT_PRECISION_SINGLE)
+# define NFFT__FGS__ "g"
+# define NFFT__FES__ "E"
+# define NFFT__FE__ "% 12.8E"
+# define NFFT__FI__ "%f"
+# define NFFT__FIS__ "f"
+# define NFFT__FR__ "%a"
+#elif defined(NFFT_PRECISION_DOUBLE)
+# define NFFT__FGS__ "lg"
+# define NFFT__FES__ "lE"
+# define NFFT__FE__ "% 20.16lE"
+# define NFFT__FI__ "%lf"
+# define NFFT__FIS__ "lf"
+# define NFFT__FR__ "%la"
+#else
+#error Either define macro NFFT_PRECISION_SINGLE, NFFT_PRECISION_DOUBLE or NFFT_PRECISION_LONG_DOUBLE for single, double or long double precision
+#endif
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif /* __cplusplus */
+
+/** Swap two vectors. */
+#define NFFT_CSWAP(x,y) {NFFT_C* NFFT_SWAP_temp__; \
+ NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
+
+#define NFFT_KPI NFFT_K(3.1415926535897932384626433832795028841971693993751)
+
+#endif /* defined(__NFFT3MP_H__) */
diff --git a/kernel/nfct/nfct.c b/kernel/nfct/nfct.c
index ab9a3f6..12a5726 100644
--- a/kernel/nfct/nfct.c
+++ b/kernel/nfct/nfct.c
@@ -1132,7 +1132,7 @@ const char* X(check)(X(plan) *ths)
for (j = 0; j < ths->d; j++)
{
if (ths->sigma[j] <= 1)
- return "nfft_check: oversampling factor too small";
+ return "Oversampling factor too small";
if(ths->N[j] - 1 <= ths->m)
return "Polynomial degree N is smaller than cut-off m";
diff --git a/kernel/nfft/nfft.c b/kernel/nfft/nfft.c
index 2b97fc2..ed54429 100644
--- a/kernel/nfft/nfft.c
+++ b/kernel/nfft/nfft.c
@@ -5893,7 +5893,7 @@ const char* X(check)(X(plan) *ths)
for (j = 0; j < ths->d; j++)
{
if (ths->sigma[j] <= 1)
- return "nfft_check: oversampling factor too small";
+ return "Oversampling factor too small";
if(ths->N[j] <= ths->m)
return "Polynomial degree N is smaller than cut-off m";
diff --git a/kernel/nfsft/nfsft.c b/kernel/nfsft/nfsft.c
index c6595a0..f02b633 100644
--- a/kernel/nfsft/nfsft.c
+++ b/kernel/nfsft/nfsft.c
@@ -263,7 +263,7 @@ void nfsft_init_advanced(nfsft_plan* plan, int N, int M,
unsigned int flags)
{
/* Call nfsft_init_guru with the flags and default NFFT cut-off. */
- nfsft_init_guru(plan, N, M, flags, PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
+ nfsft_init_guru(plan, N, M, flags, PRE_PHI_HUT | PRE_PSI | FFTW_INIT | NFFT_OMP_BLOCKWISE_ADJOINT |
FFT_OUT_OF_PLACE, NFSFT_DEFAULT_NFFT_CUTOFF);
}
diff --git a/kernel/nfsoft/nfsoft.c b/kernel/nfsoft/nfsoft.c
index 022ca88..2644f14 100755
--- a/kernel/nfsoft/nfsoft.c
+++ b/kernel/nfsoft/nfsoft.c
@@ -45,7 +45,7 @@ void nfsoft_init(nfsoft_plan *plan, int N, int M)
void nfsoft_init_advanced(nfsoft_plan *plan, int N, int M,
unsigned int nfsoft_flags)
{
- nfsoft_init_guru(plan, N, M, nfsoft_flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X
+ nfsoft_init_guru(plan, N, M, nfsoft_flags, PRE_PHI_HUT | PRE_PSI | MALLOC_X | NFFT_OMP_BLOCKWISE_ADJOINT
| MALLOC_F_HAT | MALLOC_F | FFTW_INIT | FFT_OUT_OF_PLACE,
DEFAULT_NFFT_CUTOFF, FPT_THRESHOLD);
}
diff --git a/kernel/nfst/nfst.c b/kernel/nfst/nfst.c
index 8d682bd..6551957 100644
--- a/kernel/nfst/nfst.c
+++ b/kernel/nfst/nfst.c
@@ -1131,7 +1131,7 @@ const char* X(check)(X(plan) *ths)
for (j = 0; j < ths->d; j++)
{
if (ths->sigma[j] <= 1)
- return "nfft_check: oversampling factor too small";
+ return "Oversampling factor too small";
if(ths->N[j] - 1 <= ths->m)
return "Polynomial degree N is smaller than cut-off m";
diff --git a/kernel/nnfft/nnfft.c b/kernel/nnfft/nnfft.c
index de69e9d..6269bbf 100644
--- a/kernel/nnfft/nnfft.c
+++ b/kernel/nnfft/nnfft.c
@@ -590,7 +590,7 @@ void nnfft_init_guru(nnfft_plan *ths, int d, int N_total, int M_total, int *N, i
ths->m= m;
ths->nnfft_flags= nnfft_flags;
fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
- nfft_flags= PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE;
+ nfft_flags= PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE| NFFT_OMP_BLOCKWISE_ADJOINT;
if(ths->nnfft_flags & PRE_PSI)
nfft_flags = nfft_flags | PRE_PSI;
@@ -643,7 +643,7 @@ ths->m=WINDOW_HELP_ESTIMATE_m;
}
ths->nnfft_flags=PRE_PSI| PRE_PHI_HUT| MALLOC_X| MALLOC_V| MALLOC_F_HAT| MALLOC_F;
- nfft_flags= PRE_PSI| PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE;
+ nfft_flags= PRE_PSI| PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE| NFFT_OMP_BLOCKWISE_ADJOINT;
fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
nnfft_init_help(ths,ths->m,nfft_flags,fftw_flags);
diff --git a/kernel/util/time.c b/kernel/util/time.c
index 72e563c..1876180 100644
--- a/kernel/util/time.c
+++ b/kernel/util/time.c
@@ -20,9 +20,30 @@
#include "infft.h"
+#ifdef HAVE_TIME_H
+#include <time.h>
+#endif
+
R Y(elapsed_seconds)(ticks t1, ticks t0)
{
UNUSED(t1);
UNUSED(t0);
return (R)(elapsed(t1,t0)) / (R)(TICKS_PER_SECOND);
}
+
+R Y(clock_gettime_seconds)(void)
+{
+#if defined(HAVE_CLOCK_GETTIME)
+ struct timespec tp;
+ if (clock_gettime(CLOCK_REALTIME, &tp) != 0)
+ {
+ tp.tv_sec = 0;
+ tp.tv_nsec = 0;
+ }
+
+ return tp.tv_sec+(R)tp.tv_nsec/K(1e9);
+#else
+ return K(0.0);
+#endif
+}
+
diff --git a/matlab/nfft/Makefile.am b/matlab/nfft/Makefile.am
index 26162a3..98e9952 100644
--- a/matlab/nfft/Makefile.am
+++ b/matlab/nfft/Makefile.am
@@ -23,7 +23,8 @@ dist_nfftmatlab_DATA = FFT_OUT_OF_PLACE.m FFTW_ESTIMATE.m FFTW_MEASURE.m FG_PSI.
nfft_adjoint.m nfft_finalize.m nfft_get_f.m nfft_get_f_hat.m nfft_get_x.m nfft_init_1d.m nfft_init_2d.m \
nfft_init_3d.m nfft_init_guru.m nfft_precompute_psi.m nfft_set_f.m nfft_set_f_hat.m nfft_set_x.m nfft_trafo.m \
PRE_FG_PSI.m PRE_FULL_PSI.m PRE_LIN_PSI.m PRE_PHI_HUT.m PRE_PSI.m simple_test.m \
- nfft_get_num_threads.m nfft.m test_nfft1d.m test_nfft2d.m test_nfft3d.m
+ nfft_get_num_threads.m nfft.m test_nfft1d.m test_nfft2d.m test_nfft3d.m \
+ NFFT_OMP_BLOCKWISE_ADJOINT.m
# target all-am builds .libs/libnfft at matlab_mexext@
nfftmex at matlab_mexext@: all-am
diff --git a/matlab/nnfft/FG_PSI.m b/matlab/nfft/NFFT_OMP_BLOCKWISE_ADJOINT.m
similarity index 68%
rename from matlab/nnfft/FG_PSI.m
rename to matlab/nfft/NFFT_OMP_BLOCKWISE_ADJOINT.m
index 4104144..b95d6e5 100644
--- a/matlab/nnfft/FG_PSI.m
+++ b/matlab/nfft/NFFT_OMP_BLOCKWISE_ADJOINT.m
@@ -1,8 +1,6 @@
-%FG_PSI Precompuation flag
-% If this flag is set, the convolution step (the multiplication with the
-% sparse matrix B) uses particular properties of the Gaussian window function
-% to trade multiplications for direct calls to exponential function.
-%
+%NFFT_OMP_BLOCKWISE_ADJOINT Flag which results in usage of (probably) faster algorithm in adjoint NFFT with OpenMP support.
+% If this flag is set, additionally the nonequispaced nodes are sorted
+% which may also accelerate the standard NFFT trafo due to cache effects.
% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
@@ -21,7 +19,7 @@
% this program; if not, write to the Free Software Foundation, Inc., 51
% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
%
-% $Id: FG_PSI.m 3776 2012-06-03 13:29:25Z keiner $
-function f = FG_PSI()
+% $Id: PRE_PSI.m 4272 2015-01-10 19:39:35Z keiner $
+function f = NFFT_OMP_BLOCKWISE_ADJOINT()
-f = bitshift(1, 1);
+f = bitshift(1, 12);
diff --git a/matlab/nfft/nfft.m b/matlab/nfft/nfft.m
index 484573a..b1870df 100644
--- a/matlab/nfft/nfft.m
+++ b/matlab/nfft/nfft.m
@@ -62,9 +62,9 @@ function h=nfft(d,N,M,varargin)
%
% h=nfft(d,N,M,varargin) for use of nfft_init_guru
% For example
-% h=nfft(1,N,M,n,7,'PRE_PHI_HUT','FFTW_MEASURE') for d=1
-% h=nfft(2,N,M,n,n,7,'PRE_PHI_HUT','FFTW_MEASURE') for d=2
-% h=nfft(3,N,M,n,n,n,7,'PRE_PHI_HUT','FFTW_MEASURE') for d=3
+% h=nfft(1,N,M,n,7,bitor(PRE_PHI_HUT,PRE_PSI),FFTW_MEASURE) for d=1, m=7
+% h=nfft(2,N,M,n,n,7,bitor(PRE_PHI_HUT,bitor(PRE_PSI,NFFT_OMP_BLOCKWISE_ADJOINT)),FFTW_MEASURE) for d=2, m=7
+% h=nfft(3,N,M,n,n,n,7,bitor(PRE_PHI_HUT,bitor(PRE_PSI,NFFT_OMP_BLOCKWISE_ADJOINT)),FFTW_MEASURE) for d=3, m=7
% with n=2^(ceil(log(max(N))/log(2))+1)
% Be careful: There is no error handling with using nfft_init_guru.
% Incorrect inputs can cause a Matlab crash!
diff --git a/matlab/nfft/simple_test.m b/matlab/nfft/simple_test.m
index bc17382..b8edb0f 100644
--- a/matlab/nfft/simple_test.m
+++ b/matlab/nfft/simple_test.m
@@ -70,7 +70,7 @@ for logN=3:10
N=2^logN;
M=N^2;
x=rand(2,M)-0.5;
- plan = nfft_init_guru(d,N,N,M,2*N,2*N,2,bitor(PRE_PHI_HUT,PRE_PSI),FFTW_MEASURE);
+ plan = nfft_init_guru(d,N,N,M,2*N,2*N,2,bitor(PRE_PHI_HUT,bitor(PRE_PSI,NFFT_OMP_BLOCKWISE_ADJOINT)),FFTW_MEASURE);
nfft_set_x(plan,x);
nfft_precompute_psi(plan);
f_hat = rand(N,N)+i*rand(N,N);
@@ -103,7 +103,7 @@ for logN=3:10
N3=18; n3=32;
M=N1*N2*N3;
x=rand(3,M)-0.5;
- plan = nfft_init_guru(d,N1,N2,N3,M,n1,n2,n3,2,bitor(PRE_PHI_HUT,PRE_PSI),FFTW_MEASURE);
+ plan = nfft_init_guru(d,N1,N2,N3,M,n1,n2,n3,2,bitor(PRE_PHI_HUT,bitor(PRE_PSI,NFFT_OMP_BLOCKWISE_ADJOINT)),FFTW_MEASURE);
nfft_set_x(plan,x);
nfft_precompute_psi(plan);
f_hat = rand(N1,N2,N3)+i*rand(N1,N2,N3);
diff --git a/matlab/nfft/test_nfft1d.m b/matlab/nfft/test_nfft1d.m
index 22093cd..d93f371 100644
--- a/matlab/nfft/test_nfft1d.m
+++ b/matlab/nfft/test_nfft1d.m
@@ -26,7 +26,7 @@ x=rand(M,1)-0.5; %nodes
% Initialisation
plan=nfft(1,N,M); % create plan of class type nfft
%n=2^(ceil(log(N)/log(2))+1);
-%plan=nfft(1,N,M,n,7,'PRE_PHI_HUT','FFTW_MEASURE'); % use of nfft_init_guru
+%plan=nfft(1,N,M,n,7,bitor(PRE_PHI_HUT,PRE_PSI),FFTW_MEASURE); % use of nfft_init_guru
plan.x=x; % set nodes in plan
nfft_precompute_psi(plan); % precomputations
diff --git a/matlab/nfft/test_nfft2d.m b/matlab/nfft/test_nfft2d.m
index 6285fb0..e8543c4 100644
--- a/matlab/nfft/test_nfft2d.m
+++ b/matlab/nfft/test_nfft2d.m
@@ -28,7 +28,7 @@ x=rand(M,2)-0.5; %nodes
% Initialisation
plan=nfft(2,N,M); % create plan of class type nfft
%n=2^(ceil(log(max(N))/log(2))+1);
-%plan=nfft(2,N,M,n,n,7,'PRE_PHI_HUT','FFTW_MEASURE'); % use of nfft_init_guru
+%plan=nfft(2,N,M,n,n,7,bitor(PRE_PHI_HUT,bitor(PRE_PSI,NFFT_OMP_BLOCKWISE_ADJOINT)),FFTW_MEASURE); % use of nfft_init_guru
plan.x=x; % set nodes in plan
nfft_precompute_psi(plan); % precomputations
diff --git a/matlab/nfft/test_nfft3d.m b/matlab/nfft/test_nfft3d.m
index 03bac61..f43a105 100644
--- a/matlab/nfft/test_nfft3d.m
+++ b/matlab/nfft/test_nfft3d.m
@@ -29,7 +29,7 @@ x=rand(M,3)-0.5; %nodes
% Initialisation
plan=nfft(3,N,M); % create plan of class type nfft
%n=2^(ceil(log(max(N))/log(2))+1);
-%plan=nfft(3,N,M,n,n,n,7,'PRE_PHI_HUT','FFTW_MEASURE'); % use of nfft_init_guru
+%plan=nfft(3,N,M,n,n,n,7,bitor(PRE_PHI_HUT,bitor(PRE_PSI,NFFT_OMP_BLOCKWISE_ADJOINT)),FFTW_MEASURE); % use of nfft_init_guru
plan.x=x; % set nodes in plan
nfft_precompute_psi(plan); % precomputations
diff --git a/matlab/nnfft/FFTW_ESTIMATE.m b/matlab/nnfft/FFTW_ESTIMATE.m
deleted file mode 100644
index 4892047..0000000
--- a/matlab/nnfft/FFTW_ESTIMATE.m
+++ /dev/null
@@ -1,25 +0,0 @@
-% FFTW_ESTIMATE FFT flag
-% Valid for FFTW3
-%
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program 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, write to the Free Software Foundation, Inc., 51
-% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-%
-% $Id: FFTW_ESTIMATE.m 3776 2012-06-03 13:29:25Z keiner $
-function f = FFTW_ESTIMATE()
-
-f = bitshift(1, 6);
diff --git a/matlab/nnfft/FFTW_MEASURE.m b/matlab/nnfft/FFTW_MEASURE.m
deleted file mode 100644
index 5553de5..0000000
--- a/matlab/nnfft/FFTW_MEASURE.m
+++ /dev/null
@@ -1,25 +0,0 @@
-%FFTW_MEASURE FFT flag
-% Valid for FFTW3
-%
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program 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, write to the Free Software Foundation, Inc., 51
-% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-%
-% $Id: FFTW_MEASURE.m 3776 2012-06-03 13:29:25Z keiner $
-function f = FFTW_MEASURE()
-
-f = 0;
diff --git a/matlab/nnfft/FFT_OUT_OF_PLACE.m b/matlab/nnfft/FFT_OUT_OF_PLACE.m
deleted file mode 100644
index dac4ab6..0000000
--- a/matlab/nnfft/FFT_OUT_OF_PLACE.m
+++ /dev/null
@@ -1,25 +0,0 @@
-%FFT_OUT_OF_PLACE FFT flag
-% If this flag is set, FFTW uses disjoint input/output vectors.
-%
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program 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, write to the Free Software Foundation, Inc., 51
-% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-%
-% $Id: FFT_OUT_OF_PLACE.m 3776 2012-06-03 13:29:25Z keiner $
-function f = FFT_OUT_OF_PLACE()
-
-f = bitshift(1, 9);
diff --git a/matlab/nnfft/Makefile.am b/matlab/nnfft/Makefile.am
index 898af3e..c662f07 100644
--- a/matlab/nnfft/Makefile.am
+++ b/matlab/nnfft/Makefile.am
@@ -4,7 +4,7 @@
AM_CPPFLAGS = -I$(top_srcdir)/include -I$(top_srcdir)/matlab $(matlab_CPPFLAGS)
# matlab wrapper directory
-nnfftmatlabdir = $(datadir)/nnfft/matlab/nnfft
+nnfftmatlabdir = $(datadir)/nfft/matlab/nnfft
# library
lib_LTLIBRARIES = libnnfft.la
@@ -19,11 +19,12 @@ libnnfft_la_CFLAGS = $(OPENMP_CFLAGS)
endif
-dist_nnfftmatlab_DATA = FFT_OUT_OF_PLACE.m FFTW_ESTIMATE.m FFTW_MEASURE.m FG_PSI.m Contents.m \
- nnfft_finalize.m nnfft_get_f.m nnfft_get_f_hat.m nnfft_get_x.m nnfft_init_1d.m \
- nnfft_init_guru.m nnfft_precompute_psi.m nnfft_set_f.m nnfft_set_f_hat.m nnfft_set_x.m nnfft_trafo.m \
- PRE_FG_PSI.m PRE_FULL_PSI.m PRE_LIN_PSI.m PRE_PHI_HUT.m PRE_PSI.m simple_test.m \
- nnfft_get_num_threads.m nnfft.m
+dist_nnfftmatlab_DATA = Contents.m nnfft.m nnfft_display.m nnfft_finalize.m \
+ nnfft_get_f.m nnfft_get_f_hat.m nnfft_get_x.m nnfft_init.m nnfft_init_1d.m \
+ nnfft_init_2d.m nnfft_init_3d.m nnfft_init_guru.m nnfft_precompute_psi.m \
+ nnfft_set_f.m nnfft_set_f_hat.m nnfft_set_v.m nnfft_set_x.m nnfft_trafo.m \
+ nnfft_trafo_direct.m PRE_FULL_PSI.m PRE_LIN_PSI.m PRE_PHI_HUT.m PRE_PSI.m \
+ simple_test.m test_nnfft1d.m test_nnfft2d.m nnfft_get_num_threads.m
# target all-am builds .libs/libnnfft at matlab_mexext@
nnfftmex at matlab_mexext@: all-am
diff --git a/matlab/nnfft/PRE_FG_PSI.m b/matlab/nnfft/PRE_FG_PSI.m
deleted file mode 100644
index 430dd60..0000000
--- a/matlab/nnfft/PRE_FG_PSI.m
+++ /dev/null
@@ -1,28 +0,0 @@
-%PRE_FG_PSI Precomputation flag
-% If this flag is set, the convolution step (the multiplication with the sparse
-% matrix B) uses particular properties of the Gaussian window function to trade
-% multiplications for direct calls to exponential function (the remaining 2dM
-% direct calls are precomputed).
-%
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program 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, write to the Free Software Foundation, Inc., 51
-% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-%
-% $Id: PRE_FG_PSI.m 3776 2012-06-03 13:29:25Z keiner $
-function f = PRE_FG_PSI()
-
-f = bitshift(1, 3);
diff --git a/matlab/nnfft/nnfft.m b/matlab/nnfft/nnfft.m
index 770adaf..cb3486b 100644
--- a/matlab/nnfft/nnfft.m
+++ b/matlab/nnfft/nnfft.m
@@ -24,7 +24,7 @@ classdef nnfft < handle
properties(Dependent=true)
x; % nodes in time/spatial domain (real Mxd matrix)
v; % nodes in fourier domain (real N_totalxd matrix)
- fhat; % fourier coefficients (complex column vector of length N_1, N_1*N_2 or N_1*N_2*N_3, for d=2 columnwise linearisation of N_1xN_2 matrix and for d=3 columnwise linearisation of N_1xN_2xN_3 array)
+ fhat; % fourier coefficients (complex column vector of length N_total)
f; % samples (complex column vector of length M)
end %properties
@@ -69,9 +69,9 @@ function h=nnfft(d,N_total,M,N,varargin)
%
% h=nnfft(d,N_total,M,N,varargin) for use of nnfft_init_guru
% For example
-% h=nnfft(1,N_total,M,N,N1,7,'PRE_PHI_HUT') for d=1
-% h=nnfft(2,N_total,M,N,N,N1,N1,7,'PRE_PHI_HUT') for d=2
-% h=nnfft(3,N_total,M,N,N,N,N1,N1,N1,7,'PRE_PHI_HUT') for d=3
+% h=nnfft(1,N_total,M,N,N1,7,bitor(PRE_PHI_HUT,PRE_PSI)) for d=1, m=7
+% h=nnfft(2,N_total,M,N,N1_1,N1_2,7,bitor(PRE_PHI_HUT,PRE_PSI)) for d=2, m=7
+% h=nnfft(3,N_total,M,N,N1_1,N1_2,N1_3,7,bitor(PRE_PHI_HUT,PRE_PSI)) for d=3, m=7
% with N1=sigma*N ; n=N1
% NOT IMPLEMENTED: Be careful: There is no error handling with using nfft_init_guru.
% Incorrect inputs can cause a Matlab crash!
@@ -84,7 +84,7 @@ function h=nnfft(d,N_total,M,N,varargin)
% varargin parameters for use of nnfft_init_guru (see documentation of NFFT for more details)
%
% OUTPUT
-% h object of class type nfft
+% h object of class type nnfft
h.d=d;
h.N_total=N_total;
@@ -284,38 +284,24 @@ function set.v(h,v)
end %function
function set.fhat(h,fhat)
- switch h.d
- case 1
- n=h.N_1;
- case 2
- n=h.N_1*h.N_2;
- case 3
- n=h.N_1*h.N_2*h.N_3;
- otherwise
- error('Unknown error.');
- end % switch
+ %switch h.d
+ %case 1
+ %n=h.N_1;
+ %case 2
+ % n=h.N_1*h.N_2;
+ %case 3
+ % n=h.N_1*h.N_2*h.N_3;
+ %otherwise
+ % error('Unknown error.');
+ %end % switch
+
+ n=h.N_total;
if( isempty(fhat) || ~isnumeric(fhat))
error('The Fourier coefficients fhat have to be complex numbers.');
elseif( size(fhat,1)~=(n) || size(fhat,2)~=1 )
error('The Fourier coefficients fhat have to be a column vector of length %u.',n);
else
- switch h.d
- case 1
- % Do nothing.
- case 2
- % linearization in matlab with column (:) operator is columnwise, in NFFT it is rowwise
- fhat=reshape(fhat,h.N_1,h.N_2).';
- fhat=fhat(:);
- case 3
- % linearization in matlab with column (:) operator is columnwise, in NFFT it is rowwise
- fhat=reshape(fhat,h.N_1,h.N_2,h.N_3);
- fhat=permute(fhat,[3,2,1]);
- fhat=fhat(:);
- otherwise
- error('Unknown error.');
- end %switch
-
nnfftmex('set_f_hat',h.plan,fhat);
h.fhat_is_set=true;
end %if
@@ -346,21 +332,21 @@ function fhat=get.fhat(h)
if(h.fhat_is_set)
fhat=nnfftmex('get_f_hat',h.plan);
- switch h.d
- case 1
+ %switch h.d
+ %case 1
% Do nothing.
- case 2
+ %case 2
% linearization in matlab with column (:) operator is columnwise, in NFFT it is rowwise
- fhat=reshape(fhat,h.N_2,h.N_1).';
- fhat=fhat(:);
- case 3
+ %fhat=reshape(fhat,h.N_2,h.N_1).';
+ %fhat=fhat(:);
+ %case 3
% linearization in matlab with column (:) operator is columnwise, in NFFT it is rowwise
- fhat=reshape(fhat,h.N_3,h.N_2,h.N_1);
- fhat=permute(fhat,[3,2,1]);
- fhat=fhat(:);
- otherwise
- error('Unknown error.');
- end %switch
+ %fhat=reshape(fhat,h.N_3,h.N_2,h.N_1);
+ %fhat=permute(fhat,[3,2,1]);
+ %fhat=fhat(:);
+ %otherwise
+ % error('Unknown error.');
+ %end %switch
else
fhat=[];
end %if
diff --git a/matlab/nnfft/out.txt b/matlab/nnfft/out.txt
deleted file mode 100644
index 02479dc..0000000
--- a/matlab/nnfft/out.txt
+++ /dev/null
@@ -1,32 +0,0 @@
-insgesamt 136
-drwxr-xr-x 2 sukunis dip 4096 Jun 23 14:05 .
-drwxr-xr-x 7 sukunis dip 137 Jun 23 13:14 ..
--rw-r--r-- 1 sukunis dip 1467 Apr 30 08:59 Contents.m
--rw-r--r-- 1 sukunis dip 1055 Apr 30 08:59 FFT_OUT_OF_PLACE.m
--rw-r--r-- 1 sukunis dip 1001 Apr 30 08:59 FFTW_ESTIMATE.m
--rw-r--r-- 1 sukunis dip 984 Apr 30 08:59 FFTW_MEASURE.m
--rw-r--r-- 1 sukunis dip 1195 Apr 30 08:59 FG_PSI.m
--rw-r--r-- 1 sukunis dip 1592 Apr 30 08:59 Makefile.am
--rw-r--r-- 1 sukunis dip 982 Apr 30 08:59 nnfft_finalize.m
--rw-r--r-- 1 sukunis dip 1024 Apr 30 08:59 nnfft_get_f_hat.m
--rw-r--r-- 1 sukunis dip 995 Apr 30 08:59 nnfft_get_f.m
--rw-r--r-- 1 sukunis dip 1058 Apr 30 08:59 nnfft_get_num_threads.m
--rw-r--r-- 1 sukunis dip 985 Apr 30 08:59 nnfft_get_x.m
--rw-r--r-- 1 sukunis dip 994 Apr 30 08:59 nnfft_init_1d.m
--rw-r--r-- 1 sukunis dip 1496 Apr 30 08:59 nnfft_init_guru.m
--rw-r--r-- 1 sukunis dip 10229 Apr 30 08:59 nnfft.m
--rw-r--r-- 1 sukunis dip 13874 Apr 30 08:59 nnfftmex.c
--rw-r--r-- 1 sukunis dip 951 Apr 30 08:59 nnfftmex.m
--rw-r--r-- 1 sukunis dip 1029 Apr 30 08:59 nnfft_precompute_psi.m
--rw-r--r-- 1 sukunis dip 1017 Apr 30 08:59 nnfft_set_f_hat.m
--rw-r--r-- 1 sukunis dip 988 Apr 30 08:59 nnfft_set_f.m
--rw-r--r-- 1 sukunis dip 978 Apr 30 08:59 nnfft_set_x.m
--rw-r--r-- 1 sukunis dip 995 Apr 30 08:59 nnfft_trafo.m
--rw-r--r-- 1 sukunis dip 0 Jun 23 14:05 out.txt
--rw-r--r-- 1 sukunis dip 1261 Apr 30 08:59 PRE_FG_PSI.m
--rw-r--r-- 1 sukunis dip 1207 Apr 30 08:59 PRE_FULL_PSI.m
--rw-r--r-- 1 sukunis dip 1224 Apr 30 08:59 PRE_LIN_PSI.m
--rw-r--r-- 1 sukunis dip 1157 Apr 30 08:59 PRE_PHI_HUT.m
--rw-r--r-- 1 sukunis dip 1123 Apr 30 08:59 PRE_PSI.m
--rw-r--r-- 1 sukunis dip 3189 Apr 30 08:59 simple_test.m
--rw-r--r-- 1 sukunis dip 1963 Apr 30 08:59 test_nfft1d.m
diff --git a/matlab/nnfft/test_nfft1d.m b/matlab/nnfft/test_nfft1d.m
deleted file mode 100644
index 22093cd..0000000
--- a/matlab/nnfft/test_nfft1d.m
+++ /dev/null
@@ -1,70 +0,0 @@
-
-% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
-%
-% This program is free software; you can redistribute it and/or modify it under
-% the terms of the GNU General Public License as published by the Free Software
-% Foundation; either version 2 of the License, or (at your option) any later
-% version.
-%
-% This program 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, write to the Free Software Foundation, Inc., 51
-% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-
-% Test script of class nfft for spatial dimension d=1.
-clear all;
-
-M=16; % number of nodes
-N=24; % number of Fourier coefficients in first direction
-
-x=rand(M,1)-0.5; %nodes
-
-% Initialisation
-plan=nfft(1,N,M); % create plan of class type nfft
-%n=2^(ceil(log(N)/log(2))+1);
-%plan=nfft(1,N,M,n,7,'PRE_PHI_HUT','FFTW_MEASURE'); % use of nfft_init_guru
-
-plan.x=x; % set nodes in plan
-nfft_precompute_psi(plan); % precomputations
-
-% NFFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-fhat=rand(N,1); % Fourier coefficients
-fhatv=fhat(:);
-
-% Compute samples with NFFT
-plan.fhat=fhatv; % set Fourier coefficients
-nfft_trafo(plan); % compute nonequispaced Fourier transform
-f1=plan.f; % get samples
-
-% Compute samples direct
-k1=(-N/2:N/2-1).';
-f2=zeros(M,1);
-for j=1:M
- x1j=x(j,1);
- f2(j)=sum( fhatv.*exp(-2*pi*1i*k1*x1j) );
-end %for
-
-% Compare results
-max(abs(f1-f2))
-
-% Adjoint NFFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Computation with NFFT
-nfft_adjoint(plan);
-fhat1=plan.fhat;
-
-% Direct computation
-fhat2=zeros(N,1);
-for j=1:N
- k1j=k1(j);
- fhat2(j)=sum( plan.f.*exp(2*pi*1i*k1j*x(:,1)) );
-end %for
-
-% Compare results
-max(abs(fhat1-fhat2))
-
diff --git a/matlab/nnfft/test_nnfft1d.m b/matlab/nnfft/test_nnfft1d.m
index 0076a60..667a297 100644
--- a/matlab/nnfft/test_nnfft1d.m
+++ b/matlab/nnfft/test_nnfft1d.m
@@ -19,17 +19,20 @@
clear all;
M=16; % number of nodes
-N=8; % number of Fourier coefficients in first direction
-N1=2*N;
-N_total=8;
+N=24; % number of Fourier coefficients in first direction
+N_total=N; % total number of Fourier coefficients
x=rand(M,1)-0.5; %nodes
-v=rand(N_total,1)-0.5; %nodes
+v=rand(N,1)-0.5; %nodes
-% Initialisation
-%plan=nnfft(1,N_total,M,N); % create plan of class type nfft
+% Plan initialisation simple interface
+plan=nnfft(1,N_total,M,N); % create plan of class type nfft
-plan=nnfft(1,N_total,M,N,N1,6,'PRE_PHI_HUT'); % use of nfft_init_guru
+% Plan initialisation guru interface
+%sigma=2; % oversampling factor
+%N1=sigma*N; % FFTW length, must be even natural number!
+%m=6; % window cut-off parameter
+%plan=nnfft(1,N,M,N,N1,m,bitor(PRE_PHI_HUT,PRE_PSI)); % create plan of class type nfft
plan.x=x; % set nodes in plan
plan.v=v;
@@ -43,7 +46,7 @@ fhatv=fhat(:);
% Compute samples with NFFT
plan.fhat=fhatv; % set Fourier coefficients
-nnfft_trafo(plan); % compute nonequispaced Fourier transform
+nnfft_trafo(plan); % compute nonequispaced Fourier transform (in space and time)
f1=plan.f; % get samples
% Compute samples direct
diff --git a/matlab/nnfft/test_nnfft2d.m b/matlab/nnfft/test_nnfft2d.m
index abeaa2a..7eab4be 100644
--- a/matlab/nnfft/test_nnfft2d.m
+++ b/matlab/nnfft/test_nnfft2d.m
@@ -18,22 +18,32 @@
% Test script of class nnfft for spatial dimension d=2.
clear all;
-M=16; % number of nodes
-N_1=8; % number of Fourier coefficients in first direction
-N_2=8; % number of Fourier coefficients in second direction
+%M=16; % number of nodes
+%N_1=24; % number of Fourier coefficients in first direction
+%N_2=32; % number of Fourier coefficients in second direction
+%N=[N_1;N_2];
+%N_total=N_1*N_2; % total number of Fourier coefficients
+M=4;
+N_1=2;
+N_2=2;
N=[N_1;N_2];
-N1_1=2*N_1;
-N1_2=2*N_2;
-N1=[N1_1;N1_2];
-N_total=N_1*N_2;
+N_total=6;
+
+
x=rand(M,2)-0.5; %nodes
v=rand(N_total,2)-0.5; %nodes
-% Initialisation
-%plan=nnfft(2,N_total,M,N); % create plan of class type nnfft
+% Plan initialisation simple interface
+plan=nnfft(2,N_total,M,N); % create plan of class type nnfft
-plan=nnfft(2,N_total,M,N,N1,N1,6,'PRE_PHI_HUT'); % use of nnfft_init_guru
+
+% Plan initialisation guru interface
+%sigma=2; % oversampling factor
+%N1_1=sigma*N_1; % FFTW length, must be even natural number!
+%N1_2=sigma*N_2; % FFTW length, must be even natural number!
+%m=6; % window cut-off parameter
+%plan=nnfft(2,N_total,M,N,N1_1,N1_2,m,bitor(PRE_PHI_HUT,PRE_PSI)); % create plan of class type nnfft
plan.x=x; % set nodes in plan
plan.v=v; % set nodes in plan
@@ -41,22 +51,35 @@ nnfft_precompute_psi(plan); % precomputations
% NFFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-fhat=rand(N_1,N_2); % Fourier coefficients
-fhatv=fhat(:);
+%fhat=rand(N_1,N_2); % Fourier coefficients
+%fhat=rand(N_total,1);
+%fhat=ones(N_total,1);
+
+%Test mit zweitem Einheitsvektor
+fhat=zeros(N_total,1);
+fhat(2)=1;
+fhatv=fhat;%(:);
% Compute samples with NNFFT
-plan.fhat=fhatv; % set Fourier coefficients
-nnfft_trafo(plan); % compute nonequispaced Fourier transform
-f1=plan.f; % get samples
+ plan.fhat=fhatv; % set Fourier coefficients
+ nnfft_trafo(plan); % compute nonequispaced Fourier transform
+ f1=plan.f; % get samples
% Compute samples direct
nnfft_trafo_direct(plan);
f2=plan.f;
-% Compare results
+%% Compare results
disp('NNFFT vs NNDFT');
max(abs(f1-f2))
+A=exp(-2*pi*1i*x*diag(N)*v');
+f3=A*fhatv;
+%tmpv=(v.*repmat(N',size(v,1),1)).';
+%geflipter zweiter Einheitsvektor
+%f3=exp(-2*pi*1i*(x*tmpv))*flipud(plan.fhat);
+%f3=exp(-2*pi*1i*(tmpx*v.'))*fhatv;
+max(abs(f2-f3))
diff --git a/matlab/nnfft/test_nnfft2d_N215.m b/matlab/nnfft/test_nnfft2d_N215.m
new file mode 100644
index 0000000..98deda1
--- /dev/null
+++ b/matlab/nnfft/test_nnfft2d_N215.m
@@ -0,0 +1,86 @@
+
+% Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
+%
+% This program is free software; you can redistribute it and/or modify it under
+% the terms of the GNU General Public License as published by the Free Software
+% Foundation; either version 2 of the License, or (at your option) any later
+% version.
+%
+% This program 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, write to the Free Software Foundation, Inc., 51
+% Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+% Test script of class nnfft for spatial dimension d=2.
+clear all;
+
+%M=16; % number of nodes
+%N_1=24; % number of Fourier coefficients in first direction
+%N_2=32; % number of Fourier coefficients in second direction
+%N=[N_1;N_2];
+%N_total=N_1*N_2; % total number of Fourier coefficients
+M=2^15;
+N_1=2^15;
+N_2=2^15;
+N=[N_1;N_2];
+N_total=2^15;
+
+
+
+x=rand(M,2)-0.5; %nodes
+x(1,:)=[0.5;0.5];
+v=rand(N_total,2)-0.5; %nodes
+
+% Plan initialisation simple interface
+plan=nnfft(2,N_total,M,N); % create plan of class type nnfft
+
+plan
+% Plan initialisation guru interface
+%sigma=2; % oversampling factor
+%N1_1=sigma*N_1; % FFTW length, must be even natural number!
+%N1_2=sigma*N_2; % FFTW length, must be even natural number!
+%m=6; % window cut-off parameter
+%plan=nnfft(2,N_total,M,N,N1_1,N1_2,m,bitor(PRE_PHI_HUT,PRE_PSI)); % create plan of class type nnfft
+
+plan.x=x; % set nodes in plan
+plan.v=v; % set nodes in plan
+nnfft_precompute_psi(plan); % precomputations
+
+% NFFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%fhat=rand(N_1,N_2); % Fourier coefficients
+fhatv=rand(N_total,1)-0.5;
+%fhat=ones(N_total,1);
+
+%Test mit zweitem Einheitsvektor
+%fhat=zeros(N_total,1);
+%fhat(2)=1;
+%fhatv=fhat;%(:);
+
+% Compute samples with NNFFT
+ plan.fhat=fhatv; % set Fourier coefficients
+ nnfft_trafo(plan); % compute nonequispaced Fourier transform
+ f1=plan.f; % get samples
+
+% Compute samples direct
+nnfft_trafo_direct(plan);
+f2=plan.f;
+
+%% Compare results
+disp('NNFFT vs NNDFT');
+max(abs(f1-f2))
+
+
+%A=exp(-2*pi*1i*x*diag(N)*v');
+%f3=A*fhatv;
+
+%tmpv=(v.*repmat(N',size(v,1),1)).';
+
+%geflipter zweiter Einheitsvektor
+%f3=exp(-2*pi*1i*(x*tmpv))*flipud(plan.fhat);
+%f3=exp(-2*pi*1i*(tmpx*v.'))*fhatv;
+%max(abs(f2-f3))
diff --git a/nfft3.pc.in b/nfft3.pc.in
index 3fc7bfe..2c9c89d 100644
--- a/nfft3.pc.in
+++ b/nfft3.pc.in
@@ -6,6 +6,6 @@ includedir=${prefix}/include
Name: NFFT
Description: Nonuniform fast Fourier transform library
Version: @VERSION@
-Requires: fftw3
-Libs: -L${libdir} -lnfft3
+Requires: fftw3 at PREC_SUFFIX@
+Libs: -L${libdir} -lnfft3 at PREC_SUFFIX@
Cflags: -I${includedir}
diff --git a/support/Portfile b/support/Portfile
new file mode 100644
index 0000000..c3e2502
--- /dev/null
+++ b/support/Portfile
@@ -0,0 +1,59 @@
+# -*- coding: utf-8; mode: tcl; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
+# $Id$
+
+PortSystem 1.0
+PortGroup github 1.0
+
+github.setup NFFT nfft 3.3.0.alpha4
+github.tarball_from releases
+
+name nfft-3
+categories math
+license GPL-2+
+platforms darwin
+maintainers nfft.org:jens
+homepage http://www.nfft.org/
+distname ${version}
+worksrcdir nfft-${version}
+description Fast C routines to compute the Non-equispaced Discrete Fourier Transform
+long_description NFFT3 is a software library written in C for computing nonequispaced fast Fourier \
+ and related transformations. In detail, NFFT3 implements \
+ 1) The nonequispaced fast Fourier transform (NFFT) \
+ - the forward transform (NFFT) \
+ - the adjoint transform (adjoint NFFT) \
+ 2) Generalisations of the NFFT \
+ - to arbitrary knots in time and frequency domain (NNFFT) \
+ - to the sphere S^2 (NFSFT) \
+ - to the hyperbolic cross (NSFFT) \
+ - to real-valued data, i.e. (co)sine transforms, (NFCT, NFST) \
+ - to the rotation group (NFSOFT) \
+ 3) Generalised inverses based on iterative methods, e.g. CGNR, CGNE \
+ 4) Applications in \
+ - medical imaging \
+ (i) magnetic resonance imaging \
+ (ii) computerised tomography \
+ - summation schemes \
+ (i) fast Gauss transform (FGT) \
+ (ii) singular kernels \
+ (iii) zonal kernels \
+ - polar FFT, discrete Radon transform, ridgelet transform
+homepage http://www.nfft.org
+master_sites https://github.com/NFFT/nfft/archive/
+checksums rmd160 619d92a96e3037b918de0da87da78acf5fe35ca3 \
+ sha256 dd61562cca1345871aa847b70ecc256bd65c3c2d88476a14a17fea80285c3034
+depends_lib port:fftw-3
+use_autoreconf yes
+configure.args --enable-shared --enable-static --enable-all
+use_parallel_build yes
+
+variant gaussian description {compile NFFT with Gaussian window} conflicts bspline sinc {
+ configure.args-append --with-window=gaussian
+}
+
+variant bspline description {compile NFFT with B-Spline window} conflicts gaussian sinc {
+ configure.args-append --with-window=bspline
+}
+
+variant sinc description {compile NFFT with Sinc window} conflicts gaussian bspline {
+ configure.args-append --with-window=sinc
+}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/nfft.git
More information about the debian-science-commits
mailing list