/[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 2066 - (show annotations)
Thu Nov 20 05:31:33 2008 UTC (10 years, 5 months ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 39119 byte(s)
Fixed Data::toString to look at the amount of data actually stored rather than the number of points in the domain.

Added support for GeneralTensorProduct to LazyData
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 /*
32 How does DataLazy work?
33 ~~~~~~~~~~~~~~~~~~~~~~~
34
35 Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36 denoted left and right.
37
38 A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39 This means that all "internal" nodes in the structure are instances of DataLazy.
40
41 Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42 Note that IDENITY is not considered a unary operation.
43
44 I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45 It must however form a DAG (directed acyclic graph).
46 I will refer to individual DataLazy objects with the structure as nodes.
47
48 Each node also stores:
49 - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50 evaluated.
51 - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
52 evaluate the expression.
53 - m_samplesize ~ the number of doubles stored in a sample.
54
55 When a new node is created, the above values are computed based on the values in the child nodes.
56 Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
57
58 The resolve method, which produces a DataReady from a DataLazy, does the following:
59 1) Create a DataReady to hold the new result.
60 2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
61 3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
62
63 (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
64
65 resolveSample returns a Vector* and an offset within that vector where the result is stored.
66 Normally, this would be v, but for identity nodes their internal vector is returned instead.
67
68 The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
69
70 For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
71 The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
72
73 To add a new operator you need to do the following (plus anything I might have forgotten):
74 1) Add to the ES_optype.
75 2) determine what opgroup your operation belongs to (X)
76 3) add a string for the op to the end of ES_opstrings
77 4) increase ES_opcount
78 5) add an entry (X) to opgroups
79 6) add an entry to the switch in collapseToReady
80 7) add an entry to resolveX
81 */
82
83
84 using namespace std;
85 using namespace boost;
86
87 namespace escript
88 {
89
90 namespace
91 {
92
93 enum ES_opgroup
94 {
95 G_UNKNOWN,
96 G_IDENTITY,
97 G_BINARY, // pointwise operations with two arguments
98 G_UNARY, // pointwise operations with one argument
99 G_NP1OUT, // non-pointwise op with one output
100 G_TENSORPROD // general tensor product
101 };
102
103
104
105
106 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
107 "sin","cos","tan",
108 "asin","acos","atan","sinh","cosh","tanh","erf",
109 "asinh","acosh","atanh",
110 "log10","log","sign","abs","neg","pos","exp","sqrt",
111 "1/","where>0","where<0","where>=0","where<=0",
112 "symmetric","nonsymmetric",
113 "prod"};
114 int ES_opcount=36;
115 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
116 G_UNARY,G_UNARY,G_UNARY, //10
117 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
118 G_UNARY,G_UNARY,G_UNARY, // 20
119 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
120 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 33
121 G_NP1OUT,G_NP1OUT,
122 G_TENSORPROD};
123 inline
124 ES_opgroup
125 getOpgroup(ES_optype op)
126 {
127 return opgroups[op];
128 }
129
130 // return the FunctionSpace of the result of "left op right"
131 FunctionSpace
132 resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
133 {
134 // perhaps this should call interpolate and throw or something?
135 // maybe we need an interpolate node -
136 // that way, if interpolate is required in any other op we can just throw a
137 // programming error exception.
138
139 FunctionSpace l=left->getFunctionSpace();
140 FunctionSpace r=right->getFunctionSpace();
141 if (l!=r)
142 {
143 if (r.probeInterpolation(l))
144 {
145 return l;
146 }
147 if (l.probeInterpolation(r))
148 {
149 return r;
150 }
151 throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
152 }
153 return l;
154 }
155
156 // return the shape of the result of "left op right"
157 // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
158 DataTypes::ShapeType
159 resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
160 {
161 if (left->getShape()!=right->getShape())
162 {
163 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
164 {
165 throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
166 }
167 if (left->getRank()==0) // we need to allow scalar * anything
168 {
169 return right->getShape();
170 }
171 if (right->getRank()==0)
172 {
173 return left->getShape();
174 }
175 throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
176 }
177 return left->getShape();
178 }
179
180 // determine the output shape for the general tensor product operation
181 // the additional parameters return information required later for the product
182 // the majority of this code is copy pasted from C_General_Tensor_Product
183 DataTypes::ShapeType
184 GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
185 {
186
187 // Get rank and shape of inputs
188 int rank0 = left->getRank();
189 int rank1 = right->getRank();
190 const DataTypes::ShapeType& shape0 = left->getShape();
191 const DataTypes::ShapeType& shape1 = right->getShape();
192
193 // Prepare for the loops of the product and verify compatibility of shapes
194 int start0=0, start1=0;
195 if (transpose == 0) {}
196 else if (transpose == 1) { start0 = axis_offset; }
197 else if (transpose == 2) { start1 = rank1-axis_offset; }
198 else { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
199
200
201 // Adjust the shapes for transpose
202 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
203 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
204 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
205 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
206
207 // Prepare for the loops of the product
208 SL=1, SM=1, SR=1;
209 for (int i=0; i<rank0-axis_offset; i++) {
210 SL *= tmpShape0[i];
211 }
212 for (int i=rank0-axis_offset; i<rank0; i++) {
213 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
214 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
215 }
216 SM *= tmpShape0[i];
217 }
218 for (int i=axis_offset; i<rank1; i++) {
219 SR *= tmpShape1[i];
220 }
221
222 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
223 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
224 { // block to limit the scope of out_index
225 int out_index=0;
226 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
227 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
228 }
229 return shape2;
230 }
231
232
233 // determine the number of points in the result of "left op right"
234 // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here
235 // size_t
236 // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
237 // {
238 // switch (getOpgroup(op))
239 // {
240 // case G_BINARY: return left->getLength();
241 // case G_UNARY: return left->getLength();
242 // case G_NP1OUT: return left->getLength();
243 // default:
244 // throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
245 // }
246 // }
247
248 // determine the number of samples requires to evaluate an expression combining left and right
249 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
250 // The same goes for G_TENSORPROD
251 int
252 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
253 {
254 switch(getOpgroup(op))
255 {
256 case G_IDENTITY: return 1;
257 case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
258 case G_UNARY: return max(left->getBuffsRequired(),1);
259 case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
260 case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
261 default:
262 throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
263 }
264 }
265
266
267 } // end anonymous namespace
268
269
270
271 // Return a string representing the operation
272 const std::string&
273 opToString(ES_optype op)
274 {
275 if (op<0 || op>=ES_opcount)
276 {
277 op=UNKNOWNOP;
278 }
279 return ES_opstrings[op];
280 }
281
282
283 DataLazy::DataLazy(DataAbstract_ptr p)
284 : parent(p->getFunctionSpace(),p->getShape()),
285 m_op(IDENTITY),
286 m_axis_offset(0),
287 m_transpose(0),
288 m_SL(0), m_SM(0), m_SR(0)
289 {
290 if (p->isLazy())
291 {
292 // I don't want identity of Lazy.
293 // Question: Why would that be so bad?
294 // Answer: We assume that the child of ID is something we can call getVector on
295 throw DataException("Programmer error - attempt to create identity from a DataLazy.");
296 }
297 else
298 {
299 m_id=dynamic_pointer_cast<DataReady>(p);
300 if(p->isConstant()) {m_readytype='C';}
301 else if(p->isExpanded()) {m_readytype='E';}
302 else if (p->isTagged()) {m_readytype='T';}
303 else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
304 }
305 m_buffsRequired=1;
306 m_samplesize=getNumDPPSample()*getNoValues();
307 m_maxsamplesize=m_samplesize;
308 cout << "(1)Lazy created with " << m_samplesize << endl;
309 }
310
311
312
313
314 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
315 : parent(left->getFunctionSpace(),left->getShape()),
316 m_op(op),
317 m_axis_offset(0),
318 m_transpose(0),
319 m_SL(0), m_SM(0), m_SR(0)
320 {
321 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
322 {
323 throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
324 }
325
326 DataLazy_ptr lleft;
327 if (!left->isLazy())
328 {
329 lleft=DataLazy_ptr(new DataLazy(left));
330 }
331 else
332 {
333 lleft=dynamic_pointer_cast<DataLazy>(left);
334 }
335 m_readytype=lleft->m_readytype;
336 m_left=lleft;
337 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
338 m_samplesize=getNumDPPSample()*getNoValues();
339 m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
340 }
341
342
343 // In this constructor we need to consider interpolation
344 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
345 : parent(resultFS(left,right,op), resultShape(left,right,op)),
346 m_op(op),
347 m_SL(0), m_SM(0), m_SR(0)
348 {
349 if ((getOpgroup(op)!=G_BINARY))
350 {
351 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
352 }
353
354 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
355 {
356 FunctionSpace fs=getFunctionSpace();
357 Data ltemp(left);
358 Data tmp(ltemp,fs);
359 left=tmp.borrowDataPtr();
360 }
361 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
362 {
363 Data tmp(Data(right),getFunctionSpace());
364 right=tmp.borrowDataPtr();
365 }
366 left->operandCheck(*right);
367
368 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
369 {
370 m_left=dynamic_pointer_cast<DataLazy>(left);
371 }
372 else
373 {
374 m_left=DataLazy_ptr(new DataLazy(left));
375 }
376 if (right->isLazy())
377 {
378 m_right=dynamic_pointer_cast<DataLazy>(right);
379 }
380 else
381 {
382 m_right=DataLazy_ptr(new DataLazy(right));
383 }
384 char lt=m_left->m_readytype;
385 char rt=m_right->m_readytype;
386 if (lt=='E' || rt=='E')
387 {
388 m_readytype='E';
389 }
390 else if (lt=='T' || rt=='T')
391 {
392 m_readytype='T';
393 }
394 else
395 {
396 m_readytype='C';
397 }
398 m_samplesize=getNumDPPSample()*getNoValues();
399 m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
400 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
401 cout << "(3)Lazy created with " << m_samplesize << endl;
402 }
403
404 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
405 : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
406 m_op(op),
407 m_axis_offset(axis_offset),
408 m_transpose(transpose)
409 {
410 if ((getOpgroup(op)!=G_TENSORPROD))
411 {
412 throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
413 }
414 if ((transpose>2) || (transpose<0))
415 {
416 throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
417 }
418 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
419 {
420 FunctionSpace fs=getFunctionSpace();
421 Data ltemp(left);
422 Data tmp(ltemp,fs);
423 left=tmp.borrowDataPtr();
424 }
425 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
426 {
427 Data tmp(Data(right),getFunctionSpace());
428 right=tmp.borrowDataPtr();
429 }
430 left->operandCheck(*right);
431
432 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
433 {
434 m_left=dynamic_pointer_cast<DataLazy>(left);
435 }
436 else
437 {
438 m_left=DataLazy_ptr(new DataLazy(left));
439 }
440 if (right->isLazy())
441 {
442 m_right=dynamic_pointer_cast<DataLazy>(right);
443 }
444 else
445 {
446 m_right=DataLazy_ptr(new DataLazy(right));
447 }
448 char lt=m_left->m_readytype;
449 char rt=m_right->m_readytype;
450 if (lt=='E' || rt=='E')
451 {
452 m_readytype='E';
453 }
454 else if (lt=='T' || rt=='T')
455 {
456 m_readytype='T';
457 }
458 else
459 {
460 m_readytype='C';
461 }
462 m_samplesize=getNumDPPSample()*getNoValues();
463 m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
464 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
465 cout << "(4)Lazy created with " << m_samplesize << endl;
466 }
467
468
469 DataLazy::~DataLazy()
470 {
471 }
472
473
474 int
475 DataLazy::getBuffsRequired() const
476 {
477 return m_buffsRequired;
478 }
479
480
481 size_t
482 DataLazy::getMaxSampleSize() const
483 {
484 return m_maxsamplesize;
485 }
486
487 /*
488 \brief Evaluates the expression using methods on Data.
489 This does the work for the collapse method.
490 For reasons of efficiency do not call this method on DataExpanded nodes.
491 */
492 DataReady_ptr
493 DataLazy::collapseToReady()
494 {
495 if (m_readytype=='E')
496 { // this is more an efficiency concern than anything else
497 throw DataException("Programmer Error - do not use collapse on Expanded data.");
498 }
499 if (m_op==IDENTITY)
500 {
501 return m_id;
502 }
503 DataReady_ptr pleft=m_left->collapseToReady();
504 Data left(pleft);
505 Data right;
506 if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
507 {
508 right=Data(m_right->collapseToReady());
509 }
510 Data result;
511 switch(m_op)
512 {
513 case ADD:
514 result=left+right;
515 break;
516 case SUB:
517 result=left-right;
518 break;
519 case MUL:
520 result=left*right;
521 break;
522 case DIV:
523 result=left/right;
524 break;
525 case SIN:
526 result=left.sin();
527 break;
528 case COS:
529 result=left.cos();
530 break;
531 case TAN:
532 result=left.tan();
533 break;
534 case ASIN:
535 result=left.asin();
536 break;
537 case ACOS:
538 result=left.acos();
539 break;
540 case ATAN:
541 result=left.atan();
542 break;
543 case SINH:
544 result=left.sinh();
545 break;
546 case COSH:
547 result=left.cosh();
548 break;
549 case TANH:
550 result=left.tanh();
551 break;
552 case ERF:
553 result=left.erf();
554 break;
555 case ASINH:
556 result=left.asinh();
557 break;
558 case ACOSH:
559 result=left.acosh();
560 break;
561 case ATANH:
562 result=left.atanh();
563 break;
564 case LOG10:
565 result=left.log10();
566 break;
567 case LOG:
568 result=left.log();
569 break;
570 case SIGN:
571 result=left.sign();
572 break;
573 case ABS:
574 result=left.abs();
575 break;
576 case NEG:
577 result=left.neg();
578 break;
579 case POS:
580 // it doesn't mean anything for delayed.
581 // it will just trigger a deep copy of the lazy object
582 throw DataException("Programmer error - POS not supported for lazy data.");
583 break;
584 case EXP:
585 result=left.exp();
586 break;
587 case SQRT:
588 result=left.sqrt();
589 break;
590 case RECIP:
591 result=left.oneOver();
592 break;
593 case GZ:
594 result=left.wherePositive();
595 break;
596 case LZ:
597 result=left.whereNegative();
598 break;
599 case GEZ:
600 result=left.whereNonNegative();
601 break;
602 case LEZ:
603 result=left.whereNonPositive();
604 break;
605 case SYM:
606 result=left.symmetric();
607 break;
608 case NSYM:
609 result=left.nonsymmetric();
610 break;
611 case PROD:
612 result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
613 break;
614 default:
615 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
616 }
617 return result.borrowReadyPtr();
618 }
619
620 /*
621 \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
622 This method uses the original methods on the Data class to evaluate the expressions.
623 For this reason, it should not be used on DataExpanded instances. (To do so would defeat
624 the purpose of using DataLazy in the first place).
625 */
626 void
627 DataLazy::collapse()
628 {
629 if (m_op==IDENTITY)
630 {
631 return;
632 }
633 if (m_readytype=='E')
634 { // this is more an efficiency concern than anything else
635 throw DataException("Programmer Error - do not use collapse on Expanded data.");
636 }
637 m_id=collapseToReady();
638 m_op=IDENTITY;
639 }
640
641 /*
642 \brief Compute the value of the expression (unary operation) for the given sample.
643 \return Vector which stores the value of the subexpression for the given sample.
644 \param v A vector to store intermediate results.
645 \param offset Index in v to begin storing results.
646 \param sampleNo Sample number to evaluate.
647 \param roffset (output parameter) the offset in the return vector where the result begins.
648
649 The return value will be an existing vector so do not deallocate it.
650 If the result is stored in v it should be stored at the offset given.
651 Everything from offset to the end of v should be considered available for this method to use.
652 */
653 DataTypes::ValueType*
654 DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
655 {
656 // we assume that any collapsing has been done before we get here
657 // since we only have one argument we don't need to think about only
658 // processing single points.
659 if (m_readytype!='E')
660 {
661 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
662 }
663 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
664 const double* left=&((*vleft)[roffset]);
665 double* result=&(v[offset]);
666 roffset=offset;
667 switch (m_op)
668 {
669 case SIN:
670 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
671 break;
672 case COS:
673 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
674 break;
675 case TAN:
676 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
677 break;
678 case ASIN:
679 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
680 break;
681 case ACOS:
682 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
683 break;
684 case ATAN:
685 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
686 break;
687 case SINH:
688 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
689 break;
690 case COSH:
691 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
692 break;
693 case TANH:
694 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
695 break;
696 case ERF:
697 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
698 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
699 #else
700 tensor_unary_operation(m_samplesize, left, result, ::erf);
701 break;
702 #endif
703 case ASINH:
704 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
705 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
706 #else
707 tensor_unary_operation(m_samplesize, left, result, ::asinh);
708 #endif
709 break;
710 case ACOSH:
711 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
712 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
713 #else
714 tensor_unary_operation(m_samplesize, left, result, ::acosh);
715 #endif
716 break;
717 case ATANH:
718 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
719 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
720 #else
721 tensor_unary_operation(m_samplesize, left, result, ::atanh);
722 #endif
723 break;
724 case LOG10:
725 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
726 break;
727 case LOG:
728 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
729 break;
730 case SIGN:
731 tensor_unary_operation(m_samplesize, left, result, escript::fsign);
732 break;
733 case ABS:
734 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
735 break;
736 case NEG:
737 tensor_unary_operation(m_samplesize, left, result, negate<double>());
738 break;
739 case POS:
740 // it doesn't mean anything for delayed.
741 // it will just trigger a deep copy of the lazy object
742 throw DataException("Programmer error - POS not supported for lazy data.");
743 break;
744 case EXP:
745 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
746 break;
747 case SQRT:
748 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
749 break;
750 case RECIP:
751 tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
752 break;
753 case GZ:
754 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
755 break;
756 case LZ:
757 tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
758 break;
759 case GEZ:
760 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
761 break;
762 case LEZ:
763 tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
764 break;
765
766 default:
767 throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
768 }
769 return &v;
770 }
771
772
773 /*
774 \brief Compute the value of the expression (unary operation) for the given sample.
775 \return Vector which stores the value of the subexpression for the given sample.
776 \param v A vector to store intermediate results.
777 \param offset Index in v to begin storing results.
778 \param sampleNo Sample number to evaluate.
779 \param roffset (output parameter) the offset in the return vector where the result begins.
780
781 The return value will be an existing vector so do not deallocate it.
782 If the result is stored in v it should be stored at the offset given.
783 Everything from offset to the end of v should be considered available for this method to use.
784 */
785 DataTypes::ValueType*
786 DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
787 {
788 // we assume that any collapsing has been done before we get here
789 // since we only have one argument we don't need to think about only
790 // processing single points.
791 if (m_readytype!='E')
792 {
793 throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
794 }
795 // since we can't write the result over the input, we need a result offset further along
796 size_t subroffset=roffset+m_samplesize;
797 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
798 roffset=offset;
799 switch (m_op)
800 {
801 case SYM:
802 DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
803 break;
804 case NSYM:
805 DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
806 break;
807 default:
808 throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
809 }
810 return &v;
811 }
812
813
814
815
816 #define PROC_OP(TYPE,X) \
817 for (int i=0;i<steps;++i,resultp+=resultStep) \
818 { \
819 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
820 lroffset+=leftStep; \
821 rroffset+=rightStep; \
822 }
823
824 /*
825 \brief Compute the value of the expression (binary operation) for the given sample.
826 \return Vector which stores the value of the subexpression for the given sample.
827 \param v A vector to store intermediate results.
828 \param offset Index in v to begin storing results.
829 \param sampleNo Sample number to evaluate.
830 \param roffset (output parameter) the offset in the return vector where the result begins.
831
832 The return value will be an existing vector so do not deallocate it.
833 If the result is stored in v it should be stored at the offset given.
834 Everything from offset to the end of v should be considered available for this method to use.
835 */
836 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
837 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
838 // If both children are expanded, then we can process them in a single operation (we treat
839 // the whole sample as one big datapoint.
840 // If one of the children is not expanded, then we need to treat each point in the sample
841 // individually.
842 // There is an additional complication when scalar operations are considered.
843 // For example, 2+Vector.
844 // In this case each double within the point is treated individually
845 DataTypes::ValueType*
846 DataLazy::resolveBinary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
847 {
848 cout << "Resolve binary: " << toString() << endl;
849
850 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
851 // first work out which of the children are expanded
852 bool leftExp=(m_left->m_readytype=='E');
853 bool rightExp=(m_right->m_readytype=='E');
854 bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
855 int steps=(bigloops?1:getNumDPPSample());
856 size_t chunksize=(bigloops? m_samplesize : getNoValues()); // if bigloops, pretend the whole sample is a datapoint
857 if (m_left->getRank()!=m_right->getRank()) // need to deal with scalar * ? ops
858 {
859 EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
860 steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
861 chunksize=1; // for scalar
862 }
863 int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
864 int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
865 int resultStep=max(leftStep,rightStep); // only one (at most) should be !=0
866 // Get the values of sub-expressions
867 const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
868 const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
869 // the right child starts further along.
870 double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
871 switch(m_op)
872 {
873 case ADD:
874 PROC_OP(NO_ARG,plus<double>());
875 break;
876 case SUB:
877 PROC_OP(NO_ARG,minus<double>());
878 break;
879 case MUL:
880 PROC_OP(NO_ARG,multiplies<double>());
881 break;
882 case DIV:
883 PROC_OP(NO_ARG,divides<double>());
884 break;
885 case POW:
886 PROC_OP(double (double,double),::pow);
887 break;
888 default:
889 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
890 }
891 roffset=offset;
892 return &v;
893 }
894
895
896 /*
897 \brief Compute the value of the expression (tensor product) for the given sample.
898 \return Vector which stores the value of the subexpression for the given sample.
899 \param v A vector to store intermediate results.
900 \param offset Index in v to begin storing results.
901 \param sampleNo Sample number to evaluate.
902 \param roffset (output parameter) the offset in the return vector where the result begins.
903
904 The return value will be an existing vector so do not deallocate it.
905 If the result is stored in v it should be stored at the offset given.
906 Everything from offset to the end of v should be considered available for this method to use.
907 */
908 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
909 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
910 // unlike the other resolve helpers, we must treat these datapoints separately.
911 DataTypes::ValueType*
912 DataLazy::resolveTProd(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
913 {
914 cout << "Resolve TensorProduct: " << toString() << endl;
915
916 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
917 // first work out which of the children are expanded
918 bool leftExp=(m_left->m_readytype=='E');
919 bool rightExp=(m_right->m_readytype=='E');
920 int steps=getNumDPPSample();
921 int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
922 int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
923 int resultStep=max(leftStep,rightStep); // only one (at most) should be !=0
924 // Get the values of sub-expressions (leave a gap of one sample for the result).
925 const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);
926 const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);
927 double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
928 switch(m_op)
929 {
930 case PROD:
931 for (int i=0;i<steps;++i,resultp+=resultStep)
932 {
933 const double *ptr_0 = &((*left)[lroffset]);
934 const double *ptr_1 = &((*right)[rroffset]);
935 matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
936 lroffset+=leftStep;
937 rroffset+=rightStep;
938 }
939 break;
940 default:
941 throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
942 }
943 roffset=offset;
944 return &v;
945 }
946
947
948
949 /*
950 \brief Compute the value of the expression for the given sample.
951 \return Vector which stores the value of the subexpression for the given sample.
952 \param v A vector to store intermediate results.
953 \param offset Index in v to begin storing results.
954 \param sampleNo Sample number to evaluate.
955 \param roffset (output parameter) the offset in the return vector where the result begins.
956
957 The return value will be an existing vector so do not deallocate it.
958 */
959 // the vector and the offset are a place where the method could write its data if it wishes
960 // it is not obligated to do so. For example, if it has its own storage already, it can use that.
961 // Hence the return value to indicate where the data is actually stored.
962 // Regardless, the storage should be assumed to be used, even if it isn't.
963
964 // the roffset is the offset within the returned vector where the data begins
965 const DataTypes::ValueType*
966 DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
967 {
968 cout << "Resolve sample " << toString() << endl;
969 // collapse so we have a 'E' node or an IDENTITY for some other type
970 if (m_readytype!='E' && m_op!=IDENTITY)
971 {
972 collapse();
973 }
974 if (m_op==IDENTITY)
975 {
976 const ValueType& vec=m_id->getVector();
977 if (m_readytype=='C')
978 {
979 roffset=0;
980 return &(vec);
981 }
982 roffset=m_id->getPointOffset(sampleNo, 0);
983 return &(vec);
984 }
985 if (m_readytype!='E')
986 {
987 throw DataException("Programmer Error - Collapse did not produce an expanded node.");
988 }
989 switch (getOpgroup(m_op))
990 {
991 case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
992 case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
993 case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
994 case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
995 default:
996 throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
997 }
998 }
999
1000
1001 // To simplify the memory management, all threads operate on one large vector, rather than one each.
1002 // Each sample is evaluated independently and copied into the result DataExpanded.
1003 DataReady_ptr
1004 DataLazy::resolve()
1005 {
1006
1007 cout << "Sample size=" << m_samplesize << endl;
1008 cout << "Buffers=" << m_buffsRequired << endl;
1009
1010 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1011 {
1012 collapse();
1013 }
1014 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
1015 {
1016 return m_id;
1017 }
1018 // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1019 size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1020 // storage to evaluate its expression
1021 int numthreads=1;
1022 #ifdef _OPENMP
1023 numthreads=getNumberOfThreads();
1024 int threadnum=0;
1025 #endif
1026 ValueType v(numthreads*threadbuffersize);
1027 cout << "Buffer created with size=" << v.size() << endl;
1028 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
1029 ValueType& resvec=result->getVector();
1030 DataReady_ptr resptr=DataReady_ptr(result);
1031 int sample;
1032 size_t outoffset; // offset in the output data
1033 int totalsamples=getNumSamples();
1034 const ValueType* res=0; // Vector storing the answer
1035 size_t resoffset=0; // where in the vector to find the answer
1036 #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
1037 for (sample=0;sample<totalsamples;++sample)
1038 {
1039 cout << "################################# " << sample << endl;
1040 #ifdef _OPENMP
1041 res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1042 #else
1043 res=resolveSample(v,0,sample,resoffset); // res would normally be v, but not if its a single IDENTITY op.
1044 #endif
1045 cerr << "-------------------------------- " << endl;
1046 outoffset=result->getPointOffset(sample,0);
1047 cerr << "offset=" << outoffset << endl;
1048 for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset) // copy values into the output vector
1049 {
1050 resvec[outoffset]=(*res)[resoffset];
1051 }
1052 cerr << "*********************************" << endl;
1053 }
1054 return resptr;
1055 }
1056
1057 std::string
1058 DataLazy::toString() const
1059 {
1060 ostringstream oss;
1061 oss << "Lazy Data:";
1062 intoString(oss);
1063 return oss.str();
1064 }
1065
1066
1067 void
1068 DataLazy::intoString(ostringstream& oss) const
1069 {
1070 switch (getOpgroup(m_op))
1071 {
1072 case G_IDENTITY:
1073 if (m_id->isExpanded())
1074 {
1075 oss << "E";
1076 }
1077 else if (m_id->isTagged())
1078 {
1079 oss << "T";
1080 }
1081 else if (m_id->isConstant())
1082 {
1083 oss << "C";
1084 }
1085 else
1086 {
1087 oss << "?";
1088 }
1089 oss << '@' << m_id.get();
1090 break;
1091 case G_BINARY:
1092 oss << '(';
1093 m_left->intoString(oss);
1094 oss << ' ' << opToString(m_op) << ' ';
1095 m_right->intoString(oss);
1096 oss << ')';
1097 break;
1098 case G_UNARY:
1099 case G_NP1OUT:
1100 oss << opToString(m_op) << '(';
1101 m_left->intoString(oss);
1102 oss << ')';
1103 break;
1104 case G_TENSORPROD:
1105 oss << opToString(m_op) << '(';
1106 m_left->intoString(oss);
1107 oss << ", ";
1108 m_right->intoString(oss);
1109 oss << ')';
1110 break;
1111 default:
1112 oss << "UNKNOWN";
1113 }
1114 }
1115
1116 DataAbstract*
1117 DataLazy::deepCopy()
1118 {
1119 switch (getOpgroup(m_op))
1120 {
1121 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
1122 case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1123 case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1124 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1125 case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1126 default:
1127 throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1128 }
1129 }
1130
1131
1132 // There is no single, natural interpretation of getLength on DataLazy.
1133 // Instances of DataReady can look at the size of their vectors.
1134 // For lazy though, it could be the size the data would be if it were resolved;
1135 // or it could be some function of the lengths of the DataReady instances which
1136 // form part of the expression.
1137 // Rather than have people making assumptions, I have disabled the method.
1138 DataTypes::ValueType::size_type
1139 DataLazy::getLength() const
1140 {
1141 throw DataException("getLength() does not make sense for lazy data.");
1142 }
1143
1144
1145 DataAbstract*
1146 DataLazy::getSlice(const DataTypes::RegionType& region) const
1147 {
1148 throw DataException("getSlice - not implemented for Lazy objects.");
1149 }
1150
1151
1152 // To do this we need to rely on our child nodes
1153 DataTypes::ValueType::size_type
1154 DataLazy::getPointOffset(int sampleNo,
1155 int dataPointNo)
1156 {
1157 if (m_op==IDENTITY)
1158 {
1159 return m_id->getPointOffset(sampleNo,dataPointNo);
1160 }
1161 if (m_readytype!='E')
1162 {
1163 collapse();
1164 return m_id->getPointOffset(sampleNo,dataPointNo);
1165 }
1166 // at this point we do not have an identity node and the expression will be Expanded
1167 // so we only need to know which child to ask
1168 if (m_left->m_readytype=='E')
1169 {
1170 return m_left->getPointOffset(sampleNo,dataPointNo);
1171 }
1172 else
1173 {
1174 return m_right->getPointOffset(sampleNo,dataPointNo);
1175 }
1176 }
1177
1178 // To do this we need to rely on our child nodes
1179 DataTypes::ValueType::size_type
1180 DataLazy::getPointOffset(int sampleNo,
1181 int dataPointNo) const
1182 {
1183 if (m_op==IDENTITY)
1184 {
1185 return m_id->getPointOffset(sampleNo,dataPointNo);
1186 }
1187 if (m_readytype=='E')
1188 {
1189 // at this point we do not have an identity node and the expression will be Expanded
1190 // so we only need to know which child to ask
1191 if (m_left->m_readytype=='E')
1192 {
1193 return m_left->getPointOffset(sampleNo,dataPointNo);
1194 }
1195 else
1196 {
1197 return m_right->getPointOffset(sampleNo,dataPointNo);
1198 }
1199 }
1200 if (m_readytype=='C')
1201 {
1202 return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1203 }
1204 throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1205 }
1206
1207 // It would seem that DataTagged will need to be treated differently since even after setting all tags
1208 // to zero, all the tags from all the DataTags would be in the result.
1209 // However since they all have the same value (0) whether they are there or not should not matter.
1210 // So I have decided that for all types this method will create a constant 0.
1211 // It can be promoted up as required.
1212 // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1213 // but we can deal with that if it arrises.
1214 void
1215 DataLazy::setToZero()
1216 {
1217 DataTypes::ValueType v(getNoValues(),0);
1218 m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1219 m_op=IDENTITY;
1220 m_right.reset();
1221 m_left.reset();
1222 m_readytype='C';
1223 m_buffsRequired=1;
1224 }
1225
1226 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26