Merging from nfft-3.3.

git-svn-id: https://euler.mathematik.tu-chemnitz.de/repos/nfft/lib/branches/NFFT-2@4236 c8b3e4ef-9df4-0310-9d76-b6dd82e1323a

configure.ac

@@ -556,6 +556,7 @@ AC_CONFIG_FILES(Makefile \
                 matlab/Makefile \
                 matlab/nfsft/Makefile \
                 matlab/nfft/Makefile \
+		matlab/nnfft/Makefile \
                 matlab/nfsft/@f_hat/Makefile \
                 matlab/nfsoft/Makefile \
                 doxygen/Makefile \

examples/nnfft/simple_test.c

@@ -34,7 +34,7 @@ void simple_test_nnfft_1d(void)
   nnfft_plan my_plan;                    /**< plan for the nfft                */
 
   int N[1];
-  N[0]=12;
+  N[0]=10;
 
   /** init an one dimensional plan */
   nnfft_init(&my_plan, 1, 3, 19, N);
@@ -51,18 +51,20 @@ void simple_test_nnfft_1d(void)
   }
 
   /** precompute psi, the entries of the matrix B */
-  if(my_plan.nnfft_flags & PRE_PSI)
-    nnfft_precompute_psi(&my_plan);
+/*  if(my_plan.nnfft_flags & PRE_PSI)
+        nnfft_precompute_psi(&my_plan);
 
   if(my_plan.nnfft_flags & PRE_FULL_PSI)
-    nnfft_precompute_full_psi(&my_plan);
-
+ 	nnfft_precompute_full_psi(&my_plan);
   if(my_plan.nnfft_flags & PRE_LIN_PSI)
-    nnfft_precompute_lin_psi(&my_plan);
-
+	 nnfft_precompute_lin_psi(&my_plan);
   /** precompute phi_hut, the entries of the matrix D */
-  if(my_plan.nnfft_flags & PRE_PHI_HUT)
-    nnfft_precompute_phi_hut(&my_plan);
+/*  if(my_plan.nnfft_flags & PRE_PHI_HUT)
+	  nnfft_precompute_phi_hut(&my_plan);
+*/
+
+nnfft_precompute_one_psi(&my_plan);
+
 
   /** init pseudo random Fourier coefficients and show them */
   for(k=0;k<my_plan.N_total;k++)
@@ -310,10 +312,10 @@ static void measure_time_nnfft_1d(void)
 int main(void)
 {
   system("clear");
-  printf("1) computing a one dimensional nndft, nnfft\n\n");
+  printf("1) computing a one dimensional nndft, nnfft SUSE\n\n");
   simple_test_nnfft_1d();
 
-  getc(stdin);
+  /*getc(stdin);
 
   system("clear");
   printf("1a) computing a one dimensional adjoint nndft, nnfft\n\n");
@@ -336,6 +338,6 @@ int main(void)
   system("clear");
   printf("4) computing times for one dimensional nnfft\n\n");
   measure_time_nnfft_1d();
-
+*/
   return 1;
 }

include/nfft3.h

@@ -410,6 +410,7 @@ typedef struct\
 } X(plan);\
 \
 NFFT_EXTERN void X(init)(X(plan) *ths_plan, int d, int N_total, int M_total, int *N); \
+NFFT_EXTERN void X(init_1d)(X(plan) *ths_plan, int N, int M_total); \
 NFFT_EXTERN void X(init_guru)(X(plan) *ths_plan, int d, int N_total, int M_total, \
   int *N, int *N1, int m, unsigned nnfft_flags); \
 NFFT_EXTERN void X(trafo_direct)(X(plan) *ths_plan); \
@@ -420,6 +421,7 @@ NFFT_EXTERN void X(precompute_lin_psi)(X(plan) *ths_plan); \
 NFFT_EXTERN void X(precompute_psi)(X(plan) *ths_plan); \
 NFFT_EXTERN void X(precompute_full_psi)(X(plan) *ths_plan); \
 NFFT_EXTERN void X(precompute_phi_hut)(X(plan) *ths_plan); \
+NFFT_EXTERN void X(precompute_one_psi)(X(plan) *ths);\
 NFFT_EXTERN void X(finalize)(X(plan) *ths_plan);
 
 #ifdef HAVE_NNFFT

kernel/nfft/nfft.c

@@ -1917,7 +1917,7 @@ static void nfft_adjoint_1d_compute_omp_blockwise(const C f, C *g,
 
   if (ar_u < ar_o)
   {
-    INT u = MAY(my_u0,ar_u);
+    INT u = MAX(my_u0,ar_u);
     INT o = MIN(my_o0,ar_o);
     INT offset_psij = u-ar_u;
 #ifdef OMP_ASSERT
@@ -1931,7 +1931,7 @@ static void nfft_adjoint_1d_compute_omp_blockwise(const C f, C *g,
   }
   else
   {
-    INT u = MAY(my_u0,ar_u);
+    INT u = MAX(my_u0,ar_u);
     INT o = my_o0;
     INT offset_psij = u-ar_u;
 #ifdef OMP_ASSERT
@@ -2773,7 +2773,7 @@ static void nfft_adjoint_2d_compute_omp_blockwise(const C f, C *g,
 
   if(ar_u0 < ar_o0)
   {
-    INT u0 = MAY(my_u0,ar_u0);
+    INT u0 = MAX(my_u0,ar_u0);
     INT o0 = MIN(my_o0,ar_o0);
     INT offset_psij = u0-ar_u0;
 #ifdef OMP_ASSERT
@@ -2793,7 +2793,7 @@ static void nfft_adjoint_2d_compute_omp_blockwise(const C f, C *g,
   }
   else
   {
-    INT u0 = MAY(my_u0,ar_u0);
+    INT u0 = MAX(my_u0,ar_u0);
     INT o0 = my_o0;
     INT offset_psij = u0-ar_u0;
 #ifdef OMP_ASSERT
@@ -4016,7 +4016,7 @@ static void nfft_adjoint_3d_compute_omp_blockwise(const C f, C *g,
 
   if(ar_u0<ar_o0)
   {
-    INT u0 = MAY(my_u0,ar_u0);
+    INT u0 = MAX(my_u0,ar_u0);
     INT o0 = MIN(my_o0,ar_o0);
     INT offset_psij = u0-ar_u0;
 #ifdef OMP_ASSERT
@@ -4042,7 +4042,7 @@ static void nfft_adjoint_3d_compute_omp_blockwise(const C f, C *g,
   }
   else
   {
-    INT u0 = MAY(my_u0,ar_u0);
+    INT u0 = MAX(my_u0,ar_u0);
     INT o0 = my_o0;
     INT offset_psij = u0-ar_u0;
 #ifdef OMP_ASSERT

kernel/nnfft/nnfft.c

@@ -302,6 +302,7 @@ void nnfft_trafo(nnfft_plan *ths)
     }
   }
 
+
   /* allows for external swaps of ths->f */
   ths->direct_plan->f = ths->f;
 
@@ -314,6 +315,7 @@ void nnfft_trafo(nnfft_plan *ths)
   }
 
   nnfft_D(ths);
+
 } /* nnfft_trafo */
 
 void nnfft_adjoint(nnfft_plan *ths)
@@ -478,6 +480,19 @@ void nnfft_precompute_full_psi(nnfft_plan *ths)
     } /* for(j) */
 }
 
+void nnfft_precompute_one_psi(nnfft_plan *ths)
+{
+  if(ths->nnfft_flags & PRE_PSI)
+    nnfft_precompute_psi(ths);
+  if(ths->nnfft_flags & PRE_FULL_PSI)
+    nnfft_precompute_full_psi(ths);
+  if(ths->nnfft_flags & PRE_LIN_PSI)
+    nnfft_precompute_lin_psi(ths);
+  /** precompute phi_hut, the entries of the matrix D */
+  if(ths->nnfft_flags & PRE_PHI_HUT)
+	  nnfft_precompute_phi_hut(ths);
+}
+
 static void nnfft_init_help(nnfft_plan *ths, int m2, unsigned nfft_flags, unsigned fftw_flags)
 {
   int t;                                /**< index over all dimensions       */
@@ -516,22 +531,28 @@ static void nnfft_init_help(nnfft_plan *ths, int m2, unsigned nfft_flags, unsign
 
   if(ths->nnfft_flags & MALLOC_X)
     ths->x = (double*)nfft_malloc(ths->d*ths->M_total*sizeof(double));
-  if(ths->nnfft_flags & MALLOC_F)
+  if(ths->nnfft_flags & MALLOC_F){
     ths->f=(double _Complex*)nfft_malloc(ths->M_total*sizeof(double _Complex));
-
+  }
   if(ths->nnfft_flags & MALLOC_V)
     ths->v = (double*)nfft_malloc(ths->d*ths->N_total*sizeof(double));
   if(ths->nnfft_flags & MALLOC_F_HAT)
     ths->f_hat = (double _Complex*)nfft_malloc(ths->N_total*sizeof(double _Complex));
 
+  //BUGFIX SUSE 2
+  /** precompute phi_hut, the entries of the matrix D */
+//  if(ths->nnfft_flags & PRE_PHI_HUT)
+//	  nnfft_precompute_phi_hut(ths);
+
   if(ths->nnfft_flags & PRE_LIN_PSI)
   {
     ths->K=(1U<< 10)*(ths->m+1);
     ths->psi = (double*) nfft_malloc((ths->K+1)*ths->d*sizeof(double));
   }
 
-  if(ths->nnfft_flags & PRE_PSI)
+  if(ths->nnfft_flags & PRE_PSI){
     ths->psi = (double*)nfft_malloc(ths->N_total*ths->d*(2*ths->m+2)*sizeof(double));
+  }
 
   if(ths->nnfft_flags & PRE_FULL_PSI)
   {
@@ -543,13 +564,11 @@ static void nnfft_init_help(nnfft_plan *ths, int m2, unsigned nfft_flags, unsign
       ths->psi_index_f = (int*) nfft_malloc(ths->N_total*sizeof(int));
       ths->psi_index_g = (int*) nfft_malloc(ths->N_total*lprod*sizeof(int));
   }
-
   ths->direct_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
-
   nfft_init_guru(ths->direct_plan, ths->d, ths->aN1, ths->M_total, N2, m2,
 		 nfft_flags, fftw_flags);
-
   ths->direct_plan->x = ths->x;
+
   ths->direct_plan->f = ths->f;
   ths->F = ths->direct_plan->f_hat;
 
@@ -598,26 +617,26 @@ void nnfft_init(nnfft_plan *ths, int d, int N_total, int M_total, int *N)
 
   unsigned nfft_flags;
   unsigned fftw_flags;
-
   ths->d = d;
   ths->M_total = M_total;
   ths->N_total = N_total;
-
   /* m should be greater to get the same accuracy as the nfft */
 /* Was soll dieser Ausdruck machen? Es handelt sich um eine Ganzzahl!
 
   WINDOW_HELP_ESTIMATE_m;
 */
+  //BUGFIX SUSE 1
+ths->m=WINDOW_HELP_ESTIMATE_m;
+
 
   ths->N = (int*) nfft_malloc(ths->d*sizeof(int));
   ths->N1 = (int*) nfft_malloc(ths->d*sizeof(int));
 
   for(t=0; t<d; t++) {
     ths->N[t] = N[t];
-
     /* the standard oversampling factor in the nnfft is 1.5 */
     ths->N1[t] = ceil(1.5*ths->N[t]);
-    
+
     /* N1 should be even */
     if(ths->N1[t]%2 != 0)
       ths->N1[t] = ths->N1[t] +1;
@@ -627,15 +646,17 @@ void nnfft_init(nnfft_plan *ths, int d, int N_total, int M_total, int *N)
   nfft_flags= PRE_PSI| PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE;
 
   fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
-
   nnfft_init_help(ths,ths->m,nfft_flags,fftw_flags);
 }
 
-void nnfft_finalize(nnfft_plan *ths)
+void nnfft_init_1d(nnfft_plan *ths,int N1, int M_total)
 {
+  nnfft_init(ths,1,N1,M_total,&N1);
+}
 
+void nnfft_finalize(nnfft_plan *ths)
+{
   nfft_finalize(ths->direct_plan);
-
   nfft_free(ths->direct_plan);
 
   nfft_free(ths->aN1);

matlab/Makefile.am

@@ -13,8 +13,9 @@ else
 endif
 
 DIR_NFFT=nfft
+DIR_NNFFT=nnfft
 
-SUBDIRS = . $(DIR_NFFT) $(DIR_NFSFT) $(DIR_NFSOFT)
+SUBDIRS = . $(DIR_NFFT) $(DIR_NFSFT) $(DIR_NFSOFT) $(DIR_NNFFT)
 
 AM_CPPFLAGS = -I$(top_srcdir)/include $(matlab_CPPFLAGS)
 

matlab/nnfft/Contents.m

@@ -0,0 +1,28 @@
+% NNFFT
+%
+% Files
+%   FFT_OUT_OF_PLACE    - FFT flag
+%   FFTW_ESTIMATE       - FFT flag
+%   FFTW_MEASURE        - FFT flag
+%   FG_PSI              - Precompuation flag
+%   nnfft_display	- display plan values
+%   nnfft_finalize       - Finalize plan
+%   nnfft_get_f          - Get function values from plan
+%   nnfft_get_f_hat      - Get Fourier coefficients from plan
+%   nnfft_get_x          - Get nodes from plan
+%   nnfft_init        - Initialise plans
+%   nnfft_init_guru        - Initialise plans
+%   nnfft_init_1d        - Initialise plans
+%   nnfft_set_f          - Set function values in plan
+%   nnfft_set_f_hat      - Set Fourier coefficients in plan
+%   nnfft_set_x          - Set nodes in plan
+%   nnfft_set_v          - Set nodes in plan
+%   nnfft_trafo          - nonequispaced fast Fourier transformat
+%   nnfft_trafo_direct   - nonequispaced fast diskret Fourie transformat
+%   nnfftmex             - Gateway function to NFFT module from NFFT3
+%   PRE_FG_PSI          - Precomputation flag
+%   PRE_FULL_PSI        - Precomputation flag
+%   PRE_LIN_PSI         - Precomputation flag
+%   PRE_PHI_HUT         - Precomputation flag
+%   PRE_PSI             - Precomputation flag
+%   simple_test         - Example program: Basic usage principles

matlab/nnfft/FFTW_ESTIMATE.m

@@ -0,0 +1,25 @@
+% FFTW_ESTIMATE FFT flag
+%   Valid for FFTW3
+%
+%   Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
+
+% 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: FFTW_ESTIMATE.m 3776 2012-06-03 13:29:25Z keiner $
+function f = FFTW_ESTIMATE()
+
+f = bitshift(1, 6);

matlab/nnfft/FFTW_MEASURE.m

@@ -0,0 +1,25 @@
+%FFTW_MEASURE FFT flag
+%   Valid for FFTW3
+%
+%   Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
+
+% 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: FFTW_MEASURE.m 3776 2012-06-03 13:29:25Z keiner $
+function f = FFTW_MEASURE()
+
+f = 0;

matlab/nnfft/FFT_OUT_OF_PLACE.m

@@ -0,0 +1,25 @@
+%FFT_OUT_OF_PLACE FFT flag
+%   If this flag is set, FFTW uses disjoint input/output vectors.
+%
+%   Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
+
+% 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: FFT_OUT_OF_PLACE.m 3776 2012-06-03 13:29:25Z keiner $
+function f = FFT_OUT_OF_PLACE()
+
+f = bitshift(1, 9);