00001
00002 #include "pencilfft.h"
00003
00004 #define CONST_A 1001
00005 #define CONST_B 17
00006
00010 LineFFTArray::LineFFTArray (LineFFTInfo &info, int phase): _info (info)
00011 {
00012 #if __LINEFFT_DEBUG__
00013 CkPrintf ("[%d, %d, %d] Pencil FFT array constructor \n", thisIndex.x, thisIndex.y, phase);
00014 #endif
00015
00016 _nElements = 0;
00017 _fftSize = 0;
00018 _phase = (LineFFTPhase) phase;
00019
00020 switch (phase) {
00021 case PHASE_X:
00022 _fftSize = info.sizeX;
00023 _nElements = info.sizeX * info.grainY * info.grainZ;
00024 _nFFTMessages = info.sizeX / info.grainX;
00025 break;
00026
00027 case PHASE_Y:
00028 _fftSize = info.sizeY;
00029 _nElements = info.sizeY * info.grainZ * info.grainX;
00030 _nFFTMessages = info.sizeY / info.grainY;
00031 break;
00032
00033 case PHASE_Z:
00034 _fftSize = info.sizeZ;
00035 _nElements = info.sizeZ * info.grainX * info.grainY;
00036 _nFFTMessages = info.sizeZ / info.grainZ;
00037 break;
00038 };
00039
00040 _line = new complex [_nElements];
00041
00042 LineFFTAssert ((_nFFTMessages % CONST_B) != 0);
00043
00044 #if 0 //__LINEFFT_DEBUG__
00045 for (int count = 0; count < _nElements; count++)
00046 _line[count] = complex (1.0 * count, 0.0);
00047 #endif
00048
00049 _fwdplan = fftw_create_plan(_fftSize, FFTW_FORWARD, FFTW_IN_PLACE |
00050 FFTW_MEASURE | FFTW_USE_WISDOM);
00051 _bwdplan = fftw_create_plan(_fftSize, FFTW_BACKWARD, FFTW_IN_PLACE |
00052 FFTW_MEASURE | FFTW_USE_WISDOM);
00053
00054 _curFFTMessages = 0;
00055 _curGridMessages = 0;
00056 _nGridMessages = 0;
00057 }
00058
00059
00063 void LineFFTArray::receiveFFTMessage (LineFFTMsg *msg) {
00064
00065 #if __LINEFFT_DEBUG__
00066 CkPrintf ("[%d, %d] Pencil FFT receive fft message \n",
00067 thisIndex.x, thisIndex.y);
00068 #endif
00069 int index = 0, localindex=0;
00070 int start = 0;
00071
00072 #if __LINEFFT_DEBUG__
00073 LineFFTAssert (msg->phase == _phase);
00074 #endif
00075
00076 switch (_phase) {
00077 case PHASE_X:
00078 start = msg->chunk_id * _info.grainX;
00079 for (int z = 0; z < _info.grainZ; z++)
00080 for (int y = 0; y < _info.grainY; y++)
00081 for (int x = start; x < start + _info.grainX; x++) {
00082 localindex = _info.sizeX * _info.grainY * z + _info.sizeX * y + x;
00083 #if __LINEFFT_DEBUG__
00084 LineFFTAssert (localindex < _nElements);
00085 #endif
00086 _line [localindex] = msg->data[index++];
00087 }
00088 break;
00089
00090 case PHASE_Y:
00091 start = msg->chunk_id * _info.grainY;
00092 for (int z = 0; z < _info.grainZ; z++)
00093 for (int y = start; y < start+_info.grainY; y++)
00094 for (int x = 0; x < _info.grainX; x++) {
00095 localindex = _info.sizeY * _info.grainZ * x + _info.sizeY * z + y;
00096 #if __LINEFFT_DEBUG__
00097 LineFFTAssert (localindex < _nElements);
00098 #endif
00099 _line [localindex] = msg->data[index++];
00100 }
00101 break;
00102
00103 case PHASE_Z:
00104 start = msg->chunk_id * _info.grainZ;
00105 for (int z = start; z < start+_info.grainZ; z++)
00106 for (int y = 0; y < _info.grainY; y++)
00107 for (int x = 0; x < _info.grainX; x++) {
00108 localindex = _info.sizeZ * _info.grainX * y + _info.sizeZ * x + z;
00109 #if __LINEFFT_DEBUG__
00110 LineFFTAssert (localindex < _nElements);
00111 #endif
00112 _line [localindex] = msg->data[index++];
00113 }
00114 break;
00115
00116 default:
00117 CkAbort ("Phase undefined");
00118 break;
00119 };
00120
00121 _curFFTMessages ++;
00122
00123 if (_curFFTMessages == _nFFTMessages) {
00124 _curFFTMessages = 0;
00125 startFFT (msg->direction);
00126 }
00127
00128 CkFreeMsg (msg);
00129 }
00130
00137 void LineFFTArray::sendFFTMessages (int dir) {
00138
00139 int size = _info.grainX * _info.grainY * _info.grainZ;
00140 CProxy_LineFFTArray proxy;
00141 int arridx_x = 0, arridx_y =0;
00142
00143 for (int count = 0; count < _nFFTMessages; count ++) {
00144 int start = 0, index = 0, localindex =0;
00145 int chunk_id = (CONST_A*CkMyPe() + CONST_B * count) % _nFFTMessages;
00146
00147 LineFFTMsg *msg = new (size, sizeof(int) * 8) LineFFTMsg;
00148
00149 if (dir == FFTW_FORWARD) {
00150 arridx_x = thisIndex.y;
00151 arridx_y = chunk_id;
00152 msg->chunk_id = thisIndex.x;
00153 }
00154 else {
00155 arridx_x = chunk_id;
00156 arridx_y = thisIndex.x;
00157 msg->chunk_id = thisIndex.y;
00158 }
00159
00160 switch (_phase) {
00161 case PHASE_X:
00162 start = chunk_id * _info.grainX;
00163 for (int z = 0; z < _info.grainZ; z++)
00164 for (int y = 0; y < _info.grainY; y++)
00165 for (int x = start; x < start+_info.grainX; x++) {
00166 localindex = _info.sizeX * _info.grainY * z + _info.sizeX * y + x;
00167 #if __LINEFFT_DEBUG__
00168 LineFFTAssert (localindex < _nElements);
00169 #endif
00170 msg->data[index++] = _line [localindex];
00171 }
00172
00173 LineFFTAssert (dir == FFTW_FORWARD);
00174 proxy = _info.yProxy;
00175 msg->phase = PHASE_Y;
00176
00177 break;
00178
00179 case PHASE_Y:
00180 start = chunk_id * _info.grainY;
00181 for (int z = 0; z < _info.grainZ; z++)
00182 for (int y = start; y < start+_info.grainY; y++)
00183 for (int x = 0; x < _info.grainX; x++) {
00184 localindex = _info.sizeY * _info.grainZ * x + _info.sizeY * z + y;
00185 #if __LINEFFT_DEBUG__
00186 LineFFTAssert (localindex < _nElements);
00187 #endif
00188 msg->data[index++] = _line [localindex];
00189 }
00190
00191 if (dir == FFTW_FORWARD) {
00192 proxy = _info.zProxy;
00193 msg->phase = PHASE_Z;
00194 }
00195 else {
00196 proxy = _info.xProxy;
00197 msg->phase = PHASE_X;
00198 }
00199
00200 break;
00201
00202 case PHASE_Z:
00203 start = chunk_id * _info.grainZ;
00204 for (int z = start; z < start+_info.grainZ; z++)
00205 for (int y = 0; y < _info.grainY; y++)
00206 for (int x = 0; x < _info.grainX; x++) {
00207 localindex = _info.sizeZ * _info.grainX * y + _info.sizeZ * x + z;
00208 #if __LINEFFT_DEBUG__
00209 LineFFTAssert (localindex < _nElements);
00210 #endif
00211 msg->data[index++] = _line [localindex];
00212 }
00213
00214 LineFFTAssert (dir == FFTW_BACKWARD);
00215 proxy = _info.yProxy;
00216 msg->phase = PHASE_Y;
00217
00218 break;
00219
00220 default:
00221 CkAbort ("Phase undefined");
00222 break;
00223 };
00224
00225 msg->direction = dir;
00226
00227
00228 proxy(arridx_x, arridx_y).receiveFFTMessage (msg);
00229 }
00230
00231 }
00232
00237 void LineFFTArray::startFFT (int dir) {
00238
00239 #if __LINEFFT_DEBUG__
00240 CkPrintf ("[%d, %d, %d, %d] Pencil FFT start fft %d,%d\n",
00241 thisIndex.x, thisIndex.y, _phase, dir, _nElements, _fftSize);
00242 #endif
00243
00244 fftw_plan plan = NULL;
00245
00246 if (dir == FFTW_FORWARD)
00247 plan = _fwdplan;
00248 else
00249 plan = _bwdplan;
00250
00251 for (int count = 0; count < _nElements; count += _fftSize) {
00252 CmiNetworkProgress();
00253 fftw_one (plan, (fftw_complex *) (_line + count), NULL);
00254 }
00255
00256 if (dir == FFTW_FORWARD && _phase == PHASE_Z) {
00257 if (_info.normalize) {
00258 fftw_real scale = 1.0 / (_info.sizeX * _info.sizeY * _info.sizeZ);
00259 for (int count = 0; count < _nElements; count ++)
00260 _line [count] *= scale;
00261 }
00262
00263
00264 startFFT (FFTW_BACKWARD);
00265 return;
00266 }
00267
00268 if (dir == FFTW_BACKWARD && _phase == PHASE_X) {
00269 #if __LINEFFT_DEBUG__
00270 for (int count = 0; count < _nElements; count++) {
00271 if (!((_line[count].re <= 1.001 * count)
00272 && (_line[count].re >= 0.999 * count)))
00273 CkPrintf ("[%d, %d] val = (%5.3f, %5.3f) \n", thisIndex.x, thisIndex.y,
00274 _line[count].re, _line[count].im);
00275 }
00276 #endif
00277
00278 call_donecallback ();
00279 return;
00280 }
00281
00282 sendFFTMessages (dir);
00283 }
00284
00285
00290 void LineFFTArray::receiveGridMessage (LineFFTGridMsg *msg) {
00291
00292 LineFFTAssert (_phase == PHASE_X);
00293
00294 int my_start_y = thisIndex.x * _info.grainY;
00295 int my_start_z = thisIndex.y * _info.grainZ;
00296
00297 LineFFTGrid &grid = msg->grid;
00298
00299 int index = 0;
00300 for (int z = grid.z0; z < grid.z0 + grid.zsize; z++){
00301 LineFFTAssert ((z >= my_start_z) && (z < my_start_z + _info.grainZ));
00302 for (int y = grid.y0; y < grid.y0 + grid.ysize; y++) {
00303 LineFFTAssert ((y >= my_start_y) && (y < my_start_y + _info.grainY));
00304 for (int x = grid.x0; x < grid.x0 + grid.xsize; x++) {
00305 int rel_x = x;
00306 if (rel_x < 0)
00307 rel_x = x + _info.sizeX;
00308 if (rel_x >= _info.sizeX)
00309 rel_x -= _info.sizeX;
00310
00311 int local_index = (z - my_start_z) * _info.grainY * _info.sizeX +
00312 (y - my_start_y) * _info.sizeX + rel_x;
00313
00314 LineFFTAssert (msg->data [index] == 2.0);
00315
00316 LineFFTAssert (local_index < _nElements);
00317 _line [local_index].re += msg->data [index++];
00318 }
00319 }
00320 }
00321
00322 _gridList [_curGridMessages] = msg->grid;
00323
00324 CkFreeMsg (msg);
00325 _curGridMessages ++;
00326
00327 if (_curGridMessages == _nGridMessages) {
00328 _curGridMessages = 0;
00329 startFFT (FFTW_FORWARD);
00330 }
00331 }
00332
00333
00338 void LineFFTArray::sendGridMessages () {
00339
00340 int my_start_y = thisIndex.x * _info.grainY;
00341 int my_start_z = thisIndex.y * _info.grainZ;
00342
00343 for (int count = 0; count < _nGridMessages; count ++) {
00344 LineFFTGrid &grid = _gridList[count];
00345
00346 LineFFTGridMsg *msg =
00347 new (grid.xsize * grid.ysize * grid.zsize, sizeof(int) * 4) LineFFTGridMsg;
00348
00349 int index = 0;
00350 for (int z = grid.z0; z < grid.z0 + grid.zsize; z++){
00351 LineFFTAssert ((z >= my_start_z) && (z < my_start_z + _info.grainZ));
00352 for (int y = grid.y0; y < grid.y0 + grid.ysize; y++) {
00353 LineFFTAssert ((y >= my_start_y) && (y < my_start_y + _info.grainY));
00354 for (int x = grid.x0; x < grid.x0 + grid.xsize; x++) {
00355 int rel_x = x;
00356 if (rel_x < 0)
00357 rel_x = x + _info.sizeX;
00358 if (rel_x >= _info.sizeX)
00359 rel_x -= _info.sizeX;
00360
00361 int local_index = (z - my_start_z) * _info.grainY * _info.sizeX +
00362 (y - my_start_y) * _info.sizeX + rel_x;
00363
00364 LineFFTAssert (local_index < _nElements);
00365 msg->data [index++] = _line [local_index].re;
00366 }
00367 }
00368 }
00369
00370 grid.cb_done.send (msg);
00371 }
00372 }
00373
00374 #include "PencilFFT.def.h"
00375
1.5.5