/[escript]/branches/schroedinger/escript/src/DataLazy.cpp
ViewVC logotype

Annotation of /branches/schroedinger/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1884 - (hide annotations)
Tue Oct 14 04:54:59 2008 UTC (10 years, 5 months ago) by jfenwick
File size: 9398 byte(s)
Branch commit
Not crashing.

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     #include "FunctionSpace.h"
23     #include "DataTypes.h"
24 jfenwick 1868 #include "Data.h"
25 jfenwick 1865
26     using namespace std;
27 jfenwick 1868 using namespace boost;
28 jfenwick 1865
29     namespace escript
30     {
31    
32     const std::string&
33     opToString(ES_optype op);
34    
35     namespace
36     {
37    
38     // return the FunctionSpace of the result of "left op right"
39     FunctionSpace
40     resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
41     {
42     // perhaps this should call interpolate and throw or something?
43     // maybe we need an interpolate node -
44     // that way, if interpolate is required in any other op we can just throw a
45     // programming error exception.
46 jfenwick 1879
47    
48     if (left->getFunctionSpace()!=right->getFunctionSpace())
49     {
50     throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
51     }
52 jfenwick 1865 return left->getFunctionSpace();
53     }
54    
55     // return the shape of the result of "left op right"
56     DataTypes::ShapeType
57     resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
58     {
59 jfenwick 1879 if (left->getShape()!=right->getShape())
60     {
61     throw DataException("Shapes not the same - shapes must match for lazy data.");
62     }
63     return left->getShape();
64 jfenwick 1865 }
65    
66     size_t
67     resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
68     {
69     switch(op)
70     {
71 jfenwick 1879 // case IDENTITY: return left->getLength();
72 jfenwick 1868 case ADD: // the length is preserved in these ops
73     case SUB:
74     case MUL:
75     case DIV: return left->getLength();
76 jfenwick 1865 default:
77     throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
78    
79     }
80     }
81    
82 jfenwick 1879 int
83     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
84     {
85     switch(op)
86     {
87     case IDENTITY: return 0;
88     case ADD: // the length is preserved in these ops
89     case SUB:
90     case MUL:
91     case DIV: return max(left->getBuffsRequired(),right->getBuffsRequired());
92     default:
93     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
94     }
95     }
96    
97 jfenwick 1868 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/"};
98     int ES_opcount=5;
99 jfenwick 1865
100 jfenwick 1879 // void performOp(ValueType& v, int startOffset, ES_optype op, int m_samplesize)
101     // {
102     // switch(op)
103     // {
104     // case ADD: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
105     // startOffset+m_samplesize,plus<double>());
106     // break;
107     // case SUB: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
108     // startOffset+m_samplesize,minus<double>());
109     // break;
110     // case MUL: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
111     // startOffset+m_samplesize,multiplies<double>());
112     // break;
113     // case DIV: DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),
114     // startOffset+m_samplesize,divides<double>());
115     // break;
116     // default:
117     // throw DataException("Programmer Error - attempt to performOp() for operator "+opToString(op)+".");
118     // }
119     //
120     // }
121    
122 jfenwick 1865 } // end anonymous namespace
123    
124    
125     const std::string&
126     opToString(ES_optype op)
127     {
128     if (op<0 || op>=ES_opcount)
129     {
130     op=UNKNOWNOP;
131     }
132     return ES_opstrings[op];
133     }
134    
135    
136     DataLazy::DataLazy(DataAbstract_ptr p)
137     : parent(p->getFunctionSpace(),p->getShape()),
138     m_op(IDENTITY)
139     {
140 jfenwick 1879 if (p->isLazy())
141     {
142     // TODO: fix this. We could make the new node a copy of p?
143     // I don't want identity of Lazy.
144     // Question: Why would that be so bad?
145     // Answer: We assume that the child of ID is something we can call getVector on
146     throw DataException("Programmer error - attempt to create identity from a DataLazy.");
147     }
148     else
149     {
150     m_id=dynamic_pointer_cast<DataReady>(p);
151     }
152     m_length=p->getLength();
153     m_buffsRequired=0;
154     m_samplesize=getNumDPPSample()*getNoValues();
155     cout << "(1)Lazy created with " << m_samplesize << endl;
156 jfenwick 1865 }
157    
158 jfenwick 1879 DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
159 jfenwick 1865 : parent(resultFS(left,right,op), resultShape(left,right,op)),
160     m_left(left),
161     m_right(right),
162     m_op(op)
163     {
164 jfenwick 1879 m_length=resultLength(m_left,m_right,m_op);
165     m_samplesize=getNumDPPSample()*getNoValues();
166     m_buffsRequired=calcBuffs(m_left, m_right, m_op);
167     cout << "(2)Lazy created with " << m_samplesize << endl;
168 jfenwick 1865 }
169    
170 jfenwick 1879 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
171     : parent(resultFS(left,right,op), resultShape(left,right,op)),
172     m_op(op)
173     {
174     if (left->isLazy())
175     {
176     m_left=dynamic_pointer_cast<DataLazy>(left);
177     }
178     else
179     {
180     m_left=DataLazy_ptr(new DataLazy(left));
181     }
182     if (right->isLazy())
183     {
184     m_right=dynamic_pointer_cast<DataLazy>(right);
185     }
186     else
187     {
188     m_right=DataLazy_ptr(new DataLazy(right));
189     }
190    
191     m_length=resultLength(m_left,m_right,m_op);
192     m_samplesize=getNumDPPSample()*getNoValues();
193     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
194     cout << "(3)Lazy created with " << m_samplesize << endl;
195     }
196    
197    
198 jfenwick 1865 DataLazy::~DataLazy()
199     {
200     }
201    
202 jfenwick 1879
203     int
204     DataLazy::getBuffsRequired() const
205 jfenwick 1865 {
206 jfenwick 1879 return m_buffsRequired;
207     }
208    
209 jfenwick 1884
210     // the vector and the offset are a place where the method could write its data if it wishes
211     // it is not obligated to do so. For example, if it has its own storage already, it can use that.
212     // hence the return value to indicate where the data is actually stored.
213     // regardless, the storage should be assumed to be used, even if it isn't.
214     const double*
215 jfenwick 1879 DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset ) const
216     {
217     if (m_op==IDENTITY) // copy the contents into the vector
218 jfenwick 1865 {
219 jfenwick 1879 cout << "Begin ID" << endl;
220     cout << "dpps=" << getNumDPPSample() << " novals=" << getNoValues() << endl;
221     const ValueType& vec=m_id->getVector();
222 jfenwick 1884 // size_t srcOffset=m_id->getPointOffset(sampleNo, 0);
223     // cout << "v.size()=" << v.size() << " vec=" << vec.size() << endl;
224     // for (size_t i=0;i<m_samplesize;++i,++srcOffset,++offset)
225     // {
226     // cout << "Trying offset=" << offset << " srcOffset=" << srcOffset << endl;
227     // v[offset]=vec[srcOffset];
228     // }
229     cout << "End ID - returning offset " << m_id->getPointOffset(sampleNo, 0) << " of vector@" << &vec<<endl;
230     return &(vec[m_id->getPointOffset(sampleNo, 0)]);
231     // return;
232 jfenwick 1865 }
233 jfenwick 1884 cout << "Begin op";
234 jfenwick 1879 size_t rightoffset=offset+m_samplesize;
235 jfenwick 1884 const double* left=m_left->resolveSample(v,sampleNo,offset);
236     const double* right=m_right->resolveSample(v,sampleNo,rightoffset);
237     double* result=&(v[offset]);
238     cout << "left=" << left << " right=" << right << " result=" << result << endl;
239 jfenwick 1879 // for (int i=0;i<getNumDPPSample();++i)
240 jfenwick 1865 {
241 jfenwick 1879 switch(m_op)
242     {
243     case ADD: // since these are pointwise ops, pretend each sample is one point
244 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
245 jfenwick 1879 break;
246     case SUB:
247 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
248 jfenwick 1879 break;
249     case MUL:
250 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
251 jfenwick 1879 break;
252     case DIV:
253 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
254 jfenwick 1879 break;
255 jfenwick 1868 default:
256     throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
257 jfenwick 1879 }
258 jfenwick 1868 }
259 jfenwick 1884 return result;
260 jfenwick 1865 }
261    
262 jfenwick 1879 DataReady_ptr
263     DataLazy::resolve()
264     {
265     // This is broken! We need to have a buffer per thread!
266     // so the allocation of v needs to move inside the loop somehow
267    
268     cout << "Sample size=" << m_samplesize << endl;
269     cout << "Buffers=" << m_buffsRequired << endl;
270    
271 jfenwick 1884 ValueType v(m_samplesize*(max(1,m_buffsRequired)+1)); // the +1 comes from the fact that I want to have a safe
272     // space for the RHS of ops to write to even if they don't
273     // need it.
274 jfenwick 1879 cout << "Buffer created with size=" << v.size() << endl;
275     ValueType dummy(getNoValues());
276     DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);
277     ValueType& resvec=result->getVector();
278     DataReady_ptr resptr=DataReady_ptr(result);
279     int sample;
280     #pragma omp parallel for private(sample) schedule(static)
281     for (sample=0;sample<getNumSamples();++sample)
282     {
283 jfenwick 1884 cout << "Processing sample#" << sample << endl;
284 jfenwick 1879 resolveSample(v,sample,0);
285 jfenwick 1884 cout << "Copying#" << sample << endl;
286 jfenwick 1879 for (int i=0;i<m_samplesize;++i) // copy values into the output vector
287     {
288     resvec[i]=v[i];
289     }
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