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     // allocating the data for sending to destination side
00061     lineSize = fftinfo.srcSize[1]/2+1;
00062     complex *sendData = new complex[fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize];
00063     complex *temp;
00064     // i <> pe if destPlanesPerSlab is not 1
00065     int pe, i, j;
00066     for(i = 0, pe = 0; i < fftinfo.srcSize[0]; i += fftinfo.destPlanesPerSlab, pe++) {
00067         int ti;
00068         temp = sendData;
00069 // Copy the data in (y-z-x) order
00070         for(ti = i; ti < i + fftinfo.destPlanesPerSlab; ti++) // x direction
00071             for(p = 0; p < fftinfo.srcPlanesPerSlab; p++) {   // z direction
00072                 memcpy(temp,
00073                        outDataPtr2 + p * cplaneSize + ti * lineSize,
00074                        sizeof(complex) * lineSize);
00075                 temp += lineSize;
00076             }
00077 #ifdef VERBOSE
00078         CkPrintf("[%d] sendFFTData to %d, size %d \n", thisIndex, pe, 
00079          lineSize * fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab);
00080 #endif
00081         ((CProxy_NormalSlabArray)destProxy)(pe).acceptDataForFFT(lineSize * fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab, sendData, thisIndex, dst_id);
00082     }
00083     delete [] sendData;
00084     delete [] outData;
00085     delete [] outData2;
00086 }
00087 
00088 /*
00089  * This is on the "destination" side.
00090  */
00091 void
00092 NormalRealSlabArray::acceptDataForFFT(int numPoints, complex *points, int posn, int info_id)
00093 {
00094     NormalFFTinfo &fftinfo = (infoVec[info_id]->info);
00095     CkAssert(fftinfo.transformType == COMPLEX_TO_REAL);
00096 
00097     complex *dataPtr = (complex*)fftinfo.dataPtr;
00098     int lineSize = fftinfo.destSize[1]/2+1;
00099     
00100 #if CAREFUL
00101     CkAssert(numPoints == fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize);
00102 #endif
00103     
00104     infoVec[info_id]->count++;
00105     int planeSize = fftinfo.destSize[0] * (fftinfo.destSize[1]/2+1);
00106     int p;
00107     for (p = 0; p < fftinfo.destPlanesPerSlab; p++) {
00108     memcpy(dataPtr + posn * fftinfo.srcPlanesPerSlab * lineSize + p * planeSize,
00109            points, 
00110            sizeof(complex) * lineSize * fftinfo.srcPlanesPerSlab);
00111     points += lineSize * fftinfo.srcPlanesPerSlab;
00112     }
00113 
00114     if (infoVec[info_id]->count == fftinfo.destSize[0] / fftinfo.srcPlanesPerSlab) {
00115     infoVec[info_id]->count = 0;
00116     CkAssert(fwd1DZPlan != NULL);
00117     for(p = 0; p < fftinfo.destPlanesPerSlab; p++) {
00118         fftw(fwd1DZPlan, 
00119              lineSize,
00120              (fftw_complex*)dataPtr + p * planeSize,
00121              lineSize, 1, //stride, nextFFT
00122              NULL, 0, 0);
00123     }
00124     doneFFT(info_id);
00125     }
00126 }
00127 
00128 /*
00129  * This is on the "destination" side.
00130  *    There are destSize[0] rows of destSize[1] each, in row major order, and srcSize[1] == destSize[1]
00131  */
00132 void
00133 NormalRealSlabArray::doIFFT(int src_id, int dst_id)
00134 {
00135     NormalFFTinfo &fftinfo = (infoVec[src_id]->info);
00136     CkAssert(fftinfo.transformType == COMPLEX_TO_REAL);
00137 
00138     complex *dataPtr = (complex*)fftinfo.dataPtr;
00139     int planeSize = fftinfo.destSize[0] * (fftinfo.destSize[1]/2+1);
00140     int lineSize = fftinfo.destSize[1]/2+1;
00141     
00142     CmiAssert(bwd1DZPlan!=NULL && bwd1DYPlan!=NULL);
00143     int p;
00144 
00145     for(p = 0; p < fftinfo.destPlanesPerSlab; p++){
00146     // Z direction FFT
00147     fftw(bwd1DZPlan, lineSize, 
00148          (fftw_complex*)dataPtr+p*planeSize, lineSize, 1, 
00149          NULL, 0, 0);
00150 
00151 /*  if(thisIndex==0 && src_id==0){
00152         for(int i=0;i<(fftinfo.srcSize[1]/2+1)*fftinfo.srcSize[0]; i++)
00153         CkPrintf("%d %g %g\n", i, dataPtr[i+p*planeSize].re, dataPtr[i+p*planeSize].im);
00154     }
00155 */
00156     }
00157       
00158     complex *sendData = new complex[fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize];
00159     complex *temp;
00160     int i, pe;
00161     for (i = 0, pe = 0; i < fftinfo.destSize[0]; i += fftinfo.srcPlanesPerSlab, pe++) {
00162     int ti;
00163     temp = sendData;
00164     // Sending in (y-x-z) order
00165     for (ti = i; ti < i + fftinfo.srcPlanesPerSlab; ti++)
00166         for (p = 0; p < fftinfo.destPlanesPerSlab; p++) {
00167         memcpy(temp,
00168                dataPtr + p * planeSize + ti * lineSize,
00169                sizeof(complex) * lineSize);
00170         temp += lineSize;
00171         }
00172 
00173     ((CProxy_NormalSlabArray)srcProxy)(pe).acceptDataForIFFT(lineSize * fftinfo.destPlanesPerSlab * fftinfo.srcPlanesPerSlab, sendData, thisIndex, dst_id);
00174     }
00175     delete [] sendData;
00176 }
00177 
00178 void
00179 NormalRealSlabArray::acceptDataForIFFT(int numPoints, complex *points, int posn, int info_id)
00180 {
00181     NormalFFTinfo &fftinfo = (infoVec[info_id]->info);
00182     CkAssert(fftinfo.transformType == REAL_TO_COMPLEX);
00183 
00184     int rplaneSize = fftinfo.srcSize[0] * fftinfo.srcSize[1];
00185     int cplaneSize = fftinfo.srcSize[0] * (fftinfo.srcSize[1]/2+1);
00186     int planeSize = fftinfo.destSize[0] * (fftinfo.destSize[1]/2+1);
00187     int lineSize = fftinfo.destSize[1]/2+1;
00188     if(tempdataPtr==NULL)
00189     tempdataPtr = new complex[fftinfo.srcPlanesPerSlab * cplaneSize];
00190     complex *inDataPtr = tempdataPtr;
00191 #if CAREFUL
00192     CkAssert(numPoints == fftinfo.srcPlanesPerSlab * fftinfo.destPlanesPerSlab * lineSize);
00193 #endif
00194 
00195     infoVec[info_id]->count++;
00196     int p;
00197     complex *pointPtr = points;
00198     inDataPtr = tempdataPtr + posn * lineSize * fftinfo.destPlanesPerSlab;
00199     for(p = 0; p < fftinfo.srcPlanesPerSlab; p++) {
00200     memcpy(inDataPtr, pointPtr, lineSize* fftinfo.destPlanesPerSlab*sizeof(complex));
00201     inDataPtr += planeSize;
00202     pointPtr += lineSize*fftinfo.destPlanesPerSlab;
00203     }
00204     
00205     int expectedCount = fftinfo.srcSize[0]/fftinfo.destPlanesPerSlab;
00206     if (infoVec[info_id]->count == expectedCount) {
00207     infoVec[info_id]->count = 0;
00208     CmiAssert(rbwd1DXPlan!=NULL);
00209     double *dataPtr = (double*)fftinfo.dataPtr;
00210     for(p = 0; p < fftinfo.srcPlanesPerSlab; p++) {
00211         rfftwnd_one_complex_to_real(rbwd2DXYPlan, 
00212                     (fftw_complex*)(tempdataPtr+p*cplaneSize), 
00213                     (fftw_real*)(dataPtr+p*rplaneSize));
00214 /*
00215         if(thisIndex==0 && info_id==0){
00216         for(int i=0;i<fftinfo.srcSize[1]*fftinfo.srcSize[0]; i++)
00217             CkPrintf("%d %g\n", i, dataPtr[i]);
00218         }
00219         
00220 
00221         // Doing X direction first
00222         fftw(bwd1DYPlan, lineSize, 
00223              (fftw_complex*)tempdataPtr+p*planeSize, lineSize, 1,
00224              NULL, 0, 0);
00225         
00226         // complex_to_real transform at Y direction
00227         rfftw(rbwd1DXPlan,
00228               fftinfo.srcSize[0],
00229               (fftw_real*)(tempdataPtr + p * cplaneSize),
00230               1, lineSize*2, //stride, nextFFT
00231               (fftw_real*)dataPtr + p * rplaneSize, 
00232               1, fftinfo.srcSize[1]);
00233 */
00234     }
00235     doneIFFT(info_id);
00236     }
00237 }
00238 void NormalRealSlabArray::createPlans(NormalFFTinfo &info)
00239 {
00240     if (info.isSrcSlab) {
00241     int size[]={info.srcSize[0], info.srcSize[1]};
00242     rfwd2DXYPlan = rfftw2d_create_plan(info.srcSize[0], info.srcSize[1], FFTW_REAL_TO_COMPLEX, FFTW_OUT_OF_PLACE);
00243     rbwd2DXYPlan = rfftw2d_create_plan(info.srcSize[0], info.srcSize[1], FFTW_COMPLEX_TO_REAL, FFTW_OUT_OF_PLACE);
00244 
00245     rfwd1DXPlan = rfftw_create_plan(info.srcSize[0], FFTW_REAL_TO_COMPLEX, FFTW_OUT_OF_PLACE);
00246     //fwd1DYPlan = fftw_create_plan(info.srcSize[1], FFTW_BACKWARD, FFTW_MEASURE|FFTW_IN_PLACE); 
00247     fwd1DYPlan = fftw_create_plan(info.srcSize[1], FFTW_BACKWARD, FFTW_OUT_OF_PLACE); 
00248     rbwd1DXPlan = rfftw_create_plan(info.srcSize[0], FFTW_COMPLEX_TO_REAL, FFTW_OUT_OF_PLACE);
00249     bwd1DYPlan = fftw_create_plan(info.destSize[1], FFTW_BACKWARD, FFTW_IN_PLACE);
00250     }
00251     else {
00252 
00253     bwd1DZPlan = fftw_create_plan(info.destSize[0], FFTW_BACKWARD, FFTW_IN_PLACE);
00254     bwd1DYPlan = fftw_create_plan(info.destSize[1], FFTW_BACKWARD, FFTW_IN_PLACE);
00255     fwd1DZPlan = fftw_create_plan(info.destSize[0], FFTW_FORWARD, FFTW_IN_PLACE);
00256     }
00257 }
00258 
00259 void NormalRealSlabArray::setup(NormalFFTinfo &info,
00260                 CProxy_NormalRealSlabArray src, 
00261                 CProxy_NormalRealSlabArray dest)
00262 {
00263     SlabArrayInfo *slabinfo=new SlabArrayInfo;
00264     slabinfo->info = info;
00265     slabinfo->count = 0;
00266     infoVec.insert(infoVec.size(), slabinfo);
00267 
00268     srcProxy = src;
00269     destProxy = dest;
00270     rfwd1DXPlan = rbwd1DXPlan = (rfftw_plan) NULL;
00271     fwd1DYPlan = bwd1DYPlan = (fftw_plan) NULL;
00272     fwd1DZPlan = bwd1DZPlan = (fftw_plan) NULL;
00273 
00274     createPlans(info);
00275 }
00276 
00277 NormalRealSlabArray::NormalRealSlabArray(NormalFFTinfo &info,
00278                      CProxy_NormalRealSlabArray src, 
00279                      CProxy_NormalRealSlabArray dest)
00280 { 
00281     setup(info, src, dest);
00282 }
00283 
00284 NormalRealSlabArray::~NormalRealSlabArray() 
00285 {
00286     if (rfwd1DXPlan)
00287     rfftw_destroy_plan(rfwd1DXPlan);
00288     if (rbwd1DXPlan)
00289     rfftw_destroy_plan(rbwd1DXPlan);
00290     if (fwd1DYPlan)
00291     fftw_destroy_plan(fwd1DYPlan);
00292     if (bwd1DYPlan)
00293     fftw_destroy_plan(bwd1DYPlan);
00294     if (fwd1DZPlan)
00295     fftw_destroy_plan(fwd1DZPlan);
00296     if (bwd1DZPlan)
00297     fftw_destroy_plan(bwd1DZPlan);
00298 
00299     if(rfwd2DXYPlan)
00300     fftwnd_destroy_plan(rfwd2DXYPlan);
00301     if(rbwd2DXYPlan)
00302     fftwnd_destroy_plan(rbwd2DXYPlan);
00303 
00304     infoVec.removeAll();
00305     if(tempdataPtr != NULL) 
00306     delete [] tempdataPtr;
00307 }
00308 
00309 
00310 void NormalRealSlabArray::pup(PUP::er &p)
00311 {
00312     int i;
00313 
00314 /*    for (i = 0; i < MAX_FFTS; i++) {
00315     if (p.isUnpacking()) {
00316           fftinfos[i] = NULL;
00317     }
00318     int val = fftinfos[i]? 1:0;
00319     p | val; 
00320     if (val) {                // if this entry is not NULL
00321         if (p.isUnpacking()){
00322         fftinfos[i] = new NormalFFTinfo;
00323         }  
00324         fftinfos[i]->pup(p);
00325         
00326             if (p.isUnpacking()){
00327         createPlans(*fftinfos[i]);
00328         }       
00329     }
00330     }
00331 */  
00332 }
00333 
00334
libs/ck-libs/fftlib/rfftlib-normal.C
