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

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

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