NFFT  3.3.1
nsfft.c
00001 /*
00002  * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 #include "config.h"
00019 
00020 #include <stdio.h>
00021 #include <math.h>
00022 #include <string.h>
00023 #include <stdlib.h>
00024 #ifdef HAVE_COMPLEX_H
00025 #include <complex.h>
00026 #endif
00027 #include "nfft3.h"
00028 #include "infft.h"
00029 
00030 #define NSFTT_DISABLE_TEST
00031 
00032 /* computes a 2d ndft by 1d nfft along the dimension 1 times
00033    1d ndft along dimension 0
00034 */
00035 static void short_nfft_trafo_2d(nfft_plan* ths, nfft_plan* plan_1d)
00036 {
00037   int j,k0;
00038   double omega;
00039 
00040   for(j=0;j<ths->M_total;j++)
00041     {
00042       ths->f[j]= 0;
00043       plan_1d->x[j] = ths->x[ths->d * j + 1];
00044     }
00045 
00046   for(k0=0;k0<ths->N[0];k0++) /* for shorties */
00047     {
00048       plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
00049 
00050       nfft_trafo(plan_1d);
00051 
00052       for(j=0;j<ths->M_total;j++)
00053   {
00054     omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00055           ths->f[j] += plan_1d->f[j] * cexp( - I*2*KPI*omega);
00056   }
00057     }
00058 }
00059 
00060 static void short_nfft_adjoint_2d(nfft_plan* ths, nfft_plan* plan_1d)
00061 {
00062   int j,k0;
00063   double omega;
00064 
00065   for(j=0;j<ths->M_total;j++)
00066     plan_1d->x[j] = ths->x[ths->d * j + 1];
00067 
00068   for(k0=0;k0<ths->N[0];k0++) /* for shorties */
00069     {
00070       for(j=0;j<ths->M_total;j++)
00071   {
00072     omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00073           plan_1d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
00074   }
00075 
00076       plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
00077 
00078       nfft_adjoint(plan_1d);
00079     }
00080 }
00081 
00082 /* computes a 3d ndft by 1d nfft along the dimension 2 times
00083    2d ndft along dimension 0,1
00084 */
00085 static void short_nfft_trafo_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
00086 {
00087   int j,k0,k1;
00088   double omega;
00089 
00090   for(j=0;j<ths->M_total;j++)
00091     {
00092       ths->f[j] = 0;
00093       plan_1d->x[j] = ths->x[ths->d * j + 2];
00094     }
00095 
00096   for(k0=0;k0<ths->N[0];k0++) /* for shorties */
00097     for(k1=0;k1<ths->N[1];k1++)
00098       {
00099   plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
00100 
00101   nfft_trafo(plan_1d);
00102 
00103   for(j=0;j<ths->M_total;j++)
00104     {
00105       omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
00106         +     ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
00107             ths->f[j] += plan_1d->f[j] * cexp( - I*2*KPI*omega);
00108     }
00109       }
00110 }
00111 
00112 static void short_nfft_adjoint_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
00113 {
00114   int j,k0,k1;
00115   double omega;
00116 
00117   for(j=0;j<ths->M_total;j++)
00118     plan_1d->x[j] = ths->x[ths->d * j + 2];
00119 
00120   for(k0=0;k0<ths->N[0];k0++) /* for shorties */
00121     for(k1=0;k1<ths->N[1];k1++)
00122       {
00123   for(j=0;j<ths->M_total;j++)
00124     {
00125       omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
00126         +     ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
00127             plan_1d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
00128     }
00129 
00130   plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
00131 
00132   nfft_adjoint(plan_1d);
00133       }
00134 }
00135 
00136 /* computes a 3d ndft by 2d nfft along the dimension 1,2 times
00137    1d ndft along dimension 0
00138 */
00139 static void short_nfft_trafo_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
00140 {
00141   int j,k0;
00142   double omega;
00143 
00144   for(j=0;j<ths->M_total;j++)
00145     {
00146       ths->f[j] = 0;
00147       plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
00148       plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
00149     }
00150 
00151   for(k0=0;k0<ths->N[0];k0++) /* for shorties */
00152     {
00153       plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
00154 
00155       nfft_trafo(plan_2d);
00156 
00157       for(j=0;j<ths->M_total;j++)
00158   {
00159     omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00160     ths->f[j] += plan_2d->f[j] * cexp( - I*2*KPI*omega);
00161   }
00162     }
00163 }
00164 
00165 static void short_nfft_adjoint_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
00166 {
00167   int j,k0;
00168   double omega;
00169 
00170   for(j=0;j<ths->M_total;j++)
00171     {
00172       plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
00173       plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
00174     }
00175 
00176   for(k0=0;k0<ths->N[0];k0++) /* for shorties */
00177     {
00178       for(j=0;j<ths->M_total;j++)
00179   {
00180     omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00181     plan_2d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
00182   }
00183 
00184       plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
00185 
00186       nfft_adjoint(plan_2d);
00187     }
00188 }
00189 
00190 /*---------------------------------------------------------------------------*/
00191 
00192 #ifdef GAUSSIAN
00193 static int index_sparse_to_full_direct_2d(int J, int k)
00194 {
00195     int N=X(exp2i)(J+2);               /* number of full coeffs             */
00196     int N_B=X(exp2i)(J);               /* number in each sparse block       */
00197 
00198     int j=k/N_B;                        /* consecutive number of Block       */
00199     int r=j/4;                          /* level of block                    */
00200 
00201     int i, o, a, b,s,l,m1,m2;
00202     int k1,k2;
00203 
00204     if (k>=(J+4)*X(exp2i)(J+1))
00205       {
00206   printf("Fehler!\n");
00207   return(-1);
00208       }
00209     else
00210       {
00211   if (r>(J+1)/2)                  /* center block                      */
00212     {
00213       i=k-4*((J+1)/2+1)*N_B;
00214       a=X(exp2i)(J/2+1);
00215       m1=i/a;
00216       m2=i%a;
00217       k1=N/2-a/2+m1;
00218       k2=N/2-a/2+m2;
00219     }
00220   else                            /* no center block                   */
00221     {
00222       i=k-j*N_B;                  /* index in specific block           */
00223       o=j%4;                      /* kind of specific block            */
00224       a=X(exp2i)(r);
00225       b=X(exp2i)(J-r);
00226       l=MAX(a,b);                 /* long dimension of block           */
00227       s=MIN(a,b);                 /* short dimension of block          */
00228       m1=i/l;
00229       m2=i%l;
00230 
00231       switch(o)
00232         {
00233         case 0:
00234     {
00235       k1=N/2-a/2 ;
00236       k2=N/2+ b  ;
00237 
00238       if (b>=a)
00239         {
00240           k1+=m1;
00241           k2+=m2;
00242         }
00243       else
00244         {
00245           k1+=m2;
00246           k2+=m1;
00247         }
00248       break;
00249     }
00250         case 1:
00251     {
00252       k1=N/2+ b  ;
00253       k2=N/2-a/2 ;
00254 
00255       if (b>a)
00256         {
00257           k1+=m2;
00258           k2+=m1;
00259         }
00260       else
00261         {
00262           k1+=m1;
00263           k2+=m2;
00264         }
00265       break;
00266     }
00267         case 2:
00268     {
00269       k1=N/2-a/2 ;
00270       k2=N/2-2*b ;
00271 
00272       if (b>=a)
00273         {
00274           k1+=m1;
00275           k2+=m2;
00276         }
00277       else
00278         {
00279           k1+=m2;
00280           k2+=m1;
00281         }
00282       break;
00283     }
00284         case 3:
00285     {
00286       k1=N/2-2*b ;
00287       k2=N/2-a/2 ;
00288 
00289       if (b>a)
00290         {
00291           k1+=m2;
00292           k2+=m1;
00293         }
00294       else
00295         {
00296           k1+=m1;
00297           k2+=m2;
00298         }
00299       break;
00300     }
00301         default:
00302     {
00303       k1=-1;
00304       k2=-1;
00305     }
00306         }
00307     }
00308   //printf("m1=%d, m2=%d\n",m1,m2);
00309   return(k1*N+k2);
00310       }
00311 }
00312 #endif
00313 
00314 static inline int index_sparse_to_full_2d(nsfft_plan *ths, int k)
00315 {
00316   /* only by lookup table */
00317   if( k < ths->N_total)
00318     return ths->index_sparse_to_full[k];
00319   else
00320     return -1;
00321 }
00322 
00323 #ifndef NSFTT_DISABLE_TEST
00324 static int index_full_to_sparse_2d(int J, int k)
00325 {
00326     int N=X(exp2i)(J+2);               /* number of full coeffs       */
00327     int N_B=X(exp2i)(J);               /* number in each sparse block */
00328 
00329     int k1=k/N-N/2;                     /* coordinates in the full grid */
00330     int k2=k%N-N/2;                     /* k1: row, k2: column          */
00331 
00332     int r,a,b;
00333 
00334     a=X(exp2i)(J/2+1);
00335 
00336     if ( (k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) )
00337       {
00338   return(4*((J+1)/2+1)*N_B+(k1+a/2)*a+(k2+a/2));
00339       }
00340 
00341     for (r=0; r<=(J+1)/2; r++)
00342       {
00343   b=X(exp2i)(r);
00344   a=X(exp2i)(J-r);
00345   if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=a) && (k2<2*a) )
00346     {
00347             if (a>=b)
00348         return((4*r+0)*N_B+(k1+b/2)*a+(k2-a));
00349       else
00350         return((4*r+0)*N_B+(k2-a)*b+(k1+b/2));
00351     }
00352   else if ( (k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
00353     {
00354             if (a>b)
00355         return((4*r+1)*N_B+(k2+b/2)*a+(k1-a));
00356       else
00357         return((4*r+1)*N_B+(k1-a)*b+(k2+b/2));
00358     }
00359   else if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=-2*a) && (k2<-a) )
00360     {
00361             if (a>=b)
00362         return((4*r+2)*N_B+(k1+b/2)*a+(k2+2*a));
00363       else
00364         return((4*r+2)*N_B+(k2+2*a)*b+(k1+b/2));
00365     }
00366   else if ( (k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
00367     {
00368             if (a>b)
00369         return((4*r+3)*N_B+(k2+b/2)*a+(k1+2*a));
00370       else
00371         return((4*r+3)*N_B+(k1+2*a)*b+(k2+b/2));
00372     }
00373       }
00374 
00375     return(-1);
00376 }
00377 #endif
00378 
00379 #ifdef GAUSSIAN
00380 static void init_index_sparse_to_full_2d(nsfft_plan *ths)
00381 {
00382   int k_S;
00383 
00384   for (k_S=0; k_S<ths->N_total; k_S++)
00385     ths->index_sparse_to_full[k_S]=index_sparse_to_full_direct_2d(ths->J, k_S);
00386 }
00387 #endif
00388 
00389 #ifdef GAUSSIAN
00390 static inline int index_sparse_to_full_3d(nsfft_plan *ths, int k)
00391 {
00392   /* only by lookup table */
00393   if( k < ths->N_total)
00394     return ths->index_sparse_to_full[k];
00395   else
00396     return -1;
00397 }
00398 #endif
00399 
00400 #ifndef NSFTT_DISABLE_TEST
00401 static int index_full_to_sparse_3d(int J, int k)
00402 {
00403   int N=X(exp2i)(J+2);                 /* length of the full grid           */
00404   int N_B_r;                            /* size of a sparse block in level r */
00405   int sum_N_B_less_r;                   /* sum N_B_r                         */
00406 
00407   int r,a,b;
00408 
00409   int k3=(k%N)-N/2;                     /* coordinates in the full grid      */
00410   int k2=((k/N)%N)-N/2;
00411   int k1=k/(N*N)-N/2;
00412 
00413   a=X(exp2i)(J/2+1);                   /* length of center block            */
00414 
00415   if((k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) && (k3>=(-a/2)) &&
00416      (k3<a/2))
00417     {
00418       return(6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+((k1+a/2)*a+(k2+a/2))*a+
00419              (k3+a/2));
00420     }
00421 
00422   sum_N_B_less_r=0;
00423   for (r=0; r<=(J+1)/2; r++)
00424     {
00425       a=X(exp2i)(J-r);
00426       b=X(exp2i)(r);
00427 
00428       N_B_r=a*b*b;
00429 
00430       /* right - rear - top - left - front - bottom */
00431       if ((k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
00432           (k3>=-(b/2)) && (k3<(b+1)/2)) /* right */
00433   {
00434     if(a>b)
00435       return sum_N_B_less_r+N_B_r*0 + ((k2+b/2)*b+k3+b/2)*a + (k1-a);
00436     else
00437       return sum_N_B_less_r+N_B_r*0 + ((k1-a)*b+(k2+b/2))*b + (k3+b/2);
00438   }
00439       else if ((k2>=a) && (k2<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00440                (k3>=-(b/2)) && (k3<(b+1)/2)) /* rear */
00441   {
00442     if(a>b)
00443       return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+k3+b/2)*a + (k2-a);
00444     else if (a==b)
00445       return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+(k2-a))*a + (k3+b/2);
00446     else
00447       return sum_N_B_less_r+N_B_r*1 + ((k2-a)*b+(k1+b/2))*b + (k3+b/2);
00448   }
00449        else if ((k3>=a) && (k3<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00450                 (k2>=-(b/2)) && (k2<(b+1)/2)) /* top */
00451   {
00452     if(a>=b)
00453       return sum_N_B_less_r+N_B_r*2 + ((k1+b/2)*b+k2+b/2)*a + (k3-a);
00454     else
00455       return sum_N_B_less_r+N_B_r*2 + ((k3-a)*b+(k1+b/2))*b + (k2+b/2);
00456   }
00457 
00458       else if ((k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
00459                (k3>=-(b/2)) && (k3<(b+1)/2)) /* left */
00460   {
00461     if(a>b)
00462       return sum_N_B_less_r+N_B_r*3 + ((k2+b/2)*b+k3+b/2)*a + (k1+2*a);
00463     else
00464       return sum_N_B_less_r+N_B_r*3 + ((k1+2*a)*b+(k2+b/2))*b + (k3+b/2);
00465   }
00466       else if ((k2>=-2*a) && (k2<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00467                (k3>=-(b/2)) && (k3<(b+1)/2)) /* front */
00468   {
00469     if(a>b)
00470       return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+k3+b/2)*a + (k2+2*a);
00471     else if (a==b)
00472       return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+(k2+2*a))*a + (k3+b/2);
00473     else
00474       return sum_N_B_less_r+N_B_r*4 + ((k2+2*a)*b+(k1+b/2))*b + (k3+b/2);
00475   }
00476        else if ((k3>=-2*a) && (k3<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00477                 (k2>=-(b/2)) && (k2<(b+1)/2)) /* bottom */
00478   {
00479     if(a>=b)
00480       return sum_N_B_less_r+N_B_r*5 + ((k1+b/2)*b+k2+b/2)*a + (k3+2*a);
00481     else
00482       return sum_N_B_less_r+N_B_r*5 + ((k3+2*a)*b+(k1+b/2))*b + (k2+b/2);
00483   }
00484 
00485       sum_N_B_less_r+=6*N_B_r;
00486     } /* for(r) */
00487 
00488   return(-1);
00489 }
00490 #endif
00491 
00492 #ifdef GAUSSIAN
00493 static void init_index_sparse_to_full_3d(nsfft_plan *ths)
00494 {
00495   int k1,k2,k3,k_s,r;
00496   int a,b;
00497   int N=X(exp2i)(ths->J+2);            /* length of the full grid           */
00498   int Nc=ths->center_nfft_plan->N[0];   /* length of the center block        */
00499 
00500   for (k_s=0, r=0; r<=(ths->J+1)/2; r++)
00501     {
00502       a=X(exp2i)(ths->J-r);
00503       b=X(exp2i)(r);
00504 
00505       /* right - rear - top - left - front - bottom */
00506 
00507       /* right */
00508       if(a>b)
00509   for(k2=-b/2;k2<(b+1)/2;k2++)
00510     for(k3=-b/2;k3<(b+1)/2;k3++)
00511       for(k1=a; k1<2*a; k1++,k_s++)
00512         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00513       else
00514   for(k1=a; k1<2*a; k1++)
00515     for(k2=-b/2;k2<(b+1)/2;k2++)
00516       for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00517         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00518 
00519       /* rear */
00520       if(a>b)
00521   for(k1=-b/2;k1<(b+1)/2;k1++)
00522     for(k3=-b/2;k3<(b+1)/2;k3++)
00523       for(k2=a; k2<2*a; k2++,k_s++)
00524         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00525       else if(a==b)
00526   for(k1=-b/2;k1<(b+1)/2;k1++)
00527     for(k2=a; k2<2*a; k2++)
00528       for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00529         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00530       else
00531   for(k2=a; k2<2*a; k2++)
00532     for(k1=-b/2;k1<(b+1)/2;k1++)
00533       for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00534         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00535 
00536       /* top */
00537       if(a>=b)
00538   for(k1=-b/2;k1<(b+1)/2;k1++)
00539     for(k2=-b/2;k2<(b+1)/2;k2++)
00540       for(k3=a; k3<2*a; k3++,k_s++)
00541         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00542       else
00543   for(k3=a; k3<2*a; k3++)
00544     for(k1=-b/2;k1<(b+1)/2;k1++)
00545       for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
00546         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00547 
00548       /* left */
00549       if(a>b)
00550   for(k2=-b/2;k2<(b+1)/2;k2++)
00551     for(k3=-b/2;k3<(b+1)/2;k3++)
00552       for(k1=-2*a; k1<-a; k1++,k_s++)
00553         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00554       else
00555   for(k1=-2*a; k1<-a; k1++)
00556     for(k2=-b/2;k2<(b+1)/2;k2++)
00557       for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00558         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00559 
00560       /* front */
00561       if(a>b)
00562   for(k1=-b/2;k1<(b+1)/2;k1++)
00563     for(k3=-b/2;k3<(b+1)/2;k3++)
00564       for(k2=-2*a; k2<-a; k2++,k_s++)
00565         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00566       else if(a==b)
00567   for(k1=-b/2;k1<(b+1)/2;k1++)
00568     for(k2=-2*a; k2<-a; k2++)
00569       for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00570         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00571       else
00572   for(k2=-2*a; k2<-a; k2++)
00573     for(k1=-b/2;k1<(b+1)/2;k1++)
00574       for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00575         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00576 
00577       /* top */
00578       if(a>=b)
00579   for(k1=-b/2;k1<(b+1)/2;k1++)
00580     for(k2=-b/2;k2<(b+1)/2;k2++)
00581       for(k3=-2*a; k3<-a; k3++,k_s++)
00582         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00583       else
00584   for(k3=-2*a; k3<-a; k3++)
00585     for(k1=-b/2;k1<(b+1)/2;k1++)
00586       for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
00587         ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00588     }
00589 
00590   /* center */
00591   for(k1=-Nc/2;k1<Nc/2;k1++)
00592     for(k2=-Nc/2;k2<Nc/2;k2++)
00593       for(k3=-Nc/2; k3<Nc/2; k3++,k_s++)
00594   ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00595 }
00596 #endif
00597 
00598 /* copies ths->f_hat to ths_plan->f_hat */
00599 void nsfft_cp(nsfft_plan *ths, nfft_plan *ths_full_plan)
00600 {
00601   int k;
00602 
00603   /* initialize f_hat with zero values */
00604   memset(ths_full_plan->f_hat, 0, ths_full_plan->N_total*sizeof(double _Complex));
00605 
00606    /* copy values at hyperbolic grid points */
00607   for(k=0;k<ths->N_total;k++)
00608     ths_full_plan->f_hat[ths->index_sparse_to_full[k]]=ths->f_hat[k];
00609 
00610   /* copy nodes */
00611   memcpy(ths_full_plan->x,ths->act_nfft_plan->x,ths->M_total*ths->d*sizeof(double));
00612 }
00613 
00614 #ifndef NSFTT_DISABLE_TEST
00615 /* test copy_sparse_to_full */
00616 static void test_copy_sparse_to_full_2d(nsfft_plan *ths, nfft_plan *ths_full_plan)
00617 {
00618   int r;
00619   int k1, k2;
00620   int a,b;
00621   const int J=ths->J;   /* N=2^J                  */
00622   const int N=ths_full_plan->N[0];  /* size of full NFFT      */
00623   const int N_B=X(exp2i)(J);        /* size of small blocks   */
00624 
00625   /* copy sparse plan to full plan */
00626   nsfft_cp(ths, ths_full_plan);
00627 
00628   /* show blockwise f_hat */
00629   printf("f_hat blockwise\n");
00630   for (r=0; r<=(J+1)/2; r++){
00631     a=X(exp2i)(J-r); b=X(exp2i)(r);
00632 
00633     printf("top\n");
00634     for (k1=0; k1<a; k1++){
00635       for (k2=0; k2<b; k2++){
00636         printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]),
00637                            cimag(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]));
00638       }
00639       printf("\n");
00640     }
00641 
00642     printf("bottom\n");
00643     for (k1=0; k1<a; k1++){
00644       for (k2=0; k2<b; k2++){
00645         printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]),
00646                                  cimag(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]));
00647       }
00648       printf("\n");
00649     }
00650 
00651     printf("right\n");
00652     for (k2=0; k2<b; k2++){
00653       for (k1=0; k1<a; k1++){
00654         printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]),
00655                                  cimag(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]));
00656       }
00657       printf("\n");
00658     }
00659 
00660     printf("left\n");
00661     for (k2=0; k2<b; k2++){
00662       for (k1=0; k1<a; k1++){
00663         printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]),
00664                            cimag(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]));
00665       }
00666       printf("\n");
00667     }
00668   }
00669 
00670   return;
00671   /* show full f_hat */
00672   printf("full f_hat\n");
00673   for (k1=0;k1<N;k1++){
00674     for (k2=0;k2<N;k2++){
00675       printf("(%1.1f,%1.1f) ", creal(ths_full_plan->f_hat[k1*N+k2]),
00676                                cimag(ths_full_plan->f_hat[k1*N+k2]));
00677     }
00678     printf("\n");
00679   }
00680 }
00681 #endif
00682 
00683 #ifndef NSFTT_DISABLE_TEST
00684 static void test_sparse_to_full_2d(nsfft_plan* ths)
00685 {
00686   int k_S,k1,k2;
00687   int N=X(exp2i)(ths->J+2);
00688 
00689   printf("N=%d\n\n",N);
00690 
00691   for(k1=0;k1<N;k1++)
00692     for(k2=0;k2<N;k2++)
00693       {
00694         k_S=index_full_to_sparse_2d(ths->J, k1*N+k2);
00695   if(k_S!=-1)
00696     printf("(%+d, %+d)\t= %+d \t= %+d = %+d \n",k1-N/2,k2-N/2,
00697                  k1*N+k2, k_S, ths->index_sparse_to_full[k_S]);
00698       }
00699 }
00700 #endif
00701 
00702 #ifndef NSFTT_DISABLE_TEST
00703 static void test_sparse_to_full_3d(nsfft_plan* ths)
00704 {
00705   int k_S,k1,k2,k3;
00706   int N=X(exp2i)(ths->J+2);
00707 
00708   printf("N=%d\n\n",N);
00709 
00710   for(k1=0;k1<N;k1++)
00711     for(k2=0;k2<N;k2++)
00712         for(k3=0;k3<N;k3++)
00713     {
00714       k_S=index_full_to_sparse_3d(ths->J, (k1*N+k2)*N+k3);
00715       if(k_S!=-1)
00716         printf("(%d, %d, %d)\t= %d \t= %d = %d \n",k1-N/2,k2-N/2,k3-N/2,
00717                      (k1*N+k2)*N+k3,k_S, ths->index_sparse_to_full[k_S]);
00718     }
00719 }
00720 #endif
00721 
00722 
00723 void nsfft_init_random_nodes_coeffs(nsfft_plan *ths)
00724 {
00725   int j;
00726 
00727   /* init frequencies */
00728   nfft_vrand_unit_complex(ths->f_hat, ths->N_total);
00729 
00730   /* init nodes */
00731   nfft_vrand_shifted_unit_double(ths->act_nfft_plan->x, ths->d * ths->M_total);
00732 
00733   if(ths->d==2)
00734     for(j=0;j<ths->M_total;j++)
00735       {
00736         ths->x_transposed[2*j+0]=ths->act_nfft_plan->x[2*j+1];
00737         ths->x_transposed[2*j+1]=ths->act_nfft_plan->x[2*j+0];
00738       }
00739   else /* this->d==3 */
00740     for(j=0;j<ths->M_total;j++)
00741       {
00742         ths->x_102[3*j+0]=ths->act_nfft_plan->x[3*j+1];
00743         ths->x_102[3*j+1]=ths->act_nfft_plan->x[3*j+0];
00744         ths->x_102[3*j+2]=ths->act_nfft_plan->x[3*j+2];
00745 
00746         ths->x_201[3*j+0]=ths->act_nfft_plan->x[3*j+2];
00747         ths->x_201[3*j+1]=ths->act_nfft_plan->x[3*j+0];
00748         ths->x_201[3*j+2]=ths->act_nfft_plan->x[3*j+1];
00749 
00750         ths->x_120[3*j+0]=ths->act_nfft_plan->x[3*j+1];
00751         ths->x_120[3*j+1]=ths->act_nfft_plan->x[3*j+2];
00752         ths->x_120[3*j+2]=ths->act_nfft_plan->x[3*j+0];
00753 
00754         ths->x_021[3*j+0]=ths->act_nfft_plan->x[3*j+0];
00755         ths->x_021[3*j+1]=ths->act_nfft_plan->x[3*j+2];
00756         ths->x_021[3*j+2]=ths->act_nfft_plan->x[3*j+1];
00757       }
00758 }
00759 
00760 static void nsdft_trafo_2d(nsfft_plan *ths)
00761 {
00762   int j,k_S,k_L,k0,k1;
00763   double omega;
00764   int N=X(exp2i)(ths->J+2);
00765 
00766   memset(ths->f,0,ths->M_total*sizeof(double _Complex));
00767 
00768   for(k_S=0;k_S<ths->N_total;k_S++)
00769     {
00770       k_L=ths->index_sparse_to_full[k_S];
00771       k0=k_L / N;
00772       k1=k_L % N;
00773 
00774       for(j=0;j<ths->M_total;j++)
00775   {
00776     omega =
00777       ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
00778       ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
00779           ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*KPI*omega);
00780   }
00781     }
00782 } /* void nsdft_trafo_2d */
00783 
00784 static void nsdft_trafo_3d(nsfft_plan *ths)
00785 {
00786   int j,k_S,k0,k1,k2;
00787   double omega;
00788   int N=X(exp2i)(ths->J+2);
00789   int k_L;
00790 
00791   memset(ths->f,0,ths->M_total*sizeof(double _Complex));
00792 
00793   for(k_S=0;k_S<ths->N_total;k_S++)
00794     {
00795       k_L=ths->index_sparse_to_full[k_S];
00796 
00797       k0=k_L/(N*N);
00798       k1=(k_L/N)%N;
00799       k2=k_L%N;
00800 
00801       for(j=0;j<ths->M_total;j++)
00802   {
00803     omega =
00804       ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
00805       ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
00806       ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
00807           ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*KPI*omega);
00808   }
00809     }
00810 } /* void nsdft_trafo_3d */
00811 
00812 void nsfft_trafo_direct(nsfft_plan *ths)
00813 {
00814   if(ths->d==2)
00815     nsdft_trafo_2d(ths);
00816   else
00817     nsdft_trafo_3d(ths);
00818 }
00819 
00820 static void nsdft_adjoint_2d(nsfft_plan *ths)
00821 {
00822   int j,k_S,k_L,k0,k1;
00823   double omega;
00824   int N=X(exp2i)(ths->J+2);
00825 
00826   memset(ths->f_hat,0,ths->N_total*sizeof(double _Complex));
00827 
00828   for(k_S=0;k_S<ths->N_total;k_S++)
00829     {
00830       k_L=ths->index_sparse_to_full[k_S];
00831       k0=k_L / N;
00832       k1=k_L % N;
00833 
00834       for(j=0;j<ths->M_total;j++)
00835   {
00836     omega =
00837       ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
00838       ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
00839           ths->f_hat[k_S] += ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
00840   }
00841     }
00842 } /* void nsdft_adjoint_2d */
00843 
00844 static void nsdft_adjoint_3d(nsfft_plan *ths)
00845 {
00846   int j,k_S,k0,k1,k2;
00847   double omega;
00848   int N=X(exp2i)(ths->J+2);
00849   int k_L;
00850 
00851   memset(ths->f_hat,0,ths->N_total*sizeof(double _Complex));
00852 
00853   for(k_S=0;k_S<ths->N_total;k_S++)
00854     {
00855       k_L=ths->index_sparse_to_full[k_S];
00856 
00857       k0=k_L/(N*N);
00858       k1=(k_L/N)%N;
00859       k2=k_L%N;
00860 
00861       for(j=0;j<ths->M_total;j++)
00862   {
00863     omega =
00864       ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
00865       ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
00866       ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
00867           ths->f_hat[k_S] += ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
00868   }
00869     }
00870 } /* void nsdft_adjoint_3d */
00871 
00872 void nsfft_adjoint_direct(nsfft_plan *ths)
00873 {
00874   if(ths->d==2)
00875     nsdft_adjoint_2d(ths);
00876   else
00877     nsdft_adjoint_3d(ths);
00878 }
00879 
00880 static void nsfft_trafo_2d(nsfft_plan *ths)
00881 {
00882   int r,rr,j;
00883   double temp;
00884 
00885   int M=ths->M_total;
00886   int J=ths->J;
00887 
00888   /* center */
00889   ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*X(exp2i)(J);
00890 
00891   if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
00892     nfft_trafo_direct(ths->center_nfft_plan);
00893   else
00894     nfft_trafo(ths->center_nfft_plan);
00895 
00896   for (j=0; j<M; j++)
00897     ths->f[j] = ths->center_nfft_plan->f[j];
00898 
00899   for(rr=0;rr<=(J+1)/2;rr++)
00900     {
00901       r=MIN(rr,J-rr);
00902       ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[r];
00903       ths->act_nfft_plan->N[0]=X(exp2i)(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
00904       ths->act_nfft_plan->N[1]=X(exp2i)(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
00905 
00906       /*printf("%d x %d\n",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1]);*/
00907 
00908       temp=-3.0*KPI*X(exp2i)(J-rr);
00909 
00910       /* right */
00911       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*X(exp2i)(J);
00912 
00913       if(r<rr)
00914   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00915 
00916       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00917   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00918     nfft_trafo_direct(ths->act_nfft_plan);
00919   else
00920     short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00921       else
00922   nfft_trafo(ths->act_nfft_plan);
00923 
00924       if(r<rr)
00925   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00926 
00927       for (j=0; j<M; j++)
00928         ths->f[j] +=  ths->act_nfft_plan->f[j] *
00929                       cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
00930 
00931       /* top */
00932       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*X(exp2i)(J);
00933 
00934       if((r==rr)&&(J-rr!=rr))
00935   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00936 
00937       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00938   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00939     nfft_trafo_direct(ths->act_nfft_plan);
00940   else
00941     short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00942       else
00943   nfft_trafo(ths->act_nfft_plan);
00944 
00945       if((r==rr)&&(J-rr!=rr))
00946   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00947 
00948       for (j=0; j<M; j++)
00949         ths->f[j] += ths->act_nfft_plan->f[j] *
00950                      cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
00951 
00952       /* left */
00953       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*X(exp2i)(J);
00954 
00955       if(r<rr)
00956   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00957 
00958       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00959   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00960     nfft_trafo_direct(ths->act_nfft_plan);
00961   else
00962     short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00963       else
00964   nfft_trafo(ths->act_nfft_plan);
00965 
00966       if(r<rr)
00967   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00968 
00969       for (j=0; j<M; j++)
00970         ths->f[j] += ths->act_nfft_plan->f[j] *
00971                      cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
00972 
00973       /* bottom */
00974       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*X(exp2i)(J);
00975 
00976       if((r==rr)&&(J-rr!=rr))
00977   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00978 
00979       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00980   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00981     nfft_trafo_direct(ths->act_nfft_plan);
00982   else
00983     short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00984       else
00985   nfft_trafo(ths->act_nfft_plan);
00986 
00987       if((r==rr)&&(J-rr!=rr))
00988   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
00989 
00990       for (j=0; j<M; j++)
00991         ths->f[j] += ths->act_nfft_plan->f[j] *
00992                      cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
00993     } /* for(rr) */
00994 } /* void nsfft_trafo_2d */
00995 
00996 static void nsfft_adjoint_2d(nsfft_plan *ths)
00997 {
00998   int r,rr,j;
00999   double temp;
01000 
01001   int M=ths->M_total;
01002   int J=ths->J;
01003 
01004   /* center */
01005   for (j=0; j<M; j++)
01006     ths->center_nfft_plan->f[j] = ths->f[j];
01007 
01008   ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*X(exp2i)(J);
01009 
01010   if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
01011     nfft_adjoint_direct(ths->center_nfft_plan);
01012   else
01013     nfft_adjoint(ths->center_nfft_plan);
01014 
01015   for(rr=0;rr<=(J+1)/2;rr++)
01016     {
01017       r=MIN(rr,J-rr);
01018       ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[r];
01019       ths->act_nfft_plan->N[0]=X(exp2i)(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
01020       ths->act_nfft_plan->N[1]=X(exp2i)(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
01021 
01022       /*printf("%d x %d\n",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1]);*/
01023 
01024       temp=-3.0*KPI*X(exp2i)(J-rr);
01025 
01026       /* right */
01027       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*X(exp2i)(J);
01028 
01029       for (j=0; j<M; j++)
01030         ths->act_nfft_plan->f[j]= ths->f[j] *
01031                                   cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
01032 
01033       if(r<rr)
01034   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01035 
01036       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01037   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01038     nfft_adjoint_direct(ths->act_nfft_plan);
01039   else
01040     short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01041       else
01042   nfft_adjoint(ths->act_nfft_plan);
01043 
01044       if(r<rr)
01045   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01046 
01047       /* top */
01048       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*X(exp2i)(J);
01049 
01050       for (j=0; j<M; j++)
01051         ths->act_nfft_plan->f[j]= ths->f[j] *
01052                                   cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
01053 
01054       if((r==rr)&&(J-rr!=rr))
01055   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01056 
01057       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01058   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01059     nfft_adjoint_direct(ths->act_nfft_plan);
01060   else
01061     short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01062       else
01063   nfft_adjoint(ths->act_nfft_plan);
01064 
01065       if((r==rr)&&(J-rr!=rr))
01066   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01067 
01068       /* left */
01069       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*X(exp2i)(J);
01070 
01071       for (j=0; j<M; j++)
01072         ths->act_nfft_plan->f[j]= ths->f[j] *
01073                                   cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
01074 
01075       if(r<rr)
01076   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01077 
01078       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01079   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01080     nfft_adjoint_direct(ths->act_nfft_plan);
01081   else
01082     short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01083       else
01084   nfft_adjoint(ths->act_nfft_plan);
01085 
01086       if(r<rr)
01087   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01088 
01089       /* bottom */
01090       ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*X(exp2i)(J);
01091 
01092       for (j=0; j<M; j++)
01093         ths->act_nfft_plan->f[j]= ths->f[j] *
01094                                   cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
01095 
01096       if((r==rr)&&(J-rr!=rr))
01097   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01098 
01099       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01100   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01101     nfft_adjoint_direct(ths->act_nfft_plan);
01102   else
01103     short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01104       else
01105   nfft_adjoint(ths->act_nfft_plan);
01106 
01107       if((r==rr)&&(J-rr!=rr))
01108   RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
01109     } /* for(rr) */
01110 } /* void nsfft_adjoint_2d */
01111 
01112 static void nsfft_trafo_3d(nsfft_plan *ths)
01113 {
01114   int r,rr,j;
01115   double temp;
01116   int sum_N_B_less_r,N_B_r,a,b;
01117 
01118   int M=ths->M_total;
01119   int J=ths->J;
01120 
01121   /* center */
01122   ths->center_nfft_plan->f_hat=ths->f_hat+6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1);
01123 
01124   if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
01125     nfft_trafo_direct(ths->center_nfft_plan);
01126   else
01127     nfft_trafo(ths->center_nfft_plan);
01128 
01129   for (j=0; j<M; j++)
01130     ths->f[j] = ths->center_nfft_plan->f[j];
01131 
01132   sum_N_B_less_r=0;
01133   for(rr=0;rr<=(J+1)/2;rr++)
01134     {
01135       a=X(exp2i)(J-rr);
01136       b=X(exp2i)(rr);
01137 
01138       N_B_r=a*b*b;
01139 
01140       r=MIN(rr,J-rr);
01141       ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
01142 
01143       ths->act_nfft_plan->N[0]=X(exp2i)(r);
01144       if(a<b)
01145   ths->act_nfft_plan->N[1]=X(exp2i)(J-r);
01146       else
01147   ths->act_nfft_plan->N[1]=X(exp2i)(r);
01148       ths->act_nfft_plan->N[2]=X(exp2i)(J-r);
01149 
01150       /*printf("\n\n%d x %d x %d:\t",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1],ths->act_nfft_plan->N[2]); fflush(stdout);*/
01151 
01152       ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
01153       ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
01154       ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
01155       ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
01156       ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
01157 
01158       /* only for right - rear - top */
01159       if((J==0)||((J==1)&&(rr==1)))
01160   temp=-2.0*KPI;
01161       else
01162   temp=-3.0*KPI*X(exp2i)(J-rr);
01163 
01164       /* right */
01165       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
01166 
01167       if(a>b)
01168   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01169 
01170       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01171   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01172     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01173       nfft_trafo_direct(ths->act_nfft_plan);
01174     else
01175       short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01176   else
01177     short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01178       else
01179   nfft_trafo(ths->act_nfft_plan);
01180 
01181       if(a>b)
01182   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01183 
01184       for (j=0; j<M; j++)
01185         ths->f[j] += ths->act_nfft_plan->f[j] *
01186                      cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
01187 
01188       /* rear */
01189       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
01190 
01191       if(a>b)
01192   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01193       if(a<b)
01194   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01195 
01196       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01197   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01198     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01199       nfft_trafo_direct(ths->act_nfft_plan);
01200     else
01201       short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01202   else
01203     short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01204       else
01205   nfft_trafo(ths->act_nfft_plan);
01206 
01207       if(a>b)
01208   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01209       if(a<b)
01210   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01211 
01212       for (j=0; j<M; j++)
01213         ths->f[j] += ths->act_nfft_plan->f[j] *
01214                      cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
01215 
01216       /* top */
01217       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
01218 
01219       if(a<b)
01220   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01221 
01222       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01223   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01224     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01225       nfft_trafo_direct(ths->act_nfft_plan);
01226     else
01227       short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01228   else
01229     short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01230       else
01231   nfft_trafo(ths->act_nfft_plan);
01232 
01233       if(a<b)
01234   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01235 
01236       for (j=0; j<M; j++)
01237         ths->f[j] += ths->act_nfft_plan->f[j] *
01238                      cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
01239 
01240       /* only for left - front - bottom */
01241       if((J==0)||((J==1)&&(rr==1)))
01242   temp=-4.0*KPI;
01243       else
01244   temp=-3.0*KPI*X(exp2i)(J-rr);
01245 
01246       /* left */
01247       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
01248 
01249       if(a>b)
01250   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01251 
01252       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01253   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01254     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01255       nfft_trafo_direct(ths->act_nfft_plan);
01256     else
01257       short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01258   else
01259     short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01260       else
01261   nfft_trafo(ths->act_nfft_plan);
01262 
01263       if(a>b)
01264   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01265 
01266       for (j=0; j<M; j++)
01267         ths->f[j] += ths->act_nfft_plan->f[j] *
01268                      cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
01269 
01270       /* front */
01271       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
01272 
01273       if(a>b)
01274   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01275       if(a<b)
01276   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01277 
01278       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01279   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01280     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01281       nfft_trafo_direct(ths->act_nfft_plan);
01282     else
01283       short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01284   else
01285     short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01286       else
01287   nfft_trafo(ths->act_nfft_plan);
01288 
01289       if(a>b)
01290   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01291       if(a<b)
01292   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01293 
01294       for (j=0; j<M; j++)
01295         ths->f[j] += ths->act_nfft_plan->f[j] *
01296                      cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
01297 
01298       /* bottom */
01299       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
01300 
01301       if(a<b)
01302   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01303 
01304       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01305   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01306     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01307       nfft_trafo_direct(ths->act_nfft_plan);
01308     else
01309       short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01310   else
01311     short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01312       else
01313   nfft_trafo(ths->act_nfft_plan);
01314 
01315       if(a<b)
01316   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01317 
01318       for (j=0; j<M; j++)
01319         ths->f[j] += ths->act_nfft_plan->f[j] *
01320                      cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
01321 
01322       sum_N_B_less_r+=6*N_B_r;
01323     } /* for(rr) */
01324 } /* void nsfft_trafo_3d */
01325 
01326 static void nsfft_adjoint_3d(nsfft_plan *ths)
01327 {
01328   int r,rr,j;
01329   double temp;
01330   int sum_N_B_less_r,N_B_r,a,b;
01331 
01332   int M=ths->M_total;
01333   int J=ths->J;
01334 
01335   /* center */
01336   for (j=0; j<M; j++)
01337     ths->center_nfft_plan->f[j] = ths->f[j];
01338 
01339   ths->center_nfft_plan->f_hat=ths->f_hat+6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1);
01340 
01341   if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
01342     nfft_adjoint_direct(ths->center_nfft_plan);
01343   else
01344     nfft_adjoint(ths->center_nfft_plan);
01345 
01346   sum_N_B_less_r=0;
01347   for(rr=0;rr<=(J+1)/2;rr++)
01348     {
01349       a=X(exp2i)(J-rr);
01350       b=X(exp2i)(rr);
01351 
01352       N_B_r=a*b*b;
01353 
01354       r=MIN(rr,J-rr);
01355       ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
01356       ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[rr];
01357 
01358       ths->act_nfft_plan->N[0]=X(exp2i)(r);
01359       if(a<b)
01360   ths->act_nfft_plan->N[1]=X(exp2i)(J-r);
01361       else
01362   ths->act_nfft_plan->N[1]=X(exp2i)(r);
01363       ths->act_nfft_plan->N[2]=X(exp2i)(J-r);
01364 
01365       /*printf("\n\n%d x %d x %d:\t",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1],ths->act_nfft_plan->N[2]); fflush(stdout);*/
01366 
01367       ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
01368       ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
01369       ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
01370       ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
01371       ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
01372 
01373       /* only for right - rear - top */
01374       if((J==0)||((J==1)&&(rr==1)))
01375   temp=-2.0*KPI;
01376       else
01377   temp=-3.0*KPI*X(exp2i)(J-rr);
01378 
01379       /* right */
01380       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
01381 
01382       for (j=0; j<M; j++)
01383         ths->act_nfft_plan->f[j]= ths->f[j] *
01384                                   cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
01385 
01386       if(a>b)
01387   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01388 
01389       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01390   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01391     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01392       nfft_adjoint_direct(ths->act_nfft_plan);
01393     else
01394       short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01395   else
01396     short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01397       else
01398   nfft_adjoint(ths->act_nfft_plan);
01399 
01400       if(a>b)
01401   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01402 
01403       /* rear */
01404       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
01405 
01406       for (j=0; j<M; j++)
01407         ths->act_nfft_plan->f[j]= ths->f[j] *
01408                                   cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
01409 
01410       if(a>b)
01411   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01412       if(a<b)
01413   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01414 
01415       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01416   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01417     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01418       nfft_adjoint_direct(ths->act_nfft_plan);
01419     else
01420       short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01421   else
01422     short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01423       else
01424   nfft_adjoint(ths->act_nfft_plan);
01425 
01426       if(a>b)
01427   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01428       if(a<b)
01429   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01430 
01431       /* top */
01432       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
01433 
01434       for (j=0; j<M; j++)
01435         ths->act_nfft_plan->f[j]= ths->f[j] *
01436                                   cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
01437 
01438       if(a<b)
01439   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01440 
01441       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01442   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01443     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01444       nfft_adjoint_direct(ths->act_nfft_plan);
01445     else
01446       short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01447   else
01448     short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01449       else
01450   nfft_adjoint(ths->act_nfft_plan);
01451 
01452       if(a<b)
01453   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01454 
01455       /* only for left - front - bottom */
01456       if((J==0)||((J==1)&&(rr==1)))
01457   temp=-4.0*KPI;
01458       else
01459   temp=-3.0*KPI*X(exp2i)(J-rr);
01460 
01461       /* left */
01462       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
01463 
01464       for (j=0; j<M; j++)
01465         ths->act_nfft_plan->f[j]= ths->f[j] *
01466                                   cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
01467 
01468       if(a>b)
01469   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01470 
01471       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01472   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01473     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01474       nfft_adjoint_direct(ths->act_nfft_plan);
01475     else
01476       short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01477   else
01478     short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01479       else
01480   nfft_adjoint(ths->act_nfft_plan);
01481 
01482       if(a>b)
01483   RSWAP(ths->act_nfft_plan->x,ths->x_120);
01484 
01485       /* front */
01486       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
01487 
01488       for (j=0; j<M; j++)
01489         ths->act_nfft_plan->f[j]= ths->f[j] *
01490                                   cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
01491 
01492       if(a>b)
01493   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01494       if(a<b)
01495   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01496 
01497       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01498   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01499     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01500       nfft_adjoint_direct(ths->act_nfft_plan);
01501     else
01502       short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01503   else
01504     short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01505       else
01506   nfft_adjoint(ths->act_nfft_plan);
01507 
01508       if(a>b)
01509   RSWAP(ths->act_nfft_plan->x,ths->x_021);
01510       if(a<b)
01511   RSWAP(ths->act_nfft_plan->x,ths->x_102);
01512 
01513       /* bottom */
01514       ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
01515 
01516       for (j=0; j<M; j++)
01517         ths->act_nfft_plan->f[j]= ths->f[j] *
01518                                   cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
01519 
01520       if(a<b)
01521   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01522 
01523       if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01524   if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01525     if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01526       nfft_adjoint_direct(ths->act_nfft_plan);
01527     else
01528       short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01529   else
01530     short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01531       else
01532   nfft_adjoint(ths->act_nfft_plan);
01533 
01534       if(a<b)
01535   RSWAP(ths->act_nfft_plan->x,ths->x_201);
01536 
01537       sum_N_B_less_r+=6*N_B_r;
01538     } /* for(rr) */
01539 } /* void nsfft_adjoint_3d */
01540 
01541 void nsfft_trafo(nsfft_plan *ths)
01542 {
01543   if(ths->d==2)
01544     nsfft_trafo_2d(ths);
01545   else
01546     nsfft_trafo_3d(ths);
01547 }
01548 
01549 void nsfft_adjoint(nsfft_plan *ths)
01550 {
01551   if(ths->d==2)
01552     nsfft_adjoint_2d(ths);
01553   else
01554     nsfft_adjoint_3d(ths);
01555 }
01556 
01557 
01558 /*========================================================*/
01559 /* J >1, no precomputation at all!! */
01560 #ifdef GAUSSIAN
01561 static void nsfft_init_2d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
01562 {
01563   int r;
01564   int N[2];
01565   int n[2];
01566 
01567   ths->flags=snfft_flags;
01568   ths->sigma=2;
01569   ths->J=J;
01570   ths->M_total=M;
01571   ths->N_total=(J+4)*X(exp2i)(J+1);
01572 
01573   /* memory allocation */
01574   ths->f = (double _Complex *)nfft_malloc(M*sizeof(double _Complex));
01575   ths->f_hat = (double _Complex *)nfft_malloc(ths->N_total*sizeof(double _Complex));
01576   ths->x_transposed= (double*)nfft_malloc(2*M*sizeof(double));
01577 
01578   ths->act_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
01579   ths->center_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
01580 
01581   ths->set_fftw_plan1=(fftw_plan*) nfft_malloc((J/2+1)*sizeof(fftw_plan));
01582   ths->set_fftw_plan2=(fftw_plan*) nfft_malloc((J/2+1)*sizeof(fftw_plan));
01583 
01584   ths->set_nfft_plan_1d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
01585 
01586   /* planning the small nffts */
01587   /* r=0 */
01588   N[0]=1;            n[0]=ths->sigma*N[0];
01589   N[1]=X(exp2i)(J); n[1]=ths->sigma*N[1];
01590 
01591   nfft_init_guru(ths->act_nfft_plan,2,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01592 
01593   if(ths->act_nfft_plan->flags & PRE_ONE_PSI)
01594     nfft_precompute_one_psi(ths->act_nfft_plan);
01595 
01596   ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
01597   ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
01598 
01599   for(r=1;r<=J/2;r++)
01600     {
01601       N[0]=X(exp2i)(r);   n[0]=ths->sigma*N[0];
01602       N[1]=X(exp2i)(J-r); n[1]=ths->sigma*N[1];
01603       ths->set_fftw_plan1[r] =
01604   fftw_plan_dft(2, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
01605           FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
01606 
01607       ths->set_fftw_plan2[r] =
01608   fftw_plan_dft(2, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
01609           FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
01610     }
01611 
01612   /* planning the 1d nffts */
01613   for(r=0;r<=X(log2i)(m);r++)
01614     {
01615       N[0]=X(exp2i)(J-r); n[0]=ths->sigma*N[0]; /* ==N[1] of the 2 dimensional plan */
01616 
01617       nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01618       ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags | FG_PSI;
01619       ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
01620       ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
01621     }
01622 
01623   /* center plan */
01624   /* J/2 == floor(((double)J) / 2.0) */
01625   N[0]=X(exp2i)(J/2+1); n[0]=ths->sigma*N[0];
01626   N[1]=X(exp2i)(J/2+1); n[1]=ths->sigma*N[1];
01627   nfft_init_guru(ths->center_nfft_plan,2,N,M,n, m, MALLOC_F| FFTW_INIT,
01628                      FFTW_MEASURE);
01629   ths->center_nfft_plan->x= ths->act_nfft_plan->x;
01630   ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
01631                                       FG_PSI;
01632   ths->center_nfft_plan->K=ths->act_nfft_plan->K;
01633   ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
01634 
01635   if(ths->flags & NSDFT)
01636     {
01637       ths->index_sparse_to_full=(int*)nfft_malloc(ths->N_total*sizeof(int));
01638       init_index_sparse_to_full_2d(ths);
01639     }
01640 }
01641 #endif
01642 
01643 /*========================================================*/
01644 /* J >1, no precomputation at all!! */
01645 #ifdef GAUSSIAN
01646 static void nsfft_init_3d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
01647 {
01648   int r,rr,a,b;
01649   int N[3];
01650   int n[3];
01651 
01652   ths->flags=snfft_flags;
01653   ths->sigma=2;
01654   ths->J=J;
01655   ths->M_total=M;
01656   ths->N_total=6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+X(exp2i)(3*(J/2+1));
01657 
01658   /* memory allocation */
01659   ths->f =     (double _Complex *)nfft_malloc(M*sizeof(double _Complex));
01660   ths->f_hat = (double _Complex *)nfft_malloc(ths->N_total*sizeof(double _Complex));
01661 
01662   ths->x_102= (double*)nfft_malloc(3*M*sizeof(double));
01663   ths->x_201= (double*)nfft_malloc(3*M*sizeof(double));
01664   ths->x_120= (double*)nfft_malloc(3*M*sizeof(double));
01665   ths->x_021= (double*)nfft_malloc(3*M*sizeof(double));
01666 
01667   ths->act_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
01668   ths->center_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
01669 
01670   ths->set_fftw_plan1=(fftw_plan*) nfft_malloc(((J+1)/2+1)*sizeof(fftw_plan));
01671   ths->set_fftw_plan2=(fftw_plan*) nfft_malloc(((J+1)/2+1)*sizeof(fftw_plan));
01672 
01673   ths->set_nfft_plan_1d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
01674   ths->set_nfft_plan_2d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
01675 
01676   /* planning the small nffts */
01677   /* r=0 */
01678   N[0]=1;            n[0]=ths->sigma*N[0];
01679   N[1]=1;            n[1]=ths->sigma*N[1];
01680   N[2]=X(exp2i)(J); n[2]=ths->sigma*N[2];
01681 
01682   nfft_init_guru(ths->act_nfft_plan,3,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F, FFTW_MEASURE);
01683 
01684   if(ths->act_nfft_plan->flags & PRE_ONE_PSI)
01685     nfft_precompute_one_psi(ths->act_nfft_plan);
01686 
01687   /* malloc g1, g2 for maximal size */
01688   ths->act_nfft_plan->g1 = nfft_malloc(ths->sigma*ths->sigma*ths->sigma*X(exp2i)(J+(J+1)/2)*sizeof(double _Complex));
01689   ths->act_nfft_plan->g2 = nfft_malloc(ths->sigma*ths->sigma*ths->sigma*X(exp2i)(J+(J+1)/2)*sizeof(double _Complex));
01690 
01691   ths->act_nfft_plan->my_fftw_plan1 =
01692     fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
01693       FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
01694   ths->act_nfft_plan->my_fftw_plan2 =
01695     fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
01696       FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
01697 
01698   ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
01699   ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
01700 
01701   for(rr=1;rr<=(J+1)/2;rr++)
01702     {
01703       a=X(exp2i)(J-rr);
01704       b=X(exp2i)(rr);
01705 
01706       r=MIN(rr,J-rr);
01707 
01708       n[0]=ths->sigma*X(exp2i)(r);
01709       if(a<b)
01710   n[1]=ths->sigma*X(exp2i)(J-r);
01711       else
01712   n[1]=ths->sigma*X(exp2i)(r);
01713       n[2]=ths->sigma*X(exp2i)(J-r);
01714 
01715       ths->set_fftw_plan1[rr] =
01716   fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
01717           FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
01718       ths->set_fftw_plan2[rr] =
01719   fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
01720           FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
01721     }
01722 
01723   /* planning the 1d nffts */
01724   for(r=0;r<=X(log2i)(m);r++)
01725     {
01726       N[0]=X(exp2i)(J-r); n[0]=ths->sigma*N[0];
01727       N[1]=X(exp2i)(J-r); n[1]=ths->sigma*N[1];
01728 
01729       if(N[0]>m)
01730   {
01731     nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01732     ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags | FG_PSI;
01733     ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
01734     ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
01735     nfft_init_guru(&(ths->set_nfft_plan_2d[r]),2,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01736     ths->set_nfft_plan_2d[r].flags = ths->set_nfft_plan_2d[r].flags | FG_PSI;
01737     ths->set_nfft_plan_2d[r].K=ths->act_nfft_plan->K;
01738     ths->set_nfft_plan_2d[r].psi=ths->act_nfft_plan->psi;
01739   }
01740     }
01741 
01742   /* center plan */
01743   /* J/2 == floor(((double)J) / 2.0) */
01744   N[0]=X(exp2i)(J/2+1); n[0]=ths->sigma*N[0];
01745   N[1]=X(exp2i)(J/2+1); n[1]=ths->sigma*N[1];
01746   N[2]=X(exp2i)(J/2+1); n[2]=ths->sigma*N[2];
01747   nfft_init_guru(ths->center_nfft_plan,3,N,M,n, m, MALLOC_F| FFTW_INIT,
01748                      FFTW_MEASURE);
01749   ths->center_nfft_plan->x= ths->act_nfft_plan->x;
01750   ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
01751                                       FG_PSI;
01752   ths->center_nfft_plan->K=ths->act_nfft_plan->K;
01753   ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
01754 
01755   if(ths->flags & NSDFT)
01756     {
01757       ths->index_sparse_to_full=(int*)nfft_malloc(ths->N_total*sizeof(int));
01758       init_index_sparse_to_full_3d(ths);
01759     }
01760 }
01761 #endif
01762 
01763 #ifdef GAUSSIAN
01764 void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
01765 {
01766   ths->d=d;
01767 
01768   if(ths->d==2)
01769     nsfft_init_2d(ths, J, M, m, flags);
01770   else
01771     nsfft_init_3d(ths, J, M, m, flags);
01772 
01773 
01774   ths->mv_trafo = (void (*) (void* ))nsfft_trafo;
01775   ths->mv_adjoint = (void (*) (void* ))nsfft_adjoint;
01776 }
01777 #else
01778 void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
01779 {
01780   UNUSED(ths);
01781   UNUSED(d);
01782   UNUSED(J);
01783   UNUSED(M);
01784   UNUSED(m);
01785   UNUSED(flags);
01786   fprintf(stderr,
01787     "\nError in kernel/nsfft_init: require GAUSSIAN window function\n");
01788 }
01789 #endif
01790 
01791 static void nsfft_finalize_2d(nsfft_plan *ths)
01792 {
01793   int r;
01794 
01795   if(ths->flags & NSDFT)
01796     nfft_free(ths->index_sparse_to_full);
01797 
01798   /* center plan */
01799   ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^ FG_PSI;
01800   nfft_finalize(ths->center_nfft_plan);
01801 
01802   /* the 1d nffts */
01803   for(r=0;r<=X(log2i)(ths->act_nfft_plan->m);r++)
01804     {
01805       ths->set_nfft_plan_1d[r].flags =
01806         ths->set_nfft_plan_1d[r].flags ^ FG_PSI;
01807       nfft_finalize(&(ths->set_nfft_plan_1d[r]));
01808     }
01809 
01810   /* finalize the small nffts */
01811   ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
01812   ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
01813 
01814   for(r=1;r<=ths->J/2;r++)
01815     {
01816       fftw_destroy_plan(ths->set_fftw_plan2[r]);
01817       fftw_destroy_plan(ths->set_fftw_plan1[r]);
01818     }
01819 
01820   /* r=0 */
01821   nfft_finalize(ths->act_nfft_plan);
01822 
01823   nfft_free(ths->set_nfft_plan_1d);
01824 
01825   nfft_free(ths->set_fftw_plan2);
01826   nfft_free(ths->set_fftw_plan1);
01827 
01828   nfft_free(ths->x_transposed);
01829 
01830   nfft_free(ths->f_hat);
01831   nfft_free(ths->f);
01832 }
01833 
01834 static void nsfft_finalize_3d(nsfft_plan *ths)
01835 {
01836   int r;
01837 
01838   if(ths->flags & NSDFT)
01839     nfft_free(ths->index_sparse_to_full);
01840 
01841   /* center plan */
01842   ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^ FG_PSI;
01843   nfft_finalize(ths->center_nfft_plan);
01844 
01845   /* the 1d and 2d nffts */
01846   for(r=0;r<=X(log2i)(ths->act_nfft_plan->m);r++)
01847     {
01848       if(X(exp2i)(ths->J-r)>ths->act_nfft_plan->m)
01849   {
01850     ths->set_nfft_plan_2d[r].flags = ths->set_nfft_plan_2d[r].flags ^ FG_PSI;
01851     nfft_finalize(&(ths->set_nfft_plan_2d[r]));
01852     ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags ^ FG_PSI;
01853     nfft_finalize(&(ths->set_nfft_plan_1d[r]));
01854   }
01855     }
01856 
01857   /* finalize the small nffts */
01858   ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
01859   ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
01860 
01861   for(r=1;r<=(ths->J+1)/2;r++)
01862     {
01863       fftw_destroy_plan(ths->set_fftw_plan2[r]);
01864       fftw_destroy_plan(ths->set_fftw_plan1[r]);
01865     }
01866 
01867   /* r=0 */
01868   nfft_finalize(ths->act_nfft_plan);
01869 
01870   nfft_free(ths->set_nfft_plan_1d);
01871   nfft_free(ths->set_nfft_plan_2d);
01872 
01873   nfft_free(ths->set_fftw_plan2);
01874   nfft_free(ths->set_fftw_plan1);
01875 
01876   nfft_free(ths->x_102);
01877   nfft_free(ths->x_201);
01878   nfft_free(ths->x_120);
01879   nfft_free(ths->x_021);
01880 
01881   nfft_free(ths->f_hat);
01882   nfft_free(ths->f);
01883 }
01884 
01885 void nsfft_finalize(nsfft_plan *ths)
01886 {
01887   if(ths->d==2)
01888     nsfft_finalize_2d(ths);
01889   else
01890     nsfft_finalize_3d(ths);
01891 }