Zoltan2
Zoltan2_Metric.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //   Zoltan2: A package of combinatorial algorithms for scientific computing
00006 //                  Copyright 2012 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Karen Devine      (kddevin@sandia.gov)
00039 //                    Erik Boman        (egboman@sandia.gov)
00040 //                    Siva Rajamanickam (srajama@sandia.gov)
00041 //
00042 // ***********************************************************************
00043 //
00044 // @HEADER
00045 
00051 #ifndef ZOLTAN2_METRIC_HPP
00052 #define ZOLTAN2_METRIC_HPP
00053 
00054 #include <Zoltan2_StridedData.hpp>
00055 #include <Zoltan2_PartitioningSolution.hpp>
00056 
00057 #include <Epetra_SerialDenseVector.h>
00058 
00059 #include <cmath>
00060 #include <iomanip>
00061 #include <iostream>
00062 
00063 namespace Zoltan2{
00064 
00066 // Classes
00068 
00072 template <typename scalar_t>
00073   class MetricValues{
00074 
00075 private:
00076   void resetValues(){
00077     scalar_t *tmp = new scalar_t [evalNumMetrics];
00078     memset(tmp, 0, sizeof(scalar_t) * evalNumMetrics);
00079     values_ = arcp(tmp, 0, evalNumMetrics, true);
00080   }
00081   ArrayRCP<scalar_t> values_;
00082   std::string metricName_;
00083   multiCriteriaNorm mcnorm_;   // store "actualNorm + 1"
00084 
00085 public:
00086 
00094 enum metricOffset{
00095   evalLocalSum,    
00096   evalGlobalSum,   
00097   evalGlobalMin,   
00098   evalGlobalMax,   
00099   evalGlobalAvg,   
00100   evalMinImbalance, 
00101   evalAvgImbalance, 
00102   evalMaxImbalance, 
00103   evalNumMetrics    
00104 };
00105 
00113 static void printHeader(std::ostream &os);
00114 
00116 void printLine(std::ostream &os) const;
00117 
00119 MetricValues(std::string mname) : 
00120   values_(), metricName_(mname), mcnorm_(multiCriteriaNorm(0)) { 
00121     resetValues();}
00122 
00124 MetricValues() : 
00125   values_(), metricName_("unset"), mcnorm_(multiCriteriaNorm(0)) { 
00126     resetValues();}
00127 
00129 void setName(std::string name) { metricName_ = name;}
00130 
00132 void setNorm(multiCriteriaNorm normVal) { 
00133   mcnorm_ = multiCriteriaNorm(normVal+1);}
00134 
00136 void setLocalSum(scalar_t x) { values_[evalLocalSum] = x;}
00137 
00139 void setGlobalSum(scalar_t x) { values_[evalGlobalSum] = x;}
00140 
00142 void setGlobalMin(scalar_t x) { values_[evalGlobalMin] = x;}
00143 
00145 void setGlobalMax(scalar_t x) { values_[evalGlobalMax] = x;}
00146 
00148 void setGlobalAvg(scalar_t x) { values_[evalGlobalAvg] = x;}
00149 
00151 void setMinImbalance(scalar_t x) { values_[evalMinImbalance] = x;}
00152 
00156 void setMaxImbalance(scalar_t x) { values_[evalMaxImbalance] = x;}
00157 
00159 void setAvgImbalance(scalar_t x) { values_[evalAvgImbalance] = x;}
00160 
00162 const std::string &getName() const { return metricName_; }
00163 
00165 multiCriteriaNorm getNorm() { return multiCriteriaNorm(mcnorm_-1);}
00166 
00168 scalar_t getLocalSum() const { return values_[evalLocalSum];}
00169 
00171 scalar_t getGlobalSum() const { return values_[evalGlobalSum];}
00172 
00174 scalar_t getGlobalMin() const { return values_[evalGlobalMin];}
00175 
00177 scalar_t getGlobalMax() const { return values_[evalGlobalMax];}
00178 
00180 scalar_t getGlobalAvg() const { return values_[evalGlobalAvg];}
00181 
00183 scalar_t getMinImbalance() const { return values_[evalMinImbalance];}
00184 
00188 scalar_t getMaxImbalance() const { return values_[evalMaxImbalance];}
00189 
00191 scalar_t getAvgImbalance() const { return values_[evalAvgImbalance];}
00192 
00193 };  // end class
00194 
00195 
00196 template <typename scalar_t>
00197   void MetricValues<scalar_t>::printLine(std::ostream &os) const
00198 {
00199   std::string label(metricName_);
00200   if (mcnorm_ > 0){
00201     multiCriteriaNorm realNorm = multiCriteriaNorm(mcnorm_ - 1);
00202     std::ostringstream oss;
00203     switch (realNorm){
00204       case normMinimizeTotalWeight:   // 1-norm = Manhattan norm 
00205         oss << metricName_ << " (1)";
00206         break;
00207       case normBalanceTotalMaximum:   // 2-norm = sqrt of sum of squares
00208         oss << metricName_ << " (2)";
00209         break;
00210       case normMinimizeMaximumWeight: // inf-norm = maximum norm 
00211         oss << metricName_ << " (inf)";
00212         break;
00213       default:
00214         oss << metricName_ << " (?)";
00215         break;
00216     }
00217 
00218     label = oss.str();
00219   }
00220 
00221   os << std::setw(20) << label;
00222   os << std::setw(12) << std::setprecision(4) << values_[evalGlobalMin];
00223   os << std::setw(12) << std::setprecision(4) << values_[evalGlobalMax];
00224   os << std::setw(12) << std::setprecision(4) << values_[evalGlobalAvg];
00225   os << std::setw(2) << " ";
00226   os << std::setw(6) << std::setprecision(4) << values_[evalMinImbalance];
00227   os << std::setw(6) << std::setprecision(4) << values_[evalMaxImbalance];
00228   os << std::setw(6) << std::setprecision(4) << values_[evalAvgImbalance];
00229   os << std::endl;
00230 }
00231 
00232 template <typename scalar_t>
00233   void MetricValues<scalar_t>::printHeader(std::ostream &os)
00234 {
00235   os << std::setw(20) << " ";
00236   os << std::setw(36) << "------------SUM PER PART-----------";
00237   os << std::setw(2) << " ";
00238   os << std::setw(18) << "IMBALANCE PER PART";
00239   os << std::endl;
00240 
00241   os << std::setw(20) << " ";
00242   os << std::setw(12) << "min" << std::setw(12) << "max" << std::setw(12) << "avg";
00243   os << std::setw(2) << " ";
00244   os << std::setw(6) << "min" << std::setw(6) << "max" << std::setw(6) << "avg";
00245   os << std::endl;
00246 }
00247 
00249 // Namespace methods to compute metric values
00251 
00261 template <typename scalar_t>
00262  void getStridedStats(const ArrayView<scalar_t> &v, int stride, 
00263    int offset, scalar_t &min, scalar_t &max, scalar_t &sum)
00264 {
00265   if (v.size() < 1) return;
00266   min = max = sum = v[offset];
00267 
00268   for (int i=offset+stride; i < v.size(); i += stride){
00269     if (v[i] < min) min = v[i];
00270     else if (v[i] > max) max = v[i];
00271     sum += v[i];
00272   }
00273 }
00274 
00293 template <typename scalar_t, typename pnum_t, typename lno_t>
00294   void normedPartWeights(
00295     const RCP<const Environment> &env,
00296     int numberOfParts,
00297     const ArrayView<const pnum_t> &parts,
00298     const ArrayView<StridedData<lno_t, scalar_t> > &vwgts,
00299     multiCriteriaNorm mcNorm,
00300     scalar_t *weights)
00301 {
00302   env->localInputAssertion(__FILE__, __LINE__, "parts or weights", 
00303     numberOfParts > 0 && vwgts.size() > 0, BASIC_ASSERTION);
00304 
00305   int numObjects = parts.size();
00306   int vwgtDim = vwgts.size();
00307 
00308   memset(weights, 0, sizeof(scalar_t) * numberOfParts);
00309 
00310   if (numObjects == 0)
00311     return;
00312 
00313   if (vwgtDim == 0) {
00314     for (lno_t i=0; i < numObjects; i++){
00315       weights[parts[i]]++;
00316     }
00317   }
00318   else if (vwgtDim == 1){
00319     for (lno_t i=0; i < numObjects; i++){
00320       weights[parts[i]] += vwgts[0][i];
00321     }
00322   }
00323   else{
00324     switch (mcNorm){
00325       case normMinimizeTotalWeight:   
00326         for (int wdim=0; wdim < vwgtDim; wdim++){
00327           for (lno_t i=0; i < numObjects; i++){
00328             weights[parts[i]] += vwgts[wdim][i];
00329           }
00330         }  // next weight index
00331         break;
00332        
00333       case normBalanceTotalMaximum:   
00334         for (lno_t i=0; i < numObjects; i++){
00335           scalar_t ssum = 0;
00336           for (int wdim=0; wdim < vwgtDim; wdim++)
00337             ssum += (vwgts[wdim][i] * vwgts[wdim][i]);
00338           weights[parts[i]] += sqrt(ssum);
00339         }
00340         break;
00341 
00342       case normMinimizeMaximumWeight: 
00343         for (lno_t i=0; i < numObjects; i++){
00344           scalar_t max = 0;
00345           for (int wdim=0; wdim < vwgtDim; wdim++)
00346             if (vwgts[wdim][i] > max)
00347               max = vwgts[wdim][i];
00348           weights[parts[i]] += max;
00349         }
00350         break;
00351 
00352       default:
00353         env->localBugAssertion(__FILE__, __LINE__, "invalid norm", false,
00354           BASIC_ASSERTION);
00355         break;
00356     } 
00357   } 
00358 }
00359 
00401 template <typename scalar_t, typename pnum_t, typename lno_t>
00402   void globalSumsByPart( 
00403     const RCP<const Environment> &env,
00404     const RCP<const Comm<int> > &comm, 
00405     const ArrayView<const pnum_t> &part, 
00406     int vwgtDim,
00407     const ArrayView<StridedData<lno_t, scalar_t> > &vwgts,
00408     multiCriteriaNorm mcNorm,
00409     partId_t &numParts, 
00410     partId_t &numNonemptyParts,
00411     ArrayRCP<MetricValues<scalar_t> > &metrics,
00412     ArrayRCP<scalar_t> &globalSums)
00413 {
00414   env->debug(DETAILED_STATUS, "Entering globalSumsByPart");
00416   // Initialize return values
00417 
00418   numParts = numNonemptyParts = 0;
00419 
00420   int numMetrics = 1;                       // "object count"
00421   if (vwgtDim) numMetrics++;                // "normed weight" or "weight 1"
00422   if (vwgtDim > 1) numMetrics += vwgtDim;   // "weight n"
00423 
00424   typedef MetricValues<scalar_t> mv_t;
00425   mv_t *newMetrics = new mv_t [numMetrics];
00426   env->localMemoryAssertion(__FILE__, __LINE__, numMetrics, newMetrics); 
00427   ArrayRCP<mv_t> metricArray(newMetrics, 0, numMetrics, true);
00428 
00429   metrics = metricArray;
00430 
00432   // Figure out the global number of parts in use.
00433   // Verify number of vertex weights is the same everywhere.
00434 
00435   lno_t localNumObj = part.size();
00436   partId_t localNum[2], globalNum[2];
00437   localNum[0] = static_cast<partId_t>(vwgtDim);  
00438   localNum[1] = 0;
00439 
00440   for (lno_t i=0; i < localNumObj; i++)
00441     if (part[i] > localNum[1]) localNum[1] = part[i];
00442 
00443   try{
00444     reduceAll<int, partId_t>(*comm, Teuchos::REDUCE_MAX, 2, 
00445       localNum, globalNum);
00446   }
00447   Z2_THROW_OUTSIDE_ERROR(*env)
00448 
00449   env->globalBugAssertion(__FILE__,__LINE__,
00450     "inconsistent number of vertex weights",
00451     globalNum[0] == localNum[0], DEBUG_MODE_ASSERTION, comm);
00452 
00453   partId_t nparts = globalNum[1] + 1;
00454 
00455   int globalSumSize = nparts * numMetrics;
00456   scalar_t * sumBuf = new scalar_t [globalSumSize];
00457   env->localMemoryAssertion(__FILE__, __LINE__, globalSumSize, sumBuf);
00458   globalSums = arcp(sumBuf, 0, globalSumSize);
00459 
00461   // Calculate the local totals by part.
00462 
00463   scalar_t *localBuf = new scalar_t [globalSumSize];
00464   env->localMemoryAssertion(__FILE__, __LINE__, globalSumSize, localBuf);
00465   memset(localBuf, 0, sizeof(scalar_t) * globalSumSize);
00466 
00467   scalar_t *obj = localBuf;              // # of objects
00468 
00469   for (lno_t i=0; i < localNumObj; i++)
00470     obj[part[i]]++;
00471 
00472   if (numMetrics > 1){
00473 
00474     scalar_t *wgt = localBuf + nparts; // single normed weight
00475     try{
00476       normedPartWeights<scalar_t, pnum_t, lno_t>(env, nparts, 
00477         part, vwgts, mcNorm, wgt);
00478     }
00479     Z2_FORWARD_EXCEPTIONS
00480   
00481 //KDDKDD TODO  This code assumes the solution has the part ordered the
00482 //KDDKDD TODO  same way as the user input.  That assumption is not
00483 //KDDKDD TODO  currently true, although we plan to make it true.
00484 //KDDKDD TODO  As a results, currently the weight metrics may be wrong.
00485 //KDDKDD TODO  See bug 5891.  April 5, 2013
00486     if (vwgtDim > 1){
00487       wgt += nparts;         // individual weights
00488       for (int vdim = 0; vdim < vwgtDim; vdim++){
00489         for (lno_t i=0; i < localNumObj; i++)
00490           wgt[part[i]] += vwgts[vdim][i];
00491         wgt += nparts;
00492       }
00493     }
00494   }
00495 
00496   // Metric: local sums on process
00497 
00498   int next = 0;
00499 
00500   metrics[next].setName("object count");
00501   metrics[next].setLocalSum(localNumObj);
00502   next++;
00503 
00504   if (numMetrics > 1){
00505     scalar_t *wgt = localBuf + nparts; // single normed weight
00506     scalar_t total = 0.0;
00507   
00508     for (int p=0; p < nparts; p++){
00509       total += wgt[p];
00510     }
00511 
00512     if (vwgtDim == 1)
00513       metrics[next].setName("weight 1");
00514     else
00515       metrics[next].setName("normed weight");
00516 
00517     metrics[next].setLocalSum(total);
00518     next++;
00519   
00520     if (vwgtDim > 1){
00521       for (int vdim = 0; vdim < vwgtDim; vdim++){
00522         wgt += nparts;
00523         total = 0.0;
00524         for (int p=0; p < nparts; p++){
00525           total += wgt[p];
00526         }
00527 
00528         std::ostringstream oss;
00529         oss << "weight " << vdim+1;
00530 
00531         metrics[next].setName(oss.str());
00532         metrics[next].setLocalSum(total);
00533         next++;
00534       }
00535     }
00536   }
00537 
00539   // Obtain global totals by part.
00540 
00541   try{
00542     reduceAll<int, scalar_t>(*comm, Teuchos::REDUCE_SUM, globalSumSize,
00543       localBuf, sumBuf);
00544   }
00545   Z2_THROW_OUTSIDE_ERROR(*env);
00546 
00547   delete [] localBuf;
00548 
00550   // Global sum, min, max, and average over all parts
00551 
00552   obj = sumBuf;                     // # of objects
00553   scalar_t min=0, max=0, sum=0;
00554   next = 0;
00555 
00556   ArrayView<scalar_t> objVec(obj, nparts);
00557   getStridedStats<scalar_t>(objVec, 1, 0, min, max, sum);
00558 
00559   metrics[next].setGlobalMin(min);
00560   metrics[next].setGlobalMax(max);
00561   metrics[next].setGlobalSum(sum);
00562   metrics[next].setGlobalAvg(sum / nparts);
00563   next++;
00564 
00565   if (numMetrics > 1){
00566     scalar_t *wgt = sumBuf + nparts;        // single normed weight
00567   
00568     ArrayView<scalar_t> normedWVec(wgt, nparts);
00569     getStridedStats<scalar_t>(normedWVec, 1, 0, min, max, sum);
00570 
00571     if (vwgtDim > 1)
00572       metrics[next].setNorm(multiCriteriaNorm(mcNorm));
00573 
00574     metrics[next].setGlobalMin(min);
00575     metrics[next].setGlobalMax(max);
00576     metrics[next].setGlobalSum(sum);
00577     metrics[next].setGlobalAvg(sum / nparts);
00578     next++;
00579   
00580     if (vwgtDim > 1){
00581       for (int vdim=0; vdim < vwgtDim; vdim++){
00582         wgt += nparts;       // individual weights
00583         ArrayView<scalar_t> fromVec(wgt, nparts);
00584         getStridedStats<scalar_t>(fromVec, 1, 0, min, max, sum);
00585 
00586         metrics[next].setGlobalMin(min);
00587         metrics[next].setGlobalMax(max);
00588         metrics[next].setGlobalSum(sum);
00589         metrics[next].setGlobalAvg(sum / nparts);
00590         next++;
00591       }
00592     }
00593   }
00594 
00596   // How many parts do we actually have.
00597 
00598   numParts = nparts;
00599   obj = sumBuf;               // # of objects
00600 
00601   for (partId_t p=nparts-1; p > 0; p--){
00602     if (obj[p] > 0) break;
00603     numParts--;
00604   }
00605 
00606   numNonemptyParts = numParts; 
00607 
00608   for (partId_t p=0; p < numParts; p++)
00609     if (obj[p] == 0) numNonemptyParts--;
00610 
00611   env->debug(DETAILED_STATUS, "Exiting globalSumsByPart");
00612 }
00613 
00644 template <typename scalar_t>
00645   void computeImbalances(partId_t numParts, partId_t targetNumParts,
00646     const scalar_t *psizes, scalar_t sumVals , const scalar_t *vals, 
00647     scalar_t &min, scalar_t &max, scalar_t &avg)
00648 {
00649   min = sumVals;
00650   max = avg = 0;
00651 
00652   if (sumVals <= 0 || targetNumParts < 1 || numParts < 1)
00653     return;
00654 
00655   if (targetNumParts==1 || numParts==1){
00656     min = max = avg = 0;  // 0 imbalance
00657     return;
00658   }
00659 
00660   if (!psizes){
00661     scalar_t target = sumVals / targetNumParts;
00662     for (partId_t p=0; p < numParts; p++){
00663       scalar_t diff = abs(vals[p] - target);
00664       scalar_t tmp = diff / target;
00665       avg += tmp;
00666       if (tmp > max) max = tmp;
00667       if (tmp < min) min = tmp;
00668     }
00669     partId_t emptyParts = targetNumParts - numParts;  
00670     if (emptyParts > 0){
00671       if (max < 1.0)
00672         max = 1.0;       // target divided by target
00673       avg += emptyParts;
00674     }
00675   }
00676   else{
00677     for (partId_t p=0; p < targetNumParts; p++){
00678       if (psizes[p] > 0){
00679         if (p < numParts){
00680           scalar_t target = sumVals * psizes[p];
00681           scalar_t diff = abs(vals[p] - target);
00682           scalar_t tmp = diff / target;
00683           avg += tmp;
00684           if (tmp > max) max = tmp;
00685           if (tmp < min) min = tmp;
00686         }
00687         else{
00688           if (max < 1.0)
00689             max = 1.0;  // target divided by target
00690           avg += 1.0;
00691         }
00692       }
00693     }
00694   }
00695 
00696   avg /= targetNumParts;
00697 }
00698 
00732 template <typename scalar_t>
00733  void computeImbalances(partId_t numParts, partId_t targetNumParts,
00734    int numSizes, ArrayView<ArrayRCP<scalar_t> > psizes,
00735    scalar_t sumVals , const scalar_t *vals, 
00736    scalar_t &min, scalar_t &max, scalar_t &avg)
00737 {
00738   min = sumVals;
00739   max = avg = 0;
00740 
00741   if (sumVals <= 0 || targetNumParts < 1 || numParts < 1)
00742     return;
00743 
00744   if (targetNumParts==1 && numParts==1){
00745     min = max = avg = 0;  // 0 imbalance
00746     return;
00747   }
00748 
00749   bool allUniformParts = true;
00750   for (int i=0; i < numSizes; i++){
00751     if (psizes[i].size() != 0){
00752       allUniformParts = false;
00753       break;
00754     }
00755   }
00756 
00757   if (allUniformParts){
00758     computeImbalances<scalar_t>(numParts, targetNumParts, NULL,
00759       sumVals, vals, min, max, avg);
00760     return;
00761   }
00762 
00763   double uniformSize = 1.0 / targetNumParts;
00764   ArrayRCP<double> sizeVec(new double [numSizes], 0, numSizes, true);
00765   for (int i=0; i < numSizes; i++){
00766     sizeVec[i] = uniformSize;
00767   }
00768 
00769   for (partId_t p=0; p < numParts; p++){
00770 
00771     // If we have objects in parts that should have 0 objects,
00772     // we don't compute an imbalance.  It means that other
00773     // parts have too few objects, and the imbalance will be
00774     // picked up there.
00775 
00776     if (p >= targetNumParts)
00777       break;
00778 
00779     // Vector of target amounts: T
00780 
00781     for (int i=0; i < numSizes; i++)
00782       if (psizes[i].size() > 0)
00783         sizeVec[i] = psizes[i][p];
00784 
00785     Epetra_SerialDenseVector target(View, sizeVec.getRawPtr(), numSizes);
00786     target.Scale(sumVals);
00787     double targetNorm = target.Norm2();
00788 
00789     // If part is supposed to be empty, we don't compute an
00790     // imbalance.  Same argument as above.
00791 
00792     if (targetNorm > 0){
00793 
00794       // Vector of actual amounts: A
00795 
00796       Epetra_SerialDenseVector actual(numSizes);
00797       for (int i=0; i < numSizes; i++)
00798         actual[i] = vals[p] * -1.0;
00799       
00800       actual += target;
00801 
00802       //  |A - T| / |T|
00803 
00804       scalar_t imbalance = actual.Norm2() / targetNorm;
00805 
00806       if (imbalance < min)
00807         min = imbalance;
00808       else if (imbalance > max)
00809         max = imbalance;
00810       avg += imbalance; 
00811     }
00812   }
00813 
00814   partId_t numEmptyParts = 0;
00815 
00816   for (partId_t p=numParts; p < targetNumParts; p++){
00817     bool nonEmptyPart = false;
00818     for (int i=0; !nonEmptyPart && i < numSizes; i++)
00819       if (psizes[i].size() > 0 && psizes[i][p] > 0.0)
00820         nonEmptyPart = true;
00821 
00822     if (nonEmptyPart){
00823       // The partitioning has no objects for this part, which
00824       // is supposed to be non-empty.
00825       numEmptyParts++;
00826     }
00827   }
00828 
00829   if (numEmptyParts > 0){
00830     avg += numEmptyParts;
00831     if (max < 1.0)
00832       max = 1.0;       // target divided by target
00833   }
00834 
00835   avg /= targetNumParts;
00836 }
00837 
00867 template <typename Adapter>
00868   void objectMetrics(
00869     const RCP<const Environment> &env,
00870     const RCP<const Comm<int> > &comm,
00871     multiCriteriaNorm mcNorm,
00872     const RCP<const Adapter> &ia,
00873     const RCP<const PartitioningSolution<Adapter> > &solution,
00874     partId_t &numParts,
00875     partId_t &numNonemptyParts,
00876     ArrayRCP<MetricValues<typename Adapter::scalar_t> > &metrics)
00877 {
00878   env->debug(DETAILED_STATUS, "Entering objectMetrics");
00879 
00880   typedef typename Adapter::scalar_t scalar_t;
00881   typedef typename Adapter::lno_t lno_t;
00882   typedef StridedData<lno_t, scalar_t> sdata_t;
00883 
00884   // Local number of objects.
00885 
00886   size_t numLocalObjects = ia->getLocalNumIDs();
00887 
00888   // Parts to which objects are assigned.
00889 
00890   const partId_t *parts = solution->getPartList();
00891   env->localInputAssertion(__FILE__, __LINE__, "parts not set", 
00892     parts, BASIC_ASSERTION);
00893   ArrayView<const partId_t> partArray(parts, numLocalObjects);
00894 
00895   // Weights, if any, for each object.
00896 
00897   int nWeights = ia->getNumWeightsPerID();
00898   int numCriteria = (nWeights > 0 ? nWeights : 1);
00899   Array<sdata_t> weights(numCriteria);
00900 
00901   if (nWeights == 0){
00902     // One set of uniform weights is implied.
00903     // StridedData default constructor creates length 0 strided array.
00904     weights[0] = sdata_t();
00905   }
00906   else{
00907     for (int i=0; i < nWeights; i++){
00908       int stride;
00909       const scalar_t *wgt;
00910       ia->getWeightsView(wgt, stride, i); 
00911       ArrayRCP<const scalar_t> wgtArray(wgt, 0, stride*numLocalObjects, false);
00912       weights[i] = sdata_t(wgtArray, stride);
00913     }
00914   }
00915 
00916   // Relative part sizes, if any, assigned to the parts.
00917 
00918   partId_t targetNumParts = solution->getTargetGlobalNumberOfParts();
00919   scalar_t *psizes = NULL;
00920 
00921   ArrayRCP<ArrayRCP<scalar_t> > partSizes(numCriteria);
00922   for (int dim=0; dim < numCriteria; dim++){
00923     if (solution->criteriaHasUniformPartSizes(dim) != true){
00924       psizes = new scalar_t [targetNumParts];
00925       env->localMemoryAssertion(__FILE__, __LINE__, numParts, psizes);
00926       for (partId_t i=0; i < targetNumParts; i++){
00927         psizes[i] = solution->getCriteriaPartSize(dim, i);
00928       }
00929       partSizes[dim] = arcp(psizes, 0, targetNumParts, true);
00930     }
00931   }
00932 
00934   // Get number of parts, and the number that are non-empty.
00935   // Get sums per part of objects, individual weights, and normed weight sums.
00936 
00937   ArrayRCP<scalar_t> globalSums;
00938 
00939   try{
00940     globalSumsByPart<scalar_t, partId_t, lno_t>(env, comm, 
00941       partArray, nWeights, weights.view(0, numCriteria), mcNorm,
00942       numParts, numNonemptyParts, metrics, globalSums);
00943   }
00944   Z2_FORWARD_EXCEPTIONS
00945 
00947   // Compute imbalances for the object count.  
00948   // (Use first index of part sizes.) 
00949 
00950   scalar_t *objCount  = globalSums.getRawPtr();
00951   scalar_t min, max, avg;
00952   psizes=NULL;
00953 
00954   if (partSizes[0].size() > 0)
00955     psizes = partSizes[0].getRawPtr();
00956 
00957   computeImbalances<scalar_t>(numParts, targetNumParts, psizes,
00958       metrics[0].getGlobalSum(), objCount, 
00959       min, max, avg);
00960 
00961   metrics[0].setMinImbalance(1.0 + min);
00962   metrics[0].setMaxImbalance(1.0 + max);
00963   metrics[0].setAvgImbalance(1.0 + avg);
00964 
00966   // Compute imbalances for the normed weight sum.
00967 
00968   scalar_t *wgts = globalSums.getRawPtr() + numParts;
00969 
00970   if (metrics.size() > 1){
00971   
00972     computeImbalances<scalar_t>(numParts, targetNumParts, 
00973       numCriteria, partSizes.view(0, numCriteria),
00974       metrics[1].getGlobalSum(), wgts,
00975       min, max, avg);
00976   
00977     metrics[1].setMinImbalance(1.0 + min);
00978     metrics[1].setMaxImbalance(1.0 + max);
00979     metrics[1].setAvgImbalance(1.0 + avg);
00980 
00981     if (metrics.size() > 2){
00982 
00984     // Compute imbalances for each individual weight.
00985 
00986       int next = 2;
00987   
00988       for (int vdim=0; vdim < numCriteria; vdim++){
00989         wgts += numParts;
00990         psizes = NULL;
00991   
00992         if (partSizes[vdim].size() > 0)
00993            psizes = partSizes[vdim].getRawPtr();
00994          
00995         computeImbalances<scalar_t>(numParts, targetNumParts, psizes,
00996           metrics[next].getGlobalSum(), wgts, min, max, avg);
00997   
00998         metrics[next].setMinImbalance(1.0 + min);
00999         metrics[next].setMaxImbalance(1.0 + max);
01000         metrics[next].setAvgImbalance(1.0 + avg);
01001         next++;
01002       }
01003     }
01004 
01005   }
01006   env->debug(DETAILED_STATUS, "Exiting objectMetrics");
01007 }
01008 
01012 template <typename scalar_t>
01013   void printMetrics( std::ostream &os,
01014     partId_t targetNumParts, partId_t numParts, partId_t numNonemptyParts, 
01015     const ArrayView<MetricValues<scalar_t> > &infoList)
01016 {
01017   os << "NUMBER OF PARTS IS " << numParts;
01018   if (numNonemptyParts < numParts)
01019     os << " (" << numNonemptyParts << " of which are non-empty)";
01020   os << std::endl;
01021   if (targetNumParts != numParts)
01022     os << "TARGET NUMBER OF PARTS IS " << targetNumParts << std::endl;
01023 
01024   std::string unset("unset");
01025 
01026   MetricValues<scalar_t>::printHeader(os);
01027 
01028   for (int i=0; i < infoList.size(); i++)
01029     if (infoList[i].getName() != unset)
01030       infoList[i].printLine(os);
01031 
01032   os << std::endl;
01033 }
01034 
01038 template <typename scalar_t>
01039   void printMetrics( std::ostream &os,
01040     partId_t targetNumParts, partId_t numParts, partId_t numNonemptyParts, 
01041     const MetricValues<scalar_t> &info)
01042 {
01043   ArrayView<MetricValues<scalar_t> > infoList(&info, 1);
01044   printMetrics( os, targetNumParts, numParts, numNonemptyParts, infoList);
01045 }
01046 
01047 
01050 template <typename scalar_t>
01051   scalar_t normedWeight(ArrayView <scalar_t> weights, 
01052     multiCriteriaNorm norm)
01053 {
01054   size_t dim = weights.size();
01055   if (dim == 0)
01056     return 0.0;
01057   else if (dim == 1)
01058     return weights[0];
01059  
01060   scalar_t nweight = 0;
01061 
01062   switch (norm){
01063     case normMinimizeTotalWeight:   
01065       for (size_t i=0; i <dim; i++)
01066         nweight += weights[i];
01067 
01068       break;
01069      
01070     case normBalanceTotalMaximum:   
01071       for (size_t i=0; i <dim; i++)
01072         nweight += (weights[i] * weights[i]);
01073 
01074       nweight = sqrt(nweight);
01075 
01076       break;
01077 
01078     case normMinimizeMaximumWeight: 
01080       nweight = weights[0];
01081       for (size_t i=1; i <dim; i++)
01082         if (weights[i] > nweight)
01083           nweight = weights[i];
01084 
01085       break;
01086 
01087     default:
01088       Environment env;  // default environment
01089       env.localBugAssertion(__FILE__, __LINE__, "invalid norm", false,
01090         BASIC_ASSERTION);
01091       break;
01092   } 
01093 
01094   return nweight;
01095 }
01096 
01099 template<typename lno_t, typename scalar_t>
01100   scalar_t normedWeight(ArrayView<StridedData<lno_t, scalar_t> > weights, 
01101      lno_t idx, multiCriteriaNorm norm)
01102 {
01103   size_t dim = weights.size();
01104   if (dim < 1)
01105     return 0;
01106 
01107   Array<scalar_t> vec(dim, 1.0);
01108   for (size_t i=0; i < dim; i++)
01109     if (weights[i].size() > 0)
01110        vec[dim] = weights[i][idx];
01111 
01112   return normedWeight(vec, norm);
01113 }
01114 
01115 } //namespace Zoltan2
01116 #endif