00001
00002
00003
00004
00005 #ifndef PLAYA_VECTOROPSIMPL_HPP
00006 #define PLAYA_VECTOROPSIMPL_HPP
00007
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaVectorImpl.hpp"
00010 #include "PlayaVectorOpsDecl.hpp"
00011 #include "PlayaMPIComm.hpp"
00012 #include "Teuchos_RCP.hpp"
00013
00014 namespace PlayaFunctors
00015 {
00016 template <class Scalar> class BoundedMinLocFunctor;
00017
00018 template <class Scalar> class BoundedMaxLocFunctor;
00019
00020 template <class Scalar> class Norm2Dist;
00021
00022 template <class Scalar> class Norm1Dist;
00023
00024 template <class Scalar> class NormInfDist;
00025 }
00026
00027 namespace Playa
00028 {
00029
00030 using namespace PlayaFunctors;
00031
00032
00033 template <class Scalar> inline
00034 Scalar minloc(const Vector<Scalar>& x, int& gni)
00035 {
00036 return minlocWithBound(-HUGE_VAL, x, gni);
00037 }
00038
00039
00040 template <class Scalar> inline
00041 Scalar maxloc(const Vector<Scalar>& x, int& gni)
00042 {
00043 return maxlocWithBound(-HUGE_VAL, x, gni);
00044 }
00045
00046
00047 template <class Scalar> inline
00048 Scalar minlocWithBound(const Scalar& lowerBound,
00049 const Vector<Scalar>& x, int& gni)
00050 {
00051 IndexedValue<Scalar> y =
00052 x.applyUnaryReductionFunctor(BoundedMinLocFunctor<Scalar>(x.comm(), lowerBound, x.space().baseGlobalNaturalIndex()));
00053 gni = y.where;
00054 return y.what;
00055 }
00056
00057
00058 template <class Scalar> inline
00059 Scalar maxlocWithBound(const Scalar& upperBound,
00060 const Vector<Scalar>& x, int& gni)
00061 {
00062 IndexedValue<Scalar> y =
00063 x.applyUnaryReductionFunctor(BoundedMaxLocFunctor<Scalar>(x.comm(), upperBound, x.space().baseGlobalNaturalIndex()));
00064 gni = y.where;
00065 return y.what;
00066 }
00067
00068
00069
00070 template <class Scalar>
00071 Scalar norm2Dist(const Vector<Scalar>& x, const Vector<Scalar>& y)
00072 {
00073 return x.applyBinaryReductionFunctor(
00074 PlayaFunctors::Norm2Dist<Scalar>(x.comm()), y);
00075 }
00076
00077
00078 template <class Scalar>
00079 Scalar norm1Dist(const Vector<Scalar>& x, const Vector<Scalar>& y)
00080 {
00081 return x.applyBinaryReductionFunctor(
00082 PlayaFunctors::Norm1Dist<Scalar>(x.comm()), y);
00083 }
00084
00085
00086 template <class Scalar>
00087 Scalar normInfDist(const Vector<Scalar>& x, const Vector<Scalar>& y)
00088 {
00089 return x.applyBinaryReductionFunctor(
00090 PlayaFunctors::NormInfDist<Scalar>(x.comm()), y);
00091 }
00092
00093
00094 template <class Scalar>
00095 Scalar norm2(const Vector<Scalar>& x) {return x.norm2();}
00096
00097
00098 template <class Scalar>
00099 Scalar norm1(const Vector<Scalar>& x) {return x.norm1();}
00100
00101
00102 template <class Scalar>
00103 Scalar normInf(const Vector<Scalar>& x) {return x.normInf();}
00104
00105 }
00106
00107 namespace PlayaFunctors
00108 {
00109
00110
00111 using namespace Playa;
00112
00113
00114
00115 template <class Scalar>
00116 class BoundedMinLocFunctor : public ReductionFunctorBase<Scalar>
00117 {
00118 public:
00119
00120 BoundedMinLocFunctor(const MPIComm& comm, const Scalar& bound,
00121 int baseGNI)
00122 : ReductionFunctorBase<Scalar>(comm), min_(),
00123 bound_(bound), baseGNI_(baseGNI)
00124 {
00125 min_.what = HUGE_VAL;
00126 min_.where = -1;
00127 }
00128
00129
00130 void step(int i, const Scalar& x) const
00131 {
00132 if (x < min_.what && x > bound_)
00133 {
00134 min_.what = x;
00135 min_.where = i;
00136 }
00137 }
00138
00139
00140 void postProc() const
00141 {
00142 min_.where += baseGNI_;
00143
00144 IndexedValue<Scalar> out = min_;
00145
00146 this->comm().allReduce(&min_, &out, 1, MPIDataType::doubleIntPairType(),
00147 MPIOp::minlocOp());
00148 min_ = out;
00149 }
00150
00151
00152 IndexedValue<Scalar> result() const
00153 {
00154 return min_;
00155 }
00156
00157
00158 std::string description() const {return "MinLoc()";}
00159
00160 private:
00161 MPIComm comm_;
00162 mutable IndexedValue<Scalar> min_;
00163 Scalar bound_;
00164 int baseGNI_;
00165 };
00166
00167
00168
00169
00170 template <class Scalar>
00171 class BoundedMaxLocFunctor : public ReductionFunctorBase<Scalar>
00172 {
00173 public:
00174
00175 BoundedMaxLocFunctor(const MPIComm& comm, const Scalar& bound,
00176 int baseGNI)
00177 : ReductionFunctorBase<Scalar>(comm), max_(),
00178 bound_(bound), baseGNI_(baseGNI)
00179 {
00180 max_.what = -HUGE_VAL;
00181 max_.where = -1;
00182 }
00183
00184
00185 void step(int i, const Scalar& x) const
00186 {
00187 if (x > max_.what && x < bound_)
00188 {
00189 max_.what = x;
00190 max_.where = i;
00191 }
00192 }
00193
00194
00195 void postProc() const
00196 {
00197 max_.where += baseGNI_;
00198
00199 IndexedValue<Scalar> out = max_;
00200
00201 this->comm().allReduce(&max_, &out, 1, MPIDataType::doubleIntPairType(),
00202 MPIOp::minlocOp());
00203 max_ = out;
00204 }
00205
00206
00207 IndexedValue<Scalar> result() const
00208 {
00209 return max_;
00210 }
00211
00212
00213 std::string description() const {return "MinLoc()";}
00214
00215 private:
00216 MPIComm comm_;
00217 mutable IndexedValue<Scalar> max_;
00218 Scalar bound_;
00219 int baseGNI_;
00220 };
00221
00222
00223
00224
00225 template <class Scalar>
00226 class Norm2Dist : public ReductionFunctorBase<Scalar>
00227 {
00228 public:
00229 Norm2Dist(const MPIComm& comm)
00230 : ReductionFunctorBase<Scalar>(comm), val_(0.0) {}
00231
00232 void step(int i, const Scalar& x, const Scalar& y) const
00233 {
00234 Scalar d = x-y;
00235 val_ += d*d;
00236 }
00237
00238 void postProc() const
00239 {
00240 Scalar final = val_;
00241 this->comm().allReduce(&val_, &final, 1, MPIDataType::doubleType(), MPIOp::sumOp());
00242 val_ = final;
00243 }
00244
00245 Scalar result() const
00246 {
00247 return ::sqrt(val_);
00248 }
00249
00250
00251 std::string description() const {return "Norm2Dist()";}
00252
00253 private:
00254 mutable Scalar val_;
00255 };
00256
00257
00258
00259
00260 template <class Scalar>
00261 class Norm1Dist : public ReductionFunctorBase<Scalar>
00262 {
00263 public:
00264 Norm1Dist(const MPIComm& comm)
00265 : ReductionFunctorBase<Scalar>(comm), val_(0.0) {}
00266
00267 void step(int i, const Scalar& x, const Scalar& y) const
00268 {
00269 Scalar d = x-y;
00270 val_ += ::fabs(d);
00271 }
00272
00273 void postProc() const
00274 {
00275 Scalar final = val_;
00276 this->comm().allReduce(&val_, &final, 1, MPIDataType::doubleType(), MPIOp::sumOp());
00277 val_ = final;
00278 }
00279
00280 Scalar result() const
00281 {
00282 return val_;
00283 }
00284
00285
00286 std::string description() const {return "Norm1Dist()";}
00287
00288 private:
00289 mutable Scalar val_;
00290 };
00291
00292
00293 template <class Scalar>
00294 class NormInfDist : public ReductionFunctorBase<Scalar>
00295 {
00296 public:
00297 NormInfDist(const MPIComm& comm)
00298 : ReductionFunctorBase<Scalar>(comm), val_(0.0) {}
00299
00300 void step(int i, const Scalar& x, const Scalar& y) const
00301 {
00302 Scalar d = ::fabs(x-y);
00303 if (d > val_) val_ = d;
00304 }
00305
00306 void postProc() const
00307 {
00308 Scalar final = val_;
00309 this->comm().allReduce(&val_, &final, 1, MPIDataType::doubleType(), MPIOp::maxOp());
00310 val_ = final;
00311 }
00312
00313 Scalar result() const
00314 {
00315 return val_;
00316 }
00317
00318
00319 std::string description() const {return "NormInfDist()";}
00320
00321 private:
00322 mutable Scalar val_;
00323 };
00324
00325
00326
00327 template <class Scalar>
00328 class VectorFunctorTraits<Scalar, BoundedMinLocFunctor<Scalar> >
00329 {
00330 public:
00331 typedef IndexedValue<Scalar> ReturnType ;
00332 };
00333
00334
00335
00336 template <class Scalar>
00337 class VectorFunctorTraits<Scalar, BoundedMaxLocFunctor<Scalar> >
00338 {
00339 public:
00340 typedef IndexedValue<Scalar> ReturnType ;
00341 };
00342
00343
00344 }
00345
00346 #endif