NFFT Logo 3.2.2
quadratureS2.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 
00019 /* $Id: quadratureS2.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00026 #include "config.h"
00027 
00028 /* Include standard C headers. */
00029 #include <math.h>
00030 #include <stdlib.h>
00031 #include <stdio.h>
00032 #include <string.h>
00033 #include <time.h>
00034 #ifdef HAVE_COMPLEX_H
00035 #include <complex.h>
00036 #endif
00037 
00038 /* Include NFFT 3 utilities headers. */
00039 #include "nfft3util.h"
00040 
00041 /* Include NFFT 3 library header. */
00042 #include "nfft3.h"
00043 
00044 #include "infft.h"
00045 
00047 enum boolean {NO = 0, YES = 1};
00048 
00050 enum testtype {ERROR = 0, TIMING = 1};
00051 
00053 enum gridtype {GRID_GAUSS_LEGENDRE = 0, GRID_CLENSHAW_CURTIS = 1,
00054   GRID_HEALPIX = 2, GRID_EQUIDISTRIBUTION = 3, GRID_EQUIDISTRIBUTION_UNIFORM = 4};
00055 
00057 enum functiontype {FUNCTION_RANDOM_BANDLIMITED = 0, FUNCTION_F1 = 1,
00058   FUNCTION_F2 = 2, FUNCTION_F3 = 3, FUNCTION_F4 = 4, FUNCTION_F5 = 5,
00059   FUNCTION_F6 = 6};
00060 
00062 enum modes {USE_GRID = 0, RANDOM = 1};
00063 
00072 int main (int argc, char **argv)
00073 {
00074   int tc;                      
00075   int tc_max;                  
00077   int *NQ;                     
00079   int NQ_max;                  
00081   int *SQ;                     
00083   int SQ_max;                  
00084   int *RQ;                     
00086   int iNQ;                     
00087   int iNQ_max;                 
00088   int testfunction;            
00089   int N;                       
00091   int use_nfsft;               
00092   int use_nfft;                
00093   int use_fpt;                 
00094   int cutoff;                  
00095   double threshold;            
00097   int gridtype;                
00098   int repetitions;             
00099   int mode;                    
00101   double *w;                   
00102   double *x_grid;              
00103   double *x_compare;           
00104   double _Complex *f_grid;             
00105   double _Complex *f_compare;          
00106   double _Complex *f;                  
00107   double _Complex *f_hat_gen;         
00109   double _Complex *f_hat;              
00111   nfsft_plan plan_adjoint;     
00112   nfsft_plan plan;             
00113   nfsft_plan plan_gen;         
00115   double t_avg;                
00116   double err_infty_avg;        
00117   double err_2_avg;            
00119   int i;                       
00120   int k;                       
00121   int n;                       
00122   int d;                       
00124   int m_theta;                 
00126   int m_phi;                   
00128   int m_total;                 
00129   double *theta;               
00131   double *phi;                 
00133   fftw_plan fplan;             
00135   //int nside;                   /**< The size parameter for the HEALPix grid   */
00136   int d2;
00137   int M;
00138   double theta_s;
00139   double x1,x2,x3,temp;
00140   int m_compare;
00141   nfsft_plan *plan_adjoint_ptr;
00142   nfsft_plan *plan_ptr;
00143   double *w_temp;
00144   int testmode;
00145   ticks t0, t1;
00146 
00147   /* Read the number of testcases. */
00148   fscanf(stdin,"testcases=%d\n",&tc_max);
00149   fprintf(stdout,"%d\n",tc_max);
00150 
00151   /* Process each testcase. */
00152   for (tc = 0; tc < tc_max; tc++)
00153   {
00154     /* Check if the fast transform shall be used. */
00155     fscanf(stdin,"nfsft=%d\n",&use_nfsft);
00156     fprintf(stdout,"%d\n",use_nfsft);
00157     if (use_nfsft != NO)
00158     {
00159       /* Check if the NFFT shall be used. */
00160       fscanf(stdin,"nfft=%d\n",&use_nfft);
00161       fprintf(stdout,"%d\n",use_nfsft);
00162       if (use_nfft != NO)
00163       {
00164         /* Read the cut-off parameter. */
00165         fscanf(stdin,"cutoff=%d\n",&cutoff);
00166         fprintf(stdout,"%d\n",cutoff);
00167       }
00168       else
00169       {
00170         /* TODO remove this */
00171         /* Initialize unused variable with dummy value. */
00172         cutoff = 1;
00173       }
00174       /* Check if the fast polynomial transform shall be used. */
00175       fscanf(stdin,"fpt=%d\n",&use_fpt);
00176       fprintf(stdout,"%d\n",use_fpt);
00177       if (use_fpt != NO)
00178       {
00179         /* Read the NFSFT threshold parameter. */
00180         fscanf(stdin,"threshold=%lf\n",&threshold);
00181         fprintf(stdout,"%lf\n",threshold);
00182       }
00183       else
00184       {
00185         /* TODO remove this */
00186         /* Initialize unused variable with dummy value. */
00187         threshold = 1000.0;
00188       }
00189     }
00190     else
00191     {
00192       /* TODO remove this */
00193       /* Set dummy values. */
00194       use_nfft = NO;
00195       use_fpt = NO;
00196       cutoff = 3;
00197       threshold = 1000.0;
00198     }
00199 
00200     /* Read the testmode type. */
00201     fscanf(stdin,"testmode=%d\n",&testmode);
00202     fprintf(stdout,"%d\n",testmode);
00203 
00204     if (testmode == ERROR)
00205     {
00206       /* Read the quadrature grid type. */
00207       fscanf(stdin,"gridtype=%d\n",&gridtype);
00208       fprintf(stdout,"%d\n",gridtype);
00209 
00210       /* Read the test function. */
00211       fscanf(stdin,"testfunction=%d\n",&testfunction);
00212       fprintf(stdout,"%d\n",testfunction);
00213 
00214       /* Check if random bandlimited function has been chosen. */
00215       if (testfunction == FUNCTION_RANDOM_BANDLIMITED)
00216       {
00217         /* Read the bandwidht. */
00218         fscanf(stdin,"bandlimit=%d\n",&N);
00219         fprintf(stdout,"%d\n",N);
00220       }
00221       else
00222       {
00223         N = 1;
00224       }
00225 
00226       /* Read the number of repetitions. */
00227       fscanf(stdin,"repetitions=%d\n",&repetitions);
00228       fprintf(stdout,"%d\n",repetitions);
00229 
00230       fscanf(stdin,"mode=%d\n",&mode);
00231       fprintf(stdout,"%d\n",mode);
00232 
00233       if (mode == RANDOM)
00234       {
00235         /* Read the bandwidht. */
00236         fscanf(stdin,"points=%d\n",&m_compare);
00237         fprintf(stdout,"%d\n",m_compare);
00238         x_compare = (double*) nfft_malloc(2*m_compare*sizeof(double));
00239         d = 0;
00240         while (d < m_compare)
00241         {
00242           x1 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
00243           x2 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
00244           x3 = 2.0*(((double)rand())/RAND_MAX) - 1.0;
00245           temp = sqrt(x1*x1+x2*x2+x3*x3);
00246           if (temp <= 1)
00247           {
00248             x_compare[2*d+1] = acos(x3);
00249             if (x_compare[2*d+1] == 0 || x_compare[2*d+1] == PI)
00250             {
00251               x_compare[2*d] = 0.0;
00252             }
00253             else
00254             {
00255               x_compare[2*d] = atan2(x2/sin(x_compare[2*d+1]),x1/sin(x_compare[2*d+1]));
00256             }
00257             x_compare[2*d] *= 1.0/(2.0*PI);
00258             x_compare[2*d+1] *= 1.0/(2.0*PI);
00259             d++;
00260           }
00261         }
00262         f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
00263         f = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
00264       }
00265     }
00266 
00267     /* Initialize maximum cut-off degree and grid size parameter. */
00268     NQ_max = 0;
00269     SQ_max = 0;
00270 
00271     /* Read the number of cut-off degrees. */
00272     fscanf(stdin,"bandwidths=%d\n",&iNQ_max);
00273     fprintf(stdout,"%d\n",iNQ_max);
00274 
00275     /* Allocate memory for the cut-off degrees and grid size parameters. */
00276     NQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
00277     SQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
00278     if (testmode == TIMING)
00279     {
00280       RQ = (int*) nfft_malloc(iNQ_max*sizeof(int));
00281     }
00282 
00283     /* Read the cut-off degrees and grid size parameters. */
00284     for (iNQ = 0; iNQ < iNQ_max; iNQ++)
00285     {
00286       if (testmode == TIMING)
00287       {
00288         /* Read cut-off degree and grid size parameter. */
00289         fscanf(stdin,"%d %d %d\n",&NQ[iNQ],&SQ[iNQ],&RQ[iNQ]);
00290         fprintf(stdout,"%d %d %d\n",NQ[iNQ],SQ[iNQ],RQ[iNQ]);
00291         NQ_max = NFFT_MAX(NQ_max,NQ[iNQ]);
00292         SQ_max = NFFT_MAX(SQ_max,SQ[iNQ]);
00293       }
00294       else
00295       {
00296         /* Read cut-off degree and grid size parameter. */
00297         fscanf(stdin,"%d %d\n",&NQ[iNQ],&SQ[iNQ]);
00298         fprintf(stdout,"%d %d\n",NQ[iNQ],SQ[iNQ]);
00299         NQ_max = NFFT_MAX(NQ_max,NQ[iNQ]);
00300         SQ_max = NFFT_MAX(SQ_max,SQ[iNQ]);
00301       }
00302     }
00303 
00304     /* Do precomputation. */
00305     //fprintf(stderr,"NFSFT Precomputation\n");
00306     //fflush(stderr);
00307     nfsft_precompute(NQ_max, threshold,
00308       ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U);
00309 
00310     if (testmode == TIMING)
00311     {
00312       /* Allocate data structures. */
00313       f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ_max)*sizeof(double _Complex));
00314       f = (double _Complex*) nfft_malloc(SQ_max*sizeof(double _Complex));
00315       x_grid = (double*) nfft_malloc(2*SQ_max*sizeof(double));
00316       for (d = 0; d < SQ_max; d++)
00317       {
00318         f[d] = (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5);
00319         x_grid[2*d] = (((double)rand())/RAND_MAX) - 0.5;
00320         x_grid[2*d+1] = (((double)rand())/RAND_MAX) * 0.5;
00321       }
00322     }
00323 
00324     //fprintf(stderr,"Entering loop\n");
00325     //fflush(stderr);
00326     /* Process all cut-off bandwidths. */
00327     for (iNQ = 0; iNQ < iNQ_max; iNQ++)
00328     {
00329       if (testmode == TIMING)
00330       {
00331         nfsft_init_guru(&plan,NQ[iNQ],SQ[iNQ], NFSFT_NORMALIZED |
00332           ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00333           ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00334           PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFTW_MEASURE | FFT_OUT_OF_PLACE,
00335           cutoff);
00336 
00337         plan.f_hat = f_hat;
00338         plan.x = x_grid;
00339         plan.f = f;
00340 
00341         nfsft_precompute_x(&plan);
00342 
00343         t_avg = 0.0;
00344 
00345         for (i = 0; i < RQ[iNQ]; i++)
00346         {
00347           t0 = getticks();
00348 
00349           if (use_nfsft != NO)
00350           {
00351             /* Execute the adjoint NFSFT transformation. */
00352             nfsft_adjoint(&plan);
00353           }
00354           else
00355           {
00356             /* Execute the adjoint direct NDSFT transformation. */
00357             nfsft_adjoint_direct(&plan);
00358           }
00359 
00360           t1 = getticks();
00361           t_avg += nfft_elapsed_seconds(t1,t0);
00362         }
00363 
00364         t_avg = t_avg/((double)RQ[iNQ]);
00365 
00366         nfsft_finalize(&plan);
00367 
00368         fprintf(stdout,"%+le\n", t_avg);
00369         fprintf(stderr,"%d: %4d %4d %+le\n", tc, NQ[iNQ], SQ[iNQ], t_avg);
00370       }
00371       else
00372       {
00373         /* Determine the maximum number of nodes. */
00374         switch (gridtype)
00375         {
00376           case GRID_GAUSS_LEGENDRE:
00377             /* Calculate grid dimensions. */
00378             m_theta = SQ[iNQ] + 1;
00379             m_phi = 2*SQ[iNQ] + 2;
00380             m_total = m_theta*m_phi;
00381             break;
00382           case GRID_CLENSHAW_CURTIS:
00383             /* Calculate grid dimensions. */
00384             m_theta = 2*SQ[iNQ] + 1;
00385             m_phi = 2*SQ[iNQ] + 2;
00386             m_total = m_theta*m_phi;
00387             break;
00388           case GRID_HEALPIX:
00389             m_theta = 1;
00390             m_phi = 12*SQ[iNQ]*SQ[iNQ];
00391             m_total = m_theta * m_phi;
00392             //fprintf("HEALPix: SQ = %d, m_theta = %d, m_phi= %d, m");
00393             break;
00394           case GRID_EQUIDISTRIBUTION:
00395           case GRID_EQUIDISTRIBUTION_UNIFORM:
00396             m_theta = 2;
00397             //fprintf(stderr,"ed: m_theta = %d\n",m_theta);
00398             for (k = 1; k < SQ[iNQ]; k++)
00399             {
00400               m_theta += (int)floor((2*PI)/acos((cos(PI/(double)SQ[iNQ])-
00401                 cos(k*PI/(double)SQ[iNQ])*cos(k*PI/(double)SQ[iNQ]))/
00402                 (sin(k*PI/(double)SQ[iNQ])*sin(k*PI/(double)SQ[iNQ]))));
00403               //fprintf(stderr,"ed: m_theta = %d\n",m_theta);
00404             }
00405             //fprintf(stderr,"ed: m_theta final = %d\n",m_theta);
00406             m_phi = 1;
00407             m_total = m_theta * m_phi;
00408             break;
00409         }
00410 
00411         /* Allocate memory for data structures. */
00412         w = (double*) nfft_malloc(m_theta*sizeof(double));
00413         x_grid = (double*) nfft_malloc(2*m_total*sizeof(double));
00414 
00415         //fprintf(stderr,"NQ = %d\n",NQ[iNQ]);
00416         //fflush(stderr);
00417         switch (gridtype)
00418         {
00419           case GRID_GAUSS_LEGENDRE:
00420             //fprintf(stderr,"Generating grid for NQ = %d, SQ = %d\n",NQ[iNQ],SQ[iNQ]);
00421             //fflush(stderr);
00422 
00423             /* Read quadrature weights. */
00424             for (k = 0; k < m_theta; k++)
00425             {
00426               fscanf(stdin,"%le\n",&w[k]);
00427               w[k] *= (2.0*PI)/((double)m_phi);
00428             }
00429 
00430             //fprintf(stderr,"Allocating theta and phi\n");
00431             //fflush(stderr);
00432             /* Allocate memory to store the grid's angles. */
00433             theta = (double*) nfft_malloc(m_theta*sizeof(double));
00434             phi = (double*) nfft_malloc(m_phi*sizeof(double));
00435 
00436             //if (theta == NULL || phi == NULL)
00437             //{
00438               //fprintf(stderr,"Couldn't allocate theta and phi\n");
00439               //fflush(stderr);
00440             //}
00441 
00442 
00443             /* Read angles theta. */
00444             for (k = 0; k < m_theta; k++)
00445             {
00446               fscanf(stdin,"%le\n",&theta[k]);
00447             }
00448 
00449             /* Generate the grid angles phi. */
00450             for (n = 0; n < m_phi; n++)
00451             {
00452               phi[n] = n/((double)m_phi);
00453               phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0));
00454             }
00455 
00456             //fprintf(stderr,"Generating grid nodes\n");
00457             //fflush(stderr);
00458 
00459             /* Generate the grid's nodes. */
00460             d = 0;
00461             for (k = 0; k < m_theta; k++)
00462             {
00463               for (n = 0; n < m_phi; n++)
00464               {
00465                 x_grid[2*d] = phi[n];
00466                 x_grid[2*d+1] = theta[k];
00467                 d++;
00468               }
00469             }
00470 
00471             //fprintf(stderr,"Freeing theta and phi\n");
00472             //fflush(stderr);
00473             /* Free the arrays for the grid's angles. */
00474             nfft_free(theta);
00475             nfft_free(phi);
00476 
00477             break;
00478 
00479           case GRID_CLENSHAW_CURTIS:
00480 
00481             /* Allocate memory to store the grid's angles. */
00482             theta = (double*) nfft_malloc(m_theta*sizeof(double));
00483             phi = (double*) nfft_malloc(m_phi*sizeof(double));
00484 
00485             /* Generate the grid angles theta. */
00486             for (k = 0; k < m_theta; k++)
00487             {
00488               theta[k] = k/((double)2*(m_theta-1));
00489             }
00490 
00491             /* Generate the grid angles phi. */
00492             for (n = 0; n < m_phi; n++)
00493             {
00494               phi[n] = n/((double)m_phi);
00495               phi[n] -= ((phi[n]>=0.5)?(1.0):(0.0));
00496             }
00497 
00498             /* Generate quadrature weights. */
00499             fplan = fftw_plan_r2r_1d(SQ[iNQ]+1, w, w, FFTW_REDFT00, 0U);
00500             for (k = 0; k < SQ[iNQ]+1; k++)
00501             {
00502               w[k] = -2.0/(4*k*k-1);
00503             }
00504             fftw_execute(fplan);
00505             w[0] *= 0.5;
00506 
00507             for (k = 0; k < SQ[iNQ]+1; k++)
00508             {
00509               w[k] *= (2.0*PI)/((double)(m_theta-1)*m_phi);
00510               w[m_theta-1-k] = w[k];
00511             }
00512             fftw_destroy_plan(fplan);
00513 
00514             /* Generate the grid's nodes. */
00515             d = 0;
00516             for (k = 0; k < m_theta; k++)
00517             {
00518               for (n = 0; n < m_phi; n++)
00519               {
00520                 x_grid[2*d] = phi[n];
00521                 x_grid[2*d+1] = theta[k];
00522                 d++;
00523               }
00524             }
00525 
00526             /* Free the arrays for the grid's angles. */
00527             nfft_free(theta);
00528             nfft_free(phi);
00529 
00530             break;
00531 
00532           case GRID_HEALPIX:
00533 
00534             d = 0;
00535             for (k = 1; k <= SQ[iNQ]-1; k++)
00536             {
00537               for (n = 0; n <= 4*k-1; n++)
00538               {
00539                 x_grid[2*d+1] = 1 - (k*k)/((double)(3.0*SQ[iNQ]*SQ[iNQ]));
00540                 x_grid[2*d] =  ((n+0.5)/(4*k));
00541                 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
00542                 d++;
00543               }
00544             }
00545 
00546             d2 = d-1;
00547 
00548             for (k = SQ[iNQ]; k <= 3*SQ[iNQ]; k++)
00549             {
00550               for (n = 0; n <= 4*SQ[iNQ]-1; n++)
00551               {
00552                 x_grid[2*d+1] = 2.0/(3*SQ[iNQ])*(2*SQ[iNQ]-k);
00553                 x_grid[2*d] = (n+((k%2==0)?(0.5):(0.0)))/(4*SQ[iNQ]);
00554                 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
00555                 d++;
00556               }
00557             }
00558 
00559             for (k = 1; k <= SQ[iNQ]-1; k++)
00560             {
00561               for (n = 0; n <= 4*k-1; n++)
00562               {
00563                 x_grid[2*d+1] = -x_grid[2*d2+1];
00564                 x_grid[2*d] =  x_grid[2*d2];
00565                 d++;
00566                 d2--;
00567               }
00568             }
00569 
00570             for (d = 0; d < m_total; d++)
00571             {
00572               x_grid[2*d+1] = acos(x_grid[2*d+1])/(2.0*PI);
00573             }
00574 
00575             w[0] = (4.0*PI)/(m_total);
00576             break;
00577 
00578           case GRID_EQUIDISTRIBUTION:
00579           case GRID_EQUIDISTRIBUTION_UNIFORM:
00580             /* TODO Compute the weights. */
00581 
00582             if (gridtype == GRID_EQUIDISTRIBUTION)
00583             {
00584               w_temp = (double*) nfft_malloc((SQ[iNQ]+1)*sizeof(double));
00585               fplan = fftw_plan_r2r_1d(SQ[iNQ]/2+1, w_temp, w_temp, FFTW_REDFT00, 0U);
00586               for (k = 0; k < SQ[iNQ]/2+1; k++)
00587               {
00588                 w_temp[k] = -2.0/(4*k*k-1);
00589               }
00590               fftw_execute(fplan);
00591               w_temp[0] *= 0.5;
00592 
00593               for (k = 0; k < SQ[iNQ]/2+1; k++)
00594               {
00595                 w_temp[k] *= (2.0*PI)/((double)(SQ[iNQ]));
00596                 w_temp[SQ[iNQ]-k] = w_temp[k];
00597               }
00598               fftw_destroy_plan(fplan);
00599             }
00600 
00601             d = 0;
00602             x_grid[2*d] = -0.5;
00603             x_grid[2*d+1] = 0.0;
00604             if (gridtype == GRID_EQUIDISTRIBUTION)
00605             {
00606               w[d] = w_temp[0];
00607             }
00608             else
00609             {
00610               w[d] = (4.0*PI)/(m_total);
00611             }
00612             d = 1;
00613             x_grid[2*d] = -0.5;
00614             x_grid[2*d+1] = 0.5;
00615             if (gridtype == GRID_EQUIDISTRIBUTION)
00616             {
00617               w[d] = w_temp[SQ[iNQ]];
00618             }
00619             else
00620             {
00621               w[d] = (4.0*PI)/(m_total);
00622             }
00623             d = 2;
00624 
00625             for (k = 1; k < SQ[iNQ]; k++)
00626             {
00627               theta_s = (double)k*PI/(double)SQ[iNQ];
00628               M = (int)floor((2.0*PI)/acos((cos(PI/(double)SQ[iNQ])-
00629                 cos(theta_s)*cos(theta_s))/(sin(theta_s)*sin(theta_s))));
00630 
00631               for (n = 0; n < M; n++)
00632               {
00633                 x_grid[2*d] = (n + 0.5)/M;
00634                 x_grid[2*d] -= (x_grid[2*d]>=0.5)?(1.0):(0.0);
00635                 x_grid[2*d+1] = theta_s/(2.0*PI);
00636                 if (gridtype == GRID_EQUIDISTRIBUTION)
00637                 {
00638                   w[d] = w_temp[k]/((double)(M));
00639                 }
00640                 else
00641                 {
00642                   w[d] = (4.0*PI)/(m_total);
00643                 }
00644                 d++;
00645               }
00646             }
00647 
00648             if (gridtype == GRID_EQUIDISTRIBUTION)
00649             {
00650               nfft_free(w_temp);
00651             }
00652             break;
00653 
00654           default:
00655             break;
00656         }
00657 
00658         /* Allocate memory for grid values. */
00659         f_grid = (double _Complex*) nfft_malloc(m_total*sizeof(double _Complex));
00660 
00661         if (mode == RANDOM)
00662         {
00663         }
00664         else
00665         {
00666           m_compare = m_total;
00667           f_compare = (double _Complex*) nfft_malloc(m_compare*sizeof(double _Complex));
00668           x_compare = x_grid;
00669           f = f_grid;
00670         }
00671 
00672         //fprintf(stderr,"Generating test function\n");
00673         //fflush(stderr);
00674         switch (testfunction)
00675         {
00676           case FUNCTION_RANDOM_BANDLIMITED:
00677             f_hat_gen = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(N)*sizeof(double _Complex));
00678             //fprintf(stderr,"Generating random test function\n");
00679             //fflush(stderr);
00680             /* Generate random function samples by sampling a bandlimited
00681              * function. */
00682             nfsft_init_guru(&plan_gen,N,m_total, NFSFT_NORMALIZED |
00683               ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00684               ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00685               ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00686               FFT_OUT_OF_PLACE, cutoff);
00687 
00688             plan_gen.f_hat = f_hat_gen;
00689             plan_gen.x = x_grid;
00690             plan_gen.f = f_grid;
00691 
00692             nfsft_precompute_x(&plan_gen);
00693 
00694             for (k = 0; k < plan_gen.N_total; k++)
00695             {
00696               f_hat_gen[k] = 0.0;
00697             }
00698 
00699             for (k = 0; k <= N; k++)
00700             {
00701               for (n = -k; n <= k; n++)
00702               {
00703                 f_hat_gen[NFSFT_INDEX(k,n,&plan_gen)] =
00704                 (((double)rand())/RAND_MAX)-0.5 + _Complex_I*((((double)rand())/RAND_MAX)-0.5);
00705               }
00706             }
00707 
00708             if (use_nfsft != NO)
00709             {
00710               /* Execute the NFSFT transformation. */
00711               nfsft_trafo(&plan_gen);
00712             }
00713             else
00714             {
00715               /* Execute the direct NDSFT transformation. */
00716               nfsft_trafo_direct(&plan_gen);
00717             }
00718 
00719             nfsft_finalize(&plan_gen);
00720 
00721             if (mode == RANDOM)
00722             {
00723               nfsft_init_guru(&plan_gen,N,m_compare, NFSFT_NORMALIZED |
00724                 ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00725                 ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00726                 ((N>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00727                 FFT_OUT_OF_PLACE, cutoff);
00728 
00729               plan_gen.f_hat = f_hat_gen;
00730               plan_gen.x = x_compare;
00731               plan_gen.f = f_compare;
00732 
00733               nfsft_precompute_x(&plan_gen);
00734 
00735               if (use_nfsft != NO)
00736               {
00737                 /* Execute the NFSFT transformation. */
00738                 nfsft_trafo(&plan_gen);
00739               }
00740               else
00741               {
00742                 /* Execute the direct NDSFT transformation. */
00743                 nfsft_trafo_direct(&plan_gen);
00744               }
00745 
00746               nfsft_finalize(&plan_gen);
00747             }
00748             else
00749             {
00750               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00751             }
00752 
00753             nfft_free(f_hat_gen);
00754 
00755             break;
00756 
00757           case FUNCTION_F1:
00758             for (d = 0; d < m_total; d++)
00759             {
00760               x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00761               x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00762               x3 = cos(x_grid[2*d+1]*2.0*PI);
00763               f_grid[d] = x1*x2*x3;
00764             }
00765             if (mode == RANDOM)
00766             {
00767               for (d = 0; d < m_compare; d++)
00768               {
00769                 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00770                 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00771                 x3 = cos(x_compare[2*d+1]*2.0*PI);
00772                 f_compare[d] = x1*x2*x3;
00773               }
00774             }
00775             else
00776             {
00777               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00778             }
00779             break;
00780           case FUNCTION_F2:
00781             for (d = 0; d < m_total; d++)
00782             {
00783               x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00784               x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00785               x3 = cos(x_grid[2*d+1]*2.0*PI);
00786               f_grid[d] = 0.1*exp(x1+x2+x3);
00787             }
00788             if (mode == RANDOM)
00789             {
00790               for (d = 0; d < m_compare; d++)
00791               {
00792                 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00793                 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00794                 x3 = cos(x_compare[2*d+1]*2.0*PI);
00795                 f_compare[d] = 0.1*exp(x1+x2+x3);
00796               }
00797             }
00798             else
00799             {
00800               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00801             }
00802             break;
00803           case FUNCTION_F3:
00804             for (d = 0; d < m_total; d++)
00805             {
00806               x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00807               x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00808               x3 = cos(x_grid[2*d+1]*2.0*PI);
00809               temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00810               f_grid[d] = 0.1*temp;
00811             }
00812             if (mode == RANDOM)
00813             {
00814               for (d = 0; d < m_compare; d++)
00815               {
00816                 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00817                 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00818                 x3 = cos(x_compare[2*d+1]*2.0*PI);
00819                 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00820                 f_compare[d] = 0.1*temp;
00821               }
00822             }
00823             else
00824             {
00825               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00826             }
00827             break;
00828           case FUNCTION_F4:
00829             for (d = 0; d < m_total; d++)
00830             {
00831               x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00832               x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00833               x3 = cos(x_grid[2*d+1]*2.0*PI);
00834               temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00835               f_grid[d] = 1.0/(temp);
00836             }
00837             if (mode == RANDOM)
00838             {
00839               for (d = 0; d < m_compare; d++)
00840               {
00841                 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00842                 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00843                 x3 = cos(x_compare[2*d+1]*2.0*PI);
00844                 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00845                 f_compare[d] = 1.0/(temp);
00846               }
00847             }
00848             else
00849             {
00850               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00851             }
00852             break;
00853           case FUNCTION_F5:
00854             for (d = 0; d < m_total; d++)
00855             {
00856               x1 = sin(x_grid[2*d+1]*2.0*PI)*cos(x_grid[2*d]*2.0*PI);
00857               x2 = sin(x_grid[2*d+1]*2.0*PI)*sin(x_grid[2*d]*2.0*PI);
00858               x3 = cos(x_grid[2*d+1]*2.0*PI);
00859               temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00860               f_grid[d] = 0.1*sin(1+temp)*sin(1+temp);
00861             }
00862             if (mode == RANDOM)
00863             {
00864               for (d = 0; d < m_compare; d++)
00865               {
00866                 x1 = sin(x_compare[2*d+1]*2.0*PI)*cos(x_compare[2*d]*2.0*PI);
00867                 x2 = sin(x_compare[2*d+1]*2.0*PI)*sin(x_compare[2*d]*2.0*PI);
00868                 x3 = cos(x_compare[2*d+1]*2.0*PI);
00869                 temp = sqrt(x1*x1)+sqrt(x2*x2)+sqrt(x3*x3);
00870                 f_compare[d] = 0.1*sin(1+temp)*sin(1+temp);
00871               }
00872             }
00873             else
00874             {
00875               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00876             }
00877             break;
00878           case FUNCTION_F6:
00879             for (d = 0; d < m_total; d++)
00880             {
00881               if (x_grid[2*d+1] <= 0.25)
00882               {
00883                 f_grid[d] = 1.0;
00884               }
00885               else
00886               {
00887                 f_grid[d] = 1.0/(sqrt(1+3*cos(2.0*PI*x_grid[2*d+1])*cos(2.0*PI*x_grid[2*d+1])));
00888               }
00889             }
00890             if (mode == RANDOM)
00891             {
00892               for (d = 0; d < m_compare; d++)
00893               {
00894                 if (x_compare[2*d+1] <= 0.25)
00895                 {
00896                   f_compare[d] = 1.0;
00897                 }
00898                 else
00899                 {
00900                   f_compare[d] = 1.0/(sqrt(1+3*cos(2.0*PI*x_compare[2*d+1])*cos(2.0*PI*x_compare[2*d+1])));
00901                 }
00902               }
00903             }
00904             else
00905             {
00906               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00907             }
00908             break;
00909           default:
00910             //fprintf(stderr,"Generating one function\n");
00911             //fflush(stderr);
00912             for (d = 0; d < m_total; d++)
00913             {
00914               f_grid[d] = 1.0;
00915             }
00916             if (mode == RANDOM)
00917             {
00918               for (d = 0; d < m_compare; d++)
00919               {
00920                 f_compare[d] = 1.0;
00921               }
00922             }
00923             else
00924             {
00925               memcpy(f_compare,f_grid,m_total*sizeof(double _Complex));
00926             }
00927             break;
00928         }
00929 
00930         //fprintf(stderr,"Initializing trafo\n");
00931         //fflush(stderr);
00932         /* Init transform plan. */
00933         nfsft_init_guru(&plan_adjoint,NQ[iNQ],m_total, NFSFT_NORMALIZED |
00934           ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00935           ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00936           ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00937           FFT_OUT_OF_PLACE, cutoff);
00938 
00939         plan_adjoint_ptr = &plan_adjoint;
00940 
00941         if (mode == RANDOM)
00942         {
00943           nfsft_init_guru(&plan,NQ[iNQ],m_compare, NFSFT_NORMALIZED |
00944             ((use_nfft!=NO)?(0U):(NFSFT_USE_NDFT)) |
00945             ((use_fpt!=NO)?(0U):(NFSFT_USE_DPT)),
00946             ((NQ[iNQ]>512)?(0U):(PRE_PHI_HUT | PRE_PSI)) | FFTW_INIT |
00947             FFT_OUT_OF_PLACE, cutoff);
00948           plan_ptr = &plan;
00949         }
00950         else
00951         {
00952           plan_ptr = &plan_adjoint;
00953         }
00954 
00955         f_hat = (double _Complex*) nfft_malloc(NFSFT_F_HAT_SIZE(NQ[iNQ])*sizeof(double _Complex));
00956 
00957         plan_adjoint_ptr->f_hat = f_hat;
00958         plan_adjoint_ptr->x = x_grid;
00959         plan_adjoint_ptr->f = f_grid;
00960 
00961         plan_ptr->f_hat = f_hat;
00962         plan_ptr->x = x_compare;
00963         plan_ptr->f = f;
00964 
00965         //fprintf(stderr,"Precomputing for x\n");
00966         //fflush(stderr);
00967         nfsft_precompute_x(plan_adjoint_ptr);
00968         if (plan_adjoint_ptr != plan_ptr)
00969         {
00970           nfsft_precompute_x(plan_ptr);
00971         }
00972 
00973         /* Initialize cumulative time variable. */
00974         t_avg = 0.0;
00975         err_infty_avg = 0.0;
00976         err_2_avg = 0.0;
00977 
00978         /* Cycle through all runs. */
00979         for (i = 0; i < 1/*repetitions*/; i++)
00980         {
00981           //fprintf(stderr,"Copying original values\n");
00982           //fflush(stderr);
00983           /* Copy exact funtion values to working array. */
00984           //memcpy(f,f_grid,m_total*sizeof(double _Complex));
00985 
00986           /* Initialize time measurement. */
00987           t0 = getticks();
00988 
00989           //fprintf(stderr,"Multiplying with quadrature weights\n");
00990           //fflush(stderr);
00991           /* Multiplication with the quadrature weights. */
00992           /*fprintf(stderr,"\n");*/
00993           d = 0;
00994           for (k = 0; k < m_theta; k++)
00995           {
00996             for (n = 0; n < m_phi; n++)
00997             {
00998               /*fprintf(stderr,"f_ref[%d] = %le + I*%le,\t f[%d] = %le + I*%le,  \t w[%d] = %le\n",
00999               d,creal(f_ref[d]),cimag(f_ref[d]),d,creal(f[d]),cimag(f[d]),k,
01000               w[k]);*/
01001               f_grid[d] *= w[k];
01002               d++;
01003             }
01004           }
01005 
01006           t1 = getticks();
01007           t_avg += nfft_elapsed_seconds(t1,t0);
01008 
01009           nfft_free(w);
01010 
01011           t0 = getticks();
01012 
01013           /*fprintf(stderr,"\n");
01014           d = 0;
01015           for (d = 0; d < grid_total; d++)
01016           {
01017             fprintf(stderr,"f[%d] = %le + I*%le, theta[%d] = %le, phi[%d] = %le\n",
01018                     d,creal(f[d]),cimag(f[d]),d,x[2*d+1],d,x[2*d]);
01019           }*/
01020 
01021           //fprintf(stderr,"Executing adjoint\n");
01022           //fflush(stderr);
01023           /* Check if the fast NFSFT algorithm shall be tested. */
01024           if (use_nfsft != NO)
01025           {
01026             /* Execute the adjoint NFSFT transformation. */
01027             nfsft_adjoint(plan_adjoint_ptr);
01028           }
01029           else
01030           {
01031             /* Execute the adjoint direct NDSFT transformation. */
01032             nfsft_adjoint_direct(plan_adjoint_ptr);
01033           }
01034 
01035           /* Multiplication with the Fourier-Legendre coefficients. */
01036           /*for (k = 0; k <= m[im]; k++)
01037           {
01038             for (n = -k; n <= k; n++)
01039             {
01040               fprintf(stderr,"f_hat[%d,%d] = %le\t + I*%le\n",k,n,
01041                       creal(f_hat[NFSFT_INDEX(k,n,&plan_adjoint)]),
01042                       cimag(f_hat[NFSFT_INDEX(k,n,&plan_adjoint)]));
01043             }
01044           }*/
01045 
01046           //fprintf(stderr,"Executing trafo\n");
01047           //fflush(stderr);
01048           if (use_nfsft != NO)
01049           {
01050             /* Execute the NFSFT transformation. */
01051             nfsft_trafo(plan_ptr);
01052           }
01053           else
01054           {
01055             /* Execute the direct NDSFT transformation. */
01056             nfsft_trafo_direct(plan_ptr);
01057           }
01058 
01059           t1 = getticks();
01060           t_avg += nfft_elapsed_seconds(t1,t0);
01061 
01062           //fprintf(stderr,"Finalizing\n");
01063           //fflush(stderr);
01064           /* Finalize the NFSFT plans */
01065           nfsft_finalize(plan_adjoint_ptr);
01066           if (plan_ptr != plan_adjoint_ptr)
01067           {
01068             nfsft_finalize(plan_ptr);
01069           }
01070 
01071           /* Free data arrays. */
01072           nfft_free(f_hat);
01073           nfft_free(x_grid);
01074 
01075           err_infty_avg += X(error_l_infty_complex)(f, f_compare, m_compare);
01076           err_2_avg += X(error_l_2_complex)(f, f_compare, m_compare);
01077 
01078           nfft_free(f_grid);
01079 
01080           if (mode == RANDOM)
01081           {
01082           }
01083           else
01084           {
01085             nfft_free(f_compare);
01086           }
01087 
01088           /*for (d = 0; d < m_total; d++)
01089           {
01090             fprintf(stderr,"f_ref[%d] = %le + I*%le,\t f[%d] = %le + I*%le\n",
01091               d,creal(f_ref[d]),cimag(f_ref[d]),d,creal(f[d]),cimag(f[d]));
01092           }*/
01093         }
01094 
01095         //fprintf(stderr,"Calculating the error\n");
01096         //fflush(stderr);
01097         /* Calculate average time needed. */
01098         t_avg = t_avg/((double)repetitions);
01099 
01100         /* Calculate the average error. */
01101         err_infty_avg = err_infty_avg/((double)repetitions);
01102 
01103         /* Calculate the average error. */
01104         err_2_avg = err_2_avg/((double)repetitions);
01105 
01106         /* Print out the error measurements. */
01107         fprintf(stdout,"%+le %+le %+le\n", t_avg, err_infty_avg, err_2_avg);
01108         fprintf(stderr,"%d: %4d %4d %+le %+le %+le\n", tc, NQ[iNQ], SQ[iNQ],
01109           t_avg, err_infty_avg, err_2_avg);
01110       }
01111     } /* for (im = 0; im < im_max; im++) - Process all cut-off
01112        * bandwidths.*/
01113     fprintf(stderr,"\n");
01114 
01115     /* Delete precomputed data. */
01116     nfsft_forget();
01117 
01118     /* Free memory for cut-off bandwidths and grid size parameters. */
01119     nfft_free(NQ);
01120     nfft_free(SQ);
01121     if (testmode == TIMING)
01122     {
01123       nfft_free(RQ);
01124     }
01125 
01126     if (mode == RANDOM)
01127     {
01128       nfft_free(x_compare);
01129       nfft_free(f_compare);
01130       nfft_free(f);
01131     }
01132 
01133     if (testmode == TIMING)
01134     {
01135       /* Allocate data structures. */
01136       nfft_free(f_hat);
01137       nfft_free(f);
01138       nfft_free(x_grid);
01139     }
01140 
01141   } /* for (tc = 0; tc < tc_max; tc++) - Process each testcase. */
01142 
01143   /* Return exit code for successful run. */
01144   return EXIT_SUCCESS;
01145 }
01146 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409