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

File:
1 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 
Note: See TracChangeset for help on using the changeset viewer.