Source
277
277
Float maxconv=abs(screen(0,0));
278
278
for (Int chunkId=nchunk-1; chunkId >= 0; --chunkId){
279
279
//cerr << "chunkId " << chunkId << endl;
280
280
Cube<Complex> convFuncSect( maxConvSize/2-1, maxConvSize/2-1, chunksize(chunkId));
281
281
convFuncSect.set(0.0);
282
282
Bool convFuncStor=false;
283
283
Complex *convFuncPtr=convFuncSect.getStorage(convFuncStor);
284
284
//////openmp like to share reference param ...but i don't like to share
285
285
Int cpConvSize=maxConvSize;
286
286
//cerr << "orig convsize " << convSize << endl;
287
-
Int cpWConvSize=wConvSize;
287
+
// Int cpWConvSize=wConvSize;
288
288
Double cpWscale=wScale;
289
289
Int wstart=planesPerChunk*chunkId;
290
290
Int wend=wstart+chunksize(chunkId)-1;
291
291
#ifdef _OPENMP
292
292
omp_set_nested(0);
293
293
#endif
294
-
#pragma omp parallel for default(none) firstprivate(cpWConvSize, cpConvSize, convFuncPtr, s0, s1, wsaveptr, lsav, cor, inner, cpWscale, wstart, wend)
294
+
#pragma omp parallel for default(none) firstprivate(/*cpWConvSize,*/ cpConvSize, convFuncPtr, s0, s1, wsaveptr, lsav, cor, inner, cpWscale, wstart, wend)
295
295
for (Int iw=wstart; iw < (wend+1) ; ++iw) {
296
296
Matrix<Complex> screen1(cpConvSize, cpConvSize);
297
297
makeGWplane(screen1, iw, s0, s1, wsaveptr, lsav, inner, cor, cpWscale);
298
298
ooLong offset=ooLong(ooLong(iw-wstart)*ooLong(cpConvSize/2-1)*ooLong(cpConvSize/2-1));
299
299
for (ooLong y=0; y< ooLong(cpConvSize/2)-1; ++y){
300
300
for (ooLong x=0; x< ooLong(cpConvSize/2)-1; ++x){
301
301
////////convFuncPtr[offset+y*(convSize/2-1)+x] = quarter(x,y);
302
302
convFuncPtr[offset+y*ooLong(cpConvSize/2-1)+x] = screen1(x,y);
303
303
}
304
304
}
305
305
306
306
}//iw
307
307
convFuncSect.putStorage(convFuncPtr, convFuncStor);
308
308
for (uInt iw=0; iw< uInt(chunksize(chunkId)); ++iw){
309
309
convFuncSect.xyPlane(iw)=convFuncSect.xyPlane(iw)/(maxconv);
310
310
}
311
311
convFuncPtr=convFuncSect.getStorage(convFuncStor);
312
312
Int thischunkSize=chunksize(chunkId);
313
313
//cerr << "chunkId* planesPerChunk " << chunkId* planesPerChunk << " chunksize " << thischunkSize << endl;
314
-
#pragma omp parallel for default(none) firstprivate(suppstor, cpConvSize, cpWConvSize, convFuncPtr, maxConvSize, thischunkSize, wstart, planesPerChunk, chunkId) reduction(+: warner)
314
+
#pragma omp parallel for default(none) firstprivate(suppstor, cpConvSize, /*cpWConvSize,*/ convFuncPtr, maxConvSize, thischunkSize, wstart, planesPerChunk, chunkId) reduction(+: warner)
315
315
for (Int iw=0; iw<thischunkSize; iw++) {
316
316
Bool found=false;
317
317
Int trial=0;
318
318
ooLong ploffset=(ooLong)(cpConvSize/2-1)*(ooLong)(cpConvSize/2-1)*(ooLong)iw;
319
319
for (trial=cpConvSize/2-2;trial>0;trial--) {
320
320
// if((abs(convFunc(trial,0,iw))>1e-3)||(abs(convFunc(0,trial,iw))>1e-3) ) {
321
321
if((abs(convFuncPtr[(ooLong)(trial)+ploffset])>1e-3)||(abs(convFuncPtr[(ooLong)(trial*(cpConvSize/2-1))+ploffset])>1e-3) ) {
322
322
//cout <<"iw " << iw << " x " << abs(convFunc(trial,0,iw)) << " y "
323
323
// <<abs(convFunc(0,trial,iw)) << endl;
324
324
found=true;