Zoltan2 Version of the Day
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   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(ostream &os);
00114 
00116 void printLine(ostream &os) const;
00117 
00119 MetricValues(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(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 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(ostream &os) const
00198 {
00199   string label(metricName_);
00200   if (mcnorm_ > 0){
00201     multiCriteriaNorm realNorm = multiCriteriaNorm(mcnorm_ - 1);
00202     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 << setw(20) << label;
00222   os << setw(12) << setprecision(4) << values_[evalGlobalMin];
00223   os << setw(12) << setprecision(4) << values_[evalGlobalMax];
00224   os << setw(12) << setprecision(4) << values_[evalGlobalAvg];
00225   os << setw(2) << " ";
00226   os << setw(6) << setprecision(4) << values_[evalMinImbalance];
00227   os << setw(6) << setprecision(4) << values_[evalMaxImbalance];
00228   os << setw(6) << setprecision(4) << values_[evalAvgImbalance];
00229   os << endl;
00230 }
00231 
00232 template <typename scalar_t>
00233   void MetricValues<scalar_t>::printHeader(ostream &os)
00234 {
00235   os << setw(20) << " ";
00236   os << setw(36) << "------------SUM PER PART-----------";
00237   os << setw(2) << " ";
00238   os << setw(18) << "IMBALANCE PER PART";
00239   os << endl;
00240 
00241   os << setw(20) << " ";
00242   os << setw(12) << "min" << setw(12) << "max" << setw(12) << "avg";
00243   os << setw(2) << " ";
00244   os << setw(6) << "min" << setw(6) << "max" << setw(6) << "avg";
00245   os << 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   bool haveUniform = false;
00314   bool haveNonUniform = false;
00315 
00316   for (int wdim=0; wdim < vwgtDim; wdim++){
00317     if (vwgts[wdim].size() > 0)
00318       haveNonUniform = true;
00319     if (vwgts[wdim].size() == 0)
00320       haveUniform = true;
00321   }
00322 
00323   if (!haveNonUniform){
00324     for (lno_t i=0; i < numObjects; i++){
00325       weights[parts[i]]++;
00326     }
00327   }
00328   else if (vwgtDim == 1){
00329     for (lno_t i=0; i < numObjects; i++){
00330       weights[parts[i]] += vwgts[0][i];
00331     }
00332   }
00333   else{
00334     switch (mcNorm){
00335       case normMinimizeTotalWeight:   
00337         for (int wdim=0; wdim < vwgtDim; wdim++){
00338           if (vwgts[wdim].size() == 0){
00339             for (lno_t i=0; i < numObjects; i++){
00340               weights[parts[i]]++;
00341             }
00342           }
00343           else{
00344             for (lno_t i=0; i < numObjects; i++){
00345               weights[parts[i]] += vwgts[wdim][i];
00346             }
00347           }
00348         }  // next weight dimension
00349         break;
00350        
00351       case normBalanceTotalMaximum:   
00352         if (!haveUniform){
00353           for (lno_t i=0; i < numObjects; i++){
00354             scalar_t ssum = 0;
00355             for (int wdim=0; wdim < vwgtDim; wdim++)
00356               ssum += (vwgts[wdim][i] * vwgts[wdim][i]);
00357             weights[parts[i]] += sqrt(ssum);
00358           }
00359         }
00360         else{
00361           scalar_t uniformFactor = 0;
00362           for (int wdim=0; wdim < vwgtDim; wdim++)
00363             if (vwgts[wdim].size() == 0)
00364               uniformFactor++;
00365               
00366           for (lno_t i=0; i < numObjects; i++){
00367             scalar_t ssum = 0;
00368             for (int wdim=0; wdim < vwgtDim; wdim++){
00369               if (vwgts[wdim].size() > 0)
00370                 ssum += (vwgts[wdim][i] * vwgts[wdim][i]);
00371             }
00372             ssum += uniformFactor;
00373             weights[parts[i]] += sqrt(ssum);
00374           }
00375         }
00376         break;
00377 
00378       case normMinimizeMaximumWeight: 
00380         if (!haveUniform){
00381 
00382           for (lno_t i=0; i < numObjects; i++){
00383             scalar_t max = 0;
00384             for (int wdim=0; wdim < vwgtDim; wdim++)
00385               if (vwgts[wdim][i] > max)
00386                 max = vwgts[wdim][i];
00387             weights[parts[i]] += max;
00388           }
00389         }
00390         else{
00391 
00392           for (lno_t i=0; i < numObjects; i++){
00393             scalar_t max = 1.0;
00394             for (int wdim=0; wdim < vwgtDim; wdim++)
00395               if (vwgts[wdim].size() > 0 && vwgts[wdim][i] > max)
00396                 max = vwgts[wdim][i];
00397             weights[parts[i]] += max;
00398           }
00399         }
00400         break;
00401 
00402       default:
00403         env->localBugAssertion(__FILE__, __LINE__, "invalid norm", false,
00404           BASIC_ASSERTION);
00405         break;
00406     } 
00407   } 
00408 }
00409 
00453 template <typename scalar_t, typename pnum_t, typename lno_t>
00454   void globalSumsByPart( 
00455     const RCP<const Environment> &env,
00456     const RCP<const Comm<int> > &comm, 
00457     const ArrayView<const pnum_t> &part, 
00458     const ArrayView<StridedData<lno_t, scalar_t> > &vwgts,
00459     multiCriteriaNorm mcNorm,
00460     partId_t &numParts, 
00461     partId_t &numNonemptyParts,
00462     ArrayRCP<MetricValues<scalar_t> > &metrics,
00463     ArrayRCP<scalar_t> &globalSums)
00464 {
00465   env->debug(DETAILED_STATUS, "Entering globalSumsByPart");
00467   // Initialize return values
00468 
00469   numParts = numNonemptyParts = 0;
00470 
00471   int vwgtDim = vwgts.size();
00472 
00473   int numMetrics = 1;        // "object count"
00474   if (vwgts[0].size() > 0)
00475     numMetrics++;            // "normed weight" or "weight 1"
00476   if (vwgtDim > 1)
00477     numMetrics += vwgtDim;   // "weight n"
00478 
00479   typedef MetricValues<scalar_t> mv_t;
00480   mv_t *newMetrics = new mv_t [numMetrics];
00481   env->localMemoryAssertion(__FILE__, __LINE__, numMetrics, newMetrics); 
00482   ArrayRCP<mv_t> metricArray(newMetrics, 0, numMetrics, true);
00483 
00484   metrics = metricArray;
00485 
00487   // Figure out the global number of parts in use.
00488   // Verify vertex weight dim is the same everywhere.
00489 
00490   lno_t localNumObj = part.size();
00491   partId_t localNum[2], globalNum[2];
00492   localNum[0] = static_cast<partId_t>(vwgtDim);  
00493   localNum[1] = 0;
00494 
00495   for (lno_t i=0; i < localNumObj; i++)
00496     if (part[i] > localNum[1]) localNum[1] = part[i];
00497 
00498   try{
00499     reduceAll<int, partId_t>(*comm, Teuchos::REDUCE_MAX, 2, 
00500       localNum, globalNum);
00501   }
00502   Z2_THROW_OUTSIDE_ERROR(*env)
00503 
00504   env->globalBugAssertion(__FILE__, __LINE__, "inconsistent vertex dimension",
00505     globalNum[0] > 0 && globalNum[0] == localNum[0], 
00506     DEBUG_MODE_ASSERTION, comm);
00507 
00508   partId_t nparts = globalNum[1] + 1;
00509 
00510   int globalSumSize = nparts * numMetrics;
00511   scalar_t * sumBuf = new scalar_t [globalSumSize];
00512   env->localMemoryAssertion(__FILE__, __LINE__, globalSumSize, sumBuf);
00513   globalSums = arcp(sumBuf, 0, globalSumSize);
00514 
00516   // Calculate the local totals by part.
00517 
00518   scalar_t *localBuf = new scalar_t [globalSumSize];
00519   env->localMemoryAssertion(__FILE__, __LINE__, globalSumSize, localBuf);
00520   memset(localBuf, 0, sizeof(scalar_t) * globalSumSize);
00521 
00522   scalar_t *obj = localBuf;              // # of objects
00523 
00524   for (lno_t i=0; i < localNumObj; i++)
00525     obj[part[i]]++;
00526 
00527   if (numMetrics > 1){
00528 
00529     scalar_t *wgt = localBuf + nparts; // single normed weight
00530     try{
00531       normedPartWeights<scalar_t, pnum_t, lno_t>(env, nparts, 
00532         part, vwgts, mcNorm, wgt);
00533     }
00534     Z2_FORWARD_EXCEPTIONS
00535   
00536     if (vwgtDim > 1){
00537       wgt += nparts;         // individual weights
00538       for (int vdim = 0; vdim < vwgtDim; vdim++){
00539         if (vwgts[vdim].size()){
00540          for (lno_t i=0; i < localNumObj; i++)
00541            wgt[part[i]] += vwgts[vdim][i];
00542         }
00543         else{  // uniform weights
00544           for (int p=0; p < nparts; p++)
00545             wgt[p] = obj[p];
00546         }
00547         wgt += nparts;
00548       }
00549     }
00550   }
00551 
00552   // Metric: local sums on process
00553 
00554   int next = 0;
00555 
00556   metrics[next].setName("object count");
00557   metrics[next].setLocalSum(localNumObj);
00558   next++;
00559 
00560   if (numMetrics > 1){
00561     scalar_t *wgt = localBuf + nparts; // single normed weight
00562     scalar_t total = 0.0;
00563   
00564     for (int p=0; p < nparts; p++){
00565       total += wgt[p];
00566     }
00567 
00568     if (vwgtDim == 1)
00569       metrics[next].setName("weight 1");
00570     else
00571       metrics[next].setName("normed weight");
00572 
00573     metrics[next].setLocalSum(total);
00574     next++;
00575   
00576     if (vwgtDim > 1){
00577       for (int vdim = 0; vdim < vwgtDim; vdim++){
00578         wgt += nparts;
00579         total = 0.0;
00580         for (int p=0; p < nparts; p++){
00581           total += wgt[p];
00582         }
00583 
00584         ostringstream oss;
00585         oss << "weight " << vdim+1;
00586 
00587         metrics[next].setName(oss.str());
00588         metrics[next].setLocalSum(total);
00589         next++;
00590       }
00591     }
00592   }
00593 
00595   // Obtain global totals by part.
00596 
00597   try{
00598     reduceAll<int, scalar_t>(*comm, Teuchos::REDUCE_SUM, globalSumSize,
00599       localBuf, sumBuf);
00600   }
00601   Z2_THROW_OUTSIDE_ERROR(*env);
00602 
00603   delete [] localBuf;
00604 
00606   // Global sum, min, max, and average over all parts
00607 
00608   obj = sumBuf;                     // # of objects
00609   scalar_t min=0, max=0, sum=0;
00610   next = 0;
00611 
00612   ArrayView<scalar_t> objVec(obj, nparts);
00613   getStridedStats<scalar_t>(objVec, 1, 0, min, max, sum);
00614 
00615   metrics[next].setGlobalMin(min);
00616   metrics[next].setGlobalMax(max);
00617   metrics[next].setGlobalSum(sum);
00618   metrics[next].setGlobalAvg(sum / nparts);
00619   next++;
00620 
00621   if (numMetrics > 1){
00622     scalar_t *wgt = sumBuf + nparts;        // single normed weight
00623   
00624     ArrayView<scalar_t> normedWVec(wgt, nparts);
00625     getStridedStats<scalar_t>(normedWVec, 1, 0, min, max, sum);
00626 
00627     if (vwgtDim > 1)
00628       metrics[next].setNorm(multiCriteriaNorm(mcNorm));
00629 
00630     metrics[next].setGlobalMin(min);
00631     metrics[next].setGlobalMax(max);
00632     metrics[next].setGlobalSum(sum);
00633     metrics[next].setGlobalAvg(sum / nparts);
00634     next++;
00635   
00636     if (vwgtDim > 1){
00637       for (int vdim=0; vdim < vwgtDim; vdim++){
00638         wgt += nparts;       // individual weights
00639         ArrayView<scalar_t> fromVec(wgt, nparts);
00640         getStridedStats<scalar_t>(fromVec, 1, 0, min, max, sum);
00641 
00642         metrics[next].setGlobalMin(min);
00643         metrics[next].setGlobalMax(max);
00644         metrics[next].setGlobalSum(sum);
00645         metrics[next].setGlobalAvg(sum / nparts);
00646         next++;
00647       }
00648     }
00649   }
00650 
00652   // How many parts do we actually have.
00653 
00654   numParts = nparts;
00655   obj = sumBuf;               // # of objects
00656 
00657   for (partId_t p=nparts-1; p > 0; p--){
00658     if (obj[p] > 0) break;
00659     numParts--;
00660   }
00661 
00662   numNonemptyParts = numParts; 
00663 
00664   for (partId_t p=0; p < numParts; p++)
00665     if (obj[p] == 0) numNonemptyParts--;
00666 
00667   env->debug(DETAILED_STATUS, "Exiting globalSumsByPart");
00668 }
00669 
00700 template <typename scalar_t>
00701   void computeImbalances(partId_t numParts, partId_t targetNumParts,
00702     const scalar_t *psizes, scalar_t sumVals , const scalar_t *vals, 
00703     scalar_t &min, scalar_t &max, scalar_t &avg)
00704 {
00705   min = sumVals;
00706   max = avg = 0;
00707 
00708   if (sumVals <= 0 || targetNumParts < 1 || numParts < 1)
00709     return;
00710 
00711   if (targetNumParts==1 || numParts==1){
00712     min = max = avg = 0;  // 0 imbalance
00713     return;
00714   }
00715 
00716   if (!psizes){
00717     scalar_t target = sumVals / targetNumParts;
00718     for (partId_t p=0; p < numParts; p++){
00719       scalar_t diff = abs(vals[p] - target);
00720       scalar_t tmp = diff / target;
00721       avg += tmp;
00722       if (tmp > max) max = tmp;
00723       if (tmp < min) min = tmp;
00724     }
00725     partId_t emptyParts = targetNumParts - numParts;  
00726     if (emptyParts > 0){
00727       if (max < 1.0)
00728         max = 1.0;       // target divided by target
00729       avg += emptyParts;
00730     }
00731   }
00732   else{
00733     for (partId_t p=0; p < targetNumParts; p++){
00734       if (psizes[p] > 0){
00735         if (p < numParts){
00736           scalar_t target = sumVals * psizes[p];
00737           scalar_t diff = abs(vals[p] - target);
00738           scalar_t tmp = diff / target;
00739           avg += tmp;
00740           if (tmp > max) max = tmp;
00741           if (tmp < min) min = tmp;
00742         }
00743         else{
00744           if (max < 1.0)
00745             max = 1.0;  // target divided by target
00746           avg += 1.0;
00747         }
00748       }
00749     }
00750   }
00751 
00752   avg /= targetNumParts;
00753 }
00754 
00788 template <typename scalar_t>
00789  void computeImbalances(partId_t numParts, partId_t targetNumParts,
00790    int numSizes, ArrayView<ArrayRCP<scalar_t> > psizes,
00791    scalar_t sumVals , const scalar_t *vals, 
00792    scalar_t &min, scalar_t &max, scalar_t &avg)
00793 {
00794   min = sumVals;
00795   max = avg = 0;
00796 
00797   if (sumVals <= 0 || targetNumParts < 1 || numParts < 1)
00798     return;
00799 
00800   if (targetNumParts==1 && numParts==1){
00801     min = max = avg = 0;  // 0 imbalance
00802     return;
00803   }
00804 
00805   bool allUniformParts = true;
00806   for (int i=0; i < numSizes; i++){
00807     if (psizes[i].size() != 0){
00808       allUniformParts = false;
00809       break;
00810     }
00811   }
00812 
00813   if (allUniformParts){
00814     computeImbalances<scalar_t>(numParts, targetNumParts, NULL,
00815       sumVals, vals, min, max, avg);
00816     return;
00817   }
00818 
00819   double uniformSize = 1.0 / targetNumParts;
00820   ArrayRCP<double> sizeVec(new double [numSizes], 0, numSizes, true);
00821   for (int i=0; i < numSizes; i++){
00822     sizeVec[i] = uniformSize;
00823   }
00824 
00825   for (partId_t p=0; p < numParts; p++){
00826 
00827     // If we have objects in parts that should have 0 objects,
00828     // we don't compute an imbalance.  It means that other
00829     // parts have too few objects, and the imbalance will be
00830     // picked up there.
00831 
00832     if (p >= targetNumParts)
00833       break;
00834 
00835     // Vector of target amounts: T
00836 
00837     for (int i=0; i < numSizes; i++)
00838       if (psizes[i].size() > 0)
00839         sizeVec[i] = psizes[i][p];
00840 
00841     Epetra_SerialDenseVector target(View, sizeVec.getRawPtr(), numSizes);
00842     target.Scale(sumVals);
00843     double targetNorm = target.Norm2();
00844 
00845     // If part is supposed to be empty, we don't compute an
00846     // imbalance.  Same argument as above.
00847 
00848     if (targetNorm > 0){
00849 
00850       // Vector of actual amounts: A
00851 
00852       Epetra_SerialDenseVector actual(numSizes);
00853       for (int i=0; i < numSizes; i++)
00854         actual[i] = vals[p] * -1.0;
00855       
00856       actual += target;
00857 
00858       //  |A - T| / |T|
00859 
00860       scalar_t imbalance = actual.Norm2() / targetNorm;
00861 
00862       if (imbalance < min)
00863         min = imbalance;
00864       else if (imbalance > max)
00865         max = imbalance;
00866       avg += imbalance; 
00867     }
00868   }
00869 
00870   partId_t numEmptyParts = 0;
00871 
00872   for (partId_t p=numParts; p < targetNumParts; p++){
00873     bool nonEmptyPart = false;
00874     for (int i=0; !nonEmptyPart && i < numSizes; i++)
00875       if (psizes[i].size() > 0 && psizes[i][p] > 0.0)
00876         nonEmptyPart = true;
00877 
00878     if (nonEmptyPart){
00879       // The partitioning has no objects for this part, which
00880       // is supposed to be non-empty.
00881       numEmptyParts++;
00882     }
00883   }
00884 
00885   if (numEmptyParts > 0){
00886     avg += numEmptyParts;
00887     if (max < 1.0)
00888       max = 1.0;       // target divided by target
00889   }
00890 
00891   avg /= targetNumParts;
00892 }
00893 
00923 template <typename Adapter>
00924   void objectMetrics(
00925     const RCP<const Environment> &env,
00926     const RCP<const Comm<int> > &comm,
00927     multiCriteriaNorm mcNorm,
00928     const RCP<const Adapter> &ia,
00929     const RCP<const PartitioningSolution<Adapter> > &solution,
00930     partId_t &numParts,
00931     partId_t &numNonemptyParts,
00932     ArrayRCP<MetricValues<typename Adapter::scalar_t> > &metrics)
00933 {
00934   env->debug(DETAILED_STATUS, "Entering objectMetrics");
00935 
00936   typedef typename Adapter::scalar_t scalar_t;
00937   typedef typename Adapter::lno_t lno_t;
00938   typedef StridedData<lno_t, scalar_t> sdata_t;
00939 
00940   // Local number of objects.
00941 
00942   size_t numLocalObjects = ia->getLocalNumberOfObjects();
00943 
00944   // Parts to which objects are assigned.
00945 
00946   const partId_t *parts = solution->getPartList();
00947   env->localInputAssertion(__FILE__, __LINE__, "parts not set", 
00948     parts, BASIC_ASSERTION);
00949   ArrayView<const partId_t> partArray(parts, numLocalObjects);
00950 
00951   // Weights, if any, for each object.
00952 
00953   int weightDim = ia->getNumberOfWeightsPerObject();
00954   int numCriteria = (weightDim > 0 ? weightDim : 1);
00955   Array<sdata_t> weights(numCriteria);
00956 
00957   if (weightDim == 0){
00958     // One dimension of uniform weights is implied.
00959     // StridedData default constructor creates length 0 strided array.
00960     weights[0] = sdata_t();
00961   }
00962   else{
00963     for (int i=0; i < weightDim; i++){
00964       int stride;
00965       const scalar_t *wgt;
00966       size_t len = ia->getObjectWeights(i, wgt, stride); 
00967       ArrayRCP<const scalar_t> wgtArray(wgt, 0, len, false);
00968       weights[i] = sdata_t(wgtArray, stride);
00969     }
00970   }
00971 
00972   // Relative part sizes, if any, assigned to the parts.
00973 
00974   partId_t targetNumParts = solution->getTargetGlobalNumberOfParts();
00975   scalar_t *psizes = NULL;
00976 
00977   ArrayRCP<ArrayRCP<scalar_t> > partSizes(numCriteria);
00978   for (int dim=0; dim < numCriteria; dim++){
00979     if (solution->criteriaHasUniformPartSizes(dim) != true){
00980       psizes = new scalar_t [targetNumParts];
00981       env->localMemoryAssertion(__FILE__, __LINE__, numParts, psizes);
00982       for (partId_t i=0; i < targetNumParts; i++){
00983         psizes[i] = solution->getCriteriaPartSize(dim, i);
00984       }
00985       partSizes[dim] = arcp(psizes, 0, targetNumParts, true);
00986     }
00987   }
00988 
00990   // Get number of parts, and the number that are non-empty.
00991   // Get sums per part of objects, individual weights, and normed weight sums.
00992 
00993   ArrayRCP<scalar_t> globalSums;
00994 
00995   try{
00996     globalSumsByPart<scalar_t, partId_t, lno_t>(env, comm, 
00997       partArray, weights.view(0, numCriteria), mcNorm,
00998       numParts, numNonemptyParts, metrics, globalSums);
00999   }
01000   Z2_FORWARD_EXCEPTIONS
01001 
01003   // Compute imbalances for the object count.  
01004   // (Use first dimension of part sizes.) 
01005 
01006   scalar_t *objCount  = globalSums.getRawPtr();
01007   scalar_t min, max, avg;
01008   psizes=NULL;
01009 
01010   if (partSizes[0].size() > 0)
01011     psizes = partSizes[0].getRawPtr();
01012 
01013   computeImbalances<scalar_t>(numParts, targetNumParts, psizes,
01014       metrics[0].getGlobalSum(), objCount, 
01015       min, max, avg);
01016 
01017   metrics[0].setMinImbalance(1.0 + min);
01018   metrics[0].setMaxImbalance(1.0 + max);
01019   metrics[0].setAvgImbalance(1.0 + avg);
01020 
01022   // Compute imbalances for the normed weight sum.
01023 
01024   scalar_t *wgts = globalSums.getRawPtr() + numParts;
01025 
01026   if (metrics.size() > 1){
01027   
01028     computeImbalances<scalar_t>(numParts, targetNumParts, 
01029       numCriteria, partSizes.view(0, numCriteria),
01030       metrics[1].getGlobalSum(), wgts,
01031       min, max, avg);
01032   
01033     metrics[1].setMinImbalance(1.0 + min);
01034     metrics[1].setMaxImbalance(1.0 + max);
01035     metrics[1].setAvgImbalance(1.0 + avg);
01036 
01037     if (metrics.size() > 2){
01038 
01040     // Compute imbalances for each individual weight dimension.
01041 
01042       int next = 2;
01043   
01044       for (int vdim=0; vdim < numCriteria; vdim++){
01045         wgts += numParts;
01046         psizes = NULL;
01047   
01048         if (partSizes[vdim].size() > 0)
01049            psizes = partSizes[vdim].getRawPtr();
01050          
01051         computeImbalances<scalar_t>(numParts, targetNumParts, psizes,
01052           metrics[next].getGlobalSum(), wgts, min, max, avg);
01053   
01054         metrics[next].setMinImbalance(1.0 + min);
01055         metrics[next].setMaxImbalance(1.0 + max);
01056         metrics[next].setAvgImbalance(1.0 + avg);
01057         next++;
01058       }
01059     }
01060 
01061   }
01062   env->debug(DETAILED_STATUS, "Exiting objectMetrics");
01063 }
01064 
01068 template <typename scalar_t>
01069   void printMetrics( ostream &os,
01070     partId_t targetNumParts, partId_t numParts, partId_t numNonemptyParts, 
01071     const ArrayView<MetricValues<scalar_t> > &infoList)
01072 {
01073   os << "NUMBER OF PARTS IS " << numParts;
01074   if (numNonemptyParts < numParts)
01075     os << " (" << numNonemptyParts << " of which are non-empty)";
01076   os << endl;
01077   if (targetNumParts != numParts)
01078     os << "TARGET NUMBER OF PARTS IS " << targetNumParts << std::endl;
01079 
01080   string unset("unset");
01081 
01082   MetricValues<scalar_t>::printHeader(os);
01083 
01084   for (int i=0; i < infoList.size(); i++)
01085     if (infoList[i].getName() != unset)
01086       infoList[i].printLine(os);
01087 
01088   os << endl;
01089 }
01090 
01094 template <typename scalar_t>
01095   void printMetrics( ostream &os,
01096     partId_t targetNumParts, partId_t numParts, partId_t numNonemptyParts, 
01097     const MetricValues<scalar_t> &info)
01098 {
01099   ArrayView<MetricValues<scalar_t> > infoList(&info, 1);
01100   printMetrics( os, targetNumParts, numParts, numNonemptyParts, infoList);
01101 }
01102 
01103 
01106 template <typename scalar_t>
01107   scalar_t normedWeight(ArrayView <scalar_t> weights, 
01108     multiCriteriaNorm norm)
01109 {
01110   size_t dim = weights.size();
01111   if (dim == 0)
01112     return 0.0;
01113   else if (dim == 1)
01114     return weights[0];
01115  
01116   scalar_t nweight = 0;
01117 
01118   switch (norm){
01119     case normMinimizeTotalWeight:   
01121       for (size_t i=0; i <dim; i++)
01122         nweight += weights[i];
01123 
01124       break;
01125      
01126     case normBalanceTotalMaximum:   
01127       for (size_t i=0; i <dim; i++)
01128         nweight += (weights[i] * weights[i]);
01129 
01130       nweight = sqrt(nweight);
01131 
01132       break;
01133 
01134     case normMinimizeMaximumWeight: 
01136       nweight = weights[0];
01137       for (size_t i=1; i <dim; i++)
01138         if (weights[i] > nweight)
01139           nweight = weights[i];
01140 
01141       break;
01142 
01143     default:
01144       Environment env;  // default environment
01145       env.localBugAssertion(__FILE__, __LINE__, "invalid norm", false,
01146         BASIC_ASSERTION);
01147       break;
01148   } 
01149 
01150   return nweight;
01151 }
01152 
01155 template<typename lno_t, typename scalar_t>
01156   scalar_t normedWeight(ArrayView<StridedData<lno_t, scalar_t> > weights, 
01157      lno_t idx, multiCriteriaNorm norm)
01158 {
01159   size_t dim = weights.size();
01160   if (dim < 1)
01161     return 0;
01162 
01163   Array<scalar_t> vec(dim, 1.0);
01164   for (size_t i=0; i < dim; i++)
01165     if (weights[i].size() > 0)
01166        vec[dim] = weights[i][idx];
01167 
01168   return normedWeight(vec, norm);
01169 }
01170 
01171 } //namespace Zoltan2
01172 #endif