Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
57 // this utlity method is used to allocate more metrics in the metric array
58 // the array is composed of an array of ptrs to BaseClassMetric
59 // but each ptr is allocated to the specific derived metric type
60 // So the array can access the polymorphic hierarchy
61 //
62 // Note this is currently only relevant to EvaluatePartition and the
63 // GraphMetrics and ImbalanceMetrics calculations
64 template <typename metric_t, typename scalar_t>
65 RCP<metric_t> addNewMetric(const RCP<const Environment> &env,
66  ArrayRCP<RCP<BaseClassMetrics<scalar_t> > > &metrics)
67 {
68  metrics.resize(metrics.size() + 1); // increase array size by 1
69  metric_t * newMetric = new metric_t; // allocate
70  env->localMemoryAssertion(__FILE__,__LINE__,1,newMetric); // check errors
71  RCP<metric_t> newRCP = rcp(newMetric); // rcp of derived class
72  metrics[metrics.size()-1] = newRCP; // create the new rcp
73  return newRCP;
74 }
75 
77 // Namespace methods to compute metric values
79 
89 template <typename scalar_t>
90  void getStridedStats(const ArrayView<scalar_t> &v, int stride,
91  int offset, scalar_t &min, scalar_t &max, scalar_t &sum)
92 {
93  if (v.size() < 1) return;
94  min = max = sum = v[offset];
95 
96  for (int i=offset+stride; i < v.size(); i += stride){
97  if (v[i] < min) min = v[i];
98  else if (v[i] > max) max = v[i];
99  sum += v[i];
100  }
101 }
102 
111 template <typename scalar_t>
112  void getStridedStats(const ArrayView<scalar_t> &v, int stride,
113  int offset, scalar_t &max, scalar_t &sum)
114 {
115  if (v.size() < 1) return;
116  max = sum = v[offset];
117 
118  for (int i=offset+stride; i < v.size(); i += stride){
119  if (v[i] > max) max = v[i];
120  sum += v[i];
121  }
122 }
123 
142 template <typename scalar_t, typename lno_t, typename part_t>
144  const RCP<const Environment> &env,
145  part_t numberOfParts,
146  const ArrayView<const part_t> &parts,
147  const ArrayView<StridedData<lno_t, scalar_t> > &vwgts,
148  multiCriteriaNorm mcNorm,
149  scalar_t *weights)
150 {
151  env->localInputAssertion(__FILE__, __LINE__, "parts or weights",
152  numberOfParts > 0 && vwgts.size() > 0, BASIC_ASSERTION);
153 
154  int numObjects = parts.size();
155  int vwgtDim = vwgts.size();
156 
157  memset(weights, 0, sizeof(scalar_t) * numberOfParts);
158 
159  if (numObjects == 0)
160  return;
161 
162  if (vwgtDim == 0) {
163  for (lno_t i=0; i < numObjects; i++){
164  weights[parts[i]]++;
165  }
166  }
167  else if (vwgtDim == 1){
168  for (lno_t i=0; i < numObjects; i++){
169  weights[parts[i]] += vwgts[0][i];
170  }
171  }
172  else{
173  switch (mcNorm){
175  for (int wdim=0; wdim < vwgtDim; wdim++){
176  for (lno_t i=0; i < numObjects; i++){
177  weights[parts[i]] += vwgts[wdim][i];
178  }
179  } // next weight index
180  break;
181 
183  for (lno_t i=0; i < numObjects; i++){
184  scalar_t ssum = 0;
185  for (int wdim=0; wdim < vwgtDim; wdim++)
186  ssum += (vwgts[wdim][i] * vwgts[wdim][i]);
187  weights[parts[i]] += sqrt(ssum);
188  }
189  break;
190 
192  for (lno_t i=0; i < numObjects; i++){
193  scalar_t max = 0;
194  for (int wdim=0; wdim < vwgtDim; wdim++)
195  if (vwgts[wdim][i] > max)
196  max = vwgts[wdim][i];
197  weights[parts[i]] += max;
198  }
199  break;
200 
201  default:
202  env->localBugAssertion(__FILE__, __LINE__, "invalid norm", false,
204  break;
205  }
206  }
207 }
208 
239 template <typename scalar_t, typename part_t>
241  part_t numExistingParts, // Max Part ID + 1
242  part_t targetNumParts, // comm.size() or requested global # parts from soln
243  const scalar_t *psizes,
244  scalar_t sumVals,
245  const scalar_t *vals,
246  scalar_t &min,
247  scalar_t &max,
248  scalar_t &avg)
249 {
250  min = sumVals;
251  max = avg = 0;
252 
253  if (sumVals <= 0 || targetNumParts < 1 || numExistingParts < 1)
254  return;
255 
256  if (targetNumParts==1) {
257  min = max = avg = 0; // 0 imbalance
258  return;
259  }
260 
261  if (!psizes){
262  scalar_t target = sumVals / targetNumParts;
263  for (part_t p=0; p < numExistingParts; p++){
264  scalar_t diff = vals[p] - target;
265  scalar_t adiff = (diff >= 0 ? diff : -diff);
266  scalar_t tmp = diff / target;
267  scalar_t atmp = adiff / target;
268  avg += atmp;
269  if (tmp > max) max = tmp;
270  if (tmp < min) min = tmp;
271  }
272  part_t emptyParts = targetNumParts - numExistingParts;
273  if (emptyParts > 0){
274  if (max < 1.0)
275  max = 1.0; // target divided by target
276  avg += emptyParts;
277  }
278  }
279  else{
280  for (part_t p=0; p < targetNumParts; p++){
281  if (psizes[p] > 0){
282  if (p < numExistingParts){
283  scalar_t target = sumVals * psizes[p];
284  scalar_t diff = vals[p] - target;
285  scalar_t adiff = (diff >= 0 ? diff : -diff);
286  scalar_t tmp = diff / target;
287  scalar_t atmp = adiff / target;
288  avg += atmp;
289  if (tmp > max) max = tmp;
290  if (tmp < min) min = tmp;
291  }
292  else{
293  if (max < 1.0)
294  max = 1.0; // target divided by target
295  avg += 1.0;
296  }
297  }
298  }
299  }
300 
301  avg /= targetNumParts;
302 }
303 
337 template <typename scalar_t, typename part_t>
339  part_t numExistingParts,
340  part_t targetNumParts,
341  int numSizes,
342  ArrayView<ArrayRCP<scalar_t> > psizes,
343  scalar_t sumVals,
344  const scalar_t *vals,
345  scalar_t &min,
346  scalar_t &max,
347  scalar_t &avg)
348 {
349  min = sumVals;
350  max = avg = 0;
351 
352  if (sumVals <= 0 || targetNumParts < 1 || numExistingParts < 1)
353  return;
354 
355  if (targetNumParts==1) {
356  min = max = avg = 0; // 0 imbalance
357  return;
358  }
359 
360  bool allUniformParts = true;
361  for (int i=0; i < numSizes; i++){
362  if (psizes[i].size() != 0){
363  allUniformParts = false;
364  break;
365  }
366  }
367 
368  if (allUniformParts){
369  computeImbalances<scalar_t, part_t>(numExistingParts, targetNumParts, NULL,
370  sumVals, vals, min, max, avg);
371  return;
372  }
373 
374  double uniformSize = 1.0 / targetNumParts;
375  std::vector<double> sizeVec(numSizes);
376  for (int i=0; i < numSizes; i++){
377  sizeVec[i] = uniformSize;
378  }
379 
380  for (part_t p=0; p < numExistingParts; p++){
381 
382  // If we have objects in parts that should have 0 objects,
383  // we don't compute an imbalance. It means that other
384  // parts have too few objects, and the imbalance will be
385  // picked up there.
386 
387  if (p >= targetNumParts)
388  break;
389 
390  // Vector of target amounts: T
391 
392  double targetNorm = 0;
393  for (int i=0; i < numSizes; i++) {
394  if (psizes[i].size() > 0)
395  sizeVec[i] = psizes[i][p];
396  sizeVec[i] *= sumVals;
397  targetNorm += (sizeVec[i] * sizeVec[i]);
398  }
399  targetNorm = sqrt(targetNorm);
400 
401  // If part is supposed to be empty, we don't compute an
402  // imbalance. Same argument as above.
403 
404  if (targetNorm > 0){
405 
406  // Vector of actual amounts: A
407 
408  std::vector<double> actual(numSizes);
409  double actualNorm = 0.;
410  for (int i=0; i < numSizes; i++) {
411  actual[i] = vals[p] * -1.0;
412  actual[i] += sizeVec[i];
413  actualNorm += (actual[i] * actual[i]);
414  }
415  actualNorm = sqrt(actualNorm);
416 
417  // |A - T| / |T|
418 
419  scalar_t imbalance = actualNorm / targetNorm;
420 
421  if (imbalance < min)
422  min = imbalance;
423  else if (imbalance > max)
424  max = imbalance;
425  avg += imbalance;
426  }
427  }
428 
429  part_t numEmptyParts = 0;
430 
431  for (part_t p=numExistingParts; p < targetNumParts; p++){
432  bool nonEmptyPart = false;
433  for (int i=0; !nonEmptyPart && i < numSizes; i++)
434  if (psizes[i].size() > 0 && psizes[i][p] > 0.0)
435  nonEmptyPart = true;
436 
437  if (nonEmptyPart){
438  // The partition has no objects for this part, which
439  // is supposed to be non-empty.
440  numEmptyParts++;
441  }
442  }
443 
444  if (numEmptyParts > 0){
445  avg += numEmptyParts;
446  if (max < 1.0)
447  max = 1.0; // target divided by target
448  }
449 
450  avg /= targetNumParts;
451 }
452 
455 template <typename scalar_t>
456  scalar_t normedWeight(ArrayView <scalar_t> weights,
457  multiCriteriaNorm norm)
458 {
459  size_t dim = weights.size();
460  if (dim == 0)
461  return 0.0;
462  else if (dim == 1)
463  return weights[0];
464 
465  scalar_t nweight = 0;
466 
467  switch (norm){
470  for (size_t i=0; i <dim; i++)
471  nweight += weights[i];
472 
473  break;
474 
476  for (size_t i=0; i <dim; i++)
477  nweight += (weights[i] * weights[i]);
478 
479  nweight = sqrt(nweight);
480 
481  break;
482 
485  nweight = weights[0];
486  for (size_t i=1; i <dim; i++)
487  if (weights[i] > nweight)
488  nweight = weights[i];
489 
490  break;
491 
492  default:
493  std::ostringstream emsg;
494  emsg << __FILE__ << ":" << __LINE__ << std::endl;
495  emsg << "bug: " << "invalid norm" << std::endl;
496  throw std::logic_error(emsg.str());
497  }
498 
499  return nweight;
500 }
501 
504 template<typename lno_t, typename scalar_t>
506  lno_t idx, multiCriteriaNorm norm)
507 {
508  size_t dim = weights.size();
509  if (dim < 1)
510  return 0;
511 
512  Array<scalar_t> vec(dim, 1.0);
513  for (size_t i=0; i < dim; i++)
514  if (weights[i].size() > 0)
515  vec[dim] = weights[i][idx];
516 
517  return normedWeight(vec, norm);
518 }
519 
520 } //namespace Zoltan2
521 
522 #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.
RCP< metric_t > addNewMetric(const RCP< const Environment > &env, ArrayRCP< RCP< BaseClassMetrics< scalar_t > > > &metrics)
static ArrayRCP< ArrayRCP< zscalar_t > > weights
SparseMatrixAdapter_t::part_t part_t
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
The StridedData class manages lists of weights or coordinates.
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.