/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Contents of /branches/clazy/escriptcore/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

1
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 #ifdef _OPENMP
23 #include <omp.h>
24 #endif
25 #include "FunctionSpace.h"
26 #include "DataTypes.h"
27 #include "Data.h"
28
29 using namespace std;
30 using namespace boost;
31
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
50
51 if (left->getFunctionSpace()!=right->getFunctionSpace())
52 {
53 throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
54 }
55 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 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 }
68
69 size_t
70 resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
71 {
72 switch(op)
73 {
74 // case IDENTITY: return left->getLength();
75 case ADD: // the length is preserved in these ops
76 case SUB:
77 case MUL:
78 case DIV: return left->getLength();
79 default:
80 throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
81
82 }
83 }
84
85 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 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/"};
101 int ES_opcount=5;
102
103 // 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 } // 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 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 }
160
161 DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
162 : parent(resultFS(left,right,op), resultShape(left,right,op)),
163 m_left(left),
164 m_right(right),
165 m_op(op)
166 {
167 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 }
172
173 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 DataLazy::~DataLazy()
202 {
203 }
204
205
206 int
207 DataLazy::getBuffsRequired() const
208 {
209 return m_buffsRequired;
210 }
211
212
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 // 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 const double*
218 DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset ) const
219 {
220 if (m_op==IDENTITY)
221 {
222 const ValueType& vec=m_id->getVector();
223 return &(vec[m_id->getPointOffset(sampleNo, 0)]);
224 }
225 size_t rightoffset=offset+m_samplesize;
226 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 {
230 switch(m_op)
231 {
232 case ADD: // since these are pointwise ops, pretend each sample is one point
233 tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
234 break;
235 case SUB:
236 tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
237 break;
238 case MUL:
239 tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
240 break;
241 case DIV:
242 tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
243 break;
244 default:
245 throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
246 }
247 }
248 return result;
249 }
250
251 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 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 // space for the RHS of ops to write to even if they don't
268 // need it.
269 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 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 {
280 #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 {
288 resvec[resoffset]=res[i];
289 }
290 }
291 return resptr;
292 }
293
294 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 return m_length;
317 }
318
319
320 DataAbstract*
321 DataLazy::getSlice(const DataTypes::RegionType& region) const
322 {
323 throw DataException("getSlice - not implemented for Lazy objects.");
324 }
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