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

revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC revision 122 by jgs, Thu Jun 9 05:38:05 2005 UTC
# Line 1  Line 1
1    // \$Id\$
2  /*  /*
3   ******************************************************************************   ******************************************************************************
4   *                                                                            *   *                                                                            *
# Line 12  Line 13
13   ******************************************************************************   ******************************************************************************
14  */  */
15
16  #if !defined  escript_DataAlgorithm_20040714_H  #if !defined escript_DataAlgorithm_20040714_H
17  #define escript_DataAlgorithm_20040714_H  #define escript_DataAlgorithm_20040714_H
18
19  #include "escript/Data/DataExpanded.h"  #include "escript/Data/DataExpanded.h"
# Line 24  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 43  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 58  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 70  struct AbsMax : public std::binary_funct Line 107  struct AbsMax : public std::binary_funct
107      return std::max(fabs(x),fabs(y));      return std::max(fabs(x),fabs(y));
108    }    }
109  };  };
110
111  /**  /**
112     \brief     \brief
113     Adapt algorithms so they may be used by Data.     Return the absolute minimum value of the two given values.
114    */
115    struct AbsMin : public std::binary_function<double,double,double>
116    {
117      inline double operator()(double x, double y) const
118      {
119        return std::min(fabs(x),fabs(y));
120      }
121    };
122
123     Description:  /**
124     Adapt algorithms so they may be used by Data. The functor     \brief
125     maintains state, ie the currentValue retuned by the operation.     Return the length between the two given values.
126  */  */
127  template <class BinaryFunction>  struct Length : public std::binary_function<double,double,double>
129   public:    inline double operator()(double x, double y) const
131        m_currentValue(initialValue)      return std::sqrt(std::pow(x,2)+std::pow(y,2));
132      {}    }
inline void operator()(double value)
{
m_currentValue=operation(m_currentValue,value);
return;
}
inline double getResult() const
{
return m_currentValue;
}
private:
//
// the current maximum value
double m_currentValue;
//
// The operation to perform
BinaryFunction operation;
133  };  };
134
135  /**  /**
136     \brief     \brief
137     Perform the given operation upon all DataElements and return a single     Return the trace of the two given values.
138     result.  */
139    struct Trace : public std::binary_function<double,double,double>
140    {
141      inline double operator()(double x, double y) const
142      {
143        return x+y;
144      }
145    };
146
147     Description:  /**
148     Perform the given operation upon all DataElements and return a single     \brief
149     result.     Perform the given operation upon all values in all data-points in the
150       given Data object and return the final result.
151
152       Calls DataArrayView::reductionOp
153  */  */
154  template <class UnaryFunction>  template <class UnaryFunction>
155  inline double algorithm(DataExpanded& data, UnaryFunction operation)  inline
156    double
157    algorithm(DataExpanded& data,
158              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;    int resultVectorLength=numDPPSample*numSamples;
164    #pragma omp parallel private(resultLocal)    std::vector<double> resultVector(resultVectorLength);
165    {    DataArrayView dataView=data.getPointDataView();
166      #pragma omp parallel 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().algorithm(data.getPointOffset(i,j),    //
170                                operation);    // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
171      #pragma omp critical (algorithm)    // which maintains state between calls which would be corrupted by parallel execution
172      operation(resultLocal);    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
189  template <class UnaryFunction>  template <class UnaryFunction>
190  inline double algorithm(DataTagged& data, UnaryFunction operation)  inline
191    double
192    algorithm(DataTagged& data,
193              UnaryFunction operation)
194  {  {
//
// perform the operation on each tagged value including the default
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    for (i=lookup.begin();i!=lookupEnd;++i) {    std::vector<double> resultVector;
200      operation(dataView.algorithm(i->second,operation));    int resultVectorLength;
201      // perform the operation on each tagged value
202      for (i=lookup.begin();i!=lookupEnd;i++) {
203        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().algorithm(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
220  template <class UnaryFunction>  template <class UnaryFunction>
221  inline double algorithm(DataConstant& data, UnaryFunction operation)  inline
222    double
223    algorithm(DataConstant& data,
224              UnaryFunction operation)
225  {  {
226    return data.getPointDataView().algorithm(operation);    return data.getPointDataView().reductionOp(operation);
227  }  }
228
229    /**
230       \brief
231       Perform the given data-point reduction operation on all data-points
232       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
235       of data points, but where data has data points of rank n, result must have
236       data points of rank 0.
237
238       Calls DataArrayView::reductionOp
239    */
240    template <class UnaryFunction>
241    inline
242    void
243    dp_algorithm(DataExpanded& data,
244                 DataExpanded& result,
245                 UnaryFunction operation)
246    {
247      int i,j;
248      int numSamples=data.getNumSamples();
249      int numDPPSample=data.getNumDPPSample();
250      DataArrayView dataView=data.getPointDataView();
251      DataArrayView resultView=result.getPointDataView();
252      // perform the operation on each data-point and assign
253      // this to the corresponding element in result
254      //
255      // this loop cannot be prallelised as "operation" is an instance of DataAlgorithmAdapter
256      // 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    }
264
265    template <class UnaryFunction>
266    inline
267    void
268    dp_algorithm(DataTagged& data,
269                 DataTagged& result,
270                 UnaryFunction operation)
271    {
272      const DataTagged::DataMapType& lookup=data.getTagLookup();
273      DataTagged::DataMapType::const_iterator i;
274      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++) {
283        resultView.getData(i->second) =
284          dataView.reductionOp(i->second,operation);
285      }
286      // perform the operation on the default data value
287      // and assign this to the default element in result
288      resultView.getData(0) =
289        data.getDefaultValue().reductionOp(operation);
290    }
291
292    template <class UnaryFunction>
293    inline
294    void
295    dp_algorithm(DataConstant& data,
296                 DataConstant& result,
297                 UnaryFunction operation)
298    {
299      // perform the operation on the data value
300      // and assign this to the element in result
301      result.getPointDataView().getData(0) =
302        data.getPointDataView().reductionOp(operation);
303    }
304
305  } // end of namespace  } // end of namespace
306
307  #endif  #endif

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