/[escript]/trunk/escript/src/DataLazy.cpp
ViewVC logotype

Annotation of /trunk/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1885 - (hide annotations)
Tue Oct 14 05:51:59 2008 UTC (10 years, 11 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 9122 byte(s)
memory management for basic ops done.

1 jfenwick 1865
2     /*******************************************************
3     *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14    
15     #include "DataLazy.h"
16     #ifdef USE_NETCDF
17     #include <netcdfcpp.h>
18     #endif
19     #ifdef PASO_MPI
20     #include <mpi.h>
21     #endif
22 jfenwick 1885 #ifdef _OPENMP
23     #include <omp.h>
24     #endif
25 jfenwick 1865 #include "FunctionSpace.h"
26     #include "DataTypes.h"
27 jfenwick 1868 #include "Data.h"
28 jfenwick 1865
29     using namespace std;
30 jfenwick 1868 using namespace boost;
31 jfenwick 1865
32     namespace escript
33     {
34    
35     const std::string&
36     opToString(ES_optype op);
37    
38     namespace
39     {
40    
41     // return the FunctionSpace of the result of "left op right"
42     FunctionSpace
43     resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
44     {
45     // perhaps this should call interpolate and throw or something?
46     // maybe we need an interpolate node -
47     // that way, if interpolate is required in any other op we can just throw a
48     // programming error exception.
49 jfenwick 1879
50    
51     if (left->getFunctionSpace()!=right->getFunctionSpace())
52     {
53     throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
54     }
55 jfenwick 1865 return left->getFunctionSpace();
56     }
57    
58     // return the shape of the result of "left op right"
59     DataTypes::ShapeType
60     resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
61     {
62 jfenwick 1879 if (left->getShape()!=right->getShape())
63     {
64     throw DataException("Shapes not the same - shapes must match for lazy data.");
65     }
66     return left->getShape();
67 jfenwick 1865 }
68    
69     size_t
70     resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
71     {
72     switch(op)
73     {
74 jfenwick 1879 // case IDENTITY: return left->getLength();
75 jfenwick 1868 case ADD: // the length is preserved in these ops
76     case SUB:
77     case MUL:
78     case DIV: return left->getLength();
79 jfenwick 1865 default:
80     throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
81    
82     }
83     }
84    
85 jfenwick 1879 int
86     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
87     {
88     switch(op)
89     {
90     case IDENTITY: return 0;
91     case ADD: // the length is preserved in these ops
92     case SUB:
93     case MUL:
94     case DIV: return max(left->getBuffsRequired(),right->getBuffsRequired());
95     default:
96     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
97     }
98     }
99    
100 jfenwick 1868 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/"};
101     int ES_opcount=5;
102 jfenwick 1865
103 jfenwick 1879 // void performOp(ValueType& v, int startOffset, ES_optype op, int m_samplesize)
104     // {
105     // switch(op)
106     // {
107     // case ADD: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
108     // startOffset+m_samplesize,plus<double>());
109     // break;
110     // case SUB: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
111     // startOffset+m_samplesize,minus<double>());
112     // break;
113     // case MUL: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
114     // startOffset+m_samplesize,multiplies<double>());
115     // break;
116     // case DIV: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
117     // startOffset+m_samplesize,divides<double>());
118     // break;
119     // default:
120     // throw DataException("Programmer Error - attempt to performOp() for operator "+opToString(op)+".");
121     // }
122     //
123     // }
124    
125 jfenwick 1865 } // end anonymous namespace
126    
127    
128     const std::string&
129     opToString(ES_optype op)
130     {
131     if (op<0 || op>=ES_opcount)
132     {
133     op=UNKNOWNOP;
134     }
135     return ES_opstrings[op];
136     }
137    
138    
139     DataLazy::DataLazy(DataAbstract_ptr p)
140     : parent(p->getFunctionSpace(),p->getShape()),
141     m_op(IDENTITY)
142     {
143 jfenwick 1879 if (p->isLazy())
144     {
145     // TODO: fix this. We could make the new node a copy of p?
146     // I don't want identity of Lazy.
147     // Question: Why would that be so bad?
148     // Answer: We assume that the child of ID is something we can call getVector on
149     throw DataException("Programmer error - attempt to create identity from a DataLazy.");
150     }
151     else
152     {
153     m_id=dynamic_pointer_cast<DataReady>(p);
154     }
155     m_length=p->getLength();
156     m_buffsRequired=0;
157     m_samplesize=getNumDPPSample()*getNoValues();
158     cout << "(1)Lazy created with " << m_samplesize << endl;
159 jfenwick 1865 }
160    
161 jfenwick 1879 DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
162 jfenwick 1865 : parent(resultFS(left,right,op), resultShape(left,right,op)),
163     m_left(left),
164     m_right(right),
165     m_op(op)
166     {
167 jfenwick 1879 m_length=resultLength(m_left,m_right,m_op);
168     m_samplesize=getNumDPPSample()*getNoValues();
169     m_buffsRequired=calcBuffs(m_left, m_right, m_op);
170     cout << "(2)Lazy created with " << m_samplesize << endl;
171 jfenwick 1865 }
172    
173 jfenwick 1879 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
174     : parent(resultFS(left,right,op), resultShape(left,right,op)),
175     m_op(op)
176     {
177     if (left->isLazy())
178     {
179     m_left=dynamic_pointer_cast<DataLazy>(left);
180     }
181     else
182     {
183     m_left=DataLazy_ptr(new DataLazy(left));
184     }
185     if (right->isLazy())
186     {
187     m_right=dynamic_pointer_cast<DataLazy>(right);
188     }
189     else
190     {
191     m_right=DataLazy_ptr(new DataLazy(right));
192     }
193    
194     m_length=resultLength(m_left,m_right,m_op);
195     m_samplesize=getNumDPPSample()*getNoValues();
196     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
197     cout << "(3)Lazy created with " << m_samplesize << endl;
198     }
199    
200    
201 jfenwick 1865 DataLazy::~DataLazy()
202     {
203     }
204    
205 jfenwick 1879
206     int
207     DataLazy::getBuffsRequired() const
208 jfenwick 1865 {
209 jfenwick 1879 return m_buffsRequired;
210     }
211    
212 jfenwick 1884
213     // the vector and the offset are a place where the method could write its data if it wishes
214     // it is not obligated to do so. For example, if it has its own storage already, it can use that.
215 jfenwick 1885 // Hence the return value to indicate where the data is actually stored.
216     // Regardless, the storage should be assumed to be used, even if it isn't.
217 jfenwick 1884 const double*
218 jfenwick 1879 DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset ) const
219     {
220 jfenwick 1885 if (m_op==IDENTITY)
221 jfenwick 1865 {
222 jfenwick 1879 const ValueType& vec=m_id->getVector();
223 jfenwick 1884 return &(vec[m_id->getPointOffset(sampleNo, 0)]);
224 jfenwick 1865 }
225 jfenwick 1879 size_t rightoffset=offset+m_samplesize;
226 jfenwick 1884 const double* left=m_left->resolveSample(v,sampleNo,offset);
227     const double* right=m_right->resolveSample(v,sampleNo,rightoffset);
228     double* result=&(v[offset]);
229 jfenwick 1865 {
230 jfenwick 1879 switch(m_op)
231     {
232     case ADD: // since these are pointwise ops, pretend each sample is one point
233 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
234 jfenwick 1879 break;
235     case SUB:
236 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
237 jfenwick 1879 break;
238     case MUL:
239 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
240 jfenwick 1879 break;
241     case DIV:
242 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
243 jfenwick 1879 break;
244 jfenwick 1868 default:
245     throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
246 jfenwick 1879 }
247 jfenwick 1868 }
248 jfenwick 1884 return result;
249 jfenwick 1865 }
250    
251 jfenwick 1879 DataReady_ptr
252     DataLazy::resolve()
253     {
254     // This is broken! We need to have a buffer per thread!
255     // so the allocation of v needs to move inside the loop somehow
256    
257     cout << "Sample size=" << m_samplesize << endl;
258     cout << "Buffers=" << m_buffsRequired << endl;
259    
260 jfenwick 1885 size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);
261     int numthreads=1;
262     #ifdef _OPENMP
263     numthreads=omp_get_max_threads();
264     int threadnum=0;
265     #endif
266     ValueType v(numthreads*threadbuffersize); // the +1 comes from the fact that I want to have a safe
267 jfenwick 1884 // space for the RHS of ops to write to even if they don't
268     // need it.
269 jfenwick 1879 cout << "Buffer created with size=" << v.size() << endl;
270     ValueType dummy(getNoValues());
271     DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);
272     ValueType& resvec=result->getVector();
273     DataReady_ptr resptr=DataReady_ptr(result);
274     int sample;
275 jfenwick 1885 int resoffset;
276     int totalsamples=getNumSamples();
277     #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)
278     for (sample=0;sample<totalsamples;++sample)
279 jfenwick 1879 {
280 jfenwick 1885 #ifdef _OPENMP
281     const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());
282     #else
283     const double* res=resolveSample(v,sample,0); // this would normally be v, but not if its a single IDENTITY op.
284     #endif
285     resoffset=result->getPointOffset(sample,0);
286     for (int i=0;i<m_samplesize;++i,++resoffset) // copy values into the output vector
287 jfenwick 1879 {
288 jfenwick 1885 resvec[resoffset]=res[i];
289 jfenwick 1879 }
290     }
291     return resptr;
292     }
293    
294 jfenwick 1865 std::string
295     DataLazy::toString() const
296     {
297     return "Lazy evaluation object. No details available.";
298     }
299    
300     // Note that in this case, deepCopy does not make copies of the leaves.
301     // Hopefully copy on write (or whatever we end up using) will take care of this.
302     DataAbstract*
303     DataLazy::deepCopy()
304     {
305     if (m_op==IDENTITY)
306     {
307     return new DataLazy(m_left); // we don't need to copy the child here
308     }
309     return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
310     }
311    
312    
313     DataTypes::ValueType::size_type
314     DataLazy::getLength() const
315     {
316 jfenwick 1879 return m_length;
317 jfenwick 1865 }
318    
319    
320     DataAbstract*
321     DataLazy::getSlice(const DataTypes::RegionType& region) const
322     {
323 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
324 jfenwick 1865 }
325    
326     DataTypes::ValueType::size_type
327     DataLazy::getPointOffset(int sampleNo,
328     int dataPointNo) const
329     {
330     throw DataException("getPointOffset - not implemented for Lazy objects - yet.");
331     }
332    
333     } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26