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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2199 - (show annotations)
Thu Jan 8 06:10:52 2009 UTC (10 years, 10 months ago) by jfenwick
File size: 55552 byte(s)
Misc fixes:
Added some svn:ignore properties for output files that were cluttering things up.

Lazy fixes:
Fixed shape calculations for TRACE and TRANSPOSE for rank>2.
Adjusted unit test accordingly.

As a Temporary change to DataC.cpp to test for lazy data in DataC's expanded check.
This is wrong but would only affect people using lazy data.
The proper fix will come when the numarray removal code moves over from the branch.

Made tensor product AUTOLAZY capable.
Fixed some bugs resolving tensor products (incorrect offsets in buffers).
Macro'd some stray couts.

- It appears that AUTOLAZY now passes all unit tests.
- It will not be _really_ safe for general use until I can add COW. 
- (Everything's better with COW)
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 #include "UnaryFuncs.h" // for escript::fsign
29 #include "Utils.h"
30
31 #include <iomanip> // for some fancy formatting in debug
32
33 // #define LAZYDEBUG(X) if (privdebug){X;}
34 #define LAZYDEBUG(X)
35 namespace
36 {
37 bool privdebug=false;
38
39 #define ENABLEDEBUG privdebug=true;
40 #define DISABLEDEBUG privdebug=false;
41 }
42
43 #define SIZELIMIT
44 // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
45
46
47 /*
48 How does DataLazy work?
49 ~~~~~~~~~~~~~~~~~~~~~~~
50
51 Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
52 denoted left and right.
53
54 A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
55 This means that all "internal" nodes in the structure are instances of DataLazy.
56
57 Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
58 Note that IDENITY is not considered a unary operation.
59
60 I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
61 It must however form a DAG (directed acyclic graph).
62 I will refer to individual DataLazy objects with the structure as nodes.
63
64 Each node also stores:
65 - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
66 evaluated.
67 - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
68 evaluate the expression.
69 - m_samplesize ~ the number of doubles stored in a sample.
70
71 When a new node is created, the above values are computed based on the values in the child nodes.
72 Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
73
74 The resolve method, which produces a DataReady from a DataLazy, does the following:
75 1) Create a DataReady to hold the new result.
76 2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
77 3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
78
79 (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
80
81 resolveSample returns a Vector* and an offset within that vector where the result is stored.
82 Normally, this would be v, but for identity nodes their internal vector is returned instead.
83
84 The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
85
86 For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
87 The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
88
89 To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
90 1) Add to the ES_optype.
91 2) determine what opgroup your operation belongs to (X)
92 3) add a string for the op to the end of ES_opstrings
93 4) increase ES_opcount
94 5) add an entry (X) to opgroups
95 6) add an entry to the switch in collapseToReady
96 7) add an entry to resolveX
97 */
98
99
100 using namespace std;
101 using namespace boost;
102
103 namespace escript
104 {
105
106 namespace
107 {
108
109 enum ES_opgroup
110 {
111 G_UNKNOWN,
112 G_IDENTITY,
113 G_BINARY, // pointwise operations with two arguments
114 G_UNARY, // pointwise operations with one argument
115 G_UNARY_P, // pointwise operations with one argument, requiring a parameter
116 G_NP1OUT, // non-pointwise op with one output
117 G_NP1OUT_P, // non-pointwise op with one output requiring a parameter
118 G_TENSORPROD // general tensor product
119 };
120
121
122
123
124 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
125 "sin","cos","tan",
126 "asin","acos","atan","sinh","cosh","tanh","erf",
127 "asinh","acosh","atanh",
128 "log10","log","sign","abs","neg","pos","exp","sqrt",
129 "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
130 "symmetric","nonsymmetric",
131 "prod",
132 "transpose", "trace"};
133 int ES_opcount=40;
134 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
135 G_UNARY,G_UNARY,G_UNARY, //10
136 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
137 G_UNARY,G_UNARY,G_UNARY, // 20
138 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
139 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P, // 35
140 G_NP1OUT,G_NP1OUT,
141 G_TENSORPROD,
142 G_NP1OUT_P, G_NP1OUT_P};
143 inline
144 ES_opgroup
145 getOpgroup(ES_optype op)
146 {
147 return opgroups[op];
148 }
149
150 // return the FunctionSpace of the result of "left op right"
151 FunctionSpace
152 resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
153 {
154 // perhaps this should call interpolate and throw or something?
155 // maybe we need an interpolate node -
156 // that way, if interpolate is required in any other op we can just throw a
157 // programming error exception.
158
159 FunctionSpace l=left->getFunctionSpace();
160 FunctionSpace r=right->getFunctionSpace();
161 if (l!=r)
162 {
163 if (r.probeInterpolation(l))
164 {
165 return l;
166 }
167 if (l.probeInterpolation(r))
168 {
169 return r;
170 }
171 throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
172 }
173 return l;
174 }
175
176 // return the shape of the result of "left op right"
177 // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
178 DataTypes::ShapeType
179 resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
180 {
181 if (left->getShape()!=right->getShape())
182 {
183 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
184 {
185 throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
186 }
187 if (left->getRank()==0) // we need to allow scalar * anything
188 {
189 return right->getShape();
190 }
191 if (right->getRank()==0)
192 {
193 return left->getShape();
194 }
195 throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
196 }
197 return left->getShape();
198 }
199
200 // return the shape for "op left"
201
202 DataTypes::ShapeType
203 resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
204 {
205 switch(op)
206 {
207 case TRANS:
208 { // for the scoping of variables
209 const DataTypes::ShapeType& s=left->getShape();
210 DataTypes::ShapeType sh;
211 int rank=left->getRank();
212 if (axis_offset<0 || axis_offset>rank)
213 {
214 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
215 }
216 for (int i=0; i<rank; i++)
217 {
218 int index = (axis_offset+i)%rank;
219 sh.push_back(s[index]); // Append to new shape
220 }
221 return sh;
222 }
223 break;
224 case TRACE:
225 {
226 int rank=left->getRank();
227 if (rank<2)
228 {
229 throw DataException("Trace can only be computed for objects with rank 2 or greater.");
230 }
231 if ((axis_offset>rank-2) || (axis_offset<0))
232 {
233 throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
234 }
235 if (rank==2)
236 {
237 return DataTypes::scalarShape;
238 }
239 else if (rank==3)
240 {
241 DataTypes::ShapeType sh;
242 if (axis_offset==0)
243 {
244 sh.push_back(left->getShape()[2]);
245 }
246 else // offset==1
247 {
248 sh.push_back(left->getShape()[0]);
249 }
250 return sh;
251 }
252 else if (rank==4)
253 {
254 DataTypes::ShapeType sh;
255 const DataTypes::ShapeType& s=left->getShape();
256 if (axis_offset==0)
257 {
258 sh.push_back(s[2]);
259 sh.push_back(s[3]);
260 }
261 else if (axis_offset==1)
262 {
263 sh.push_back(s[0]);
264 sh.push_back(s[3]);
265 }
266 else // offset==2
267 {
268 sh.push_back(s[0]);
269 sh.push_back(s[1]);
270 }
271 return sh;
272 }
273 else // unknown rank
274 {
275 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
276 }
277 }
278 break;
279 default:
280 throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
281 }
282 }
283
284 // determine the output shape for the general tensor product operation
285 // the additional parameters return information required later for the product
286 // the majority of this code is copy pasted from C_General_Tensor_Product
287 DataTypes::ShapeType
288 GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
289 {
290
291 // Get rank and shape of inputs
292 int rank0 = left->getRank();
293 int rank1 = right->getRank();
294 const DataTypes::ShapeType& shape0 = left->getShape();
295 const DataTypes::ShapeType& shape1 = right->getShape();
296
297 // Prepare for the loops of the product and verify compatibility of shapes
298 int start0=0, start1=0;
299 if (transpose == 0) {}
300 else if (transpose == 1) { start0 = axis_offset; }
301 else if (transpose == 2) { start1 = rank1-axis_offset; }
302 else { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
303
304 if (rank0<axis_offset)
305 {
306 throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
307 }
308
309 // Adjust the shapes for transpose
310 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
311 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
312 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
313 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
314
315 // Prepare for the loops of the product
316 SL=1, SM=1, SR=1;
317 for (int i=0; i<rank0-axis_offset; i++) {
318 SL *= tmpShape0[i];
319 }
320 for (int i=rank0-axis_offset; i<rank0; i++) {
321 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
322 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
323 }
324 SM *= tmpShape0[i];
325 }
326 for (int i=axis_offset; i<rank1; i++) {
327 SR *= tmpShape1[i];
328 }
329
330 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
331 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
332 { // block to limit the scope of out_index
333 int out_index=0;
334 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
335 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
336 }
337
338 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
339 {
340 ostringstream os;
341 os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
342 throw DataException(os.str());
343 }
344
345 return shape2;
346 }
347
348 // determine the number of samples requires to evaluate an expression combining left and right
349 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
350 // The same goes for G_TENSORPROD
351 // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefths.
352 // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
353 // multiple times
354 int
355 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
356 {
357 switch(getOpgroup(op))
358 {
359 case G_IDENTITY: return 1;
360 case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
361 case G_UNARY:
362 case G_UNARY_P: return max(left->getBuffsRequired(),1);
363 case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
364 case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
365 case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
366 default:
367 throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
368 }
369 }
370
371
372 } // end anonymous namespace
373
374
375
376 // Return a string representing the operation
377 const std::string&
378 opToString(ES_optype op)
379 {
380 if (op<0 || op>=ES_opcount)
381 {
382 op=UNKNOWNOP;
383 }
384 return ES_opstrings[op];
385 }
386
387
388 DataLazy::DataLazy(DataAbstract_ptr p)
389 : parent(p->getFunctionSpace(),p->getShape())
390 {
391 if (p->isLazy())
392 {
393 // I don't want identity of Lazy.
394 // Question: Why would that be so bad?
395 // Answer: We assume that the child of ID is something we can call getVector on
396 throw DataException("Programmer error - attempt to create identity from a DataLazy.");
397 }
398 else
399 {
400 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
401 makeIdentity(dr);
402 LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
403 }
404 LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
405 }
406
407
408
409
410 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
411 : parent(left->getFunctionSpace(),left->getShape()),
412 m_op(op),
413 m_axis_offset(0),
414 m_transpose(0),
415 m_SL(0), m_SM(0), m_SR(0)
416 {
417 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
418 {
419 throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
420 }
421
422 DataLazy_ptr lleft;
423 if (!left->isLazy())
424 {
425 lleft=DataLazy_ptr(new DataLazy(left));
426 }
427 else
428 {
429 lleft=dynamic_pointer_cast<DataLazy>(left);
430 }
431 m_readytype=lleft->m_readytype;
432 m_left=lleft;
433 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
434 m_samplesize=getNumDPPSample()*getNoValues();
435 m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
436 m_children=m_left->m_children+1;
437 m_height=m_left->m_height+1;
438 SIZELIMIT
439 }
440
441
442 // In this constructor we need to consider interpolation
443 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
444 : parent(resultFS(left,right,op), resultShape(left,right,op)),
445 m_op(op),
446 m_SL(0), m_SM(0), m_SR(0)
447 {
448 LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
449 if ((getOpgroup(op)!=G_BINARY))
450 {
451 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
452 }
453
454 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
455 {
456 FunctionSpace fs=getFunctionSpace();
457 Data ltemp(left);
458 Data tmp(ltemp,fs);
459 left=tmp.borrowDataPtr();
460 }
461 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
462 {
463 Data tmp(Data(right),getFunctionSpace());
464 right=tmp.borrowDataPtr();
465 LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
466 }
467 left->operandCheck(*right);
468
469 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
470 {
471 m_left=dynamic_pointer_cast<DataLazy>(left);
472 LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
473 }
474 else
475 {
476 m_left=DataLazy_ptr(new DataLazy(left));
477 LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
478 }
479 if (right->isLazy())
480 {
481 m_right=dynamic_pointer_cast<DataLazy>(right);
482 LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
483 }
484 else
485 {
486 m_right=DataLazy_ptr(new DataLazy(right));
487 LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
488 }
489 char lt=m_left->m_readytype;
490 char rt=m_right->m_readytype;
491 if (lt=='E' || rt=='E')
492 {
493 m_readytype='E';
494 }
495 else if (lt=='T' || rt=='T')
496 {
497 m_readytype='T';
498 }
499 else
500 {
501 m_readytype='C';
502 }
503 m_samplesize=getNumDPPSample()*getNoValues();
504 m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
505 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
506 m_children=m_left->m_children+m_right->m_children+2;
507 m_height=max(m_left->m_height,m_right->m_height)+1;
508 SIZELIMIT
509 LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
510 }
511
512 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
513 : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
514 m_op(op),
515 m_axis_offset(axis_offset),
516 m_transpose(transpose)
517 {
518 if ((getOpgroup(op)!=G_TENSORPROD))
519 {
520 throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
521 }
522 if ((transpose>2) || (transpose<0))
523 {
524 throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
525 }
526 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
527 {
528 FunctionSpace fs=getFunctionSpace();
529 Data ltemp(left);
530 Data tmp(ltemp,fs);
531 left=tmp.borrowDataPtr();
532 }
533 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
534 {
535 Data tmp(Data(right),getFunctionSpace());
536 right=tmp.borrowDataPtr();
537 }
538 // left->operandCheck(*right);
539
540 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
541 {
542 m_left=dynamic_pointer_cast<DataLazy>(left);
543 }
544 else
545 {
546 m_left=DataLazy_ptr(new DataLazy(left));
547 }
548 if (right->isLazy())
549 {
550 m_right=dynamic_pointer_cast<DataLazy>(right);
551 }
552 else
553 {
554 m_right=DataLazy_ptr(new DataLazy(right));
555 }
556 char lt=m_left->m_readytype;
557 char rt=m_right->m_readytype;
558 if (lt=='E' || rt=='E')
559 {
560 m_readytype='E';
561 }
562 else if (lt=='T' || rt=='T')
563 {
564 m_readytype='T';
565 }
566 else
567 {
568 m_readytype='C';
569 }
570 m_samplesize=getNumDPPSample()*getNoValues();
571 m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
572 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
573 m_children=m_left->m_children+m_right->m_children+2;
574 m_height=max(m_left->m_height,m_right->m_height)+1;
575 SIZELIMIT
576 LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
577 }
578
579
580 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
581 : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
582 m_op(op),
583 m_axis_offset(axis_offset),
584 m_transpose(0),
585 m_tol(0)
586 {
587 if ((getOpgroup(op)!=G_NP1OUT_P))
588 {
589 throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
590 }
591 DataLazy_ptr lleft;
592 if (!left->isLazy())
593 {
594 lleft=DataLazy_ptr(new DataLazy(left));
595 }
596 else
597 {
598 lleft=dynamic_pointer_cast<DataLazy>(left);
599 }
600 m_readytype=lleft->m_readytype;
601 m_left=lleft;
602 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
603 m_samplesize=getNumDPPSample()*getNoValues();
604 m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
605 m_children=m_left->m_children+1;
606 m_height=m_left->m_height+1;
607 SIZELIMIT
608 LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
609 }
610
611 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
612 : parent(left->getFunctionSpace(), left->getShape()),
613 m_op(op),
614 m_axis_offset(0),
615 m_transpose(0),
616 m_tol(tol)
617 {
618 if ((getOpgroup(op)!=G_UNARY_P))
619 {
620 throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
621 }
622 DataLazy_ptr lleft;
623 if (!left->isLazy())
624 {
625 lleft=DataLazy_ptr(new DataLazy(left));
626 }
627 else
628 {
629 lleft=dynamic_pointer_cast<DataLazy>(left);
630 }
631 m_readytype=lleft->m_readytype;
632 m_left=lleft;
633 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
634 m_samplesize=getNumDPPSample()*getNoValues();
635 m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
636 m_children=m_left->m_children+1;
637 m_height=m_left->m_height+1;
638 SIZELIMIT
639 LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
640 }
641
642 DataLazy::~DataLazy()
643 {
644 }
645
646
647 int
648 DataLazy::getBuffsRequired() const
649 {
650 return m_buffsRequired;
651 }
652
653
654 size_t
655 DataLazy::getMaxSampleSize() const
656 {
657 return m_maxsamplesize;
658 }
659
660 /*
661 \brief Evaluates the expression using methods on Data.
662 This does the work for the collapse method.
663 For reasons of efficiency do not call this method on DataExpanded nodes.
664 */
665 DataReady_ptr
666 DataLazy::collapseToReady()
667 {
668 if (m_readytype=='E')
669 { // this is more an efficiency concern than anything else
670 throw DataException("Programmer Error - do not use collapse on Expanded data.");
671 }
672 if (m_op==IDENTITY)
673 {
674 return m_id;
675 }
676 DataReady_ptr pleft=m_left->collapseToReady();
677 Data left(pleft);
678 Data right;
679 if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
680 {
681 right=Data(m_right->collapseToReady());
682 }
683 Data result;
684 switch(m_op)
685 {
686 case ADD:
687 result=left+right;
688 break;
689 case SUB:
690 result=left-right;
691 break;
692 case MUL:
693 result=left*right;
694 break;
695 case DIV:
696 result=left/right;
697 break;
698 case SIN:
699 result=left.sin();
700 break;
701 case COS:
702 result=left.cos();
703 break;
704 case TAN:
705 result=left.tan();
706 break;
707 case ASIN:
708 result=left.asin();
709 break;
710 case ACOS:
711 result=left.acos();
712 break;
713 case ATAN:
714 result=left.atan();
715 break;
716 case SINH:
717 result=left.sinh();
718 break;
719 case COSH:
720 result=left.cosh();
721 break;
722 case TANH:
723 result=left.tanh();
724 break;
725 case ERF:
726 result=left.erf();
727 break;
728 case ASINH:
729 result=left.asinh();
730 break;
731 case ACOSH:
732 result=left.acosh();
733 break;
734 case ATANH:
735 result=left.atanh();
736 break;
737 case LOG10:
738 result=left.log10();
739 break;
740 case LOG:
741 result=left.log();
742 break;
743 case SIGN:
744 result=left.sign();
745 break;
746 case ABS:
747 result=left.abs();
748 break;
749 case NEG:
750 result=left.neg();
751 break;
752 case POS:
753 // it doesn't mean anything for delayed.
754 // it will just trigger a deep copy of the lazy object
755 throw DataException("Programmer error - POS not supported for lazy data.");
756 break;
757 case EXP:
758 result=left.exp();
759 break;
760 case SQRT:
761 result=left.sqrt();
762 break;
763 case RECIP:
764 result=left.oneOver();
765 break;
766 case GZ:
767 result=left.wherePositive();
768 break;
769 case LZ:
770 result=left.whereNegative();
771 break;
772 case GEZ:
773 result=left.whereNonNegative();
774 break;
775 case LEZ:
776 result=left.whereNonPositive();
777 break;
778 case NEZ:
779 result=left.whereNonZero(m_tol);
780 break;
781 case EZ:
782 result=left.whereZero(m_tol);
783 break;
784 case SYM:
785 result=left.symmetric();
786 break;
787 case NSYM:
788 result=left.nonsymmetric();
789 break;
790 case PROD:
791 result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
792 break;
793 case TRANS:
794 result=left.transpose(m_axis_offset);
795 break;
796 case TRACE:
797 result=left.trace(m_axis_offset);
798 break;
799 default:
800 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
801 }
802 return result.borrowReadyPtr();
803 }
804
805 /*
806 \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
807 This method uses the original methods on the Data class to evaluate the expressions.
808 For this reason, it should not be used on DataExpanded instances. (To do so would defeat
809 the purpose of using DataLazy in the first place).
810 */
811 void
812 DataLazy::collapse()
813 {
814 if (m_op==IDENTITY)
815 {
816 return;
817 }
818 if (m_readytype=='E')
819 { // this is more an efficiency concern than anything else
820 throw DataException("Programmer Error - do not use collapse on Expanded data.");
821 }
822 m_id=collapseToReady();
823 m_op=IDENTITY;
824 }
825
826 /*
827 \brief Compute the value of the expression (unary operation) for the given sample.
828 \return Vector which stores the value of the subexpression for the given sample.
829 \param v A vector to store intermediate results.
830 \param offset Index in v to begin storing results.
831 \param sampleNo Sample number to evaluate.
832 \param roffset (output parameter) the offset in the return vector where the result begins.
833
834 The return value will be an existing vector so do not deallocate it.
835 If the result is stored in v it should be stored at the offset given.
836 Everything from offset to the end of v should be considered available for this method to use.
837 */
838 DataTypes::ValueType*
839 DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
840 {
841 // we assume that any collapsing has been done before we get here
842 // since we only have one argument we don't need to think about only
843 // processing single points.
844 if (m_readytype!='E')
845 {
846 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
847 }
848 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
849 const double* left=&((*vleft)[roffset]);
850 double* result=&(v[offset]);
851 roffset=offset;
852 switch (m_op)
853 {
854 case SIN:
855 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
856 break;
857 case COS:
858 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
859 break;
860 case TAN:
861 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
862 break;
863 case ASIN:
864 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
865 break;
866 case ACOS:
867 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
868 break;
869 case ATAN:
870 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
871 break;
872 case SINH:
873 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
874 break;
875 case COSH:
876 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
877 break;
878 case TANH:
879 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
880 break;
881 case ERF:
882 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
883 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
884 #else
885 tensor_unary_operation(m_samplesize, left, result, ::erf);
886 break;
887 #endif
888 case ASINH:
889 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
890 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
891 #else
892 tensor_unary_operation(m_samplesize, left, result, ::asinh);
893 #endif
894 break;
895 case ACOSH:
896 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
897 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
898 #else
899 tensor_unary_operation(m_samplesize, left, result, ::acosh);
900 #endif
901 break;
902 case ATANH:
903 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
904 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
905 #else
906 tensor_unary_operation(m_samplesize, left, result, ::atanh);
907 #endif
908 break;
909 case LOG10:
910 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
911 break;
912 case LOG:
913 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
914 break;
915 case SIGN:
916 tensor_unary_operation(m_samplesize, left, result, escript::fsign);
917 break;
918 case ABS:
919 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
920 break;
921 case NEG:
922 tensor_unary_operation(m_samplesize, left, result, negate<double>());
923 break;
924 case POS:
925 // it doesn't mean anything for delayed.
926 // it will just trigger a deep copy of the lazy object
927 throw DataException("Programmer error - POS not supported for lazy data.");
928 break;
929 case EXP:
930 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
931 break;
932 case SQRT:
933 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
934 break;
935 case RECIP:
936 tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
937 break;
938 case GZ:
939 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
940 break;
941 case LZ:
942 tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
943 break;
944 case GEZ:
945 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
946 break;
947 case LEZ:
948 tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
949 break;
950 // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
951 case NEZ:
952 tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
953 break;
954 case EZ:
955 tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
956 break;
957
958 default:
959 throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
960 }
961 return &v;
962 }
963
964
965
966
967
968
969 /*
970 \brief Compute the value of the expression (unary operation) for the given sample.
971 \return Vector which stores the value of the subexpression for the given sample.
972 \param v A vector to store intermediate results.
973 \param offset Index in v to begin storing results.
974 \param sampleNo Sample number to evaluate.
975 \param roffset (output parameter) the offset in the return vector where the result begins.
976
977 The return value will be an existing vector so do not deallocate it.
978 If the result is stored in v it should be stored at the offset given.
979 Everything from offset to the end of v should be considered available for this method to use.
980 */
981 DataTypes::ValueType*
982 DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
983 {
984 // we assume that any collapsing has been done before we get here
985 // since we only have one argument we don't need to think about only
986 // processing single points.
987 if (m_readytype!='E')
988 {
989 throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
990 }
991 // since we can't write the result over the input, we need a result offset further along
992 size_t subroffset=roffset+m_samplesize;
993 LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
994 const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
995 roffset=offset;
996 size_t loop=0;
997 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
998 size_t step=getNoValues();
999 switch (m_op)
1000 {
1001 case SYM:
1002 for (loop=0;loop<numsteps;++loop)
1003 {
1004 DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1005 subroffset+=step;
1006 offset+=step;
1007 }
1008 break;
1009 case NSYM:
1010 for (loop=0;loop<numsteps;++loop)
1011 {
1012 DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1013 subroffset+=step;
1014 offset+=step;
1015 }
1016 break;
1017 default:
1018 throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1019 }
1020 return &v;
1021 }
1022
1023 /*
1024 \brief Compute the value of the expression (unary operation) for the given sample.
1025 \return Vector which stores the value of the subexpression for the given sample.
1026 \param v A vector to store intermediate results.
1027 \param offset Index in v to begin storing results.
1028 \param sampleNo Sample number to evaluate.
1029 \param roffset (output parameter) the offset in the return vector where the result begins.
1030
1031 The return value will be an existing vector so do not deallocate it.
1032 If the result is stored in v it should be stored at the offset given.
1033 Everything from offset to the end of v should be considered available for this method to use.
1034 */
1035 DataTypes::ValueType*
1036 DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1037 {
1038 // we assume that any collapsing has been done before we get here
1039 // since we only have one argument we don't need to think about only
1040 // processing single points.
1041 if (m_readytype!='E')
1042 {
1043 throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
1044 }
1045 // since we can't write the result over the input, we need a result offset further along
1046 size_t subroffset;
1047 const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1048 LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1049 LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1050 roffset=offset;
1051 size_t loop=0;
1052 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1053 size_t outstep=getNoValues();
1054 size_t instep=m_left->getNoValues();
1055 LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1056 switch (m_op)
1057 {
1058 case TRACE:
1059 for (loop=0;loop<numsteps;++loop)
1060 {
1061 size_t zz=sampleNo*getNumDPPSample()+loop;
1062 if (zz==5800)
1063 {
1064 LAZYDEBUG(cerr << "point=" << zz<< endl;)
1065 LAZYDEBUG(cerr << "Input to trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
1066 LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
1067 LAZYDEBUG(cerr << subroffset << endl;)
1068 LAZYDEBUG(cerr << "output=" << offset << endl;)
1069 }
1070 DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1071 if (zz==5800)
1072 {
1073 LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1074 }
1075 subroffset+=instep;
1076 offset+=outstep;
1077 }
1078 break;
1079 case TRANS:
1080 for (loop=0;loop<numsteps;++loop)
1081 {
1082 DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1083 subroffset+=instep;
1084 offset+=outstep;
1085 }
1086 break;
1087 default:
1088 throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1089 }
1090 return &v;
1091 }
1092
1093
1094 #define PROC_OP(TYPE,X) \
1095 for (int j=0;j<onumsteps;++j)\
1096 {\
1097 for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1098 { \
1099 LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1100 LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1101 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1102 LAZYDEBUG(cout << " result= " << resultp[0] << endl;) \
1103 lroffset+=leftstep; \
1104 rroffset+=rightstep; \
1105 }\
1106 lroffset+=oleftstep;\
1107 rroffset+=orightstep;\
1108 }
1109
1110 /*
1111 \brief Compute the value of the expression (binary operation) for the given sample.
1112 \return Vector which stores the value of the subexpression for the given sample.
1113 \param v A vector to store intermediate results.
1114 \param offset Index in v to begin storing results.
1115 \param sampleNo Sample number to evaluate.
1116 \param roffset (output parameter) the offset in the return vector where the result begins.
1117
1118 The return value will be an existing vector so do not deallocate it.
1119 If the result is stored in v it should be stored at the offset given.
1120 Everything from offset to the end of v should be considered available for this method to use.
1121 */
1122 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1123 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1124 // If both children are expanded, then we can process them in a single operation (we treat
1125 // the whole sample as one big datapoint.
1126 // If one of the children is not expanded, then we need to treat each point in the sample
1127 // individually.
1128 // There is an additional complication when scalar operations are considered.
1129 // For example, 2+Vector.
1130 // In this case each double within the point is treated individually
1131 DataTypes::ValueType*
1132 DataLazy::resolveBinary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1133 {
1134 LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1135
1136 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1137 // first work out which of the children are expanded
1138 bool leftExp=(m_left->m_readytype=='E');
1139 bool rightExp=(m_right->m_readytype=='E');
1140 if (!leftExp && !rightExp)
1141 {
1142 throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1143 }
1144 bool leftScalar=(m_left->getRank()==0);
1145 bool rightScalar=(m_right->getRank()==0);
1146 if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1147 {
1148 throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1149 }
1150 size_t leftsize=m_left->getNoValues();
1151 size_t rightsize=m_right->getNoValues();
1152 size_t chunksize=1; // how many doubles will be processed in one go
1153 int leftstep=0; // how far should the left offset advance after each step
1154 int rightstep=0;
1155 int numsteps=0; // total number of steps for the inner loop
1156 int oleftstep=0; // the o variables refer to the outer loop
1157 int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1158 int onumsteps=1;
1159
1160 bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1161 bool RES=(rightExp && rightScalar);
1162 bool LS=(!leftExp && leftScalar); // left is a single scalar
1163 bool RS=(!rightExp && rightScalar);
1164 bool LN=(!leftExp && !leftScalar); // left is a single non-scalar
1165 bool RN=(!rightExp && !rightScalar);
1166 bool LEN=(leftExp && !leftScalar); // left is an expanded non-scalar
1167 bool REN=(rightExp && !rightScalar);
1168
1169 if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1170 {
1171 chunksize=m_left->getNumDPPSample()*leftsize;
1172 leftstep=0;
1173 rightstep=0;
1174 numsteps=1;
1175 }
1176 else if (LES || RES)
1177 {
1178 chunksize=1;
1179 if (LES) // left is an expanded scalar
1180 {
1181 if (RS)
1182 {
1183 leftstep=1;
1184 rightstep=0;
1185 numsteps=m_left->getNumDPPSample();
1186 }
1187 else // RN or REN
1188 {
1189 leftstep=0;
1190 oleftstep=1;
1191 rightstep=1;
1192 orightstep=(RN ? -(int)rightsize : 0);
1193 numsteps=rightsize;
1194 onumsteps=m_left->getNumDPPSample();
1195 }
1196 }
1197 else // right is an expanded scalar
1198 {
1199 if (LS)
1200 {
1201 rightstep=1;
1202 leftstep=0;
1203 numsteps=m_right->getNumDPPSample();
1204 }
1205 else
1206 {
1207 rightstep=0;
1208 orightstep=1;
1209 leftstep=1;
1210 oleftstep=(LN ? -(int)leftsize : 0);
1211 numsteps=leftsize;
1212 onumsteps=m_right->getNumDPPSample();
1213 }
1214 }
1215 }
1216 else // this leaves (LEN, RS), (LEN, RN) and their transposes
1217 {
1218 if (LEN) // and Right will be a single value
1219 {
1220 chunksize=rightsize;
1221 leftstep=rightsize;
1222 rightstep=0;
1223 numsteps=m_left->getNumDPPSample();
1224 if (RS)
1225 {
1226 numsteps*=leftsize;
1227 }
1228 }
1229 else // REN
1230 {
1231 chunksize=leftsize;
1232 rightstep=leftsize;
1233 leftstep=0;
1234 numsteps=m_right->getNumDPPSample();
1235 if (LS)
1236 {
1237 numsteps*=rightsize;
1238 }
1239 }
1240 }
1241
1242 int resultStep=max(leftstep,rightstep); // only one (at most) should be !=0
1243 // Get the values of sub-expressions
1244 const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset); // see note on
1245 // calcBufss for why we can't put offset as the 2nd param above
1246 const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1247 // the right child starts further along.
1248 LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1249 LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1250 LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1251 LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1252 LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1253 LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1254 LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN << endl;)
1255
1256
1257 double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
1258 switch(m_op)
1259 {
1260 case ADD:
1261 PROC_OP(NO_ARG,plus<double>());
1262 break;
1263 case SUB:
1264 PROC_OP(NO_ARG,minus<double>());
1265 break;
1266 case MUL:
1267 PROC_OP(NO_ARG,multiplies<double>());
1268 break;
1269 case DIV:
1270 PROC_OP(NO_ARG,divides<double>());
1271 break;
1272 case POW:
1273 PROC_OP(double (double,double),::pow);
1274 break;
1275 default:
1276 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1277 }
1278 roffset=offset;
1279 return &v;
1280 }
1281
1282
1283
1284 /*
1285 \brief Compute the value of the expression (tensor product) for the given sample.
1286 \return Vector which stores the value of the subexpression for the given sample.
1287 \param v A vector to store intermediate results.
1288 \param offset Index in v to begin storing results.
1289 \param sampleNo Sample number to evaluate.
1290 \param roffset (output parameter) the offset in the return vector where the result begins.
1291
1292 The return value will be an existing vector so do not deallocate it.
1293 If the result is stored in v it should be stored at the offset given.
1294 Everything from offset to the end of v should be considered available for this method to use.
1295 */
1296 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1297 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1298 // unlike the other resolve helpers, we must treat these datapoints separately.
1299 DataTypes::ValueType*
1300 DataLazy::resolveTProd(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1301 {
1302 LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << " to offset " << offset<< endl;)
1303
1304 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1305 // first work out which of the children are expanded
1306 bool leftExp=(m_left->m_readytype=='E');
1307 bool rightExp=(m_right->m_readytype=='E');
1308 int steps=getNumDPPSample();
1309 /* int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1310 int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1311 int leftStep=(leftExp? m_left->getNoValues() : 0); // do not have scalars as input to this method
1312 int rightStep=(rightExp?m_right->getNoValues() : 0);
1313
1314 int resultStep=getNoValues();
1315 // Get the values of sub-expressions (leave a gap of one sample for the result).
1316 int gap=offset+m_samplesize;
1317
1318 LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1319
1320 const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1321 gap+=m_left->getMaxSampleSize();
1322
1323
1324 LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1325
1326
1327 const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1328
1329 LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) << endl;
1330 cout << getNoValues() << endl;)
1331 LAZYDEBUG(cerr << "Result of left=";)
1332 LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1333
1334 for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1335 {
1336 cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1337 if (i%4==0) cout << endl;
1338 })
1339 LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1340 LAZYDEBUG(
1341 for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1342 {
1343 cerr << "[" << setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1344 if (i%4==0) cout << endl;
1345 }
1346 cerr << endl;
1347 )
1348 LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1349 LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1350 LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1351 LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1352 LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1353 LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1354 LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1355
1356 double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
1357 switch(m_op)
1358 {
1359 case PROD:
1360 for (int i=0;i<steps;++i,resultp+=resultStep)
1361 {
1362
1363 LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1364 LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1365 LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1366
1367 const double *ptr_0 = &((*left)[lroffset]);
1368 const double *ptr_1 = &((*right)[rroffset]);
1369
1370 LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1371 LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1372
1373 matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1374
1375 LAZYDEBUG(cout << "Results--\n";
1376 {
1377 DataVector dv(getNoValues());
1378 for (int z=0;z<getNoValues();++z)
1379 {
1380 cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1381 if (z%4==0) cout << endl;
1382 dv[z]=resultp[z];
1383 }
1384 cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1385 cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1386 }
1387 )
1388 lroffset+=leftStep;
1389 rroffset+=rightStep;
1390 }
1391 break;
1392 default:
1393 throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1394 }
1395 roffset=offset;
1396 return &v;
1397 }
1398
1399
1400
1401 /*
1402 \brief Compute the value of the expression for the given sample.
1403 \return Vector which stores the value of the subexpression for the given sample.
1404 \param v A vector to store intermediate results.
1405 \param offset Index in v to begin storing results.
1406 \param sampleNo Sample number to evaluate.
1407 \param roffset (output parameter) the offset in the return vector where the result begins.
1408
1409 The return value will be an existing vector so do not deallocate it.
1410 */
1411 // the vector and the offset are a place where the method could write its data if it wishes
1412 // it is not obligated to do so. For example, if it has its own storage already, it can use that.
1413 // Hence the return value to indicate where the data is actually stored.
1414 // Regardless, the storage should be assumed to be used, even if it isn't.
1415
1416 // the roffset is the offset within the returned vector where the data begins
1417 const DataTypes::ValueType*
1418 DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1419 {
1420 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1421 // collapse so we have a 'E' node or an IDENTITY for some other type
1422 if (m_readytype!='E' && m_op!=IDENTITY)
1423 {
1424 collapse();
1425 }
1426 if (m_op==IDENTITY)
1427 {
1428 const ValueType& vec=m_id->getVector();
1429 if (m_readytype=='C')
1430 {
1431 roffset=0;
1432 LAZYDEBUG(cout << "Finish sample " << toString() << endl;)
1433 return &(vec);
1434 }
1435 roffset=m_id->getPointOffset(sampleNo, 0);
1436 LAZYDEBUG(cout << "Finish sample " << toString() << endl;)
1437 return &(vec);
1438 }
1439 if (m_readytype!='E')
1440 {
1441 throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1442 }
1443 switch (getOpgroup(m_op))
1444 {
1445 case G_UNARY:
1446 case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1447 case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1448 case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1449 case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1450 case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1451 default:
1452 throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1453 }
1454
1455 }
1456
1457 // This needs to do the work of the idenity constructor
1458 void
1459 DataLazy::resolveToIdentity()
1460 {
1461 if (m_op==IDENTITY)
1462 return;
1463 DataReady_ptr p=resolve();
1464 makeIdentity(p);
1465 }
1466
1467 void DataLazy::makeIdentity(const DataReady_ptr& p)
1468 {
1469 m_op=IDENTITY;
1470 m_axis_offset=0;
1471 m_transpose=0;
1472 m_SL=m_SM=m_SR=0;
1473 m_children=m_height=0;
1474 m_id=p;
1475 if(p->isConstant()) {m_readytype='C';}
1476 else if(p->isExpanded()) {m_readytype='E';}
1477 else if (p->isTagged()) {m_readytype='T';}
1478 else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1479 m_buffsRequired=1;
1480 m_samplesize=p->getNumDPPSample()*p->getNoValues();
1481 m_maxsamplesize=m_samplesize;
1482 m_left.reset();
1483 m_right.reset();
1484 }
1485
1486 // To simplify the memory management, all threads operate on one large vector, rather than one each.
1487 // Each sample is evaluated independently and copied into the result DataExpanded.
1488 DataReady_ptr
1489 DataLazy::resolve()
1490 {
1491
1492 LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1493 LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1494
1495 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1496 {
1497 collapse();
1498 }
1499 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
1500 {
1501 return m_id;
1502 }
1503 // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1504 size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1505 // storage to evaluate its expression
1506 int numthreads=1;
1507 #ifdef _OPENMP
1508 numthreads=getNumberOfThreads();
1509 #endif
1510 ValueType v(numthreads*threadbuffersize);
1511 LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1512 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
1513 ValueType& resvec=result->getVector();
1514 DataReady_ptr resptr=DataReady_ptr(result);
1515 int sample;
1516 size_t outoffset; // offset in the output data
1517 int totalsamples=getNumSamples();
1518 const ValueType* res=0; // Vector storing the answer
1519 size_t resoffset=0; // where in the vector to find the answer
1520 LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1521 #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1522 for (sample=0;sample<totalsamples;++sample)
1523 {
1524 if (sample==0) {ENABLEDEBUG}
1525
1526 // if (sample==5800/getNumDPPSample()) {ENABLEDEBUG}
1527 LAZYDEBUG(cout << "################################# " << sample << endl;)
1528 #ifdef _OPENMP
1529 res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1530 #else
1531 res=resolveSample(v,0,sample,resoffset); // res would normally be v, but not if its a single IDENTITY op.
1532 #endif
1533 LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1534 LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1535 outoffset=result->getPointOffset(sample,0);
1536 LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1537 for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset) // copy values into the output vector
1538 {
1539 LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1540 resvec[outoffset]=(*res)[resoffset];
1541 }
1542 LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1543 LAZYDEBUG(cerr << "*********************************" << endl;)
1544 DISABLEDEBUG
1545 }
1546 return resptr;
1547 }
1548
1549 std::string
1550 DataLazy::toString() const
1551 {
1552 ostringstream oss;
1553 oss << "Lazy Data:";
1554 intoString(oss);
1555 return oss.str();
1556 }
1557
1558
1559 void
1560 DataLazy::intoString(ostringstream& oss) const
1561 {
1562 // oss << "[" << m_children <<";"<<m_height <<"]";
1563 switch (getOpgroup(m_op))
1564 {
1565 case G_IDENTITY:
1566 if (m_id->isExpanded())
1567 {
1568 oss << "E";
1569 }
1570 else if (m_id->isTagged())
1571 {
1572 oss << "T";
1573 }
1574 else if (m_id->isConstant())
1575 {
1576 oss << "C";
1577 }
1578 else
1579 {
1580 oss << "?";
1581 }
1582 oss << '@' << m_id.get();
1583 break;
1584 case G_BINARY:
1585 oss << '(';
1586 m_left->intoString(oss);
1587 oss << ' ' << opToString(m_op) << ' ';
1588 m_right->intoString(oss);
1589 oss << ')';
1590 break;
1591 case G_UNARY:
1592 case G_UNARY_P:
1593 case G_NP1OUT:
1594 case G_NP1OUT_P:
1595 oss << opToString(m_op) << '(';
1596 m_left->intoString(oss);
1597 oss << ')';
1598 break;
1599 case G_TENSORPROD:
1600 oss << opToString(m_op) << '(';
1601 m_left->intoString(oss);
1602 oss << ", ";
1603 m_right->intoString(oss);
1604 oss << ')';
1605 break;
1606 default:
1607 oss << "UNKNOWN";
1608 }
1609 }
1610
1611 DataAbstract*
1612 DataLazy::deepCopy()
1613 {
1614 switch (getOpgroup(m_op))
1615 {
1616 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
1617 case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1618 case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1619 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1620 case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1621 default:
1622 throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1623 }
1624 }
1625
1626
1627 // There is no single, natural interpretation of getLength on DataLazy.
1628 // Instances of DataReady can look at the size of their vectors.
1629 // For lazy though, it could be the size the data would be if it were resolved;
1630 // or it could be some function of the lengths of the DataReady instances which
1631 // form part of the expression.
1632 // Rather than have people making assumptions, I have disabled the method.
1633 DataTypes::ValueType::size_type
1634 DataLazy::getLength() const
1635 {
1636 throw DataException("getLength() does not make sense for lazy data.");
1637 }
1638
1639
1640 DataAbstract*
1641 DataLazy::getSlice(const DataTypes::RegionType& region) const
1642 {
1643 throw DataException("getSlice - not implemented for Lazy objects.");
1644 }
1645
1646
1647 // To do this we need to rely on our child nodes
1648 DataTypes::ValueType::size_type
1649 DataLazy::getPointOffset(int sampleNo,
1650 int dataPointNo)
1651 {
1652 if (m_op==IDENTITY)
1653 {
1654 return m_id->getPointOffset(sampleNo,dataPointNo);
1655 }
1656 if (m_readytype!='E')
1657 {
1658 collapse();
1659 return m_id->getPointOffset(sampleNo,dataPointNo);
1660 }
1661 // at this point we do not have an identity node and the expression will be Expanded
1662 // so we only need to know which child to ask
1663 if (m_left->m_readytype=='E')
1664 {
1665 return m_left->getPointOffset(sampleNo,dataPointNo);
1666 }
1667 else
1668 {
1669 return m_right->getPointOffset(sampleNo,dataPointNo);
1670 }
1671 }
1672
1673 // To do this we need to rely on our child nodes
1674 DataTypes::ValueType::size_type
1675 DataLazy::getPointOffset(int sampleNo,
1676 int dataPointNo) const
1677 {
1678 if (m_op==IDENTITY)
1679 {
1680 return m_id->getPointOffset(sampleNo,dataPointNo);
1681 }
1682 if (m_readytype=='E')
1683 {
1684 // at this point we do not have an identity node and the expression will be Expanded
1685 // so we only need to know which child to ask
1686 if (m_left->m_readytype=='E')
1687 {
1688 return m_left->getPointOffset(sampleNo,dataPointNo);
1689 }
1690 else
1691 {
1692 return m_right->getPointOffset(sampleNo,dataPointNo);
1693 }
1694 }
1695 if (m_readytype=='C')
1696 {
1697 return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1698 }
1699 throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1700 }
1701
1702 // It would seem that DataTagged will need to be treated differently since even after setting all tags
1703 // to zero, all the tags from all the DataTags would be in the result.
1704 // However since they all have the same value (0) whether they are there or not should not matter.
1705 // So I have decided that for all types this method will create a constant 0.
1706 // It can be promoted up as required.
1707 // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1708 // but we can deal with that if it arrises.
1709 void
1710 DataLazy::setToZero()
1711 {
1712 DataTypes::ValueType v(getNoValues(),0);
1713 m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1714 m_op=IDENTITY;
1715 m_right.reset();
1716 m_left.reset();
1717 m_readytype='C';
1718 m_buffsRequired=1;
1719 }
1720
1721 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26