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