# Diff of /trunk/escript/src/DataAlgorithm.h

revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC revision 122 by jgs, Thu Jun 9 05:38:05 2005 UTC
# Line 25  Line 25
25  #include <algorithm>  #include <algorithm>
26  #include <math.h>  #include <math.h>
27  #include <limits>  #include <limits>
28    #include <vector>
29
30  namespace escript {  namespace escript {
31
32  /**  /**
33     \brief     \brief
34     Return the maximum value.     Adapt binary algorithms so they may be used in DataArrayView reduction operations.
35
36     Description:     Description:
37     Return the maximum value.     This functor adapts the given BinaryFunction operation by starting with the
38       given inital value applying this operation to successive values, storing the
39       rolling result in m_currentValue - which can be accessed or reset by getResult
40       and resetResult respectively.
41    */
42    template <class BinaryFunction>
44      public:
46          m_initialValue(initialValue),
47          m_currentValue(initialValue)
48        {
49        }
50        inline void operator()(double value)
51        {
52          m_currentValue=operation(m_currentValue,value);
53          return;
54        }
55        inline void resetResult()
56        {
57          m_currentValue=m_initialValue;
58        }
59        inline double getResult() const
60        {
61          return m_currentValue;
62        }
63      private:
64        //
65        // the initial operation value
66        double m_initialValue;
67        //
68        // the current operation value
69        double m_currentValue;
70        //
71        // The operation to perform
72        BinaryFunction operation;
73    };
74
75    /**
76       \brief
77       Return the maximum value of the two given values.
78  */  */
79  struct FMax : public std::binary_function<double,double,double>  struct FMax : public std::binary_function<double,double,double>
80  {  {
# Line 45  struct FMax : public std::binary_functio Line 86  struct FMax : public std::binary_functio
86
87  /**  /**
88     \brief     \brief
89     Return the minimum value.     Return the minimum value of the two given values.

Description:
Return the minimum value.
90  */  */
91  struct FMin : public std::binary_function<double,double,double>  struct FMin : public std::binary_function<double,double,double>
92  {  {
# Line 60  struct FMin : public std::binary_functio Line 98  struct FMin : public std::binary_functio
98
99  /**  /**
100     \brief     \brief
101     Return the absolute maximum value.     Return the absolute maximum value of the two given values.

Description:
Return the absolute maximum value.
102  */  */
103  struct AbsMax : public std::binary_function<double,double,double>  struct AbsMax : public std::binary_function<double,double,double>
104  {  {
# Line 75  struct AbsMax : public std::binary_funct Line 110  struct AbsMax : public std::binary_funct
110
111  /**  /**
112     \brief     \brief
113     Return the length.     Return the absolute minimum value of the two given values.

Description:
Return the length.
114  */  */
115  struct Length : public std::binary_function<double,double,double>  struct AbsMin : public std::binary_function<double,double,double>
116  {  {
117    inline double operator()(double x, double y) const    inline double operator()(double x, double y) const
118    {    {
119      return std::sqrt(std::pow(x,2)+std::pow(y,2));      return std::min(fabs(x),fabs(y));
120    }    }
121  };  };
122
123  /**  /**
124     \brief     \brief
125     Return the trace.     Return the length between the two given values.

Description:
Return the trace.
126  */  */
127  struct Trace : public std::binary_function<double,double,double>  struct Length : public std::binary_function<double,double,double>
128  {  {
129    inline double operator()(double x, double y) const    inline double operator()(double x, double y) const
130    {    {
131      return x+y;      return std::sqrt(std::pow(x,2)+std::pow(y,2));
132    }    }
133  };  };
134
135  /**  /**
136     \brief     \brief
137     Adapt algorithms so they may be used by Data.     Return the trace of the two given values.

Description:
Adapt algorithms so they may be used by Data. The functor
maintains state, the currentValue returned by the operation,
and the initial value.
138  */  */
139  template <class BinaryFunction>  struct Trace : public std::binary_function<double,double,double>
141   public:    inline double operator()(double x, double y) const
143        m_initialValue(initialValue),      return x+y;
144        m_currentValue(initialValue)    }
{}
inline void operator()(double value)
{
m_currentValue=operation(m_currentValue,value);
return;
}
inline void resetResult()
{
m_currentValue=m_initialValue;
}
inline double getResult() const
{
return m_currentValue;
}
private:
//
// the initial operation value
double m_initialValue;
//
// the current operation value
double m_currentValue;
//
// The operation to perform
BinaryFunction operation;
145  };  };
146
147  /**  /**
148     \brief     \brief
149     Perform the given operation upon all Data elements and return a single     Perform the given operation upon all values in all data-points in the
150     result.     given Data object and return the final result.
151
152     Description:     Calls DataArrayView::reductionOp
Perform the given operation upon all Data elements and return a single
result.
153  */  */
154  template <class UnaryFunction>  template <class UnaryFunction>
155  inline  inline
# Line 160  algorithm(DataExpanded& data, Line 158  algorithm(DataExpanded& data,
158            UnaryFunction operation)            UnaryFunction operation)
159  {  {
160    int i,j;    int i,j;
161    DataArrayView::ValueType::size_type numDPPSample=data.getNumDPPSample();    int numDPPSample=data.getNumDPPSample();
162    DataArrayView::ValueType::size_type numSamples=data.getNumSamples();    int numSamples=data.getNumSamples();
163    double resultLocal=0;    int resultVectorLength=numDPPSample*numSamples;
164  #pragma omp parallel private(resultLocal)    std::vector<double> resultVector(resultVectorLength);
165    {    DataArrayView dataView=data.getPointDataView();
166  #pragma omp for private(i,j) schedule(static)    // calculate the reduction operation value for each data point
167      for (i=0;i<numSamples;i++) {    // storing the result for each data-point in successive entries
168        for (j=0;j<numDPPSample;j++) {    // in resultVector
169      resultLocal=data.getPointDataView().reductionOp(data.getPointOffset(i,j), operation);    //
170  #pragma omp critical (algorithm)    // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
171      operation(resultLocal);    // which maintains state between calls which would be corrupted by parallel execution
172        }    for (i=0;i<numSamples;i++) {
173        for (j=0;j<numDPPSample;j++) {
174          resultVector[j*numSamples+i]=dataView.reductionOp(data.getPointOffset(i,j),operation);
175      }      }
176    }    }
177      // now calculate the reduction operation value across the results
178      // for each data-point
179      //
180      // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
181      // which maintains state between calls which would be corrupted by parallel execution
182      operation.resetResult();
183      for (int l=0;l<resultVectorLength;l++) {
184        operation(resultVector[l]);
185      }
186    return operation.getResult();    return operation.getResult();
187  }  }
188
# Line 183  double Line 192  double
192  algorithm(DataTagged& data,  algorithm(DataTagged& data,
193            UnaryFunction operation)            UnaryFunction operation)
194  {  {
//
// perform the operation on each tagged value
195    const DataTagged::DataMapType& lookup=data.getTagLookup();    const DataTagged::DataMapType& lookup=data.getTagLookup();
196    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
197    DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();    DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();
198    DataArrayView& dataView=data.getPointDataView();    DataArrayView& dataView=data.getPointDataView();
199      std::vector<double> resultVector;
200      int resultVectorLength;
201      // perform the operation on each tagged value
202    for (i=lookup.begin();i!=lookupEnd;i++) {    for (i=lookup.begin();i!=lookupEnd;i++) {
203      operation(dataView.reductionOp(i->second,operation));      resultVector.push_back(dataView.reductionOp(i->second,operation));
204    }    }
205      // perform the operation on the default value
206      resultVector.push_back(data.getDefaultValue().reductionOp(operation));
207      // now calculate the reduction operation value across the results
208      // for each tagged value
209    //    //
210    // finally perform the operation on the default value    // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
211    operation(data.getDefaultValue().reductionOp(operation));    // which maintains state between calls which would be corrupted by parallel execution
212      resultVectorLength=resultVector.size();
213      operation.resetResult();
214      for (int l=0;l<resultVectorLength;l++) {
215        operation(resultVector[l]);
216      }
217    return operation.getResult();    return operation.getResult();
218  }  }
219
# Line 209  algorithm(DataConstant& data, Line 228  algorithm(DataConstant& data,
228
229  /**  /**
230     \brief     \brief
231     Perform the given data point reduction operation on all data points     Perform the given data-point reduction operation on all data-points
232     in data, storing results in corresponding elements of result.     in data, storing results in corresponding data-points of result.
233
234     Objects data and result must be of the same type, and have the same number     Objects data and result must be of the same type, and have the same number
235     of data points, but where data has data points of rank n, result must have     of data points, but where data has data points of rank n, result must have
236     data points of rank 0.     data points of rank 0.
237
238     Calls DataArrayView::dp_algorithm.     Calls DataArrayView::reductionOp
239  */  */
240  template <class UnaryFunction>  template <class UnaryFunction>
241  inline  inline
# Line 225  dp_algorithm(DataExpanded& data, Line 244  dp_algorithm(DataExpanded& data,
244               DataExpanded& result,               DataExpanded& result,
245               UnaryFunction operation)               UnaryFunction operation)
246  {  {
//
// perform the operation on each data value
// and assign this to the corresponding element in result
247    int i,j;    int i,j;
248    DataArrayView::ValueType::size_type numDPPSample=data.getNumDPPSample();    int numSamples=data.getNumSamples();
249    DataArrayView::ValueType::size_type numSamples=data.getNumSamples();    int numDPPSample=data.getNumDPPSample();
250    {    DataArrayView dataView=data.getPointDataView();
251  #pragma omp for private(i,j) schedule(static)    DataArrayView resultView=result.getPointDataView();
252      for (i=0;i<numSamples;i++) {    // perform the operation on each data-point and assign
253        for (j=0;j<numDPPSample;j++) {    // this to the corresponding element in result
254  #pragma omp critical (dp_algorithm)    //
255          result.getPointDataView().getData(data.getPointOffset(i,j)) =    // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
256            data.getPointDataView().dp_reductionOp(data.getPointOffset(i,j),operation);    // which maintains state between calls which would be corrupted by parallel execution
257        }    for (i=0;i<numSamples;i++) {
258        for (j=0;j<numDPPSample;j++) {
259          resultView.getData(result.getPointOffset(i,j)) =
260            dataView.reductionOp(data.getPointOffset(i,j),operation);
261      }      }
262    }    }
263  }  }
# Line 250  dp_algorithm(DataTagged& data, Line 269  dp_algorithm(DataTagged& data,
269               DataTagged& result,               DataTagged& result,
270               UnaryFunction operation)               UnaryFunction operation)
271  {  {
//
// perform the operation on each tagged data value
// and assign this to the corresponding element in result
272    const DataTagged::DataMapType& lookup=data.getTagLookup();    const DataTagged::DataMapType& lookup=data.getTagLookup();
273    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
274    DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();    DataTagged::DataMapType::const_iterator lookupEnd=lookup.end();
275      DataArrayView dataView=data.getPointDataView();
276      DataArrayView resultView=result.getPointDataView();
277      // perform the operation on each tagged data value
278      // and assign this to the corresponding element in result
279      //
280      // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
281      // which maintains state between calls which would be corrupted by parallel execution
282    for (i=lookup.begin();i!=lookupEnd;i++) {    for (i=lookup.begin();i!=lookupEnd;i++) {
283      result.getPointDataView().getData(i->second) =      resultView.getData(i->second) =
284        data.getPointDataView().dp_reductionOp(i->second,operation);        dataView.reductionOp(i->second,operation);
285    }    }
286    //    // perform the operation on the default data value
// finally perform the operation on the default data value
287    // and assign this to the default element in result    // and assign this to the default element in result
288    result.getPointDataView().getData(0) =    resultView.getData(0) =
289      data.getDefaultValue().dp_reductionOp(operation);      data.getDefaultValue().reductionOp(operation);
290  }  }
291
292  template <class UnaryFunction>  template <class UnaryFunction>
# Line 274  dp_algorithm(DataConstant& data, Line 296  dp_algorithm(DataConstant& data,
296               DataConstant& result,               DataConstant& result,
297               UnaryFunction operation)               UnaryFunction operation)
298  {  {
299    //    // perform the operation on the data value
300    // perform the operation on the default data value    // and assign this to the element in result
// and assign this to the default element in result
301    result.getPointDataView().getData(0) =    result.getPointDataView().getData(0) =
302      data.getPointDataView().dp_reductionOp(operation);      data.getPointDataView().reductionOp(operation);
303  }  }
304
305  } // end of namespace  } // end of namespace
306
307  #endif  #endif

Legend:
 Removed from v.108 changed lines Added in v.122