escript
Revision_
|
00001 00002 /***************************************************************************** 00003 * 00004 * Copyright (c) 2003-2014 by University of Queensland 00005 * http://www.uq.edu.au 00006 * 00007 * Primary Business: Queensland, Australia 00008 * Licensed under the Open Software License version 3.0 00009 * http://www.opensource.org/licenses/osl-3.0.php 00010 * 00011 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 00012 * Development 2012-2013 by School of Earth Sciences 00013 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 00014 * 00015 *****************************************************************************/ 00016 00017 00018 #if !defined escript_DataAlgorithm_20040714_H 00019 #define escript_DataAlgorithm_20040714_H 00020 #include "system_dep.h" 00021 00022 #include "DataExpanded.h" 00023 #include "DataTagged.h" 00024 #include "DataConstant.h" 00025 00026 #include "DataMaths.h" 00027 00028 #include <iostream> 00029 #include <algorithm> 00030 #include <list> 00031 00032 namespace escript { 00033 00044 template <class BinaryFunction> 00045 class DataAlgorithmAdapter { 00046 public: 00047 DataAlgorithmAdapter(double initialValue): 00048 m_initialValue(initialValue), 00049 m_currentValue(initialValue) 00050 { 00051 } 00052 DataAlgorithmAdapter(const DataAlgorithmAdapter& other): 00053 m_initialValue(other.m_initialValue), 00054 m_currentValue(other.m_initialValue), 00055 operation(other.operation) 00056 { 00057 } 00058 inline void operator()(double value) 00059 { 00060 m_currentValue=operation(m_currentValue,value); 00061 return; 00062 } 00063 inline void resetResult() 00064 { 00065 m_currentValue=m_initialValue; 00066 } 00067 inline double getResult() const 00068 { 00069 return m_currentValue; 00070 } 00071 private: 00072 // 00073 // the initial operation value 00074 double m_initialValue; 00075 // 00076 // the current operation value 00077 double m_currentValue; 00078 // 00079 // The operation to perform 00080 BinaryFunction operation; 00081 }; 00082 00087 struct FMax : public std::binary_function<double,double,double> 00088 { 00089 inline double operator()(double x, double y) const 00090 { 00091 return std::max(x,y); 00092 } 00093 }; 00094 00099 struct FMin : public std::binary_function<double,double,double> 00100 { 00101 inline double operator()(double x, double y) const 00102 { 00103 return std::min(x,y); 00104 } 00105 }; 00106 00111 struct AbsMax : public std::binary_function<double,double,double> 00112 { 00113 inline double operator()(double x, double y) const 00114 { 00115 return std::max(fabs(x),fabs(y)); 00116 } 00117 }; 00118 00123 struct AbsMin : public std::binary_function<double,double,double> 00124 { 00125 inline double operator()(double x, double y) const 00126 { 00127 return std::min(fabs(x),fabs(y)); 00128 } 00129 }; 00130 00135 struct Length : public std::binary_function<double,double,double> 00136 { 00137 inline double operator()(double x, double y) const 00138 { 00139 return std::sqrt(std::pow(x,2)+std::pow(y,2)); 00140 } 00141 }; 00142 00147 struct Trace : public std::binary_function<double,double,double> 00148 { 00149 inline double operator()(double x, double y) const 00150 { 00151 return x+y; 00152 } 00153 }; 00154 00158 struct AbsGT : public std::binary_function<double,double,double> 00159 { 00160 inline double operator()(double x, double y) const 00161 { 00162 return fabs(x)>y; 00163 } 00164 }; 00165 00169 struct AbsLTE : public std::binary_function<double,double,double> 00170 { 00171 inline double operator()(double x, double y) const 00172 { 00173 return fabs(x)<=y; 00174 } 00175 }; 00176 00177 00183 template <class BinaryFunction> 00184 inline 00185 double 00186 algorithm(const DataExpanded& data, 00187 BinaryFunction operation, 00188 double initial_value) 00189 { 00190 int i,j; 00191 int numDPPSample=data.getNumDPPSample(); 00192 int numSamples=data.getNumSamples(); 00193 double global_current_value=initial_value; 00194 double local_current_value; 00195 // DataArrayView dataView=data.getPointDataView(); 00196 const DataTypes::ValueType& vec=data.getVectorRO(); 00197 const DataTypes::ShapeType& shape=data.getShape(); 00198 // calculate the reduction operation value for each data point 00199 // reducing the result for each data-point into the current_value variables 00200 #pragma omp parallel private(local_current_value) 00201 { 00202 local_current_value=initial_value; 00203 #pragma omp for private(i,j) schedule(static) 00204 for (i=0;i<numSamples;i++) { 00205 for (j=0;j<numDPPSample;j++) { 00206 /* local_current_value=operation(local_current_value,dataView.reductionOp(data.getPointOffset(i,j),operation,initial_value));*/ 00207 local_current_value=operation(local_current_value,DataMaths::reductionOp(vec,shape,data.getPointOffset(i,j),operation,initial_value)); 00208 00209 } 00210 } 00211 #pragma omp critical 00212 global_current_value=operation(global_current_value,local_current_value); 00213 } 00214 return global_current_value; 00215 } 00216 00217 // It is important that the algorithm only be applied to tags which are actually in use. 00218 template <class BinaryFunction> 00219 inline 00220 double 00221 algorithm(DataTagged& data, 00222 BinaryFunction operation, 00223 double initial_value) 00224 { 00225 double current_value=initial_value; 00226 00227 const DataTypes::ValueType& vec=data.getVectorRO(); 00228 const DataTypes::ShapeType& shape=data.getShape(); 00229 const DataTagged::DataMapType& lookup=data.getTagLookup(); 00230 const std::list<int> used=data.getFunctionSpace().getListOfTagsSTL(); 00231 for (std::list<int>::const_iterator i=used.begin();i!=used.end();++i) 00232 { 00233 int tag=*i; 00234 if (tag==0) // check for the default tag 00235 { 00236 current_value=operation(current_value,DataMaths::reductionOp(vec,shape,data.getDefaultOffset(),operation,initial_value)); 00237 } 00238 else 00239 { 00240 DataTagged::DataMapType::const_iterator it=lookup.find(tag); 00241 if (it!=lookup.end()) 00242 { 00243 current_value=operation(current_value,DataMaths::reductionOp(vec,shape,it->second,operation,initial_value)); 00244 } 00245 } 00246 } 00247 return current_value; 00248 } 00249 00250 template <class BinaryFunction> 00251 inline 00252 double 00253 algorithm(DataConstant& data, 00254 BinaryFunction operation, 00255 double initial_value) 00256 { 00257 return DataMaths::reductionOp(data.getVectorRO(),data.getShape(),0,operation,initial_value); 00258 } 00259 00271 template <class BinaryFunction> 00272 inline 00273 void 00274 dp_algorithm(const DataExpanded& data, 00275 DataExpanded& result, 00276 BinaryFunction operation, 00277 double initial_value) 00278 { 00279 int i,j; 00280 int numSamples=data.getNumSamples(); 00281 int numDPPSample=data.getNumDPPSample(); 00282 // DataArrayView dataView=data.getPointDataView(); 00283 // DataArrayView resultView=result.getPointDataView(); 00284 const DataTypes::ValueType& dataVec=data.getVectorRO(); 00285 const DataTypes::ShapeType& shape=data.getShape(); 00286 DataTypes::ValueType& resultVec=result.getVectorRW(); 00287 // perform the operation on each data-point and assign 00288 // this to the corresponding element in result 00289 #pragma omp parallel for private(i,j) schedule(static) 00290 for (i=0;i<numSamples;i++) { 00291 for (j=0;j<numDPPSample;j++) { 00292 /* resultView.getData(result.getPointOffset(i,j)) = 00293 dataView.reductionOp(data.getPointOffset(i,j),operation,initial_value);*/ 00294 resultVec[result.getPointOffset(i,j)] = 00295 DataMaths::reductionOp(dataVec, shape, data.getPointOffset(i,j),operation,initial_value); 00296 00297 } 00298 } 00299 } 00300 00301 template <class BinaryFunction> 00302 inline 00303 void 00304 dp_algorithm(const DataTagged& data, 00305 DataTagged& result, 00306 BinaryFunction operation, 00307 double initial_value) 00308 { 00309 // perform the operation on each tagged value in data 00310 // and assign this to the corresponding element in result 00311 const DataTypes::ShapeType& shape=data.getShape(); 00312 const DataTypes::ValueType& vec=data.getVectorRO(); 00313 const DataTagged::DataMapType& lookup=data.getTagLookup(); 00314 for (DataTagged::DataMapType::const_iterator i=lookup.begin(); i!=lookup.end(); i++) { 00315 result.getDataByTagRW(i->first,0) = 00316 DataMaths::reductionOp(vec,shape,data.getOffsetForTag(i->first),operation,initial_value); 00317 } 00318 // perform the operation on the default data value 00319 // and assign this to the default element in result 00320 result.getVectorRW()[result.getDefaultOffset()] = DataMaths::reductionOp(data.getVectorRO(),data.getShape(),data.getDefaultOffset(),operation,initial_value); 00321 } 00322 00323 template <class BinaryFunction> 00324 inline 00325 void 00326 dp_algorithm(DataConstant& data, 00327 DataConstant& result, 00328 BinaryFunction operation, 00329 double initial_value) 00330 { 00331 // perform the operation on the data value 00332 // and assign this to the element in result 00333 result.getVectorRW()[0] = 00334 DataMaths::reductionOp(data.getVectorRO(),data.getShape(),0,operation,initial_value); 00335 } 00336 00337 } // end of namespace 00338 00339 #endif