00001 #include "fftlib.h"
00002
00003
00004
00005
00006
00007 void
00008 NormalLineArray::doFirstFFT(int fftid, int direction)
00009 {
00010 LineFFTinfo &fftinfo = (infoVec[fftid]->info);
00011 int ptype = fftinfo.ptype;
00012 int pblock = fftinfo.pblock;
00013 complex *line = fftinfo.dataPtr;
00014 int sizeX = fftinfo.sizeX;
00015 int sizeZ = fftinfo.sizeZ;
00016 int *xsquare = fftinfo.xsquare;
00017 int *ysquare = fftinfo.ysquare;
00018 int *zsquare = fftinfo.zsquare;
00019
00020 #ifdef HEAVYVERBOSE
00021 {
00022 char fname[80];
00023 if(direction)
00024 snprintf(fname,80,"xline_%d.y%d.z%d.out", fftid,thisIndex.x, thisIndex.y);
00025 else
00026 snprintf(fname,80,"zline_%d.x%d.y%d.out", fftid,thisIndex.x, thisIndex.y);
00027 FILE *fp=fopen(fname,"w");
00028 for(int x = 0; x < sizeX*xsquare[0]*xsquare[1]; x++)
00029 fprintf(fp, "%d %g %g\n", x, line[x].re, line[x].im);
00030 fclose(fp);
00031 }
00032 #endif
00033
00034 if(direction && ptype==PencilType::XLINE)
00035 fftw(fwdplan, xsquare[0]*xsquare[1], (fftw_complex*)line, 1, sizeX, NULL, 0, 0);
00036 else if(!direction && ptype==PencilType::ZLINE)
00037 fftw(bwdplan, zsquare[0]*zsquare[1], (fftw_complex*)line, 1, sizeZ, NULL, 0, 0);
00038 else
00039 CkAbort("Can't do this FFT\n");
00040
00041 int x, y, z=0;
00042 #ifdef VERBOSE
00043 CkPrintf("First FFT done at [%d %d] [%d %d]\n", thisIndex.x, thisIndex.y,sizeX,sizeZ);
00044 #endif
00045 int baseX, ix, iy, iz;
00046 if(true) {
00047 if(direction)
00048 {
00049 int sendSquarethick = ysquare[1] <= xsquare[1] ? ysquare[1]:xsquare[1];
00050 int sendDataSize = ysquare[0]*xsquare[0] * sendSquarethick;
00051 int zpos = thisIndex.y;
00052 int index=0;
00053 complex *sendData = NULL;
00054
00055 for(z = 0; z < xsquare[1]; z+=sendSquarethick){
00056 for(x = 0; x < sizeX; x+=ysquare[0]) {
00057 SendFFTMsg *msg = new(sendDataSize, sizeof(int)*8) SendFFTMsg;
00058 sendData = msg->data;
00059 msg->ypos = thisIndex.x;
00060 msg->size = sendDataSize;
00061 msg->id = fftid;
00062 msg->direction = direction;
00063 msg->data = sendData;
00064 CkSetQueueing(msg, CK_QUEUEING_IFIFO);
00065 #ifdef _PRIOMSG_
00066 int prioNum = (zpos+z) + x*sizeX;
00067 *(int*)CkPriorityPtr(msg) = prioNum;
00068 #endif
00069 index = 0;
00070 for(iz = z; iz < z+sendSquarethick; iz++)
00071 for(ix = x; ix < x+ysquare[0]; ix++)
00072 for(y = 0; y < xsquare[0]; y++)
00073 sendData[index++] = line[iz*sizeX*xsquare[0]+y*sizeX+ix];
00074
00075 #ifdef VERBOSE
00076 CkPrintf(" [%d %d] sending to YLINES [ %d %d] \n", thisIndex.x, thisIndex.y, thisIndex.y, x);
00077 #endif
00078 yProxy(zpos+z, x).doSecondFFT(msg);
00079 }
00080
00081 }
00082 }
00083 else
00084 {
00085 int sendSquarewidth = ysquare[0]<=zsquare[0] ? ysquare[0]:zsquare[0];
00086 int sendDataSize = ysquare[1] * sendSquarewidth * zsquare[1];
00087 int xpos = thisIndex.x;
00088 int ypos = thisIndex.y;
00089 int index=0;
00090 complex *sendData = NULL;
00091
00092 for(x = 0; x < zsquare[0]; x+=sendSquarewidth)
00093 for(z = 0; z < sizeZ; z+=ysquare[1]){
00094 SendFFTMsg *msg = new(sendDataSize, sizeof(int)*8) SendFFTMsg;
00095 sendData = msg->data;
00096 msg->ypos = thisIndex.y;
00097 msg->size = sendDataSize;
00098 msg->id = fftid;
00099 msg->direction = direction;
00100 msg->data = sendData;
00101 CkSetQueueing(msg, CK_QUEUEING_IFIFO);
00102 #ifdef _PRIOMSG_
00103 int prioNum = (z) + (x+xpos)*sizeX;
00104 *(int*)CkPriorityPtr(msg) = prioNum;
00105 #endif
00106 index = 0;
00107 for(iz = z; iz < z+ysquare[1]; iz++)
00108 for (ix = x; ix < x+sendSquarewidth; ix++)
00109 for(iy = 0; iy < zsquare[1]; iy++)
00110 sendData[index++] = line[iz+ix*sizeZ+iy*sizeZ*zsquare[0]];
00111
00112 #ifdef VERBOSE
00113 CkPrintf(" [%d %d] sending to YLINES [%d %d] \n", thisIndex.x, thisIndex.y, z, thisIndex.x);
00114 #endif
00115 yProxy(z, xpos+x).doSecondFFT(msg);
00116 }
00117 }
00118 }
00119 }
00120
00121 void
00122 NormalLineArray::doSecondFFT(int ypos, complex *val, int datasize, int fftid, int direction)
00123 {
00124 LineFFTinfo &fftinfo = (infoVec[fftid]->info);
00125 int ptype = fftinfo.ptype;
00126 complex *line = fftinfo.dataPtr;
00127 int sizeY = fftinfo.sizeY;
00128 int *xsquare = fftinfo.xsquare;
00129 int *ysquare = fftinfo.ysquare;
00130 int *zsquare = fftinfo.zsquare;
00131 int expectSize = 0, expectMsg = 0;
00132 int x,y,z,baseY;
00133 int ix,iy,iz;
00134 if(direction){
00135 int sendSquarethick = ysquare[1]<=xsquare[1] ? ysquare[1]:xsquare[1];
00136 expectSize = ysquare[0]*xsquare[0] * sendSquarethick;
00137 expectMsg = sizeY/xsquare[0] * (ysquare[1]/sendSquarethick);
00138 CkAssert(datasize == expectSize);
00139 int idx = 0;
00140 for(z=0; z<sendSquarethick; z++)
00141 for(x=0; x<ysquare[0]; x++)
00142 for(y=0; y<xsquare[0]; y++)
00143 line[z*sizeY*ysquare[0]+x*sizeY+ypos+y] = val[idx++];
00144 }
00145 else {
00146 int sendSquarewidth = ysquare[0]<=zsquare[0] ? ysquare[0]:zsquare[0];
00147 expectSize = ysquare[1] * sendSquarewidth * zsquare[1];
00148 expectMsg = sizeY/zsquare[1] * (ysquare[0]/sendSquarewidth);
00149 CkAssert(datasize == expectSize);
00150 int idx=0;
00151 for(z=0; z<ysquare[1]; z++)
00152 for(x=0; x<sendSquarewidth; x++)
00153 for(y=0; y<zsquare[1]; y++)
00154 line[z*sizeY*ysquare[0]+x*sizeY+ypos+y] = val[idx++];
00155 }
00156 infoVec[fftid]->count++;
00157 if (infoVec[fftid]->count == expectMsg) {
00158 infoVec[fftid]->count = 0;
00159 int y;
00160
00161
00162 if(direction && ptype==PencilType::YLINE)
00163 fftw(fwdplan, ysquare[0]*ysquare[1], (fftw_complex*)line, 1, sizeY, NULL, 0, 0);
00164 else if(!direction && ptype==PencilType::YLINE)
00165 fftw(bwdplan, ysquare[0]*ysquare[1], (fftw_complex*)line, 1, sizeY, NULL, 0, 0);
00166 else
00167 CkAbort("Can't do this FFT\n");
00168
00169 #ifdef HEAVYVERBOSE
00170 {
00171 char fname[80];
00172 snprintf(fname,80,"yline_%d.x%d.z%d.out", fftid, thisIndex.y, thisIndex.x);
00173 FILE *fp=fopen(fname,"w");
00174 for(int x = 0; x < sizeY*ysquare[0]*ysquare[1]; x++)
00175 fprintf(fp, "%d %g %g\n", x, line[x].re, line[x].im);
00176 fclose(fp);
00177 }
00178 #endif
00179
00180 #ifdef VERBOSE
00181 CkPrintf("Second FFT done at [%d %d]\n", thisIndex.x, thisIndex.y);
00182 #endif
00183
00184 if(direction){
00185 int sendSquarewidth = ysquare[0]<=zsquare[0] ? ysquare[0]:zsquare[0];
00186 int sendDataSize = sendSquarewidth * ysquare[1] * zsquare[1];
00187 int xpos = thisIndex.y;
00188 int index=0;
00189 complex *sendData = NULL;
00190
00191 for(x = 0; x < ysquare[0]; x+=sendSquarewidth)
00192 for(y = 0; y < sizeY; y+=zsquare[1]) {
00193 SendFFTMsg *msg = new(sendDataSize, sizeof(int)*8) SendFFTMsg;
00194 sendData = msg->data;
00195 msg->zpos = thisIndex.x;
00196 msg->ypos = 0;
00197 msg->size = sendDataSize;
00198 msg->id = fftid;
00199 msg->direction = direction;
00200 msg->data = sendData;
00201 CkSetQueueing(msg, CK_QUEUEING_IFIFO);
00202 #ifdef _PRIOMSG_
00203 int prioNum = (xpos+x) + y*sizeY;
00204 *(int*)CkPriorityPtr(msg) = prioNum;
00205 #endif
00206 index = 0;
00207 for(iy = y; iy < y+zsquare[1]; iy++)
00208 for(ix = x; ix < x+sendSquarewidth; ix++)
00209 for(iz = 0; iz < ysquare[1]; iz++)
00210 sendData[index++] = line[iz*sizeY*ysquare[0]+ix*sizeY+iy];
00211
00212 #ifdef VERBOSE
00213 CkPrintf(" [%d %d] sending to ZLINES [ %d %d] \n", thisIndex.x, thisIndex.y, thisIndex.y, y);
00214 #endif
00215 (zProxy)(x+xpos, y).doThirdFFT(msg);
00216 }
00217 }
00218 else {
00219 int sendSquarethick = ysquare[1]<=xsquare[1] ? ysquare[1]:xsquare[1];
00220 int sendDataSize = ysquare[0]*xsquare[0] * sendSquarethick;
00221 int zpos = thisIndex.x;
00222 int index, ix, iy, iz;
00223 complex *sendData = NULL;
00224
00225 for(z = 0; z < ysquare[1]; z+=sendSquarethick)
00226 for (y = 0; y < sizeY; y+=xsquare[0]) {
00227 SendFFTMsg *msg = new(sendDataSize, sizeof(int)*8) SendFFTMsg;
00228 sendData = msg->data;
00229 msg->zpos = 0;
00230 msg->ypos = thisIndex.y;
00231 msg->size = sendDataSize;
00232 msg->id = fftid;
00233 msg->direction = direction;
00234 msg->data = sendData;
00235 CkSetQueueing(msg, CK_QUEUEING_IFIFO);
00236 #ifdef _PRIOMSG_
00237 int prioNum = (y)*sizeY + (zpos+z);
00238 *(int*)CkPriorityPtr(msg) = prioNum;
00239 #endif
00240 index = 0;
00241 for(iz = z; iz < z+sendSquarethick; iz++)
00242 for(iy = y; iy < y+xsquare[0]; iy++)
00243 for(x = 0; x < ysquare[0]; x++)
00244 sendData[index++] = line[iz*sizeY*ysquare[0]+x*sizeY+iy];
00245
00246 #ifdef VERBOSE
00247 CkPrintf(" [%d %d] sending to XLINES [%d %d] \n", thisIndex.x, thisIndex.y, y, zpos+z);
00248 #endif
00249 (xProxy)(y, zpos+z).doThirdFFT(msg);
00250 }
00251 }
00252 }
00253 }
00254
00255 void
00256 NormalLineArray::doThirdFFT(int zpos, int xpos, complex *val, int datasize, int fftid, int direction)
00257 {
00258 LineFFTinfo &fftinfo = (infoVec[fftid]->info);
00259 int ptype = fftinfo.ptype;
00260 complex *line = fftinfo.dataPtr;
00261 int sizeX = fftinfo.sizeX;
00262 int sizeZ = fftinfo.sizeZ;
00263 int *xsquare = fftinfo.xsquare;
00264 int *ysquare = fftinfo.ysquare;
00265 int *zsquare = fftinfo.zsquare;
00266 int expectSize=0, expectMsg=0, offset=0, i;
00267
00268 int x,y,z,idx;
00269 if(direction){
00270 int sendSquarewidth = ysquare[0]<=zsquare[0] ? ysquare[0]:zsquare[0];
00271 expectSize = sendSquarewidth * ysquare[1] * zsquare[1];
00272 expectMsg = sizeZ/ysquare[1] * (zsquare[0]/sendSquarewidth);
00273 CkAssert(datasize == expectSize);
00274 idx=0;
00275 for(y=0; y<zsquare[1]; y++)
00276 for(x=0; x<sendSquarewidth; x++)
00277 for(z=0; z<ysquare[1]; z++)
00278 line[z+zpos+(x+xpos)*sizeZ+y*sizeZ*zsquare[0]] = val[idx++];
00279 }
00280 else{
00281 int sendSquarethick = ysquare[1]<=xsquare[1] ? ysquare[1]:xsquare[1];
00282 expectSize = ysquare[0]*xsquare[0] * sendSquarethick;
00283 expectMsg = sizeX/ysquare[0] * (xsquare[1]/sendSquarethick);
00284 CkAssert(datasize == expectSize);
00285 int idx=0;
00286 for(z=0; z<sendSquarethick; z++)
00287 for(y=0; y<xsquare[0]; y++)
00288 for(x=0; x<ysquare[0]; x++)
00289 line[(z+zpos)*sizeX*xsquare[0]+y*sizeX+xpos+x] = val[idx++];
00290 }
00291
00292 infoVec[fftid]->count ++;
00293
00294 if (infoVec[fftid]->count == expectMsg) {
00295 infoVec[fftid]->count = 0;
00296
00297
00298 #ifdef HEAVYVERBOSE
00299 {
00300 char fname[80];
00301 if(direction)
00302 snprintf(fname,80,"zline_%d.x%d.y%d.out", fftid, thisIndex.x, thisIndex.y);
00303 else
00304 snprintf(fname,80,"xline_%d.y%d.z%d.out", fftid, thisIndex.x, thisIndex.y);
00305 FILE *fp=fopen(fname,"w");
00306 for(int x = 0; x < sizeX*xsquare[0]*xsquare[1]; x++)
00307 fprintf(fp, "%g %g\n", line[x].re, line[x].im);
00308 fclose(fp);
00309 }
00310 #endif
00311
00312
00313 if(direction && ptype==PencilType::ZLINE)
00314 fftw(fwdplan, zsquare[0]*zsquare[1], (fftw_complex*)line, 1, sizeX, NULL, 0, 0);
00315 else if(!direction && ptype==PencilType::XLINE)
00316 fftw(bwdplan, xsquare[0]*xsquare[1], (fftw_complex*)line, 1, sizeX, NULL, 0, 0);
00317 else
00318 CkAbort("Can't do this FFT\n");
00319 #ifdef VERBOSE
00320 CkPrintf("Third FFT done at [%d %d]\n", thisIndex.x, thisIndex.y);
00321 #endif
00322 doneFFT(fftid, direction);
00323
00324 }
00325 }
00326
00327 void
00328 NormalLineArray::doSecondFFT(SendFFTMsg *msg){
00329 doSecondFFT(msg->ypos, msg->data, msg->size, msg->id, msg->direction);
00330 delete msg;
00331 }
00332
00333 void
00334 NormalLineArray::doThirdFFT(SendFFTMsg *msg){
00335 doThirdFFT(msg->zpos, msg->ypos, msg->data, msg->size, msg->id, msg->direction);
00336 delete msg;
00337 }
00338
00339 void
00340 NormalLineArray::doneFFT(int id, int direction){
00341 #ifdef VERBOSE
00342 CkPrintf("FFT finished \n");
00343 #endif
00344 }
00345
00346 NormalLineArray::NormalLineArray (LineFFTinfo &info, CProxy_NormalLineArray _xProxy, CProxy_NormalLineArray _yProxy, CProxy_NormalLineArray _zProxy) {
00347 #ifdef VERBOSE
00348 CkPrintf("inserted line %d index[%d %d]\n", info.ptype, thisIndex.x, thisIndex.y);
00349 #endif
00350 setup(info, _xProxy, _yProxy, _zProxy);
00351 }
00352
00353 void
00354 NormalLineArray::setup (LineFFTinfo &info, CProxy_NormalLineArray _xProxy, CProxy_NormalLineArray _yProxy, CProxy_NormalLineArray _zProxy) {
00355 xProxy = _xProxy;
00356 yProxy = _yProxy;
00357 zProxy = _zProxy;
00358
00359 PencilArrayInfo *pencilinfo = new PencilArrayInfo();
00360 pencilinfo->info = info;
00361 pencilinfo->count = 0;
00362 infoVec.insert(infoVec.size(), pencilinfo);
00363
00364 line = NULL;
00365 fwdplan = fftw_create_plan(info.sizeX, FFTW_FORWARD, FFTW_IN_PLACE);
00366 bwdplan = fftw_create_plan(info.sizeY, FFTW_BACKWARD, FFTW_IN_PLACE);
00367 id = -1;
00368 }
1.5.5