Charm++: Charm Source Code Documentation

00001 #include "fftlib.h"
00002 
00003 /* 
00004  * This is on the "source" side.
00005  * there are srcSize[0] rows of size srcSize[1] each, in row-major order
00006  */
00007 void
00008 NormalRealSlabArray::doFFT(int src_id, int dst_id)
00009 {
00010     NormalFFTinfo &fftinfo = (infoVec[src_id]->info);
00011     
00012     if(fftinfo.transformType != REAL_TO_COMPLEX) CkPrintf("Transform Type at doFFT is %d\n", fftinfo.transformType);
00013     CkAssert(fftinfo.transformType == REAL_TO_COMPLEX);
00014     
00015     int rplaneSize = fftinfo.srcSize[0] * fftinfo.srcSize[1];
00016     int cplaneSize = (fftinfo.srcSize[0]/2+1) * fftinfo.srcSize[1];
00017     int lineSize = fftinfo.srcSize[1];
00018 
00019     double *dataPtr = (double*)fftinfo.dataPtr;
00020     complex *outData = new complex[cplaneSize * fftinfo.srcPlanesPerSlab];
00021     complex *outDataPtr = outData;
00022     complex *outData2 = new complex[cplaneSize * fftinfo.srcPlanesPerSlab];
00023     complex *outDataPtr2 = outData2;
00024 
00025         // do the 2D forward ffts
00026     CmiAssert(rfwd1DXPlan!=NULL && fwd1DYPlan!=NULL);
00027     int p;
00028 
00029     for(p = 0; p < fftinfo.srcPlanesPerSlab; p++){
00030         // real_to_complex on Y dimension
00031         rfftwnd_one_real_to_complex(rfwd2DXYPlan, (fftw_real*)dataPtr,
00032                                     (fftw_complex*)outDataPtr2);
00033 
00034 /*
00035 
00036         if(thisIndex==0 && src_id==0){
00037             for(int i=0;i<fftinfo.srcSize[1]*fftinfo.srcSize[0]; i++)
00038                 CkPrintf("%d %g\n", i, dataPtr[i]);
00039         }
00040 
00041         fftw(fwd1DYPlan, fftinfo.srcSize[0]/2+1,
00042              (fftw_complex*)outDataPtr, fftinfo.srcSize[1]/2+1, 1, 
00043              (fftw_complex*)outDataPtr2, fftinfo.srcSize[1]/2+1, 1);
00044 
00045 
00046         if(thisIndex==0 && src_id==0){
00047             for(int i=0;i<(fftinfo.srcSize[1]/2+1)*fftinfo.srcSize[0]; i++)
00048                 CkPrintf("%d %g %g\n", i, outDataPtr2[i].re, outDataPtr2[i].im);
00049         }
00050 */
00051         dataPtr += rplaneSize;
00052         outDataPtr += cplaneSize;
00053         outDataPtr2 += cplaneSize;
00054     }
00055     
00056     dataPtr = (double*)fftinfo.dataPtr;
00057     outDataPtr2 = outData2;
00058 
00059 
00060 //    CProxy_NormalSlabArray destProxy = (CProxy_NormalSlabArray)(infoVec[src_id]->destProxy);
00061 //    CProxy_NormalSlabArray srcProxy = (CProxy_NormalSlabArray)(infoVec[src_id]->srcProxy);
00062 
00063 #if FFTLIB_USE_COMLIB
00064     if (fftuseCommlib) {
00065         CProxy_NormalSlabArray destProxy_com;
00066         if(fftinfo.isSrcSlab)
00067             destProxy_com = (CProxy_NormalSlabArray)destProxy;
00068         else
00069             destProxy_com = (CProxy_NormalSlabArray)srcProxy;
00070         fftcommInstance.beginIteration();
00071         ComlibDelegateProxy(&destProxy_com);
00072     }
00073 #endif
00074 
00075     // allocating the data for sending to destination side
00076     lineSize = fftinfo.srcSize[1]/2+1;
00077     complex *sendData = new complex[fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize];
00078     complex *temp;
00079     // i <> pe if destPlanesPerSlab is not 1
00080     int pe, i, j;
00081     for(i = 0, pe = 0; i < fftinfo.srcSize[0]; i += fftinfo.destPlanesPerSlab, pe++) {
00082         int ti;
00083         temp = sendData;
00084 // Copy the data in (y-z-x) order
00085         for(ti = i; ti < i + fftinfo.destPlanesPerSlab; ti++) // x direction
00086             for(p = 0; p < fftinfo.srcPlanesPerSlab; p++) {   // z direction
00087                 memcpy(temp,
00088                        outDataPtr2 + p * cplaneSize + ti * lineSize,
00089                        sizeof(complex) * lineSize);
00090                 temp += lineSize;
00091             }
00092 #ifdef VERBOSE
00093         CkPrintf("[%d] sendFFTData to %d, size %d \n", thisIndex, pe, 
00094                  lineSize * fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab);
00095 #endif
00096         ((CProxy_NormalSlabArray)destProxy)(pe).acceptDataForFFT(lineSize * fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab, sendData, thisIndex, dst_id);
00097     }
00098     delete [] sendData;
00099     delete [] outData;
00100     delete [] outData2;
00101 }
00102 
00103 /*
00104  * This is on the "destination" side.
00105  */
00106 void
00107 NormalRealSlabArray::acceptDataForFFT(int numPoints, complex *points, int posn, int info_id)
00108 {
00109     NormalFFTinfo &fftinfo = (infoVec[info_id]->info);
00110     CkAssert(fftinfo.transformType == COMPLEX_TO_REAL);
00111 
00112     complex *dataPtr = (complex*)fftinfo.dataPtr;
00113     int lineSize = fftinfo.destSize[1]/2+1;
00114     
00115 #if CAREFUL
00116     CkAssert(numPoints == fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize);
00117 #endif
00118         
00119     infoVec[info_id]->count++;
00120     int planeSize = fftinfo.destSize[0] * (fftinfo.destSize[1]/2+1);
00121     int p;
00122     for (p = 0; p < fftinfo.destPlanesPerSlab; p++) {
00123         memcpy(dataPtr + posn * fftinfo.srcPlanesPerSlab * lineSize + p * planeSize,
00124                points, 
00125                sizeof(complex) * lineSize * fftinfo.srcPlanesPerSlab);
00126         points += lineSize * fftinfo.srcPlanesPerSlab;
00127     }
00128 
00129     if (infoVec[info_id]->count == fftinfo.destSize[0] / fftinfo.srcPlanesPerSlab) {
00130         infoVec[info_id]->count = 0;
00131         CkAssert(fwd1DZPlan != NULL);
00132         for(p = 0; p < fftinfo.destPlanesPerSlab; p++) {
00133                 fftw(fwd1DZPlan, 
00134                      lineSize,
00135                      (fftw_complex*)dataPtr + p * planeSize,
00136                      lineSize, 1, //stride, nextFFT
00137                      NULL, 0, 0);
00138         }
00139         doneFFT(info_id);
00140     }
00141 }
00142 
00143 /*
00144  * This is on the "destination" side.
00145  *    There are destSize[0] rows of destSize[1] each, in row major order, and srcSize[1] == destSize[1]
00146  */
00147 void
00148 NormalRealSlabArray::doIFFT(int src_id, int dst_id)
00149 {
00150     NormalFFTinfo &fftinfo = (infoVec[src_id]->info);
00151     CkAssert(fftinfo.transformType == COMPLEX_TO_REAL);
00152 
00153     complex *dataPtr = (complex*)fftinfo.dataPtr;
00154     int planeSize = fftinfo.destSize[0] * (fftinfo.destSize[1]/2+1);
00155     int lineSize = fftinfo.destSize[1]/2+1;
00156     
00157     CmiAssert(bwd1DZPlan!=NULL && bwd1DYPlan!=NULL);
00158     int p;
00159 
00160     for(p = 0; p < fftinfo.destPlanesPerSlab; p++){
00161         // Z direction FFT
00162         fftw(bwd1DZPlan, lineSize, 
00163              (fftw_complex*)dataPtr+p*planeSize, lineSize, 1, 
00164              NULL, 0, 0);
00165 
00166 /*      if(thisIndex==0 && src_id==0){
00167             for(int i=0;i<(fftinfo.srcSize[1]/2+1)*fftinfo.srcSize[0]; i++)
00168                 CkPrintf("%d %g %g\n", i, dataPtr[i+p*planeSize].re, dataPtr[i+p*planeSize].im);
00169         }
00170 */
00171     }
00172       
00173 //    CProxy_NormalSlabArray destProxy = (CProxy_NormalSlabArray)(infoVec[src_id]->destProxy);
00174 //    CProxy_NormalSlabArray srcProxy = (CProxy_NormalSlabArray)(infoVec[src_id]->srcProxy);
00175 
00176 #if FFTLIB_USE_COMLIB
00177     if (fftuseCommlib) {
00178         CProxy_NormalSlabArray destProxy_com;
00179         if(fftinfo.isSrcSlab)
00180             destProxy_com = (CProxy_NormalSlabArray)destProxy;
00181         else
00182             destProxy_com = (CProxy_NormalSlabArray)srcProxy;
00183         fftcommInstance.beginIteration();
00184         ComlibDelegateProxy(&destProxy_com);
00185     }
00186 #endif
00187 
00188     complex *sendData = new complex[fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize];
00189     complex *temp;
00190     int i, pe;
00191     for (i = 0, pe = 0; i < fftinfo.destSize[0]; i += fftinfo.srcPlanesPerSlab, pe++) {
00192         int ti;
00193         temp = sendData;
00194         // Sending in (y-x-z) order
00195         for (ti = i; ti < i + fftinfo.srcPlanesPerSlab; ti++)
00196             for (p = 0; p < fftinfo.destPlanesPerSlab; p++) {
00197                 memcpy(temp,
00198                        dataPtr + p * planeSize + ti * lineSize,
00199                        sizeof(complex) * lineSize);
00200                 temp += lineSize;
00201             }
00202 
00203         ((CProxy_NormalSlabArray)srcProxy)(pe).acceptDataForIFFT(lineSize * fftinfo.destPlanesPerSlab * fftinfo.srcPlanesPerSlab, sendData, thisIndex, dst_id);
00204     }
00205     delete [] sendData;
00206 }
00207 
00208 void
00209 NormalRealSlabArray::acceptDataForIFFT(int numPoints, complex *points, int posn, int info_id)
00210 {
00211     NormalFFTinfo &fftinfo = (infoVec[info_id]->info);
00212     CkAssert(fftinfo.transformType == REAL_TO_COMPLEX);
00213 
00214     int rplaneSize = fftinfo.srcSize[0] * fftinfo.srcSize[1];
00215     int cplaneSize = fftinfo.srcSize[0] * (fftinfo.srcSize[1]/2+1);
00216     int planeSize = fftinfo.destSize[0] * (fftinfo.destSize[1]/2+1);
00217     int lineSize = fftinfo.destSize[1]/2+1;
00218     if(tempdataPtr==NULL)
00219         tempdataPtr = new complex[fftinfo.srcPlanesPerSlab * cplaneSize];
00220     complex *inDataPtr = tempdataPtr;
00221 #if CAREFUL
00222     CkAssert(numPoints == fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize);
00223 #endif
00224 
00225     infoVec[info_id]->count++;
00226     int p;
00227     complex *pointPtr = points;
00228     inDataPtr = tempdataPtr + posn * lineSize * fftinfo.destPlanesPerSlab;
00229     for(p = 0; p < fftinfo.srcPlanesPerSlab; p++) {
00230         memcpy(inDataPtr, pointPtr, lineSize* fftinfo.destPlanesPerSlab*sizeof(complex));
00231         inDataPtr += planeSize;
00232         pointPtr += lineSize*fftinfo.destPlanesPerSlab;
00233     }
00234     
00235     int expectedCount = fftinfo.srcSize[0]/fftinfo.destPlanesPerSlab;
00236     if (infoVec[info_id]->count == expectedCount) {
00237         infoVec[info_id]->count = 0;
00238         CmiAssert(rbwd1DXPlan!=NULL);
00239         double *dataPtr = (double*)fftinfo.dataPtr;
00240         for(p = 0; p < fftinfo.srcPlanesPerSlab; p++) {
00241             rfftwnd_one_complex_to_real(rbwd2DXYPlan, 
00242                                         (fftw_complex*)(tempdataPtr+p*cplaneSize), 
00243                                         (fftw_real*)(dataPtr+p*rplaneSize));
00244 /*
00245             if(thisIndex==0 && info_id==0){
00246                 for(int i=0;i<fftinfo.srcSize[1]*fftinfo.srcSize[0]; i++)
00247                     CkPrintf("%d %g\n", i, dataPtr[i]);
00248             }
00249             
00250 
00251             // Doing X direction first
00252                 fftw(bwd1DYPlan, lineSize, 
00253                      (fftw_complex*)tempdataPtr+p*planeSize, lineSize, 1,
00254                      NULL, 0, 0);
00255                 
00256             // complex_to_real transform at Y direction
00257                 rfftw(rbwd1DXPlan,
00258                       fftinfo.srcSize[0],
00259                       (fftw_real*)(tempdataPtr + p * cplaneSize),
00260                       1, lineSize*2, //stride, nextFFT
00261                       (fftw_real*)dataPtr + p * rplaneSize, 
00262                       1, fftinfo.srcSize[1]);
00263 */
00264         }
00265         doneIFFT(info_id);
00266     }
00267 }
00268 void NormalRealSlabArray::createPlans(NormalFFTinfo &info)
00269 {
00270     if (info.isSrcSlab) {
00271         int size[]={info.srcSize[0], info.srcSize[1]};
00272         rfwd2DXYPlan = rfftw2d_create_plan(info.srcSize[0], info.srcSize[1], FFTW_REAL_TO_COMPLEX, FFTW_OUT_OF_PLACE);
00273         rbwd2DXYPlan = rfftw2d_create_plan(info.srcSize[0], info.srcSize[1], FFTW_COMPLEX_TO_REAL, FFTW_OUT_OF_PLACE);
00274 
00275         rfwd1DXPlan = rfftw_create_plan(info.srcSize[0], FFTW_REAL_TO_COMPLEX, FFTW_OUT_OF_PLACE);
00276         //fwd1DYPlan = fftw_create_plan(info.srcSize[1], FFTW_BACKWARD, FFTW_MEASURE|FFTW_IN_PLACE); 
00277         fwd1DYPlan = fftw_create_plan(info.srcSize[1], FFTW_BACKWARD, FFTW_OUT_OF_PLACE); 
00278         rbwd1DXPlan = rfftw_create_plan(info.srcSize[0], FFTW_COMPLEX_TO_REAL, FFTW_OUT_OF_PLACE);
00279         bwd1DYPlan = fftw_create_plan(info.destSize[1], FFTW_BACKWARD, FFTW_IN_PLACE);
00280     }
00281     else {
00282 
00283         bwd1DZPlan = fftw_create_plan(info.destSize[0], FFTW_BACKWARD, FFTW_IN_PLACE);
00284         bwd1DYPlan = fftw_create_plan(info.destSize[1], FFTW_BACKWARD, FFTW_IN_PLACE);
00285         fwd1DZPlan = fftw_create_plan(info.destSize[0], FFTW_FORWARD, FFTW_IN_PLACE);
00286     }
00287 }
00288 
00289 void NormalRealSlabArray::setup(NormalFFTinfo &info,
00290                                 CProxy_NormalRealSlabArray src, 
00291                                 CProxy_NormalRealSlabArray dest, 
00292                                 bool useCommlib, ComlibInstanceHandle inst)
00293 {
00294     SlabArrayInfo *slabinfo=new SlabArrayInfo;
00295     slabinfo->info = info;
00296     slabinfo->count = 0;
00297     infoVec.insert(infoVec.size(), slabinfo);
00298 
00299     srcProxy = src;
00300     destProxy = dest;
00301     rfwd1DXPlan = rbwd1DXPlan = (rfftw_plan) NULL;
00302     fwd1DYPlan = bwd1DYPlan = (fftw_plan) NULL;
00303     fwd1DZPlan = bwd1DZPlan = (fftw_plan) NULL;
00304 
00305     createPlans(info);
00306 
00307 #if FFTLIB_USE_COMLIB
00308     fftuseCommlib = useCommlib;
00309     fftcommInstance = ComlibInstanceHandle();
00310     if (fftuseCommlib) {        
00311         fftcommInstance = inst;
00312    }
00313 #else 
00314     fftuseCommlib = false;
00315 #endif
00316 
00317 }
00318 
00319 NormalRealSlabArray::NormalRealSlabArray(NormalFFTinfo &info,
00320                                          CProxy_NormalRealSlabArray src, 
00321                                          CProxy_NormalRealSlabArray dest, 
00322                                          bool useCommlib, ComlibInstanceHandle inst)
00323 { 
00324     setup(info, src, dest, useCommlib, inst); 
00325 }
00326 
00327 NormalRealSlabArray::~NormalRealSlabArray() 
00328 {
00329     if (rfwd1DXPlan)
00330         rfftw_destroy_plan(rfwd1DXPlan);
00331     if (rbwd1DXPlan)
00332         rfftw_destroy_plan(rbwd1DXPlan);
00333     if (fwd1DYPlan)
00334         fftw_destroy_plan(fwd1DYPlan);
00335     if (bwd1DYPlan)
00336         fftw_destroy_plan(bwd1DYPlan);
00337     if (fwd1DZPlan)
00338         fftw_destroy_plan(fwd1DZPlan);
00339     if (bwd1DZPlan)
00340         fftw_destroy_plan(bwd1DZPlan);
00341 
00342     if(rfwd2DXYPlan)
00343         fftwnd_destroy_plan(rfwd2DXYPlan);
00344     if(rbwd2DXYPlan)
00345         fftwnd_destroy_plan(rbwd2DXYPlan);
00346 
00347     infoVec.removeAll();
00348     if(tempdataPtr != NULL) 
00349         delete [] tempdataPtr;
00350 }
00351 
00352 
00353 void NormalRealSlabArray::pup(PUP::er &p)
00354 {
00355     int i;
00356 
00357 /*    for (i = 0; i < MAX_FFTS; i++) {
00358         if (p.isUnpacking()) {
00359           fftinfos[i] = NULL;
00360         }
00361         int val = fftinfos[i]? 1:0;
00362         p | val; 
00363         if (val) {                // if this entry is not NULL
00364             if (p.isUnpacking()){
00365                 fftinfos[i] = new NormalFFTinfo;
00366             }  
00367             fftinfos[i]->pup(p);
00368             
00369             if (p.isUnpacking()){
00370                 createPlans(*fftinfos[i]);
00371             }       
00372         }
00373     }
00374 */  
00375 }
00376 
00377
libs/ck-libs/fftlib/rfftlib-normal.C
