Ignore:
Timestamp:
02/01/22 15:28:48 (2 years ago)
Author:
ymipsl
Message:

Improve reduction transformation

  • make the difference between reduction over geometry or reduction between process.
  • geometrical reduction :

domain -> axis
axis -> scalar
domain -> scalar

  • reduction across processes for redondant geometrical cell :

axis -> axis
scalar -> scalar

Reduction can be local (only for the geometrical cell owned by current process) or global, using the "local" attribute (bool) over the reduction.

YM

Location:
XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/gatherer_connector.hpp

    r2267 r2291  
    99#include "distributed_view.hpp" 
    1010#include "context_client.hpp" 
     11#include "reduction_types.hpp" 
    1112 
    1213 
    … …  
    2930       
    3031      template<typename T> 
    31       void transfer(int repeat, int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut) 
     32      void transfer(int repeat, int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, EReduction op = EReduction::none) 
    3233      { 
    3334        // for future, make a specific transfer function for sizeT=1 to avoid multiplication (increasing performance) 
     35        
    3436        size_t dstSlice = dstSize_*sizeT ; 
    3537        dataOut.resize(repeat* dstSlice) ;   
    3638         
    37         for(auto& data : dataIn) 
    38         { 
     39         
     40        if (op == EReduction::none) // tranfer without reduction 
     41        { 
     42          for(auto& data : dataIn) 
     43          { 
     44            T* output = dataOut.dataFirst() ; 
     45            int rank=data.first ; 
     46            auto input  = data.second.dataFirst() ; 
     47            auto& connector=connector_[rank] ; 
     48            auto& mask=mask_[rank] ; 
     49            int size=mask.size() ; 
     50            size_t srcSlice = size * sizeT ; 
     51            for(int l=0; l<repeat; l++) 
     52            {    
     53              for(int i=0, j=0 ;i<size;i++) 
     54              { 
     55                if (mask[i])  
     56                { 
     57                  int cj = connector[j]*sizeT ; 
     58                  int ci = i*sizeT ; 
     59                  for (int k=0;k<sizeT;k++) output[cj+k] = input[ci+k] ; 
     60                  j++ ; 
     61                } 
     62              } 
     63              input+=srcSlice ; 
     64              output+=dstSlice ; 
     65           } 
     66          } 
     67        } 
     68        else // with reduction, to be optimized, see how to vectorize 
     69        { 
     70          vector<int> vcount(dataOut.size(),0) ; 
     71          int* count = vcount.data() ; 
     72          T defaultValue = std::numeric_limits<T>::quiet_NaN(); 
     73          for(auto& data : dataIn) 
     74          { 
     75            T* output = dataOut.dataFirst() ; 
     76            int rank=data.first ; 
     77            auto input  = data.second.dataFirst() ; 
     78            auto& connector=connector_[rank] ; 
     79            auto& mask=mask_[rank] ; 
     80            int size=mask.size() ; 
     81            size_t srcSlice = size * sizeT ; 
     82            for(int l=0; l<repeat; l++) 
     83            {    
     84              for(int i=0, j=0 ;i<size;i++) 
     85              { 
     86                if (mask[i])  
     87                { 
     88                  int cj = connector[j]*sizeT ; 
     89                  int ci = i*sizeT ; 
     90                  for (int k=0;k<sizeT;k++) 
     91                  {  
     92                    if (!std::isnan(input[ci+k])) // manage missing value 
     93                    { 
     94                      if (count[cj+k]==0)  output[cj+k] = input[ci+k] ; 
     95                      else 
     96                      { 
     97                        switch(op) 
     98                        { 
     99                          case EReduction::sum : 
     100                            output[cj+k] += input[ci+k] ; 
     101                            break ; 
     102                          case EReduction::min : 
     103                            output[cj+k]= std::min(output[cj+k],input[ci+k]) ; 
     104                            break ; 
     105                          case EReduction::max : 
     106                            output[cj+k]= std::max(output[cj+k],input[ci+k]) ; 
     107                            break ; 
     108                          case EReduction::average : 
     109                            output[cj+k] += input[ci+k] ; 
     110                            break ; 
     111                          default : 
     112                            ERROR("CGathererConnector::transfer", 
     113                               <<"reduction operator "<<(int)op<<" is not defined for this operation") ; 
     114                          break ; 
     115                        } 
     116                      } 
     117                      count[cj+k]++ ;  
     118                    } 
     119                  } 
     120                  j++ ; 
     121                } 
     122              } 
     123              input+=srcSlice ; 
     124              output+=dstSlice ; 
     125              count+=dstSlice ; 
     126            } 
     127          } 
     128 
    39129          T* output = dataOut.dataFirst() ; 
    40           int rank=data.first ; 
    41           auto input  = data.second.dataFirst() ; 
    42           auto& connector=connector_[rank] ; 
    43           auto& mask=mask_[rank] ; 
    44           int size=mask.size() ; 
    45           size_t srcSlice = size * sizeT ; 
    46           for(int l=0; l<repeat; l++) 
    47           {    
    48             for(int i=0, j=0 ;i<size;i++) 
    49             { 
    50               if (mask[i])  
    51               { 
    52                 int cj = connector[j]*sizeT ; 
    53                 int ci = i*sizeT ; 
    54                 for (int k=0;k<sizeT;k++) output[cj+k] = input[ci+k] ; 
    55                 j++ ; 
    56               } 
    57             } 
    58             input+=srcSlice ; 
    59             output+=dstSlice ; 
    60           } 
    61         } 
    62       } 
    63  
    64       template<typename T> 
    65       void transfer(int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut) 
    66       { 
    67         transfer(1, sizeT, dataIn, dataOut) ; 
     130          if (op==EReduction::average) 
     131            for(int i=0;  i < dataOut.size() ; i++)  
     132            { 
     133              if (count[i]>0) dataOut[i]/=count[i] ; 
     134              else dataOut[i] = defaultValue ; 
     135            } 
     136          else for(int i=0;  i < dataOut.size() ; i++) if (count[i]==0) dataOut[i] = defaultValue ; 
     137        } 
     138      } 
     139 
     140      template<typename T> 
     141      void transfer(int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, EReduction op = EReduction::none) 
     142      { 
     143        transfer(1, sizeT, dataIn, dataOut, op) ; 
    68144      } 
    69145     
    70146      template<typename T> 
    71       void transfer(map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut) 
    72       { 
    73         transfer(1,dataIn,dataOut) ; 
    74       } 
    75  
    76       template<typename T> 
    77       void transfer(int rank,  shared_ptr<CGathererConnector>* connectors, int nConnectors, const T* input, T* output) 
     147      void transfer(map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, EReduction op = EReduction::none) 
     148      { 
     149        transfer(1,dataIn,dataOut, op) ; 
     150      } 
     151 
     152      template<typename T> 
     153      void transfer(int rank,  shared_ptr<CGathererConnector>* connectors, int nConnectors, const T* input, T* output, EReduction op = EReduction::none, int* count=nullptr) 
    78154      { 
    79155        auto& connector = connector_[rank] ; // probably costly, find a better way to avoid the map 
    … …  
    83159        if (nConnectors==0) 
    84160        { 
    85           for(int i=0, j=0; i<srcSize; i++) 
    86             if (mask[i])  
    87             { 
    88               *(output+connector[j]) = *(input + i) ; 
    89               j++ ; 
    90             } 
     161          if (op == EReduction::none) 
     162          { 
     163            for(int i=0, j=0; i<srcSize; i++) 
     164              if (mask[i])  
     165              { 
     166                *(output+connector[j]) = *(input + i) ; 
     167                j++ ; 
     168              } 
     169          } 
     170          else 
     171          { 
     172            switch(op) 
     173            { 
     174              case EReduction::sum : 
     175                for(int i=0, j=0; i<srcSize; i++) 
     176                  if (mask[i])  
     177                  { 
     178                    if (!std::isnan(*(input + i))) 
     179                    { 
     180                      if (*(count+connector[j])==0) *(output+connector[j]) = *(input + i) ; 
     181                      else *(output+connector[j]) += *(input + i) ; 
     182                      (*(count+connector[j]))++ ; 
     183                    }  
     184                    j++ ; 
     185                  } 
     186                break ; 
     187              case EReduction::min : 
     188                for(int i=0, j=0; i<srcSize; i++) 
     189                  if (mask[i])  
     190                  { 
     191                    if (!std::isnan(*(input + i))) 
     192                    { 
     193                      if (*(count+connector[j])==0) *(output+connector[j]) = *(input + i) ; 
     194                      else *(output+connector[j]) = std::min(*(output+connector[j]),*(input + i)) ; 
     195                      (*(count+connector[j]))++ ; 
     196                    }  
     197                    j++ ; 
     198                  } 
     199                break ; 
     200              case EReduction::max : 
     201                for(int i=0, j=0; i<srcSize; i++) 
     202                  if (mask[i])  
     203                  { 
     204                    if (!std::isnan(*(input + i))) 
     205                    { 
     206                      if (*(count+connector[j])==0) *(output+connector[j]) = *(input + i) ; 
     207                      else *(output+connector[j]) = std::max(*(output+connector[j]),*(input + i)) ; 
     208                      (*(count+connector[j]))++ ; 
     209                    }  
     210                    j++ ; 
     211                  } 
     212                break ; 
     213              case EReduction::average : 
     214                for(int i=0, j=0; i<srcSize; i++) 
     215                  if (mask[i])  
     216                  { 
     217                    if (!std::isnan(*(input + i))) 
     218                    { 
     219                      if (*(count+connector[j])==0) *(output+connector[j]) = *(input + i) ; 
     220                      else *(output+connector[j]) += *(input + i) ; 
     221                      (*(count+connector[j]))++ ; 
     222                    }  
     223                    j++ ; 
     224                  } 
     225                break ; 
     226              default : 
     227                ERROR("CGathererConnector::transfer", 
     228                     <<"reduction operator "<<(int)op<<" is not defined for this operation") ; 
     229                break ; 
     230            } 
     231          } 
    91232 
    92233        } 
    … …  
    101242            if (mask[i])  
    102243            { 
    103               (*(connectors-1))->transfer(rank, connectors-1, nConnectors-1, in, output+connector[j]*dstSliceSize) ; // the multiplication must be avoid in further optimization 
     244              (*(connectors-1))->transfer(rank, connectors-1, nConnectors-1, in, output+connector[j]*dstSliceSize, op, count+connector[j]*dstSliceSize) ; // the multiplication must be avoid in further optimization 
    104245              j++ ; 
    105246            } 
    … …  
    149290 
    150291      template<typename T> 
    151       void transfer(map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, T missingValue) 
    152       { 
    153         transfer(1, 1, dataIn, dataOut, missingValue); 
    154       } 
    155        
    156       template<typename T> 
    157       void transfer(int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, T missingValue) 
    158       { 
    159          transfer(1, sizeT, dataIn, dataOut, missingValue) ; 
    160       } 
    161  
    162       template<typename T> 
    163       void transfer(int repeat , int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, T missingValue) 
     292      void transfer(map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, T missingValue, EReduction op = EReduction::none) 
     293      { 
     294        transfer(1, 1, dataIn, dataOut, missingValue, op); 
     295      } 
     296       
     297      template<typename T> 
     298      void transfer(int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, T missingValue, EReduction op = EReduction::none) 
     299      { 
     300         transfer(1, sizeT, dataIn, dataOut, missingValue, op) ; 
     301      } 
     302 
     303      template<typename T> 
     304      void transfer(int repeat , int sizeT, map<int, CArray<T,1>>& dataIn, CArray<T,1>& dataOut, T missingValue, EReduction op = EReduction::none) 
    164305      { 
    165306        dataOut.resize(repeat*dstSize_*sizeT) ; 
    166307        dataOut=missingValue ; 
    167         transfer(repeat, sizeT, dataIn, dataOut) ; 
    168       } 
    169  
    170       template<typename T> 
    171       void transfer(CEventServer& event, int sizeT, CArray<T,1>& dataOut) 
     308        transfer(repeat, sizeT, dataIn, dataOut, op) ; 
     309      } 
     310 
     311      template<typename T> 
     312      void transfer(CEventServer& event, int sizeT, CArray<T,1>& dataOut, EReduction op = EReduction::none) 
    172313      { 
    173314        map<int, CArray<T,1>> dataIn ; 
    … …  
    177318          (*subEvent.buffer) >> data ; 
    178319        } 
    179         transfer(1, sizeT, dataIn, dataOut) ; 
    180       } 
    181        
    182       template<typename T> 
    183       void transfer(CEventServer& event, CArray<T,1>& dataOut) 
    184       { 
    185         transfer(event, 1, dataOut) ; 
    186       } 
    187  
    188       template<typename T> 
    189       void transfer(CEventServer& event, int sizeT, CArray<T,1>& dataOut, T missingValue) 
     320        transfer(1, sizeT, dataIn, dataOut, op) ; 
     321      } 
     322       
     323      template<typename T> 
     324      void transfer(CEventServer& event, CArray<T,1>& dataOut, EReduction op = EReduction::none) 
     325      { 
     326        transfer(event, 1, dataOut, op) ; 
     327      } 
     328 
     329      template<typename T> 
     330      void transfer(CEventServer& event, int sizeT, CArray<T,1>& dataOut, T missingValue, EReduction op = EReduction::none) 
    190331      { 
    191332        map<int, CArray<T,1>> dataIn ; 
    … …  
    195336          (*subEvent.buffer) >> data ; 
    196337        } 
    197         transfer(1, sizeT, dataIn, dataOut, missingValue) ; 
    198       } 
    199  
    200       template<typename T> 
    201       void transfer(CEventServer& event, CArray<T,1>& dataOut, T missingValue) 
     338        transfer(1, sizeT, dataIn, dataOut, missingValue, op) ; 
     339      } 
     340 
     341      template<typename T> 
     342      void transfer(CEventServer& event, CArray<T,1>& dataOut, T missingValue, EReduction op = EReduction::none) 
    202343      { 
    203344        map<int, CArray<T,1>> dataIn ; 
    … …  
    207348          (*subEvent.buffer) >> data ; 
    208349        } 
    209         transfer(1, 1, dataIn, dataOut, missingValue) ; 
     350        transfer(1, 1, dataIn, dataOut, missingValue, op) ; 
    210351      } 
    211352 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/grid_gatherer_connector.hpp

    r2267 r2291  
    1010#include "context_client.hpp" 
    1111#include "gatherer_connector.hpp" 
     12#include "reduction_types.hpp" 
    1213 
    1314 
    … …  
    3031 
    3132      template<typename T>  
    32       void transfer(const map<int, CArray<T,1>>& input, CArray<T,1>& output) 
     33      void transfer(const map<int, CArray<T,1>>& input, CArray<T,1>& output, EReduction op = EReduction::none) 
    3334      { 
    3435        int n = elementsConnector_.size()-1 ; 
    35        shared_ptr<CGathererConnector>* connector = elementsConnector_.data() + n ; 
     36        shared_ptr<CGathererConnector>* connector = elementsConnector_.data() + n ; 
    3637        output.resize(dstSize_) ; 
    37         for(auto& rankDataIn : input)  
     38        
     39        if (op == EReduction::none) 
     40          for(auto& rankDataIn : input)  
     41            elementsConnector_[n]->transfer(rankDataIn.first, connector, n, rankDataIn.second.dataFirst(), output.dataFirst()) ; 
     42        else 
    3843        { 
    39           elementsConnector_[n]->transfer(rankDataIn.first, connector, n, rankDataIn.second.dataFirst(), output.dataFirst()) ; 
     44          T defaultValue = std::numeric_limits<T>::quiet_NaN(); 
     45          vector<int> count(dstSize_,0) ; 
     46          for(auto& rankDataIn : input)  
     47            elementsConnector_[n]->transfer(rankDataIn.first, connector, n, rankDataIn.second.dataFirst(), output.dataFirst(), op, count.data()) ; 
     48 
     49          for(int i=0;i<dstSize_;i++) if (count[i]==0) output(i)=defaultValue ; 
    4050        } 
    4151      }  
    … …  
    5565 
    5666      template<typename T> 
    57       void transfer(CEventServer& event, CArray<T,1>& dataOut) 
     67      void transfer(CEventServer& event, CArray<T,1>& dataOut, EReduction op = EReduction::none) 
    5868      { 
    5969        map<int, CArray<T,1>> dataIn ; 
    … …  
    6373          (*subEvent.buffer) >> data ; 
    6474        } 
    65         transfer(dataIn, dataOut) ; 
     75        transfer(dataIn, dataOut, op) ; 
    6676      } 
    6777       

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/grid_remote_connector.cpp

    r2267 r2291  
    114114  *         some redondant ranks can be removed from the elements_ map. 
    115115  */ 
    116   void CGridRemoteConnector::computeConnector(void) 
    117   { 
    118     computeViewDistribution() ; 
    119     computeConnectorMethods() ; 
    120     computeRedondantRanks() ;  
    121     for(auto& rank : rankToRemove_) 
    122       for(auto& element : elements_) element.erase(rank) ; 
    123   } 
     116  void CGridRemoteConnector::computeConnector(bool eliminateRedundant) 
     117  { 
     118    if (eliminateRedundant) 
     119    { 
     120      computeViewDistribution() ; 
     121      computeConnectorMethods() ; 
     122      computeRedondantRanks() ;  
     123      for(auto& rank : rankToRemove_) 
     124        for(auto& element : elements_) element.erase(rank) ; 
     125    } 
     126    else 
     127    { 
     128       computeViewDistribution() ; 
     129       computeConnectorRedundant() ; 
     130    } 
     131  } 
     132 
     133/** 
     134  * \brief Compute the connector, i.e. compute the \b elements_ attribute.  
     135  * \detail In this routine we don't eliminate redundant cells as it it performed in  
     136  *         computeConnectorMethods. It can be usefull to perform reduce operation over processes. 
     137            In future, some optimisation could be done considering full redondance of the  
     138            source view or the destination view. 
     139  */ 
     140  void CGridRemoteConnector::computeConnectorRedundant(void) 
     141  { 
     142    vector<shared_ptr<CLocalView>> srcView ; 
     143    vector<shared_ptr<CDistributedView>> dstView ; 
     144    vector<int> indElements ; 
     145    elements_.resize(srcView_.size()) ; 
     146     
     147    bool srcViewsNonDistributed=true ; // not usefull now but later for optimization 
     148    for(int i=0;i<srcView_.size();i++) srcViewsNonDistributed = srcViewsNonDistributed && !isSrcViewDistributed_[i]  ; 
     149     
     150    bool dstViewsNonDistributed=true ;  // not usefull now but later for optimization 
     151    for(int i=0;i<dstView_.size();i++) dstViewsNonDistributed = dstViewsNonDistributed && !isDstViewDistributed_[i] ; 
     152     
     153    for(int i=0;i<srcView_.size();i++)  
     154    { 
     155      srcView.push_back(srcView_[i]) ; 
     156      dstView.push_back(dstView_[i]) ; 
     157      indElements.push_back(i) ; 
     158    } 
     159 
     160    computeGenericMethod(srcView, dstView, indElements) ; 
     161     
     162    map<int,bool> ranks ;   
     163    for(auto& it : elements_[indElements[0]])  
     164    { 
     165      if (it.second.numElements()==0) ranks[it.first] = false ; 
     166      else  ranks[it.first] = true ; 
     167    } 
     168    
     169  } 
     170 
     171 
    124172/** 
    125173  * \brief Compute the connector, i.e. compute the \b elements_ attribute.  
    … …  
    426474  } 
    427475 
    428    
     476 
    429477 /** 
    430478  * \brief Generic method the compute the grid remote connector. Only distributed elements are specifed in the source view and remote view.  
    … …  
    650698    } 
    651699  } 
    652  
     700   
    653701} 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/grid_remote_connector.hpp

    r2267 r2291  
    2121      CGridRemoteConnector(vector<shared_ptr<CLocalView>>& srcView, vector<shared_ptr<CLocalView>>& dstView, MPI_Comm localComm, int remoteSize) ; 
    2222      void computeViewDistribution(void) ; 
    23       void computeConnector(void) ; 
     23      void computeConnector(bool eliminateRedundant=true) ; 
    2424      void computeConnectorMethods(void) ; 
     25      void computeConnectorRedundant(void) ; 
    2526      void computeGenericMethod(vector<shared_ptr<CLocalView>>& srcView, vector<shared_ptr<CDistributedView>>& dstView, vector<int>& indElements) ; 
    2627      void computeSrcDstNonDistributed(int i, map<int,bool>& ranks) ; 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/grid_transform_connector.cpp

    r2267 r2291  
    44#include "gatherer_connector.hpp" 
    55#include "grid_remote_connector.hpp" 
     6#include "grid_redondant_remote_connector.hpp" 
    67 
    78 
    89namespace xios 
    910{ 
    10   void CGridTransformConnector::computeConnector(void) 
     11  void CGridTransformConnector::computeConnector(bool eliminateRedundant) 
    1112  { 
    1213    int commSize ; 
    … …  
    1718 
    1819    auto remoteConnector = make_shared<CGridRemoteConnector>(srcViews_, remoteViews_, localComm_, commSize) ;   
    19     remoteConnector->computeConnector() ; 
     20    remoteConnector->computeConnector(eliminateRedundant) ;  
    2021     
    2122    vector<shared_ptr<CDistributedElement>> sendElements(nElements) ; 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/grid_transform_connector.hpp

    r2267 r2291  
    88#include "grid_scatterer_connector.hpp" 
    99#include "grid_gatherer_connector.hpp" 
     10#include "reduction_types.hpp" 
    1011#include "mpi.hpp" 
    1112 
    … …  
    2021      CGridTransformConnector(vector<shared_ptr<CLocalView>> srcViews, vector<shared_ptr<CLocalView>> remoteViews, MPI_Comm localComm)  
    2122                          : srcViews_(srcViews), remoteViews_(remoteViews), localComm_(localComm)  
    22                           { computeConnector();} 
     23                          { } 
    2324     
    24       void computeConnector(void) ;  
     25      void computeConnector(bool eliminateRedundant=true) ;  
    2526    protected: 
    2627     MPI_Comm localComm_ ; 
    … …  
    3637    public: 
    3738      template<typename T>  
    38       void transfer(const CArray<T,1>& dataIn, CArray<T,1>& dataOut) 
     39      void transfer(const CArray<T,1>& dataIn, CArray<T,1>& dataOut, EReduction op = EReduction::none) 
    3940      { 
    4041        map<int,CArray<T,1>> tmpArrayIn ; 
    … …  
    6263         
    6364        const double nanValue = std::numeric_limits<double>::quiet_NaN(); 
    64         gridGathererConnector_->transfer(tmpArrayOut, dataOut, nanValue) ; 
     65 
     66        if (op == EReduction::none) gridGathererConnector_->transfer(tmpArrayOut, dataOut, nanValue) ; 
     67        else gridGathererConnector_->transfer(tmpArrayOut, dataOut, op) ; 
    6568      } 
    6669   

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/reduce_transform_connector.hpp

    r2267 r2291  
    6868      { 
    6969        vector<bool> touched(repeat* dstSlice, false) ; 
    70         int pos=0 ; 
    71         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    72         { 
    73           const double* in = input; 
    74           double* out = output ; 
    75           int k=0 ; 
     70        int posr=0 ; 
     71        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     72        { 
     73          const double* in = input; 
     74          double* out = output ; 
     75          int k=0 ; 
     76          int pos = posr ; 
    7677          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    7778          { 
    … …  
    9798 
    9899        output = dataOut.dataFirst()  ; 
    99         pos=0 ; 
    100         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    101         { 
    102           const double* in = input; 
    103           double* out = output ; 
    104           int k=0 ; 
    105           int pos=0 ; 
     100        posr=0 ; 
     101        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     102        { 
     103          const double* in = input; 
     104          double* out = output ; 
     105          int k=0 ; 
     106          int pos=posr ; 
    106107          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    107108              for(int j=0 ; j<nSrc_[i] ; j++,k++) 
    … …  
    144145      { 
    145146        vector<bool> touched(repeat* dstSlice, false) ; 
    146         int pos=0 ; 
    147         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    148         { 
    149           const double* in = input; 
    150           double* out = output ; 
    151           int k=0 ; 
     147        int posr=0 ; 
     148        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     149        { 
     150          const double* in = input; 
     151          double* out = output ; 
     152          int k=0 ; 
     153          int pos=posr ; 
    152154          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    153155          { 
    … …  
    173175 
    174176        output = dataOut.dataFirst()  ; 
    175         pos=0 ; 
    176         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    177         { 
    178           const double* in = input; 
    179           double* out = output ; 
    180           int k=0 ; 
    181           int pos=0 ; 
     177        posr=0 ; 
     178        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     179        { 
     180          const double* in = input; 
     181          double* out = output ; 
     182          int k=0 ; 
     183          int pos=posr ; 
    182184          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    183185              for(int j=0 ; j<nSrc_[i] ; j++,k++) 
    … …  
    220222      { 
    221223        vector<bool> touched(repeat* dstSlice, false) ; 
    222         int pos=0 ; 
    223         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    224         { 
    225           const double* in = input; 
    226           double* out = output ; 
    227           int k=0 ; 
     224        int posr=0 ; 
     225        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     226        { 
     227          const double* in = input; 
     228          double* out = output ; 
     229          int k=0 ; 
     230          int pos=posr ; 
    228231          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    229232          { 
    … …  
    249252 
    250253        output = dataOut.dataFirst()  ; 
    251         pos=0 ; 
    252         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    253         { 
    254           const double* in = input; 
    255           double* out = output ; 
    256           int k=0 ; 
    257           int pos=0 ; 
     254        posr=0 ; 
     255        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     256        { 
     257          const double* in = input; 
     258          double* out = output ; 
     259          int k=0 ; 
     260          int pos=posr ; 
    258261          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    259262              for(int j=0 ; j<nSrc_[i] ; j++,k++) 
    … …  
    296299      { 
    297300        vector<int> touched(repeat* dstSlice, 0) ; 
    298         int pos=0 ; 
    299         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    300         { 
    301           const double* in = input; 
    302           double* out = output ; 
    303           int k=0 ; 
     301        int posr=0 ; 
     302        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     303        { 
     304          const double* in = input; 
     305          double* out = output ; 
     306          int k=0 ; 
     307          int pos=posr ; 
    304308          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    305309          { 
    … …  
    331335 
    332336        output = dataOut.dataFirst()  ; 
    333         pos=0 ; 
    334         for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, pos+=dstSlice) 
    335         { 
    336           const double* in = input; 
    337           double* out = output ; 
    338           int k=0 ; 
    339           int pos=0 ; 
     337        posr=0 ; 
     338        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice, posr+=dstSlice) 
     339        { 
     340          const double* in = input; 
     341          double* out = output ; 
     342          int k=0 ; 
     343          int pos=posr ; 
    340344          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT) 
    341345              for(int j=0 ; j<nSrc_[i] ; j++,k++) 
Note: See TracChangeset for help on using the changeset viewer.