Line data    Source code

       1             : #define _LOCAL_SINGLETHREADED 1
       2             : 
       3             : #if (_LOCAL_SINGLETHREADED)
       4             :   template <class T>
       5           0 :   void AWVisResampler::accumulateFromGrid(T& nvalue, 
       6             :      const T* __restrict__& grid, Vector<Int>& iGrdPos,
       7             :      Complex* __restrict__& convFuncV, Double& wVal, 
       8             :      Vector<Int>& scaledSupport, Vector<Float>& scaledSampling, 
       9             :      Vector<Double>& offset, Vector<Int>& convOrigin, Vector<Int>& cfShape,
      10             :      Vector<Int>& loc, Complex& phasor, Double& sinDPA, Double& cosDPA,
      11             :      Bool& finitePointingOffset, Matrix<Complex>& cached_phaseGrad_p)
      12             :   {
      13             :     (void)sinDPA; (void)cosDPA; (void)cfShape; // Keep the compiler warnings off for now
      14           0 :     Complex wt, norm=0.0;
      15           0 :     Vector<Int> iCFPos(4,0),iLoc(4,0);
      16             : 
      17             :     // cerr << scaledSupport << endl
      18             :     //   << scaledSampling << endl
      19             :     //   << offset << endl
      20             :     //   << cfShape << endl
      21             :     //   << convOrigin << endl;
      22           0 :     IPosition phaseGradOrigin_l; 
      23           0 :     phaseGradOrigin_l = cached_phaseGrad_p.shape()/2;
      24             : 
      25           0 :     for(Int iy=-scaledSupport[1]; iy <= scaledSupport[1]; iy++) 
      26             :       {
      27           0 :         iLoc[1]    = static_cast<int>(scaledSampling[1]*iy+offset[1]);
      28           0 :         iCFPos[1]  = iLoc[1] + convOrigin[1];
      29           0 :         iGrdPos[1] = loc[1]+iy;
      30             : 
      31           0 :         for(Int ix=-scaledSupport[0]; ix <= scaledSupport[0]; ix++) 
      32             :           {
      33           0 :             iLoc[0]    = static_cast<int>(scaledSampling[0]*ix+offset[0]);
      34           0 :             iCFPos[0]  = iLoc[0] + convOrigin[0];
      35           0 :             iGrdPos[0] = loc[0]+ix;
      36             : 
      37             :             {
      38           0 :               wt = getFrom4DArray((const Complex* __restrict__ &) convFuncV,
      39           0 :                                   iCFPos,cfInc_p);
      40             :               //              wt = convFuncV(iCFPos);
      41           0 :               if (wVal > 0.0) wt = conj(wt);
      42             : 
      43           0 :               norm+=(wt);
      44           0 :               if (finitePointingOffset) 
      45           0 :                 wt *= cached_phaseGrad_p(iLoc[0]+phaseGradOrigin_l[0],
      46           0 :                                          iLoc[1]+phaseGradOrigin_l[1]);
      47             :               //              nvalue+=wt*grid(iGrdPos);
      48           0 :               nvalue +=  wt * getFrom4DArray(grid, iGrdPos, gridInc_p);
      49             :             }
      50             :           }
      51             :       }
      52           0 :     nvalue *= conj(phasor)/norm;
      53           0 :   }
      54             : #else
      55             :   //
      56             :   //-----------------------------------------------------------------------------------
      57             :   //
      58             :   template <class T>
      59             :   void AWVisResampler::accumulateFromGrid(T& nvalue, 
      60             :      const T* __restrict__& grid, Vector<Int>& iGrdPos,
      61             :      Complex* __restrict__& convFuncV, Double& wVal, 
      62             :      Vector<Int>& scaledSupport, Vector<Float>& scaledSampling, 
      63             :      Vector<Double>& offset, Vector<Int>& convOrigin, Vector<Int>& cfShape,
      64             :      Vector<Int>& loc, Complex& phasor, Double& sinDPA, Double& cosDPA,
      65             :      Bool& finitePointingOffset, Matrix<Complex>& cached_phaseGrad_p)
      66             :   {
      67             :     Complex wt, norm=0.0;
      68             :     Vector<Int> iCFPos(4,0);
      69             :     Bool dummy;
      70             :   
      71             :     //---------Multi-threading related code starts----------------
      72             :     // Code to make the loops below thread safe.  CASA Array (and
      73             :     // derived classes) are not thread-safe.  So extract the pointers
      74             :     // and use them with raw indexing.
      75             :     //
      76             :     Int *iCFPosPtr=iCFPos.getStorage(dummy), 
      77             :       *scaledSupportPtr=scaledSupport.getStorage(dummy),
      78             :       *iGrdPosPtr=iGrdPos.getStorage(dummy),
      79             :       *locPtr=loc.getStorage(dummy),
      80             :       *cfInc_pPtr=cfInc_p.getStorage(dummy),
      81             :       *convOriginPtr=convOrigin.getStorage(dummy),
      82             :       *gridInc_pPtr=gridInc_p.getStorage(dummy),
      83             :       *phaseGradOriginPtr;
      84             :     Float *scaledSamplingPtr=scaledSampling.getStorage(dummy);
      85             :     Double *offsetPtr=offset.getStorage(dummy);
      86             :     Complex *cached_phaseGrad_pPtr=cached_phaseGrad_p.getStorage(dummy);
      87             :     Int phaseGradNx=cached_phaseGrad_p.shape()[0],
      88             :       phaseGradNy=cached_phaseGrad_p.shape()[1];
      89             :     Vector<Int> phaseGradOrigin_l; 
      90             :     phaseGradOrigin_l = cached_phaseGrad_p.shape().asVector()/2;
      91             :     phaseGradOriginPtr = phaseGradOrigin_l.getStorage(dummy);
      92             :     //---------Multi-threading related code ends-----------------
      93             :     Int Nth = 1;
      94             : #ifdef _OPENMP
      95             :     Nth=max(omp_get_max_threads()-2,1);
      96             : #endif
      97             :     T *iTHNValue = new T[Nth];
      98             :     Complex *iTHNorm = new Complex[Nth];
      99             :     // cerr << scaledSupport << endl
     100             :     //   << scaledSampling << endl
     101             :     //   << offset << endl
     102             :     //   << cfShape << endl
     103             :     //   << convOrigin << endl;
     104             :     
     105             :     for(Int iy=-scaledSupportPtr[1]; iy <= scaledSupportPtr[1]; iy++) 
     106             :       {
     107             :         iCFPosPtr[1]=(scaledSamplingPtr[1]*iy+offsetPtr[1])+convOriginPtr[1];
     108             :         iGrdPosPtr[1]=locPtr[1]+iy;
     109             : 
     110             :         for (Int th=0;th<Nth;th++) {iTHNValue[th]=0.0;iTHNorm[th]=0.0;}
     111             :         Int thID=0;
     112             : #pragma omp parallel default(none) firstprivate(iCFPosPtr,iGrdPosPtr) \
     113             :   private(wt,thID) \
     114             :   shared(scaledSupportPtr, scaledSamplingPtr, convOriginPtr, offsetPtr,\
     115             :          locPtr,cached_phaseGrad_pPtr, phaseGradNx, phaseGradNy,gridInc_pPtr, \
     116             :          cfInc_pPtr,iTHNValue, iTHNorm,phaseGradOriginPtr)              \
     117             :   num_threads(Nth)
     118             :         {
     119             :           //=================================================================
     120             :           //============= OMP LOOP===========================================
     121             : #pragma omp for
     122             :         for(Int ix=-scaledSupportPtr[0]; ix <= scaledSupportPtr[0]; ix++) 
     123             :           {
     124             : #ifdef _OPENMP
     125             :             thID=omp_get_thread_num();
     126             : #endif
     127             :             Int localCFPos[4];localCFPos[2]=localCFPos[3]=0;
     128             :             Int localGrdPos[4];localGrdPos[2]=iGrdPosPtr[2];localGrdPos[3]=iGrdPosPtr[3];
     129             :             localCFPos[1]=iCFPosPtr[1];
     130             :             localGrdPos[1]=iGrdPosPtr[1];
     131             : 
     132             :             // iCFPosPtr[0]=(scaledSamplingPtr[0]*ix+offsetPtr[0])+convOriginPtr[0];
     133             :             // iGrdPosPtr[0]=locPtr[0]+ix;
     134             : 
     135             :             localCFPos[0]=(scaledSamplingPtr[0]*ix+offsetPtr[0])+convOriginPtr[0];
     136             :             localGrdPos[0]=locPtr[0]+ix;
     137             :             {
     138             :               wt = getFrom4DArray((const Complex* __restrict__ &) convFuncV,
     139             :                                   localCFPos,cfInc_pPtr);
     140             :                                   //  iCFPosPtr,cfInc_pPtr);
     141             :               //              wt = convFuncV(iCFPos);
     142             :               if (wVal > 0.0) wt = conj(wt);
     143             :               //norm+=(wt);
     144             :               iTHNorm[thID]+=(wt);
     145             : 
     146             :               // if (finitePointingOffset) wt *= cached_phaseGrad_p(iCFPosPtr[0],
     147             :               //                                                         iCFPosPtr[1]);
     148             : 
     149             :               // Using raw indexing to make it thread safe
     150             :               if (finitePointingOffset) 
     151             :                 wt *= cached_phaseGrad_pPtr[(localCFPos[0]+phaseGradOriginPtr[0])*phaseGradNx+
     152             :                                             (localCFPos[1]+phaseGradOriginPtr[1])*phaseGradNy];
     153             :               //              nvalue+=wt*grid(iGrdPos);
     154             :               //              nvalue +=  wt * getFrom4DArray(grid, iGrdPosPtr, gridInc_pPtr);
     155             :          
     156             :               // nvalue +=  wt * getFrom4DArray(grid, localGrdPos, gridInc_pPtr);
     157             :               iTHNValue[thID] +=  wt * getFrom4DArray(grid, localGrdPos, gridInc_pPtr);
     158             :             }
     159             :           }
     160             :         } // End omp for
     161             :           //=================================================================
     162             :           //============= OMP LOOP===========================================
     163             :         for (Int th=0;th<Nth;th++) {nvalue += iTHNValue[th]; norm += iTHNorm[th];}
     164             :       }
     165             :     nvalue *= conj(phasor)/norm;
     166             :     delete [] iTHNorm;
     167             :     delete [] iTHNValue;
     168             :   }
     169             : #endif