source: XIOS/dev/dev_ym/XIOS_COUPLING/src/distribution/reduce_transform_connector.hpp @ 2291

Last change on this file since 2291 was 2291, checked in by ymipsl, 2 years ago

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
  • Property svn:executable set to *
File size: 12.0 KB
Line 
1#ifndef __REDUCE_TRANSFORM_CONNECTOR_HPP__
2#define __REDUCE_TRANSFORM_CONNECTOR_HPP__
3
4#include "xios_spl.hpp"
5#include "array_new.hpp"
6#include "local_view.hpp"
7#include "reduction_types.hpp"
8#include <algorithm>
9
10
11
12namespace xios
13{
14 
15  class CReduceTransformConnector
16  {
17   
18    private:
19      shared_ptr<CLocalView> srcView_;
20      shared_ptr<CLocalView> dstView_;
21
22      vector<int> connector_;  //  sizeof sum(nWeights_) 
23      vector<int> nSrc_ ;  //  sizeof dstSize_
24      int srcSize_ ;
25      int dstSize_ ;
26      bool detectMissingValue_ ;
27      EReduction type_ ;
28     
29      typedef void (CReduceTransformConnector::* transferPtr)(int, int, const CArray<double,1>& , CArray<double,1>&) ;
30      transferPtr transfer_ ;
31     
32     void computeConnector(unordered_map<int, std::vector<int>>& indexMap) ;
33     
34    public:
35
36    CReduceTransformConnector(shared_ptr<CLocalView> srcView, shared_ptr<CLocalView> dstView, EReduction op, unordered_map<int, std::vector<int>>& indexMap, 
37                              bool detectMissingValue=true) ;
38 
39    void transfer(const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
40    {
41      transfer(1,1,dataIn, dataOut) ;
42    }
43        void transfer(int sizeT, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
44    { 
45      transfer(1,sizeT, dataIn, dataOut) ;
46    }
47
48    void transfer(int repeat, int sizeT, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
49    {
50      (this->*transfer_)(repeat, sizeT, dataIn, dataOut) ;
51    }
52   
53    private :
54   
55    void transferSum(int repeat, int sizeT, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
56    {
57      double defaultValue = std::numeric_limits<double>::quiet_NaN();
58
59      int dstSlice = dstSize_*sizeT ;
60      int srcSlice = srcSize_*sizeT ;
61      dataOut.resize(repeat* dstSlice) ;
62      dataOut=0 ;
63           
64      const double* input = dataIn.dataFirst()  ;
65      double* output = dataOut.dataFirst()  ;
66
67      if (detectMissingValue_)
68      {
69        vector<bool> touched(repeat* dstSlice, false) ;
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 ;
77          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
78          {
79            if (nSrc_[i]==0) for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
80            else 
81            {
82              for(int j=0 ; j<nSrc_[i] ; j++,k++)
83                for(int l=0; l<sizeT; l++) 
84                {
85                  if (! std::isnan(in[connector_[k]*sizeT+l]) ) 
86                  {
87                    if (touched[pos+l]) out[l] += in[connector_[k]*sizeT+l] ;
88                    else 
89                    {
90                      touched[pos+l] = true ;
91                      out[l] = in[connector_[k]*sizeT+l] ;
92                    }
93                  }
94                }
95            }
96          }
97        }
98
99        output = dataOut.dataFirst()  ;
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 ;
107          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
108              for(int j=0 ; j<nSrc_[i] ; j++,k++)
109                for(int l=0; l<sizeT; l++) if (!touched[pos+l]) out[l] = defaultValue ;
110        }
111
112      }
113      else
114      {
115        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice)
116        {
117          const double* in = input;
118          double* out = output ;
119          int k=0 ;
120          for(int i=0; i<dstSize_; i++, out+=sizeT)
121            if (nSrc_[i]==0)
122              for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
123            else 
124            {
125              for(int j=0 ; j<nSrc_[i] ; j++,k++)
126                for(int l=0; l<sizeT; l++) out[l] += in[connector_[k]*sizeT+l] ;
127            }
128        }
129      }
130    }
131
132    void transferMin(int repeat, int sizeT, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
133    {
134      double defaultValue = std::numeric_limits<double>::quiet_NaN();
135
136      int dstSlice = dstSize_*sizeT ;
137      int srcSlice = srcSize_*sizeT ;
138      dataOut.resize(repeat* dstSlice) ;
139      dataOut=0 ;
140           
141      const double* input = dataIn.dataFirst()  ;
142      double* output = dataOut.dataFirst()  ;
143
144      if (detectMissingValue_)
145      {
146        vector<bool> touched(repeat* dstSlice, false) ;
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 ;
154          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
155          {
156            if (nSrc_[i]==0) for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
157            else 
158            {
159              for(int j=0 ; j<nSrc_[i] ; j++,k++)
160                for(int l=0; l<sizeT; l++) 
161                {
162                  if (! std::isnan(in[connector_[k]*sizeT+l]) ) 
163                  {
164                    if (touched[pos+l]) out[l]=min(out[l],in[connector_[k]*sizeT+l]) ;
165                    else 
166                    {
167                      touched[pos+l] = true ;
168                      out[l] = in[connector_[k]*sizeT+l] ;
169                    }
170                  }
171                }
172            }
173          }
174        }
175
176        output = dataOut.dataFirst()  ;
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 ;
184          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
185              for(int j=0 ; j<nSrc_[i] ; j++,k++)
186                for(int l=0; l<sizeT; l++) if (!touched[pos+l]) out[l] = defaultValue ;
187        }
188
189      }
190      else
191      {
192        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice)
193        {
194          const double* in = input;
195          double* out = output ;
196          int k=0 ;
197          for(int i=0; i<dstSize_; i++, out+=sizeT)
198            if (nSrc_[i]==0)
199              for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
200            else 
201            {
202              for(int j=0 ; j<nSrc_[i] ; j++,k++)
203                for(int l=0; l<sizeT; l++) out[l]=min(out[l],in[connector_[k]*sizeT+l]) ;
204            }
205        }
206      }
207    }
208
209    void transferMax(int repeat, int sizeT, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
210    {
211      double defaultValue = std::numeric_limits<double>::quiet_NaN();
212
213      int dstSlice = dstSize_*sizeT ;
214      int srcSlice = srcSize_*sizeT ;
215      dataOut.resize(repeat* dstSlice) ;
216      dataOut=0 ;
217           
218      const double* input = dataIn.dataFirst()  ;
219      double* output = dataOut.dataFirst()  ;
220
221      if (detectMissingValue_)
222      {
223        vector<bool> touched(repeat* dstSlice, false) ;
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 ;
231          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
232          {
233            if (nSrc_[i]==0) for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
234            else 
235            {
236              for(int j=0 ; j<nSrc_[i] ; j++,k++)
237                for(int l=0; l<sizeT; l++) 
238                {
239                  if (! std::isnan(in[connector_[k]*sizeT+l]) ) 
240                  {
241                    if (touched[pos+l]) out[l]=max(out[l],in[connector_[k]*sizeT+l]) ;
242                    else 
243                    {
244                      touched[pos+l] = true ;
245                      out[l] = in[connector_[k]*sizeT+l] ;
246                    }
247                  }
248                }
249            }
250          }
251        }
252
253        output = dataOut.dataFirst()  ;
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 ;
261          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
262              for(int j=0 ; j<nSrc_[i] ; j++,k++)
263                for(int l=0; l<sizeT; l++) if (!touched[pos+l]) out[l] = defaultValue ;
264        }
265
266      }
267      else
268      {
269        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice)
270        {
271          const double* in = input;
272          double* out = output ;
273          int k=0 ;
274          for(int i=0; i<dstSize_; i++, out+=sizeT)
275            if (nSrc_[i]==0)
276              for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
277            else 
278            {
279              for(int j=0 ; j<nSrc_[i] ; j++,k++)
280                for(int l=0; l<sizeT; l++) out[l]=max(out[l],in[connector_[k]*sizeT+l]) ;
281            }
282        }
283      }
284    }
285
286    void transferAverage(int repeat, int sizeT, const CArray<double,1>& dataIn, CArray<double,1>& dataOut)
287    {
288      double defaultValue = std::numeric_limits<double>::quiet_NaN();
289
290      int dstSlice = dstSize_*sizeT ;
291      int srcSlice = srcSize_*sizeT ;
292      dataOut.resize(repeat* dstSlice) ;
293      dataOut=0 ;
294           
295      const double* input = dataIn.dataFirst()  ;
296      double* output = dataOut.dataFirst()  ;
297
298      if (detectMissingValue_)
299      {
300        vector<int> touched(repeat* dstSlice, 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 ;
308          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
309          {
310            if (nSrc_[i]==0) for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
311            else 
312            {
313              for(int j=0 ; j<nSrc_[i] ; j++,k++)
314                for(int l=0; l<sizeT; l++) 
315                {
316                  if (! std::isnan(in[connector_[k]*sizeT+l]) ) 
317                  {
318                    if (touched[pos+l]) out[l] += in[connector_[k]*sizeT+l] ;
319                    else 
320                    {
321                      out[l] = in[connector_[k]*sizeT+l] ;
322                    }
323                    touched[pos+l]++ ;
324
325                  }
326                }
327              for(int l=0; l<sizeT; l++) {
328                  if (touched[pos+l]!=0)
329                      out[l] /= touched[pos+l];
330              }
331
332            }
333          }
334        }
335
336        output = dataOut.dataFirst()  ;
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 ;
344          for(int i=0; i<dstSize_; i++, out+=sizeT, pos+=sizeT)
345              for(int j=0 ; j<nSrc_[i] ; j++,k++)
346                for(int l=0; l<sizeT; l++) 
347                  if (touched[pos+l]==0) out[l] = defaultValue ;
348        }
349
350      }
351      else
352      {
353        for(int r=0;r<repeat;r++, input+=srcSlice, output+=dstSlice)
354        {
355          const double* in = input;
356          double* out = output ;
357          int k=0 ;
358          for(int i=0; i<dstSize_; i++, out+=sizeT)
359            if (nSrc_[i]==0)
360              for(int l=0; l<sizeT; l++) out[l] = defaultValue ;
361            else 
362            {
363              for(int j=0 ; j<nSrc_[i] ; j++,k++)
364                for(int l=0; l<sizeT; l++) out[l] += in[connector_[k]*sizeT+l] ;
365            }
366         
367          out = output ;
368          for(int i=0; i<dstSize_; i++, out+=sizeT)
369            if (nSrc_[i]!=0)
370              for(int l=0; l<sizeT; l++) out[l]/=nSrc_[i];
371        }
372      }
373    }
374 
375  };
376
377}
378
379#endif
Note: See TracBrowser for help on using the repository browser.