remap_store_link.cc 8.73 KB
Newer Older
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
6
7
8
9
10
11
12
13
14
15
16
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/
Oliver Heidmann's avatar
Oliver Heidmann committed
17

18
#include <algorithm>
19
#include <array>
20

21
#include "dmemory.h"
22
#include "cdo_options.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23
24
25
#include "remap.h"
#include "remap_store_link.h"

26
static bool
27
compareAdds(const Addweight &a, const Addweight &b)
28
29
30
31
32
{
  return a.add < b.add;
}

static bool
33
compareAdds4(const Addweight4 &a, const Addweight4 &b)
34
35
36
37
{
  return a.add < b.add;
}

38
static int
39
qcompareAdds(const void *a, const void *b)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
  const size_t x = ((const Addweight *) a)->add;
  const size_t y = ((const Addweight *) b)->add;
43
  return ((x > y) - (x < y)) * 2 + (x > y) - (x < y);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
44
45
}

46
static int
47
qcompareAdds4(const void *a, const void *b)
48
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
  const size_t x = ((const Addweight4 *) a)->add;
  const size_t y = ((const Addweight4 *) b)->add;
51
  return ((x > y) - (x < y)) * 2 + (x > y) - (x < y);
52
53
}

54
static void
55
sortAddweights(size_t numWeights, Addweight *addweights)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
{
57
  size_t n;
58
  for (n = 1; n < numWeights; ++n)
59
    if (addweights[n].add < addweights[n - 1].add) break;
60
  if (n == numWeights) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61

62
  std::qsort(addweights, numWeights, sizeof(Addweight), qcompareAdds);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
64
}

65
static void
66
sortAddweights4(Addweight4 *addweights)
67
{
68
69
  unsigned n;
  for (n = 1; n < 4; ++n)
70
    if (addweights[n].add < addweights[n - 1].add) break;
71
  if (n == 4) return;
72

73
  std::qsort(addweights, 4, sizeof(Addweight4), qcompareAdds4);
74
75
}

76
void
77
sort_weights_n4(size_t *src_add, double *weights)
78
79
80
81
82
83
84
85
86
87
88
89
{
  constexpr size_t numWeights = 4;
  size_t n;
  for (n = 1; n < numWeights; ++n)
    if (src_add[n] < src_add[n - 1]) break;
  if (n == numWeights) return;

  std::array<Addweight, numWeights> addweights;

  for (n = 0; n < numWeights; ++n)
    {
      addweights[n].add = src_add[n];
90
      addweights[n].weight = weights[n];
91
92
93
94
95
96
97
    }

  std::sort(addweights.begin(), addweights.end(), compareAdds);

  for (n = 0; n < numWeights; ++n)
    {
      src_add[n] = addweights[n].add;
98
      weights[n] = addweights[n].weight;
99
100
101
    }
}

102
void
103
sort_weights(size_t numWeights, size_t *src_add, double *weights)
104
{
105
  size_t n;
106
  for (n = 1; n < numWeights; ++n)
107
    if (src_add[n] < src_add[n - 1]) break;
108
  if (n == numWeights) return;
109

110
  if (numWeights > 1)
111
    {
112
      std::vector<Addweight> addweights(numWeights);
113

114
      for (n = 0; n < numWeights; ++n)
115
        {
116
          addweights[n].add = src_add[n];
117
          addweights[n].weight = weights[n];
118
        }
119

120
      std::sort(addweights.begin(), addweights.end(), compareAdds);
121

122
      for (n = 0; n < numWeights; ++n)
123
124
        {
          src_add[n] = addweights[n].add;
125
          weights[n] = addweights[n].weight;
126
        }
127
    }
128
129
}

130
void
131
sort_weights_bicubic(size_t *src_add, double (&weights)[4][4])
132
{
133
  constexpr size_t numWeights = 4;
134
  size_t n;
135
  for (n = 1; n < numWeights; ++n)
136
    if (src_add[n] < src_add[n - 1]) break;
137
  if (n == numWeights) return;
138

139
  std::array<Addweight4, numWeights> addweights;
140

141
142
143
  for (n = 0; n < numWeights; ++n)
    {
      addweights[n].add = src_add[n];
144
      for (unsigned k = 0; k < 4; ++k) addweights[n].weight[k] = weights[n][k];
145
    }
146

147
  std::sort(addweights.begin(), addweights.end(), compareAdds4);
148

149
150
151
  for (n = 0; n < numWeights; ++n)
    {
      src_add[n] = addweights[n].add;
152
      for (unsigned k = 0; k < 4; ++k) weights[n][k] = addweights[n].weight[k];
153
154
155
    }
}

156
void
157
158
store_weightlinks(int lalloc, size_t numWeights, size_t *srch_add, double *weights, size_t cell_add,
                  std::vector<WeightLinks> &weightLinks)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
{
160
161
  weightLinks[cell_add].nlinks = 0;
  weightLinks[cell_add].offset = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162

163
  if (numWeights)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
      Addweight *addweights = nullptr;
166
      if (lalloc)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
        addweights = (Addweight *) Malloc(numWeights * sizeof(Addweight));
168
      else
169
        addweights = weightLinks[cell_add].addweights;
170

171
      for (size_t n = 0; n < numWeights; ++n)
172
173
174
175
        {
          addweights[n].add = srch_add[n];
          addweights[n].weight = weights[n];
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176

177
      if (numWeights > 1) sortAddweights(numWeights, addweights);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178

179
      weightLinks[cell_add].nlinks = numWeights;
180

181
      if (lalloc) weightLinks[cell_add].addweights = addweights;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184

185
void
186
store_weightlinks_bicubic(size_t *srch_add, double (&weights)[4][4], size_t cell_add, std::vector<WeightLinks4> &weightLinks)
187
{
188
189
  weightLinks[cell_add].nlinks = 0;
  weightLinks[cell_add].offset = 0;
190

191
  Addweight4 *addweights = weightLinks[cell_add].addweights;
192

193
194
195
196
197
  for (unsigned n = 0; n < 4; ++n)
    {
      addweights[n].add = srch_add[n];
      for (unsigned k = 0; k < 4; ++k) addweights[n].weight[k] = weights[n][k];
    }
198

199
  sortAddweights4(addweights);
200

201
  weightLinks[cell_add].nlinks = 4;
202
203
}

204
void
205
weightLinksToRemapLinks(int lalloc, size_t gridSize, std::vector<WeightLinks> &weightLinks, RemapVars &rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
{
207
  size_t nlinks = 0;
208
  for (size_t i = 0; i < gridSize; ++i)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209
    {
210
      if (weightLinks[i].nlinks)
211
        {
212
213
          weightLinks[i].offset = nlinks;
          nlinks += weightLinks[i].nlinks;
214
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
216
    }

217
218
  rv.max_links = nlinks;
  rv.num_links = nlinks;
219

220
  if (nlinks)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
    {
222
      auto num_wts = rv.num_wts;
223
224
      rv.src_cell_add.resize(nlinks);
      rv.tgt_cell_add.resize(nlinks);
225
      rv.wts.resize(nlinks * num_wts);
226
227
228
      auto &src_cell_adds = rv.src_cell_add;
      auto &tgt_cell_adds = rv.tgt_cell_add;
      auto &wts = rv.wts;
229

230
#ifdef _OPENMP
231
#pragma omp parallel for schedule(static) default(none) shared(src_cell_adds, tgt_cell_adds, wts, weightLinks, gridSize, num_wts)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
#endif
233
      for (size_t i = 0; i < gridSize; ++i)
234
        {
235
          const auto num_links = weightLinks[i].nlinks;
236
237
          if (num_links)
            {
238
              const auto offset = weightLinks[i].offset;
239
              Addweight *addweights = weightLinks[i].addweights;
240
241
242
              for (size_t ilink = 0; ilink < num_links; ++ilink)
                {
                  src_cell_adds[offset + ilink] = addweights[ilink].add;
243
                  tgt_cell_adds[offset + ilink] = i;
244
                  wts[(offset + ilink) * num_wts] = addweights[ilink].weight;
245
                }
246
            }
247
        }
248

249
      if (lalloc)
250
        {
251
          for (size_t i = 0; i < gridSize; ++i)
252
            {
253
              const auto num_links = weightLinks[i].nlinks;
254
              if (num_links) Free(weightLinks[i].addweights);
255
256
257
258
            }
        }
      else
        {
259
          Free(weightLinks[0].addweights);
260
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
261
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
262
}
263

264
void
265
weightLinks4ToRemapLinks(size_t gridSize, std::vector<WeightLinks4> &weightLinks, RemapVars &rv)
266
{
267
  size_t nlinks = 0;
268
  for (size_t i = 0; i < gridSize; ++i)
269
    {
270
      if (weightLinks[i].nlinks)
271
        {
272
273
          weightLinks[i].offset = nlinks;
          nlinks += weightLinks[i].nlinks;
274
        }
275
276
    }

277
278
  rv.max_links = nlinks;
  rv.num_links = nlinks;
279
  if (nlinks)
280
    {
281
282
283
      rv.src_cell_add.resize(nlinks);
      rv.tgt_cell_add.resize(nlinks);
      rv.wts.resize(4 * nlinks);
284
285
286
      auto &src_cell_adds = rv.src_cell_add;
      auto &tgt_cell_adds = rv.tgt_cell_add;
      auto &wts = rv.wts;
287

288
#ifdef _OPENMP
289
#pragma omp parallel for default(none) shared(src_cell_adds, tgt_cell_adds, wts, weightLinks, gridSize)
290
#endif
291
      for (size_t i = 0; i < gridSize; ++i)
292
        {
293
          const auto num_links = weightLinks[i].nlinks;
294
295
          if (num_links)
            {
296
297
              const auto offset = weightLinks[i].offset;
              const auto addweights = weightLinks[i].addweights;
298
299
300
              for (size_t ilink = 0; ilink < num_links; ++ilink)
                {
                  src_cell_adds[offset + ilink] = addweights[ilink].add;
301
                  tgt_cell_adds[offset + ilink] = i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
                  for (size_t k = 0; k < 4; ++k) wts[(offset + ilink) * 4 + k] = addweights[ilink].weight[k];
303
                }
304
            }
305
        }
306

307
      Free(weightLinks[0].addweights);
308
309
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310
311

void
312
weightLinksAlloc(size_t numNeighbors, size_t gridSize, std::vector<WeightLinks> &weightLinks)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313
{
314
  weightLinks[0].addweights = (Addweight *) Malloc(numNeighbors * gridSize * sizeof(Addweight));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
  for (size_t i = 1; i < gridSize; ++i) weightLinks[i].addweights = weightLinks[0].addweights + numNeighbors * i;
316
317
318
}

void
319
weightLinks4Alloc(size_t gridSize, std::vector<WeightLinks4> &weightLinks)
320
{
321
  weightLinks[0].addweights = (Addweight4 *) Malloc(4 * gridSize * sizeof(Addweight4));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
322
  for (size_t i = 1; i < gridSize; ++i) weightLinks[i].addweights = weightLinks[0].addweights + 4 * i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
323
}