Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

emfft-debug.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  * 
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  * 
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  * 
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  * 
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  * 
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  * 
00034  * */
00035  
00036 #include <string>
00037 #include "emfft.h"
00038 #include "log.h"
00039 
00040 
00041 #include <iostream>
00042 using std::cout;
00043 using std::endl;
00044 
00045 #include "util.h"
00046 
00047 #ifdef EMAN2_USING_CUDA
00048 #include "cuda/cuda_emfft.h"
00049 #endif 
00050 
00051 #ifdef DJBFFT
00052 extern "C" {
00053         #include <fftr4.h>
00054 }
00055 #endif  //DJBFFT
00056 
00057 using namespace EMAN;
00058 
00059 namespace {
00060         int get_rank(int ny, int nz)
00061         {
00062                 int rank = 3;
00063                 if (ny == 1) {
00064                         rank = 1;
00065                 }
00066                 else if (nz == 1) {
00067                         rank = 2;
00068                 }
00069                 return rank;
00070         }
00071 }
00072 
00073 #ifdef FFTW_PLAN_CACHING
00074 // The only thing important about these constants is that they don't equal each other
00075 const int EMfft::EMAN2_REAL_2_COMPLEX = 1;
00076 const int EMfft::EMAN2_COMPLEX_2_REAL = 2;
00077 
00078 // It is important that these constants are 1 and 0, exactly as they are here
00079 const int EMfft::EMAN2_FFTW2_INPLACE = 1;
00080 const int EMfft::EMAN2_FFTW2_OUT_OF_PLACE=0;
00081 
00082 
00083 #ifdef FFTW3
00084 EMfft::EMfftw3_cache::EMfftw3_cache() :
00085                 num_plans(0)
00086 {
00087         for(int i = 0; i < EMFFTW3_CACHE_SIZE; ++i)
00088         {
00089                 rank[i] = 0;
00090                 plan_dims[i][0] = 0; plan_dims[i][1] = 0; plan_dims[i][2] = 0;
00091                 r2c[i] = -1;
00092                 ip[i] = -1;
00093                 fftwplans[i] = NULL;
00094         }
00095 }
00096 
00097 void EMfft::EMfftw3_cache::debug_plans()
00098 {
00099         for(int i = 0; i < EMFFTW3_CACHE_SIZE; ++i)
00100         {
00101                 cout << "Plan " << i << " has dims " << plan_dims[i][0] << " " 
00102                                 << plan_dims[i][1] << " " << 
00103                                 plan_dims[i][2] << ", rank " <<
00104                                 rank[i] << ", rc flag " 
00105                                 << r2c[i] << ", ip flag " << ip[i] << endl;
00106         }
00107 }
00108 
00109 EMfft::EMfftw3_cache::~EMfftw3_cache()
00110 {
00111         for(int i = 0; i < EMFFTW3_CACHE_SIZE; ++i)
00112         {
00113                 if (fftwplans[i] != NULL)
00114                 {
00115                         fftwf_destroy_plan(fftwplans[i]);
00116                         fftwplans[i] = NULL;
00117                 }
00118         }
00119 }
00120 
00121 fftwf_plan EMfft::EMfftw3_cache::get_plan(const int rank_in, const int x, const int y, const int z, const int r2c_flag, const int ip_flag, fftwf_complex* complex_data, float* real_data )
00122 {
00123 
00124         if ( rank_in > 3 || rank_in < 1 ) throw InvalidValueException(rank_in, "Error, can not get an FFTW plan using rank out of the range [1,3]");
00125         if ( r2c_flag != EMAN2_REAL_2_COMPLEX && r2c_flag != EMAN2_COMPLEX_2_REAL ) throw InvalidValueException(r2c_flag, "The real two complex flag is not supported");
00126         
00127 //      static int num_added = 0;
00128 //      cout << "Was asked for " << rank_in << " " << x << " " << y << " " << z << " " << r2c_flag << endl;
00129         
00130         int dims[3];
00131         dims[0] = z;
00132         dims[1] = y;
00133         dims[2] = x;
00134         
00135         // First check to see if we already have the plan
00136         int i;
00137         for (i=0; i<num_plans; i++) {
00138                 if (plan_dims[i][0]==x && plan_dims[i][1]==y && plan_dims[i][2]==z 
00139                                   && rank[i]==rank_in && r2c[i]==r2c_flag && ip[i]==ip_flag) return fftwplans[i];
00140         }
00141         
00142         fftwf_plan plan;
00143         // Create the plan
00144         if ( y == 1 && z == 1 )
00145         {
00146                 if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00147                         plan = fftwf_plan_dft_r2c_1d(x, real_data, complex_data, FFTW_ESTIMATE);
00148                 else // r2c_flag == EMAN2_COMPLEX_2_REAL, this is guaranteed by the error checking at the beginning of the function
00149                         plan = fftwf_plan_dft_c2r_1d(x, complex_data, real_data, FFTW_ESTIMATE);
00150         }
00151         else
00152         {
00153                 if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00154                         plan = fftwf_plan_dft_r2c(rank_in, dims + (3 - rank_in), real_data, complex_data, FFTW_ESTIMATE);
00155                 else // r2c_flag == EMAN2_COMPLEX_2_REAL, this is guaranteed by the error checking at the beginning of the function
00156                         plan = fftwf_plan_dft_c2r(rank_in, dims + (3 - rank_in), complex_data, real_data, FFTW_ESTIMATE);
00157         }
00158 
00159         if (fftwplans[EMFFTW3_CACHE_SIZE-1] != NULL )
00160         {
00161                 fftwf_destroy_plan(fftwplans[EMFFTW3_CACHE_SIZE-1]);
00162                 fftwplans[EMFFTW3_CACHE_SIZE-1] = NULL;
00163         }
00164                                 
00165         int upper_limit = num_plans;
00166         if ( upper_limit == EMFFTW3_CACHE_SIZE ) upper_limit -= 1;
00167         for (int i=upper_limit-1; i>0; i--)
00168         {
00169                 fftwplans[i]=fftwplans[i-1];
00170                 rank[i]=rank[i-1];
00171                 r2c[i]=r2c[i-1];
00172                 ip[i]=ip[i-1];
00173                 plan_dims[i][0]=plan_dims[i-1][0];
00174                 plan_dims[i][1]=plan_dims[i-1][1];
00175                 plan_dims[i][2]=plan_dims[i-1][2];
00176         }
00177                 //dimplan[0]=-1;
00178 
00179         plan_dims[0][0]=x;
00180         plan_dims[0][1]=y;
00181         plan_dims[0][2]=z;
00182         r2c[0]=r2c_flag;
00183         ip[0]=ip_flag;
00184         fftwplans[0] = plan;
00185         rank[0]=rank_in;
00186         if (num_plans<EMFFTW3_CACHE_SIZE) num_plans++;
00187 //                      debug_plans();
00188 //                      cout << "Created plan 0" << endl;
00189 //      ++num_added;
00190 //      cout << "I have created " << num_added << " plans" << endl;
00191         return fftwplans[0];
00192 
00193 }
00194 
00195 // Static init
00196 EMfft::EMfftw3_cache EMfft::plan_cache;
00197 
00198 #endif // FFTW3
00199 
00200 #ifdef FFTW2
00201 
00202 
00203 EMfft::EMfftw2_cache_nd::EMfftw2_cache_nd() :
00204                 num_plans(0)
00205 {
00206         for(int i = 0; i < EMFFTW2_ND_CACHE_SIZE; ++i)
00207         {
00208                 rank[i] = 0;
00209                 plan_dims[i][0] = 0; plan_dims[i][1] = 0; plan_dims[i][2] = 0;
00210                 r2c[i] = -1;
00211                 ip[i] = -1;
00212                 rfftwnd_plans[i] = NULL;
00213         }
00214 }
00215 
00216 void EMfft::EMfftw2_cache_nd::debug_plans()
00217 {
00218         for(int i = 0; i < EMFFTW2_ND_CACHE_SIZE; ++i)
00219         {
00220                 cout << "Plan " << i << " has dims " << plan_dims[i][0] << " "  << plan_dims[i][1] << " " <<  plan_dims[i][2] << ", rank " <<
00221                                 rank[i] << ", rc flag " << r2c[i] << ", ip flag " << ip[i] << endl;
00222         }
00223 }
00224 
00225 EMfft::EMfftw2_cache_nd::~EMfftw2_cache_nd()
00226 {
00227         for(int i = 0; i < EMFFTW2_ND_CACHE_SIZE; ++i)
00228         {
00229                 if (rfftwnd_plans[i] != NULL)
00230                 {
00231                         rfftwnd_destroy_plan(rfftwnd_plans[i]);
00232                         rfftwnd_plans[i] = NULL;
00233                 }
00234         }
00235 }
00236 
00237 rfftwnd_plan EMfft::EMfftw2_cache_nd::get_plan(const int rank_in, const int x, const int y, const int z, const int r2c_flag, const int ip_flag)
00238 {
00239         
00240         if ( rank_in > 3 || rank_in < 1 ) throw InvalidValueException(rank_in, "Error, can not get an FFTW2 plan using rank out of the range [2,3]");
00241         if ( r2c_flag != EMAN2_REAL_2_COMPLEX && r2c_flag != EMAN2_COMPLEX_2_REAL ) throw InvalidValueException(r2c_flag, "The real two complex flag is not supported");
00242         
00243         int dims[3];
00244         dims[0] = z;
00245         dims[1] = y;
00246         dims[2] = x;
00247 
00248         // First check to see if we already have the plan
00249         int i;
00250         for (i=0; i<num_plans; i++) {
00251                 if (plan_dims[i][0]==x && plan_dims[i][1]==y && plan_dims[i][2]==z 
00252                                   && rank[i]==rank_in && r2c[i]==r2c_flag && ip[i] == ip_flag) break;
00253         }
00254         
00255         // If we found the plan above, the i will not equal num_plans
00256         if (i!=num_plans) {
00257                 return rfftwnd_plans[i];
00258         }       
00259         else
00260         {
00261                 rfftwnd_plan plan;
00262                 // Create the plan
00263                 if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00264                 {
00265                         if (ip_flag == EMAN2_FFTW2_INPLACE) 
00266                                 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE|FFTW_IN_PLACE);
00267                         else 
00268                                 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00269                 }
00270                 else // r2c_flag == EMAN2_COMPLEX_2_REAL, this is guaranteed by the error checking at the beginning of the function
00271                 {
00272                         if (ip_flag == EMAN2_FFTW2_INPLACE)
00273                                 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE|FFTW_IN_PLACE);
00274                         else 
00275                                 plan = rfftwnd_create_plan(rank_in, dims + (3 - rank_in),FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE);
00276                 }
00277                 
00278                                 
00279                 if (rfftwnd_plans[EMFFTW2_ND_CACHE_SIZE-1] != NULL )
00280                 {
00281                         rfftwnd_destroy_plan(rfftwnd_plans[EMFFTW2_ND_CACHE_SIZE-1]);
00282                         rfftwnd_plans[EMFFTW2_ND_CACHE_SIZE-1] = NULL;
00283                 }
00284                 
00285                 int upper_limit = num_plans;
00286                 if ( upper_limit == EMFFTW2_ND_CACHE_SIZE ) upper_limit -= 1;
00287                 for (int i=upper_limit-1; i>0; i--)
00288                 {
00289                         rfftwnd_plans[i]=rfftwnd_plans[i-1];
00290                         rank[i]=rank[i-1];
00291                         r2c[i]=r2c[i-1];
00292                         ip[i]=ip[i-1];
00293                         plan_dims[i][0]=plan_dims[i-1][0];
00294                         plan_dims[i][1]=plan_dims[i-1][1];
00295                         plan_dims[i][2]=plan_dims[i-1][2];
00296                 }
00297                 //dimplan[0]=-1;
00298 
00299                 plan_dims[0][0]=x;
00300                 plan_dims[0][1]=y;
00301                 plan_dims[0][2]=z;
00302                 r2c[0]=r2c_flag;
00303                 rfftwnd_plans[0] = plan;
00304                 rank[0]=rank_in;
00305                 ip[0]=ip_flag;
00306                 if (num_plans<EMFFTW2_ND_CACHE_SIZE) num_plans++;
00307                 
00308                 return rfftwnd_plans[0];
00309         }
00310 }
00311 
00312 EMfft::EMfftw2_cache_1d::EMfftw2_cache_1d() :
00313                 num_plans(0)
00314 {
00315         for(int i = 0; i < EMFFTW2_1D_CACHE_SIZE; ++i)
00316         {
00317                 plan_dims[i] = 0;
00318                 r2c[i] = -1;
00319                 rfftw1d_plans[i] = NULL;
00320         }
00321 }
00322 
00323 void EMfft::EMfftw2_cache_1d::debug_plans()
00324 {
00325         for(int i = 0; i < EMFFTW2_1D_CACHE_SIZE; ++i)
00326         {
00327                 cout << "Plan " << i << " has length " << plan_dims[i] << ", rc flag " << r2c[i] << endl;
00328         }
00329 }
00330 
00331 EMfft::EMfftw2_cache_1d::~EMfftw2_cache_1d()
00332 {
00333         for(int i = 0; i < EMFFTW2_1D_CACHE_SIZE; ++i)
00334         {
00335                 if (rfftw1d_plans[i] != NULL)
00336                 {
00337                         rfftw_destroy_plan(rfftw1d_plans[i]);
00338                         rfftw1d_plans[i] = NULL;
00339                 }
00340         }
00341 }
00342 
00343 rfftw_plan EMfft::EMfftw2_cache_1d::get_plan(const int x, const int r2c_flag )
00344 {
00345         
00346         if ( r2c_flag != EMAN2_REAL_2_COMPLEX && r2c_flag != EMAN2_COMPLEX_2_REAL ) throw InvalidValueException(r2c_flag, "The real two complex flag is not supported");
00347         
00348         // First check to see if we already have the plan
00349         int i;
00350         for (i=0; i<num_plans; i++) {
00351                 if (plan_dims[i]==x && r2c[i]==r2c_flag ) return rfftw1d_plans[i];
00352         }
00353         
00354         rfftw_plan plan;
00355         // Create the plan
00356         if ( r2c_flag == EMAN2_REAL_2_COMPLEX )
00357                         plan = rfftw_create_plan(x, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00358         else // r2c_flag == EMAN2_COMPLEX_2_REAL, this is guaranteed by the error checking at the beginning of the function
00359                         plan = rfftw_create_plan(x, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00360                 
00361 
00362                         
00363         if (rfftw1d_plans[EMFFTW2_1D_CACHE_SIZE-1] != NULL )
00364         {
00365                 rfftw_destroy_plan(rfftw1d_plans[EMFFTW2_1D_CACHE_SIZE-1]);
00366                 rfftw1d_plans[EMFFTW2_1D_CACHE_SIZE-1] = NULL;
00367         }
00368                                 
00369         int upper_limit = num_plans;
00370         if ( upper_limit == EMFFTW2_1D_CACHE_SIZE ) upper_limit -= 1;
00371         for (int i=upper_limit-1; i>0; i--)
00372         {
00373                 rfftw1d_plans[i]=rfftw1d_plans[i-1];
00374                 plan_dims[i]=plan_dims[i-1];
00375                 r2c[i] = r2c[i-1];
00376         }
00377         //dimplan[0]=-1;
00378 
00379         plan_dims[0]=x;
00380         r2c[0]=r2c_flag;
00381         rfftw1d_plans[0] = plan;
00382         if (num_plans<EMFFTW2_1D_CACHE_SIZE) num_plans++;
00383                 
00384         return rfftw1d_plans[0];
00385 
00386 }
00387 
00388 // Static init
00389 EMfft::EMfftw2_cache_nd EMfft::plan_nd_cache;
00390 EMfft::EMfftw2_cache_1d EMfft::plan_1d_cache;
00391 
00392 #endif // FFTW2
00393 #endif // FFTW_PLAN_CACHING
00394 
00395 #ifdef FFTW2
00396 
00397 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00398 {
00399 #ifdef DJBFFT
00400         switch(n)
00401         {
00402                 if ( n==2 || n==4 || n==8 || n==16 || n==32 || n==64 || n==128
00403                         || n==256 || n==512 || n==1024 || n==2048 || n==4096 || n==8192 )
00404                 {
00405                         memcpy( complex_data, real_data, n * sizeof(float) );           
00406                 }
00407                 
00408                 case 2:
00409                         fftr4_2( (real4 *)complex_data );
00410                         
00411                 case 4: 
00412                         fftr4_4( (real4 *)complex_data );
00413                         
00414                 case 8:
00415                         fftr4_8( (real4 *)complex_data );
00416                          
00417                 case 16:
00418                         fftr4_16( (real4 *)complex_data ); 
00419                         
00420                 case 32:
00421                         fftr4_32( (real4 *)complex_data ); 
00422                         
00423                 case 64:
00424                         fftr4_64( (real4 *)complex_data ); 
00425                         
00426                 case 128: 
00427                         fftr4_128( (real4 *)complex_data );
00428                         
00429                 case 256:
00430                         fftr4_256( (real4 *)complex_data ); 
00431                         
00432                 case 512: 
00433                         fftr4_512( (real4 *)complex_data );
00434                         
00435                 case 1024: 
00436                         fftr4_1024( (real4 *)complex_data );
00437                         
00438                 case 2048: 
00439                         fftr4_2048( (real4 *)complex_data );
00440                         
00441                 case 4096: 
00442                         fftr4_4096( (real4 *)complex_data );
00443                         
00444                 case 8192:
00445                         fftr4_8192( (real4 *)complex_data );
00446                         
00447                 default:
00448                         const int complex_n = n + 2 - n%2;
00449                         float * fft_data = new float[n];
00450                         rfftw_plan p = rfftw_create_plan(n, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00451                         rfftw_one(p, (fftw_real *) real_data, (fftw_real *) complex_data);
00452                         rfftw_destroy_plan(p);
00453                         for(int i=0; i<complex_n; ++i) {
00454                                 if(i%2==0) {    //copy real part of complex array
00455                                         complex_data[i] = fft_data[i/2];
00456                                 }
00457                                 else {  //copy imaginary part of complex array
00458                                         if(i==1) {
00459                                                 complex_data[i] = 0.0f;
00460                                         }
00461                                         else {
00462                                                 if(n%2 == 0 && i == complex_n-1 ) {
00463                                                         complex_data[i] = 0.0f;
00464                                                 }
00465                                                 else {
00466                                                         complex_data[i] = fft_data[n-i/2];
00467                                                 }
00468                                         }
00469                                 }
00470                         }
00471         
00472                         delete [] fft_data;
00473         }       
00474 #else
00475         const int complex_n = n + 2 - n%2;
00476         float * fft_data = new float[n];
00477 #ifdef FFTW_PLAN_CACHING
00478         rfftw_plan plan_1d = plan_1d_cache.get_plan(n,EMAN2_REAL_2_COMPLEX);
00479         rfftw_one(plan_1d, (fftw_real *) real_data, (fftw_real *) fft_data);
00480 #else
00481         rfftw_plan      p = rfftw_create_plan(n, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00482         rfftw_one(p, (fftw_real *) real_data, (fftw_real *) fft_data);
00483         rfftw_destroy_plan(p);
00484 #endif //FFTW_PLAN_CACHING
00485         
00486         for(int i=0; i<complex_n; ++i) {
00487                 if(i%2==0) {    //copy real part of complex array
00488                         complex_data[i] = fft_data[i/2];
00489                 }
00490                 else {  //copy imaginary part of complex array
00491                         if(i==1) {
00492                                 complex_data[i] = 0.0f;
00493                         }
00494                         else {
00495                                 if(n%2 == 0 && i == complex_n-1 ) {
00496                                         complex_data[i] = 0.0f;
00497                                 }
00498                                 else {
00499                                         complex_data[i] = fft_data[n-i/2];
00500                                 }
00501                         }
00502                 }
00503         }
00504         
00505         delete [] fft_data;
00506 #endif  //DJBFFT        
00507 
00508         return 0;
00509 }
00510 
00511 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00512 {
00513 #ifdef DJBFFT
00514         switch(n)
00515         {
00516                 if ( n==2 || n==4 || n==8 || n==16 || n==32 || n==64 || n==128
00517                         || n==256 || n==512 || n==1024 || n==2048 || n==4096 || n==8192 )
00518                 {
00519                         memcpy( real_data, complex_data, n * sizeof(float) );           
00520                 }
00521                 
00522                 case 2: 
00523                         fftr4_scale2( (real4 *)real_data );
00524                 fftr4_un2( (real4 *)real_data );
00525                         
00526                 case 4: 
00527                         fftr4_scale4( (real4 *)real_data );
00528                 fftr4_un4( (real4 *)real_data );
00529                         
00530                 case 8:
00531                         fftr4_scale8( (real4 *)real_data );
00532                 fftr4_un8( (real4 *)real_data );
00533                          
00534                 case 16:
00535                         fftr4_scale16( (real4 *)real_data );
00536                 fftr4_un16( (real4 *)real_data );
00537                         
00538                 case 32:
00539                         fftr4_scale32( (real4 *)real_data );
00540                 fftr4_un32( (real4 *)real_data );
00541                         
00542                 case 64:
00543                         fftr4_scale64( (real4 *)real_data );
00544                 fftr4_un64( (real4 *)real_data );
00545                         
00546                 case 128: 
00547                         fftr4_scale128( (real4 *)real_data );
00548                 fftr4_un128( (real4 *)real_data );
00549                         
00550                 case 256:
00551                         fftr4_scale256( (real4 *)real_data );
00552                 fftr4_un256( (real4 *)real_data );
00553                         
00554                 case 512: 
00555                         fftr4_scale512( (real4 *)real_data );
00556                 fftr4_un512( (real4 *)real_data );
00557                         
00558                 case 1024: 
00559                         fftr4_scale1024( (real4 *)real_data );
00560                 fftr4_un1024( (real4 *)real_data );
00561                         
00562                 case 2048: 
00563                         fftr4_scale2048( (real4 *)real_data );
00564                 fftr4_un2048( (real4 *)real_data );
00565                         
00566                 case 4096: 
00567                         fftr4_scale4096( (real4 *)real_data );
00568                 fftr4_un4096( (real4 *)real_data );
00569                         
00570                 case 8192:
00571                         fftr4_scale8192( (real4 *)real_data );
00572                 fftr4_un8192( (real4 *)real_data );
00573                         
00574                 default:
00575                         const int complex_n = n + 2 - n%2;
00576                         float * fft_data = new float[n];
00577                         
00578                         for(int i=0; i<complex_n; ++i) {
00579                                 if(i%2 == 0) {  //copy real part of complex array
00580                                         fft_data[i/2] = complex_data[i];
00581                                 }
00582                                 else {  //copy imaginary part of complex array
00583                                         if(i==1) {continue;}
00584                                         if(!(n%2 == 0 && i == complex_n-1)) {
00585                                                 fft_data[n-i/2] = complex_data[i];
00586                                         }
00587                                 }
00588                         }
00589                         
00590                         rfftw_plan p = rfftw_create_plan(n, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00591                         rfftw_one(p, (fftw_real *) fft_data, (fftw_real *) real_data);
00592                         rfftw_destroy_plan(p);
00593         }
00594 #else
00595         const int complex_n = n + 2 - n%2;
00596         float * fft_data = new float[n];
00597         
00598         for(int i=0; i<complex_n; ++i) {
00599                 if(i%2 == 0) {  //copy real part of complex array
00600                         fft_data[i/2] = complex_data[i];
00601                 }
00602                 else {  //copy imaginary part of complex array
00603                         if(i==1) {continue;}
00604                         if(!(n%2 == 0 && i == complex_n-1)) {
00605                                 fft_data[n-i/2] = complex_data[i];
00606                         }
00607                 }
00608         }
00609 #ifdef FFTW_PLAN_CACHING
00610         rfftw_plan plan_1d = plan_1d_cache.get_plan(n,EMAN2_COMPLEX_2_REAL);
00611         rfftw_one(plan_1d, (fftw_real *) fft_data, (fftw_real *) real_data);
00612 #else
00613         rfftw_plan p = rfftw_create_plan(n, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00614         rfftw_one(p, (fftw_real *) fft_data, (fftw_real *) real_data);
00615         rfftw_destroy_plan(p);
00616 #endif // FFTW_PLAN_CACHING
00617         delete [] fft_data;
00618 #endif  //DJBFFT
00619         
00620         return 0;
00621         
00622 }
00623 
00624 
00625 // ming add c->c fft with fftw2 library//
00626 int EMfft::complex_to_complex_1d(float *complex_data_in, float *complex_data_out, int n)
00627 {
00628         fftw_complex *in=(fftw_complex *) complex_data_in;
00629         fftw_complex *out=(fftw_complex *) complex_data_out;
00630         fftw_plan p = fftw_create_plan(n/2, FFTW_FORWARD, FFTW_ESTIMATE);
00631         fftw_one(p,in,out);
00632         fftw_destroy_plan(p);
00633         return 0;
00634 }
00635 
00636 
00637 int EMfft::complex_to_complex_nd(float *complex_data_in, float *complex_data_out, int nx,int ny,int nz)
00638 {
00639         const int rank = get_rank(ny, nz);
00640         int dims[3];
00641         dims[0] = nz;
00642         dims[1] = ny;
00643         dims[2] = nx;
00644 
00645         switch(rank) {
00646                 case 1:
00647                         complex_to_complex_1d(complex_data_in, complex_data_out, nx);
00648                         break;
00649 
00650                 case 2:
00651                 case 3:
00652                 {
00653                         fftw_complex *in=(fftw_complex *) complex_data_in;
00654                         fftw_complex *out=(fftw_complex *) complex_data_out;
00655                         fftwnd_plan plan_nd;
00656                         if(complex_data_out == complex_data_in) {
00657                                 plan_nd = fftwnd_create_plan(rank, dims+(3-rank),FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
00658                         }
00659                         else {
00660                                 plan_nd = fftw3d_create_plan(nx/2,ny,nz,FFTW_FORWARD, FFTW_ESTIMATE);
00661                         }
00662                         fftwnd_one(plan_nd, in, out);
00663                         fftwnd_destroy_plan(plan_nd);
00664                 }
00665         }
00666         return 0;
00667 }
00668 
00669 
00670 
00671 
00672 
00673 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00674 {
00675 //      const int complex_nx = nx + 2 - nx%2;
00676         const int rank = get_rank(ny, nz);
00677         int dims[3];
00678         dims[0] = nz;
00679         dims[1] = ny;
00680         dims[2] = nx;
00681         
00682         switch(rank) {
00683                 case 1:
00684                         real_to_complex_1d(real_data, complex_data, nx);
00685                         break;
00686                 
00687                 case 2: 
00688                 case 3:
00689                 {
00690 #ifdef FFTW_PLAN_CACHING
00691                         bool ip = ( complex_data == real_data );
00692                         rfftwnd_plan plan_nd = plan_nd_cache.get_plan(rank,nx,ny,nz,EMAN2_REAL_2_COMPLEX,ip);
00693                         rfftwnd_one_real_to_complex(plan_nd, (fftw_real *) real_data, (fftw_complex *) complex_data);
00694 #else
00695                         rfftwnd_plan plan_nd;
00696                         if(real_data == complex_data) {
00697                                 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank), 
00698                                                 FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
00699                         }
00700                         else {
00701                                 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank), 
00702                                                 FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
00703                         }
00704                         rfftwnd_one_real_to_complex(plan_nd, (fftw_real *) real_data,
00705                                                                                 (fftw_complex *) complex_data);
00706                         rfftwnd_destroy_plan(plan_nd);
00707 #endif // FFTW_PLAN_CACHING
00708                 }
00709                         break;
00710                 
00711                 default:
00712                         LOGERR("Should NEVER be here!!!");
00713                         break;
00714         }
00715         
00716         return 0;
00717 }
00718 
00719 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00720 {
00721 //      const int complex_nx = nx + 2 - nx%2;
00722         const int rank = get_rank(ny, nz);
00723         int dims[3];
00724         dims[0] = nz;
00725         dims[1] = ny;
00726         dims[2] = nx;
00727         
00728         switch(rank) {
00729                 case 1:
00730                         complex_to_real_1d(complex_data, real_data, nx);
00731                         break;
00732                 
00733                 case 2:
00734                 case 3:
00735                 {
00736 #ifdef FFTW_PLAN_CACHING
00737                         bool ip = ( complex_data == real_data );
00738                         rfftwnd_plan plan_nd = plan_nd_cache.get_plan(rank,nx,ny,nz,EMAN2_COMPLEX_2_REAL,ip);
00739                         rfftwnd_one_complex_to_real(plan_nd, (fftw_complex *) complex_data,     (fftw_real *) real_data);
00740 #else
00741                         rfftwnd_plan plan_nd;
00742                         if(real_data == complex_data) {
00743                                 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank), 
00744                                                 FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
00745                         }
00746                         else {
00747                                 plan_nd = rfftwnd_create_plan(rank, dims + (3 - rank), 
00748                                                 FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
00749                         }       
00750                         rfftwnd_one_complex_to_real(plan_nd, (fftw_complex *) complex_data,     
00751                                 (fftw_real *) real_data);
00752                         rfftwnd_destroy_plan(plan_nd);
00753 #endif // FFTW_PLAN_CACHING
00754                 }               
00755                         break;
00756                 
00757                 default:
00758                         LOGERR("Should NEVER be here!!!");
00759                         break;
00760         }                       
00761 
00762         return 0;
00763 }
00764 #endif  //FFTW2
00765 
00766 //#ifdef CUDA_FFT
00767 //int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00768 //{
00769 //      return  cuda_fft_real_to_complex_1d(real_data,complex_data,n);
00770 //}
00771 
00772 //int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00773 //{
00774 //      return cuda_fft_complex_to_real_1d(complex_data,real_data,n);
00775 //}
00776 
00777 //int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00778 //{
00779 //      return cuda_fft_real_to_complex_nd(real_data,complex_data,nx,ny,nz);
00780 //}
00781 
00782 //int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00783 //{
00784 //      return cuda_fft_complex_to_real_nd(complex_data,real_data,nx,ny,nz);
00785 //}
00786 
00787 //#endif
00788 
00789 #ifdef FFTW3
00790 
00791 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00792 {//cout<<"doing fftw3"<<endl;
00793 #ifdef FFTW_PLAN_CACHING
00794         bool ip = ( complex_data == real_data );
00795         fftwf_plan plan = plan_cache.get_plan(1,n,1,1,EMAN2_REAL_2_COMPLEX,ip,(fftwf_complex *) complex_data, real_data);
00796         // According to FFTW3, this is making use of the "guru" interface - this is necessary if plans are to be reused
00797         fftwf_execute_dft_r2c(plan, real_data,(fftwf_complex *) complex_data);
00798 #else
00799         fftwf_plan plan = fftwf_plan_dft_r2c_1d(n, real_data, (fftwf_complex *) complex_data,
00800                                                                                         FFTW_ESTIMATE);
00801         fftwf_execute(plan);
00802         fftwf_destroy_plan(plan);
00803 #endif // FFTW_PLAN_CACHING
00804         return 0;
00805 };
00806 
00807 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
00808 {
00809 #ifdef FFTW_PLAN_CACHING
00810         bool ip = ( complex_data == real_data );
00811         fftwf_plan plan = plan_cache.get_plan(1,n,1,1,EMAN2_COMPLEX_2_REAL,ip,(fftwf_complex *) complex_data, real_data);
00812         // According to FFTW3, this is making use of the "guru" interface - this is necessary if plans are to be reused
00813         fftwf_execute_dft_c2r(plan, (fftwf_complex *) complex_data, real_data);
00814 #else
00815         fftwf_plan plan = fftwf_plan_dft_c2r_1d(n, (fftwf_complex *) complex_data, real_data,
00816                                                                                         FFTW_ESTIMATE);
00817         fftwf_execute(plan);
00818         fftwf_destroy_plan(plan);
00819 #endif // FFTW_PLAN_CACHING
00820         
00821         return 0;
00822 }
00823 
00824 
00825 // ming add c->c fft with fftw3 library//
00826 int EMfft::complex_to_complex_1d(float *complex_data_in, float *complex_data_out, int n)
00827 {
00828         fftwf_plan p;
00829         fftwf_complex *in=(fftwf_complex *) complex_data_in;
00830         fftwf_complex *out=(fftwf_complex *) complex_data_out;
00831         p=fftwf_plan_dft_1d(n/2,in,out, FFTW_FORWARD, FFTW_ESTIMATE);
00832         fftwf_execute(p);
00833         fftwf_destroy_plan(p);
00834         return 0;
00835 }
00836 
00837 
00838 int EMfft::complex_to_complex_nd(float *in, float *out, int nx,int ny,int nz)
00839 {
00840         const int rank = get_rank(ny, nz);
00841         int dims[3];
00842         dims[0] = nz;
00843         dims[1] = ny;
00844         dims[2] = nx;
00845 
00846         switch(rank) {
00847                 case 1:
00848                         complex_to_complex_1d(in, out, nx);
00849                         break;
00850 
00851                 case 2:
00852                 case 3:
00853                 {
00854                         fftwf_plan p;
00855 
00856                         if(out == in) {
00857                                 p=fftwf_plan_dft_3d(nx/2,ny,nz,(fftwf_complex *) in,(fftwf_complex *) out, FFTW_FORWARD, FFTW_ESTIMATE);
00858                         }
00859                         else {
00860 
00861                                 p=fftwf_plan_dft_3d(nx/2,ny,nz,(fftwf_complex *) in,(fftwf_complex *) out, FFTW_FORWARD, FFTW_ESTIMATE);
00862                         }
00863                         fftwf_execute(p);
00864                         fftwf_destroy_plan(p);
00865                 }
00866         }
00867         return 0;
00868 }
00869 // end ming
00870 
00871 
00872 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
00873 {
00874         const int rank = get_rank(ny, nz);
00875         int dims[3];
00876         dims[0] = nz;
00877         dims[1] = ny;
00878         dims[2] = nx;
00879         
00880         switch(rank) {
00881                 case 1:
00882                         real_to_complex_1d(real_data, complex_data, nx);
00883                         break;
00884                 
00885                 case 2:
00886                 case 3:
00887                 {
00888 #ifdef FFTW_PLAN_CACHING
00889                         bool ip = ( complex_data == real_data );
00890                         fftwf_plan plan = plan_cache.get_plan(rank,nx,ny,nz,EMAN2_REAL_2_COMPLEX,ip,(fftwf_complex *) complex_data, real_data);
00891                         // According to FFTW3, this is making use of the "guru" interface - this is necessary if plans are to be re-used
00892                         fftwf_execute_dft_r2c(plan, real_data,(fftwf_complex *) complex_data );
00893 #else
00894                         fftwf_plan plan = fftwf_plan_dft_r2c(rank, dims + (3 - rank), 
00895                                         real_data, (fftwf_complex *) complex_data, FFTW_ESTIMATE);
00896                         fftwf_execute(plan);
00897                         fftwf_destroy_plan(plan);
00898 #endif // FFTW_PLAN_CACHING
00899                 }
00900                         break;
00901                 
00902                 default:
00903                         LOGERR("Should NEVER be here!!!");
00904                         break;
00905         }
00906         
00907         return 0;
00908 }
00909 
00910 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
00911 {
00912         const int rank = get_rank(ny, nz);
00913         int dims[3];
00914         dims[0] = nz;
00915         dims[1] = ny;
00916         dims[2] = nx;
00917         
00918         switch(rank) {
00919                 case 1:
00920                         complex_to_real_1d(complex_data, real_data, nx);
00921                         break;
00922                 
00923                 case 2:
00924                 case 3:
00925                 {
00926 #ifdef FFTW_PLAN_CACHING
00927                         bool ip = ( complex_data == real_data );
00928                         fftwf_plan plan = plan_cache.get_plan(rank,nx,ny,nz,EMAN2_COMPLEX_2_REAL,ip,(fftwf_complex *) complex_data, real_data);
00929                         // According to FFTW3, this is making use of the "guru" interface - this is necessary if plans are to be re-used
00930                         fftwf_execute_dft_c2r(plan, (fftwf_complex *) complex_data, real_data);
00931 #else
00932                         fftwf_plan plan = fftwf_plan_dft_c2r(rank, dims + (3 - rank), 
00933                                         (fftwf_complex *) complex_data, real_data, FFTW_ESTIMATE);
00934                         fftwf_execute(plan);
00935                         fftwf_destroy_plan(plan);
00936 #endif // FFTW_PLAN_CACHING
00937                         
00938                         
00939                 }
00940                         break;
00941                         
00942                 default:
00943                         LOGERR("Should NEVER be here!!!");
00944                         break;
00945         }
00946         
00947         return 0;
00948 }
00949 
00950 #endif  //FFTW3
00951 
00952 #ifdef NATIVE_FFT
00953 #include "sparx/native_fft.h"
00954 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
00955 {
00956         //int complex_size = n + 2 - n%2;
00957         float * work = (float*) malloc((2*n+15)*sizeof(float));
00958         if (!work) {
00959                 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
00960                 LOGERR("real_to_complex_1d: failed to allocate work\n");        
00961         }
00962         
00963         rffti(n, work);
00964         memcpy(&complex_data[1], real_data, n * sizeof(float));
00965         rfftf(n, &complex_data[1], work);
00966         complex_data[0] = complex_data[1] ;
00967         complex_data[1] = 0.0f ;
00968         if (n%2 == 0)  complex_data[n+1] = 0.0f ;
00969 
00970         free(work);
00971         return 0;
00972 }
00973 /*{
00974         //int complex_size = n + 2 - n%2;
00975         
00976         //memcpy(complex_data, real_data, n * sizeof(float));
00977         //float * work = (float*) malloc((2*n+15)*sizeof(float));
00978         //if (!work) {
00979         //      fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
00980         //      LOGERR("real_to_complex_1d: failed to allocate work\n");        
00981         //}
00982         //Nativefft::fmrs_1rf(complex_data, work, n);
00983         //cout<<"doing rfftf"<<endl;
00984         //rffti(n, work);
00985         
00986         //rfftf(n, real_data, work);
00987         
00988         //cout<<"doing fftr_q"<<endl;
00989         int l=(int)(log2(n));
00990         Util::fftr_q(real_data,l);
00991 
00992         //free(work);
00993         return 0;
00994 }//
00995 {
00996         int complex_size = n + 2 - n%2;
00997         
00998         memcpy(complex_data, real_data, n * sizeof(float));
00999         float * work = (float*) malloc(complex_size*sizeof(float));
01000         if (!work) {
01001                 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
01002                 LOGERR("real_to_complex_1d: failed to allocate work\n");        
01003         }
01004         
01005         Nativefft::fmrs_1rf(complex_data, work, n);
01006         
01007         free(work);
01008         return 0;
01009 }*/
01010 
01011 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
01012 {
01013         //here, n is the "logical" size of DFT, not the array size
01014         float * work = (float*) malloc((2*n+15)*sizeof(float));
01015         if (!work) {
01016                 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
01017                 LOGERR("complex_to_real_1d: failed to allocate work\n");        
01018         }
01019         
01020         rffti(n, work);
01021         
01022         memcpy(&real_data[1], &complex_data[2], (n-1) * sizeof(float));
01023         real_data[0] = complex_data[0] ;
01024         rfftb(n, real_data, work);
01025         //  Normalize
01026         float nrm = 1.0f/float(n);
01027         for (int i = 0; i<n; i++) real_data[i] *= nrm;
01028         free(work);
01029         return 0;
01030 }
01031 /*{
01032         int complex_size = n + 2 - n%2;
01033         
01034         //here, n is the "logical" size of DFT, not the array size
01035         memcpy(real_data, complex_data, complex_size * sizeof(float));
01036         float * work = (float*)malloc(complex_size*sizeof(float));
01037         if (!work) {
01038                 fprintf(stderr,"real_to_complex_1d: failed to allocate work\n");
01039                 LOGERR("complex_to_real_1d: failed to allocate work\n");        
01040         }
01041         
01042         Nativefft::fmrs_1rb(real_data, work, n);
01043         
01044         free(work);
01045         return 0;
01046 }*/
01047 
01048 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
01049 {
01050         const int rank = get_rank(ny, nz);
01051         const int complex_nx = nx + 2 - nx%2;
01052         
01053         switch(rank) {
01054                 case 1:         //for 1D fft
01055                         real_to_complex_1d(real_data, complex_data, nx);
01056                         return 0;
01057                 case 2:         //for 2D fft
01058                 {
01059                         /*if(real_data != complex_data) {
01060                                 for (int j = 0; j < ny; j++) {
01061                                         memcpy(&complex_data[complex_nx*j], &real_data[nx*j], nx*sizeof(float));        
01062                                 }
01063                         }
01064                         float * work = (float*) malloc(complex_nx*sizeof(float));
01065                         if (!work) {
01066                                 fprintf(stderr,"real_to_complex_nd(2df): failed to allocate work\n");
01067                                 LOGERR("real_to_complex_nd(2df): failed to allocate work\n");
01068                         }
01069                         
01070                         // 2d inplace fft, overwrite y
01071                         int status = Nativefft::fmrs_2rf(complex_data, work, complex_nx, nx, ny);
01072                         if (status !=0) {
01073                         fprintf(stderr, "real_to_complex_nd(2df): status = %d\n", status);
01074                         LOGWARN("real_to_complex_nd(2df): status = %d\n", status);
01075                         }
01076                  
01077                         free(work);*/
01078                         
01079                         int status = Nativefft::ftp_2rf(real_data, complex_data, complex_nx, nx, ny);
01080                         if (status !=0) {
01081                         fprintf(stderr, "real_to_complex_nd(2df): status = %d\n", status);
01082                         LOGWARN("real_to_complex_nd(2df): status = %d\n", status);
01083                         }
01084                         return 0;
01085                 }
01086                 case 3:         //for 3D fft
01087                 {
01088                         /*if(real_data != complex_data) {
01089                                 for (int k = 0; k<nz; k++) {
01090                                 for (int j = 0; j < ny; j++) {
01091                                         memcpy(&complex_data[complex_nx*ny*k+j*complex_nx], &real_data[nx*ny*k+j*nx], nx*sizeof(float));
01092                                 }
01093                                 }
01094                         }
01095                         float * work = (float*) malloc(1*sizeof(float));//malloc(complex_nx*sizeof(float));
01096                         if (!work) {
01097                                 fprintf(stderr,"real_to_complex_nd(3df): failed to allocate work\n");
01098                                 LOGERR("real_to_complex_nd(3df): failed to allocate work\n");
01099                         }*/
01100                         
01101                         // 3d inplace fft, overwrite complex_data
01102                         int status = Nativefft::ftp_3rf(real_data, complex_data, complex_nx, nx, ny, nz);
01103                         if (status !=0) {
01104                         fprintf(stderr, "real_to_complex_nd(3df): status = %d\n", status);
01105                         LOGWARN("real_to_complex_nd(3df): status = %d\n", status);
01106                         }
01107                         
01108                         //free(work);
01109                         return 0;
01110                 }
01111                 default:
01112                         LOGERR("Never should be here...\n");
01113                         return -1;      
01114         }
01115 }
01116 
01117 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
01118 {
01119         const int rank = get_rank(ny, nz);
01120         const int complex_nx = nx + 2 - nx%2;
01121         
01122         switch(rank) {
01123                 case 1:         //for 1D ift
01124                         complex_to_real_1d(complex_data, real_data, nx);
01125                         return 0;
01126                 case 2:         //for 2D ift
01127                 /*{     
01128                         if(real_data != complex_data) {
01129                                 memcpy(real_data, complex_data, complex_nx*ny*sizeof(float));
01130                         }
01131                         
01132                         float * work = (float*) malloc(complex_nx*sizeof(float));
01133                         if (!work) {
01134                                 fprintf(stderr,"complex_to_real_nd(2db): failed to allocate work\n");
01135                                 LOGERR("complex_to_real_nd(2db): failed to allocate work\n");
01136                         }
01137                         
01138                         // 2d inplace ift, overwrite real_data
01139                         int status = Nativefft::fmrs_2rb(real_data, work, complex_nx, nx, ny);
01140                         if (status !=0) {
01141                         fprintf(stderr, "complex_to_real_nd(2db): status = %d\n", status);
01142                         LOGWARN("complex_to_real_nd(2db): status = %d\n", status);
01143                         }
01144                         
01145                         free(work);
01146                         return 0;
01147                 }*/
01148                 {       //  Only out of place!
01149                         memcpy(real_data, complex_data, complex_nx*ny*sizeof(float));
01150                         
01151                         // 2d inplace ift, overwrite real_data
01152                         int status = Nativefft::ftp_2rb(real_data, complex_nx, nx, ny);
01153 
01154                         if (status !=0) {
01155                         fprintf(stderr, "complex_to_real_nd(2db): status = %d\n", status);
01156                         LOGWARN("complex_to_real_nd(2db): status = %d\n", status);
01157                         }
01158                         return 0;
01159                 }
01160                 case 3:         //for 3D ift
01161                 {       // Only out of place!
01162                         memcpy(real_data, complex_data, complex_nx*ny*nz*sizeof(float));
01163                         
01164                         // 3d inplace fft, overwrite real_data
01165                         int status = Nativefft::ftp_3rb(real_data, complex_nx, nx, ny, nz);
01166                         if (status !=0) {
01167                         fprintf(stderr, "complex_to_real_nd(3db): status = %d\n", status);
01168                         LOGWARN("complex_to_real_nd(3db): status = %d\n", status);
01169                         }
01170                         return 0;
01171                 }
01172                 default:
01173                         LOGERR("Never should be here...\n");
01174                         return -1;
01175         }
01176 }
01177 #endif  //NATIVE_FFT
01178 
01179 #ifdef  ACML
01180 #include <iostream>
01181 using std::cout;
01182 using std::endl;
01183 
01184 int EMfft::real_to_complex_1d(float *real_data, float *complex_data, int n)
01185 {
01186         const int complex_n = n + 2 - n%2;      //the size for 1D complex array
01187         int info;
01188         
01189         /* Allocate communication work array */
01190         float * comm = (float *)malloc((3*n+100)*sizeof(float));
01191         /* Allocate work array to store ACML complex array*/
01192         float * fft_data = (float *)malloc(n * sizeof(float));
01193         
01194         //copy real_data to complex_data then apply inplace FFT on complex data
01195         memcpy(fft_data, real_data, n * sizeof(float));
01196         
01197         /* Initialize communication work array */
01198         scfft(0, n, fft_data, comm, &info);
01199         if(info != 0) {
01200                 LOGERR("Error happening in Initialize communication work array: %d", info);
01201         }
01202         
01203         /* Compute a real --> Hermitian transform */
01204         scfft(1, n, fft_data, comm, &info);
01205         if(info != 0) {
01206                 LOGERR("Error happening in Initialize communication work array: %d", info);
01207         }
01208         
01215         transform(fft_data, fft_data+n, fft_data, time_sqrt_n(n));
01216         
01217         for(int i=0; i<complex_n; ++i) {
01218                 if(i%2==0) {    //copy real part of complex array
01219                         complex_data[i] = fft_data[i/2];
01220                 }
01221                 else {  //copy imaginary part of complex array
01222                         if(i==1) {
01223                                 complex_data[i] = 0.0f;
01224                         }
01225                         else {
01226                                 if(n%2 == 0 && i == complex_n-1 ) {
01227                                         complex_data[i] = 0.0f;
01228                                 }
01229                                 else {
01230                                         complex_data[i] = fft_data[n-i/2];
01231                                 }
01232                         }
01233                 }
01234         }
01235         
01236         free(fft_data);
01237         free(comm);
01238         return 0;
01239 }
01240 
01241 int EMfft::complex_to_real_1d(float *complex_data, float *real_data, int n)
01242 {
01243         int complex_n = n + 2 - n%2;    //the size for 1D complex array
01244         int info;
01245         
01246         /* Allocate communication work array */
01247         float * comm = (float *)malloc((3*n+100)*sizeof(float));
01248                 
01249         for(int i=0; i<complex_n; ++i) {
01250                 if(i%2 == 0) {  //copy real part of complex array
01251                         real_data[i/2] = complex_data[i];
01252                 }
01253                 else {  //copy imaginary part of complex array
01254                         if(i==1) {continue;}
01255                         if(!(n%2 == 0 && i == complex_n-1)) {
01256                                 real_data[n-i/2] = complex_data[i];
01257                         }
01258                 }
01259         }
01260         transform(real_data, real_data+n, real_data, divide_sqrt_n(n));
01261         
01262         /* Initialize communication work array */
01263         csfft(0, n, real_data, comm, &info);
01264         if(info != 0) {
01265                 LOGERR("Error happening in Initialize communication work array: %d", info);
01266         }
01267         
01268         /* Conjugate the Vector X to simulate inverse transform */
01269         for (int j = n/2+1; j < n; j++) {
01270         real_data[j] = -real_data[j];
01271         }
01272         
01273         /* Compute a Hermitian --> real transform */
01274         csfft(1, n, real_data, comm, &info);
01275         if(info != 0) {
01276                 LOGERR("Error happening in Initialize communication work array: %d", info);
01277         }
01278         
01279         free(comm);
01280         return 0;
01281 }
01282 
01283 int EMfft::real_to_complex_2d(float *real_data, float *complex_data, int nx, int ny)
01284 {
01285         const int complex_nx = nx + 2 - nx%2;
01286         int info;
01287         /* Allocate communication work array */
01288         float * comm = (float *)malloc((3*nx+100)*sizeof(float));
01289         
01290         /* Allocate work array to store ACML complex array*/
01291         float * fft_data = (float *)malloc(complex_nx * ny * sizeof(float));
01292         cout << "fft_data after allocation:" << endl;
01293         for(int j=0; j<ny; ++j) {
01294                 for(int i=0; i<complex_nx; ++i) {
01295                         cout << "data[" << i << "][" << j << "] = " << fft_data[i+j*complex_nx] << "\t";
01296                 }
01297                 cout << endl;
01298         }
01299         
01300         /*copy real_data to complex_data then apply inplace FFT on complex data*/
01301         for(int i=0; i<ny; ++i) {
01302                 memcpy(fft_data+i*complex_nx, real_data+i*nx, nx*sizeof(float));
01303         }
01304         //memcpy(fft_data, real_data, nx * ny * sizeof(float));
01305         cout << endl << "real_data array: " << endl;
01306         for(int j=0; j<ny; ++j) {
01307                 for(int i=0; i<nx; ++i) {
01308                         cout << real_data[i+j*nx] << "\t";
01309                 }
01310                 cout << endl;
01311         }
01312         cout << endl << "the fft_data array: " << endl;
01313         for(int j=0; j<ny; ++j) {
01314                 for(int i=0; i<complex_nx; ++i) {
01315                         cout << fft_data[i+j*complex_nx] << "\t";
01316                 }
01317                 cout << endl;
01318         }
01319         
01320         //do a multiple 1d real-to-complex transform on x direction
01321         scfftm(ny, nx, fft_data, comm, &info);
01322         
01323         cout << endl << "the fft_data array after x dim transform: " << endl;
01324         for(int j=0; j<ny; ++j) {
01325                 for(int i=0; i<complex_nx; ++i) {
01326                         cout << fft_data[i+j*complex_nx] << "\t";
01327                 }
01328                 cout << endl;
01329         }
01330         
01331 /*      cout << "original fft array" << endl;
01332 /       cout << "complex_nx = " << complex_nx << " ny = " << n*(fft_data2+i+j*2*ny+1) << endl;
01333         for(int i=0; i<ny; ++i) {
01334                 for(int j=0; j<complex_nx; ++j) {
01335                         cout << *(fft_data+j+i*complex_nx) << "\t";
01336                 }
01337                 cout << endl;
01338         }
01339 */      
01340         //do a multiple 1d complex to complex transformation on y direction
01341         float * fft_data2 = (float *)malloc(complex_nx * ny * sizeof(float));
01342 /*      cout << "fft array rearranged in y dimension" << endl;
01343         for(int i=0; i<complex_nx; ++i) {
01344                 for(int j=0; j<ny; ++j) {
01345                         *(fft_data2+i+j*2*ny) = *(fft_data+i+j*complex_nx);
01346                         *(fft_data2+i+j*2*ny+1) = *(fft_data+i+complex_nx/2+j*complex_nx);
01347                         cout << *(fft_data2+i+j*2*ny) << "\t" << *(fft_data2+i+j*2*ny+1) << "\t";
01348                 }
01349                 cout << endl;
01350         }
01351 */      
01352         if(info != 0) {
01353                 LOGERR("Error happening in scfftm: %d", info);
01354         }
01355         
01356         return 0;       
01357 }
01358 
01359 int EMfft::complex_to_real_2d(float *complex_data, float *real_data, int nx, int ny)
01360 {
01361         return 0;
01362 }
01363 
01364 int EMfft::real_to_complex_3d(float *real_data, float *complex_data, int nx, int ny, int nz)
01365 {
01366         return 0;
01367 }
01368 
01369 int EMfft::complex_to_real_3d(float *complex_data, float *real_data, int nx, int ny, int nz)
01370 {
01371         return 0;
01372 }
01373 
01374 int EMfft::real_to_complex_nd(float *real_data, float *complex_data, int nx, int ny, int nz)
01375 {
01376         const int rank = get_rank(ny, nz);
01377         
01378         switch(rank) {
01379                 case 1:
01380                         return real_to_complex_1d(real_data, complex_data, nx);
01381                 case 2:
01382                         return real_to_complex_2d(real_data, complex_data, nx, ny);
01383                 case 3:
01384                         return real_to_complex_3d(real_data, complex_data, nx, ny, nz);         
01385                 default:
01386                         LOGERR("Never should be here...\n");
01387                         return -1;
01388         }       
01389 }
01390 
01391 int EMfft::complex_to_real_nd(float *complex_data, float *real_data, int nx, int ny, int nz)
01392 {
01393         const int rank = get_rank(ny, nz);
01394         
01395         switch(rank) {
01396                 case 1:
01397                         return complex_to_real_1d(complex_data, real_data, nx);
01398                 case 2:
01399                         return complex_to_real_2d(complex_data, real_data, nx, ny);
01400                 case 3:
01401                         return complex_to_real_3d(complex_data, real_data, nx, ny, nz);
01402                 default:
01403                         LOGERR("Never should be here...\n");
01404                         return -1;
01405         }
01406 }
01407 
01408 
01409 #endif  //ACML

Generated on Tue Jun 11 13:46:14 2013 for EMAN2 by  doxygen 1.3.9.1