Charm++: Charm Source Code Documentation

00001 #include "fftlib.h"
00002 
00003 /*
00004  * organization of data: 
00005  * On the "source" side -
00006  * there are srcSize[0] rows of size srcSize[1] each. 
00007  */
00008 #if 0
00009 void
00010 SlabArray::doFFT()
00011 {
00012     // pointer to region of data; this may be in run-length form.
00013     complex *srcData = data();
00014     FFTinfo *infoPtr = getFFTinfo();
00015     const CkIndex2D *idx = myIdx();
00016     complex *dataPtr = 0;
00017     if (info->packingMode == SPARSE_MODE) {
00018         // get the run-length information
00019         RunDescriptor *runs;
00020         int numRuns, numPoints;
00021         getRuns(&runs, &numRuns, &numPoints);
00022         dataPtr = new complex[infoPtr->srcSize[0] * infoPtr->srcSize[1] * infoPtr->srcPlanesPerSlab];
00023         int i, l, numCovered = 0;
00024         int planeSize = info->srcSize[0] * info->srcSize[1];
00025         bool *nonZero = new int[infoPtr->srcSize[1] * infoPtr->srcPlanesPerSlab];
00026         int numNonZero = 0;
00027         for (i = 0; i < infoPtr->srcSize[1] * infoPtr->srcPlanesPerSlab; i++)
00028             nonZero[i] = false;
00029         for (r = 0; r < numRuns; r++) {
00030             // the starting position of the run in the data array (dataPtr)
00031             int pos = runs[r].x * planeSize + runs[r].y * info->srcSize[1];
00032 
00033             // copying the data into the actual array where the fft will be
00034             // done.
00035             for (l = 0; l < runs[r].length; l++) {
00036                 dataPtr[pos + runs[r].offset(l)] = srcData[numCovered++];
00037                 nonZero[runs[r].x * infoPtr->srcSize[1] + runs[r].offset(l)] = true;
00038                 numNonZero++;
00039             }
00040         }
00041         // Here we do only a one-dimensional fft.
00042         int p;
00043         for(p = 0; p < srcPlanesPerSlab; p++)
00044             for (z = 0; z < info->srcSize[1]; z++)
00045                 if (nonZero[p * info->srcSize[1] + z])
00046                     fftw_one(fwd1DPlan, // the forward plan
00047                              1,
00048                              dataPtr + p * infoPtr->srcSize[0] * infoPtr->srcSize[1] + z, // the pointer to the data
00049                              infoPtr->srcSize[0], 1,
00050                              NULL, 0, 0);
00051 
00052         // Now, to send the data: we send only the non-zero parts.
00053         // This is the advantage of doing the 1D ffts on the source side.
00054         // this is done by augmenting the data with an array of integer flags
00055         // The array does not change within the context of a single plane
00056         complex *sendData = new complex[numNonZero * infoPtr->destPlanesPerSlab];
00057         int z, y, numCovered;
00058         for (y = 0; y < infoPtr->srcSize[0]; y += infoPtr->destPlanesPerSlab) {
00059             int t;
00060             numCovered = 0;
00061             for (t = y; t < y + infoPtr->destPlanesPerSlab; t++) 
00062                 for (p = 0; p < infoPtr->srcPlanesPerSlab; p++) 
00063                     for (z = 0; z < infoPtr->srcSize[1]; z++)
00064                         if (nonZero[p * infoPtr->srcSize[1] + z]) {
00065                             sendData[numCovered] = dataPtr[p * infoPtr->srcSize[0] * infoPtr->srcSize[1] + t * infoPtr->srcSize[1] + z];
00066                             numCovered++;
00067                         }
00068             
00069             destProxy(idx->x, y).acceptDataForFFT(infoPtr->srcSize[1] * infoPtr->srcPlanesPerSlab, numNonZero * infoPtr->destPlanesPerSlab, sendData, idx->y);
00070         }
00071 
00072         delete [] nonZero;
00073         delete [] sendData;
00074         return;
00075     }
00076     
00077     
00078     // The following part handles the case where packingMode equals DENSE_MODE
00079     
00080     // Do the fft in the planes
00081     int p;
00082     for(p = 0; p < srcPlanesPerSlab; p++)
00083     fftwnd_one(fwd2DPlan, 
00084            dataPtr + p * infoPtr->srcSize[0] * infoPtr->srcSize[1],
00085            NULL);
00086     
00087     // allocating the data for sending to destination side
00088     complex *sendData = new complex[srcPlanesPerSlab * destPlanesPerSlab * infoPtr->srcSize[1]];
00089     complex *temp;
00090     int i;
00091     for(i = 0; i < infoPtr->srcSize[0]; i += destPlanesPerSlab) {
00092         int ti;
00093         temp = sendData;
00094         for(ti = i; ti < i + destPlanesPerSlab; ti++)
00095             for(p = 0; p < srcPlanesPerSlab; p++) {
00096                 memcpy(temp, 
00097                        dataPtr + p * infoPtr->srcSize[0] * infoPtr->srcSize[1] + ti * infoPtr->srcSize[1],
00098                        sizeof(complex) * infoPtr->srcSize[1]);
00099                 temp += infoPtr->srcSize[1];
00100             }
00101     
00102         destProxy(idx->x, i).acceptDataForFFT(0, 0, infoPtr->srcSize[1] * srcPlanesPerSlab * destPlanesPerSlab, sendData, idx->y);
00103     }
00104     delete [] sendData;
00105     delete [] dataPtr;
00106 }
00107 
00108 void
00109 SlabArray::acceptDataForFFT(int n, int *flags, int dataSize, complex *data, int posn)
00110 {
00111     complex *dataPtr = this->data();
00112     FFTinfo *infoPtr = getFFTinfo();    
00113 
00114     count++;
00115     if (infoPtr->packingMode == SPARSE_MODE) {
00116         CkAssert (n && flags);
00117 
00118         //first, copy the flags in the right place
00119         memcpy (nonZero + posn * infoPtr->destSize[], flags, sizeof(int) * n);
00120     }
00121     else {
00122         int p;
00123         for (p = 0; p < infoPtr->destPlanesPerSlab; p++) {
00124             memcpy(dataPtr + posn * infoPtr->destSize[1] + p * infoPtr->destSize[0] * infoPtr->destSize[1],
00125                    data, 
00126                    sizeof(complex) * infoPtr->srcPlanesPerSlab * infoPtr->destSize[1]);
00127             data += destSize2 * srcPlanesPerSlab;
00128         }
00129         if (count == info->destSize[0] / srcPlanesPerSlab) {
00130             count = 0;
00131             for(p = 0; p < destPlanesPerSlab; p++)
00132                 fftw(fwd1DPlan, 
00133                      destSize2,
00134                      dataPtr + p * destSize1 * destSize2,
00135                      destSize2, 1,
00136                      NULL, 0, 0);
00137     
00138             doneFFT();
00139         }
00140     }
00141 }
00142 
00143 void
00144 SlabArray::doIFFT()
00145 {
00146     complex *dataPtr = data();
00147     const CkIndex2D *idx = myIdx();
00148     
00149     int p;
00150     for(p = 0; p < destPlanesPerSlab; p++)
00151     fftwnd_one(bwd2DPlan, 
00152            dataPtr + p * destSize1 * destSize2,
00153            NULL);
00154     
00155     complex *sendData = new complex[srcPlanesPerSlab * destPlanesPerSlab * destSize2];
00156     complex *temp;
00157     int i;
00158     for (i = 0; i < destSize1; i += srcPlanesPerSlab) {
00159     int ti;
00160     temp = sendData;
00161     for (ti = i; ti < i + srcPlanesPerSlab; ti++)
00162         for (p = 0; p < destPlanesPerSlab; p++) {
00163         memcpy(temp,
00164                dataPtr + p * destSize1 * destSize2 + ti * destSize2,
00165                sizeof(complex) * destSize2);
00166         temp += destSize2;
00167         }
00168     srcProxy(idx->x, i).acceptDataForIFFT(destSize2 * destPlanesPerSlab * srcPlanesPerSlab, sendData, idx->y);
00169     }
00170     delete [] sendData;
00171 }
00172 
00173 void
00174 SlabArray::acceptDataForIFFT(int dataSize, complex *data, int posn)
00175 {
00176     complex *dataPtr = this->data();
00177     
00178     count++;
00179     int p;
00180     for(p = 0; p < srcPlanesPerSlab; p++) {
00181         memcpy(dataPtr + p * infoPtr->srcSize[0] * infoPtr->srcSize[1] + posn * infoPtr->srcSize[1],
00182                data , 
00183                sizeof(complex) * infoPtr->srcSize[1] * destPlanesPerSlab);
00184         data += destPlanesPerSlab * infoPtr->srcSize[1];
00185     }
00186     
00187     if (count == infoPtr->srcSize[0] / destPlanesPerSlab) {
00188         count = 0;
00189         for(p = 0; p < srcPlanesPerSlab; p++)
00190             fftw(bwd1DPlan,
00191                  infoPtr->srcSize[1],
00192                  dataPtr + p * infoPtr->srcSize[0] * infoPtr->srcSize[1],
00193                  infoPtr->srcSize[1], 1,
00194                  NULL, 0, 0);
00195         doneIFFT();
00196     }
00197 }
00198 
00199 SlabArray::SlabArray(FFTinfo &info) 
00200 {
00201     count = 0;
00202     if (info.isSrcSlab) {
00203         if (info.packingMode == SPARSE_MODE) {
00204             fwd1DPlan = fftw_create_plan(info.srcSize[0], FFTW_FORWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00205             bwd1DPlan = fftw_create_plan(info.srcSize[0], FFTW_BACKWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00206             nonZero = new int[info.destSize[0] * info.destSize[1]];
00207         }
00208         else {
00209             fwd2DPlan = fftw2d_create_plan(info.srcSize[0], info.srcSize[1], FFTW_FORWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00210             bwd1DPlan = fftw_create_plan(info.srcSize[1], FFTW_BACKWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00211         }
00212     }
00213     else {
00214         if (info.packingMode == SPARSE_MODE) {
00215             fwd1DaPlan = fftw_create_plan(info.destSize[], FFTW_FORWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00216             fwd1DbPlan = fftw_create_plan(info.destSize[], FFTW_FORWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00217             bwd1DaPlan = fftw_create_plan(info.destSize[], FFTW_BACKWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00218             bwd1DbPlan = fftw_create_plan(info.destSize[], FFTW_BACKWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00219             fwd2DPlan = fftw2d_create_plan(info.destSize[0], info.destSize[1], FFTW_FORWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00220             bwd2DPlan = fftw2d_create_plan(info.destSize[0], info.destSize[1], FFTW_BACWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00221         }
00222         else {
00223             bwd2DPlan = fftw2d_create_plan(info.destSize[0], info.destSize[1], FFTW_BACKWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00224             fwd1DPlan = fftw_create_plan(info.destSize[1], FFTW_FORWARD, FFTW_USE_WISDOM|FFTW_MEASURE|FFTW_IN_PLACE);
00225         }
00226     }
00227 }
00228 
00229 #endif
libs/ck-libs/fftlib/fftlib-sparse.C
