Zoltan2
Zoltan2_MetricUtility.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 
50 #ifndef ZOLTAN2_METRICFUNCTIONS_HPP
51 #define ZOLTAN2_METRICFUNCTIONS_HPP
52 
53 #include <Zoltan2_StridedData.hpp>
54 
55 namespace Zoltan2{
56 
58 // Namespace methods to compute metric values
60 
70 template <typename scalar_t>
71  void getStridedStats(const ArrayView<scalar_t> &v, int stride,
72  int offset, scalar_t &min, scalar_t &max, scalar_t &sum)
73 {
74  if (v.size() < 1) return;
75  min = max = sum = v[offset];
76 
77  for (int i=offset+stride; i < v.size(); i += stride){
78  if (v[i] < min) min = v[i];
79  else if (v[i] > max) max = v[i];
80  sum += v[i];
81  }
82 }
83 
92 template <typename scalar_t>
93  void getStridedStats(const ArrayView<scalar_t> &v, int stride,
94  int offset, scalar_t &max, scalar_t &sum)
95 {
96  if (v.size() < 1) return;
97  max = sum = v[offset];
98 
99  for (int i=offset+stride; i < v.size(); i += stride){
100  if (v[i] > max) max = v[i];
101  sum += v[i];
102  }
103 }
104 
123 template <typename scalar_t, typename lno_t, typename part_t>
125  const RCP<const Environment> &env,
126  part_t numberOfParts,
127  const ArrayView<const part_t> &parts,
128  const ArrayView<StridedData<lno_t, scalar_t> > &vwgts,
129  multiCriteriaNorm mcNorm,
130  scalar_t *weights)
131 {
132  env->localInputAssertion(__FILE__, __LINE__, "parts or weights",
133  numberOfParts > 0 && vwgts.size() > 0, BASIC_ASSERTION);
134 
135  int numObjects = parts.size();
136  int vwgtDim = vwgts.size();
137 
138  memset(weights, 0, sizeof(scalar_t) * numberOfParts);
139 
140  if (numObjects == 0)
141  return;
142 
143  if (vwgtDim == 0) {
144  for (lno_t i=0; i < numObjects; i++){
145  weights[parts[i]]++;
146  }
147  }
148  else if (vwgtDim == 1){
149  for (lno_t i=0; i < numObjects; i++){
150  weights[parts[i]] += vwgts[0][i];
151  }
152  }
153  else{
154  switch (mcNorm){
156  for (int wdim=0; wdim < vwgtDim; wdim++){
157  for (lno_t i=0; i < numObjects; i++){
158  weights[parts[i]] += vwgts[wdim][i];
159  }
160  } // next weight index
161  break;
162 
164  for (lno_t i=0; i < numObjects; i++){
165  scalar_t ssum = 0;
166  for (int wdim=0; wdim < vwgtDim; wdim++)
167  ssum += (vwgts[wdim][i] * vwgts[wdim][i]);
168  weights[parts[i]] += sqrt(ssum);
169  }
170  break;
171 
173  for (lno_t i=0; i < numObjects; i++){
174  scalar_t max = 0;
175  for (int wdim=0; wdim < vwgtDim; wdim++)
176  if (vwgts[wdim][i] > max)
177  max = vwgts[wdim][i];
178  weights[parts[i]] += max;
179  }
180  break;
181 
182  default:
183  env->localBugAssertion(__FILE__, __LINE__, "invalid norm", false,
185  break;
186  }
187  }
188 }
189 
220 template <typename scalar_t, typename part_t>
222  part_t numExistingParts, // Max Part ID + 1
223  part_t targetNumParts, // comm.size() or requested global # parts from soln
224  const scalar_t *psizes,
225  scalar_t sumVals,
226  const scalar_t *vals,
227  scalar_t &min,
228  scalar_t &max,
229  scalar_t &avg)
230 {
231  min = sumVals;
232  max = avg = 0;
233 
234  if (sumVals <= 0 || targetNumParts < 1 || numExistingParts < 1)
235  return;
236 
237  if (targetNumParts==1) {
238  min = max = avg = 0; // 0 imbalance
239  return;
240  }
241 
242  if (!psizes){
243  scalar_t target = sumVals / targetNumParts;
244  for (part_t p=0; p < numExistingParts; p++){
245  scalar_t diff = vals[p] - target;
246  scalar_t adiff = (diff >= 0 ? diff : -diff);
247  scalar_t tmp = diff / target;
248  scalar_t atmp = adiff / target;
249  avg += atmp;
250  if (tmp > max) max = tmp;
251  if (tmp < min) min = tmp;
252  }
253  part_t emptyParts = targetNumParts - numExistingParts;
254  if (emptyParts > 0){
255  if (max < 1.0)
256  max = 1.0; // target divided by target
257  avg += emptyParts;
258  }
259  }
260  else{
261  for (part_t p=0; p < targetNumParts; p++){
262  if (psizes[p] > 0){
263  if (p < numExistingParts){
264  scalar_t target = sumVals * psizes[p];
265  scalar_t diff = vals[p] - target;
266  scalar_t adiff = (diff >= 0 ? diff : -diff);
267  scalar_t tmp = diff / target;
268  scalar_t atmp = adiff / target;
269  avg += atmp;
270  if (tmp > max) max = tmp;
271  if (tmp < min) min = tmp;
272  }
273  else{
274  if (max < 1.0)
275  max = 1.0; // target divided by target
276  avg += 1.0;
277  }
278  }
279  }
280  }
281 
282  avg /= targetNumParts;
283 }
284 
318 template <typename scalar_t, typename part_t>
320  part_t numExistingParts,
321  part_t targetNumParts,
322  int numSizes,
323  ArrayView<ArrayRCP<scalar_t> > psizes,
324  scalar_t sumVals,
325  const scalar_t *vals,
326  scalar_t &min,
327  scalar_t &max,
328  scalar_t &avg)
329 {
330  min = sumVals;
331  max = avg = 0;
332 
333  if (sumVals <= 0 || targetNumParts < 1 || numExistingParts < 1)
334  return;
335 
336  if (targetNumParts==1) {
337  min = max = avg = 0; // 0 imbalance
338  return;
339  }
340 
341  bool allUniformParts = true;
342  for (int i=0; i < numSizes; i++){
343  if (psizes[i].size() != 0){
344  allUniformParts = false;
345  break;
346  }
347  }
348 
349  if (allUniformParts){
350  computeImbalances<scalar_t, part_t>(numExistingParts, targetNumParts, NULL,
351  sumVals, vals, min, max, avg);
352  return;
353  }
354 
355  double uniformSize = 1.0 / targetNumParts;
356  std::vector<double> sizeVec(numSizes);
357  for (int i=0; i < numSizes; i++){
358  sizeVec[i] = uniformSize;
359  }
360 
361  for (part_t p=0; p < numExistingParts; p++){
362 
363  // If we have objects in parts that should have 0 objects,
364  // we don't compute an imbalance. It means that other
365  // parts have too few objects, and the imbalance will be
366  // picked up there.
367 
368  if (p >= targetNumParts)
369  break;
370 
371  // Vector of target amounts: T
372 
373  double targetNorm = 0;
374  for (int i=0; i < numSizes; i++) {
375  if (psizes[i].size() > 0)
376  sizeVec[i] = psizes[i][p];
377  sizeVec[i] *= sumVals;
378  targetNorm += (sizeVec[i] * sizeVec[i]);
379  }
380  targetNorm = sqrt(targetNorm);
381 
382  // If part is supposed to be empty, we don't compute an
383  // imbalance. Same argument as above.
384 
385  if (targetNorm > 0){
386 
387  // Vector of actual amounts: A
388 
389  std::vector<double> actual(numSizes);
390  double actualNorm = 0.;
391  for (int i=0; i < numSizes; i++) {
392  actual[i] = vals[p] * -1.0;
393  actual[i] += sizeVec[i];
394  actualNorm += (actual[i] * actual[i]);
395  }
396  actualNorm = sqrt(actualNorm);
397 
398  // |A - T| / |T|
399 
400  scalar_t imbalance = actualNorm / targetNorm;
401 
402  if (imbalance < min)
403  min = imbalance;
404  else if (imbalance > max)
405  max = imbalance;
406  avg += imbalance;
407  }
408  }
409 
410  part_t numEmptyParts = 0;
411 
412  for (part_t p=numExistingParts; p < targetNumParts; p++){
413  bool nonEmptyPart = false;
414  for (int i=0; !nonEmptyPart && i < numSizes; i++)
415  if (psizes[i].size() > 0 && psizes[i][p] > 0.0)
416  nonEmptyPart = true;
417 
418  if (nonEmptyPart){
419  // The partition has no objects for this part, which
420  // is supposed to be non-empty.
421  numEmptyParts++;
422  }
423  }
424 
425  if (numEmptyParts > 0){
426  avg += numEmptyParts;
427  if (max < 1.0)
428  max = 1.0; // target divided by target
429  }
430 
431  avg /= targetNumParts;
432 }
433 
436 template <typename scalar_t>
437  scalar_t normedWeight(ArrayView <scalar_t> weights,
438  multiCriteriaNorm norm)
439 {
440  size_t dim = weights.size();
441  if (dim == 0)
442  return 0.0;
443  else if (dim == 1)
444  return weights[0];
445 
446  scalar_t nweight = 0;
447 
448  switch (norm){
451  for (size_t i=0; i <dim; i++)
452  nweight += weights[i];
453 
454  break;
455 
457  for (size_t i=0; i <dim; i++)
458  nweight += (weights[i] * weights[i]);
459 
460  nweight = sqrt(nweight);
461 
462  break;
463 
466  nweight = weights[0];
467  for (size_t i=1; i <dim; i++)
468  if (weights[i] > nweight)
469  nweight = weights[i];
470 
471  break;
472 
473  default:
474  Environment env; // default environment
475  env.localBugAssertion(__FILE__, __LINE__, "invalid norm", false,
477  break;
478  }
479 
480  return nweight;
481 }
482 
485 template<typename lno_t, typename scalar_t>
487  lno_t idx, multiCriteriaNorm norm)
488 {
489  size_t dim = weights.size();
490  if (dim < 1)
491  return 0;
492 
493  Array<scalar_t> vec(dim, 1.0);
494  for (size_t i=0; i < dim; i++)
495  if (weights[i].size() > 0)
496  vec[dim] = weights[i][idx];
497 
498  return normedWeight(vec, norm);
499 }
500 
501 } //namespace Zoltan2
502 
503 #endif
fast typical checks for valid arguments
void normedPartWeights(const RCP< const Environment > &env, part_t numberOfParts, const ArrayView< const part_t > &parts, const ArrayView< StridedData< lno_t, scalar_t > > &vwgts, multiCriteriaNorm mcNorm, scalar_t *weights)
Compute the total weight in each part on this process.
static ArrayRCP< ArrayRCP< zscalar_t > > weights
void getStridedStats(const ArrayView< scalar_t > &v, int stride, int offset, scalar_t &min, scalar_t &max, scalar_t &sum)
Find min, max and sum of metric values.
dictionary vals
Definition: xml2dox.py:186
void localBugAssertion(const char *file, int lineNum, const char *msg, bool ok, AssertionLevel level) const
Test for valid library behavior on local process only.
The StridedData class manages lists of weights or coordinates.
The user parameters, debug, timing and memory profiling output objects, and error checking methods...
multiCriteriaNorm
Enumerator used in code for multicriteria norm choice.
void computeImbalances(part_t numExistingParts, part_t targetNumParts, const scalar_t *psizes, scalar_t sumVals, const scalar_t *vals, scalar_t &min, scalar_t &max, scalar_t &avg)
Compute the imbalance.
This file defines the StridedData class.
scalar_t normedWeight(ArrayView< scalar_t > weights, multiCriteriaNorm norm)
Compute the norm of the vector of weights.