/[escript]/branches/trilinos_from_5897/escriptcore/src/Data.h
ViewVC logotype

Contents of /branches/trilinos_from_5897/escriptcore/src/Data.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5963 - (show annotations)
Mon Feb 22 06:59:27 2016 UTC (3 years, 3 months ago) by caltinay
File MIME type: text/plain
File size: 110905 byte(s)
sync and fix.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /** \file Data.h */
19
20 #ifndef DATA_H
21 #define DATA_H
22 #include "system_dep.h"
23
24 #include "DataTypes.h"
25 #include "DataAbstract.h"
26 #include "DataAlgorithm.h"
27 #include "FunctionSpace.h"
28 #include "BinaryOp.h"
29 #include "UnaryOp.h"
30 #include "DataException.h"
31
32 #ifdef _OPENMP
33 #include <omp.h>
34 #endif
35
36 #include "esysUtils/Esys_MPI.h"
37 #include <string>
38 #include <algorithm>
39 #include <sstream>
40
41 #include <boost/shared_ptr.hpp>
42 #include <boost/python/object.hpp>
43 #include <boost/python/tuple.hpp>
44 #include <boost/math/special_functions/bessel.hpp>
45
46 namespace escript {
47
48 //
49 // Forward declaration for various implementations of Data.
50 class DataConstant;
51 class DataTagged;
52 class DataExpanded;
53 class DataLazy;
54
55 /**
56 \brief
57 Data represents a collection of datapoints.
58
59 Description:
60 Internally, the datapoints are actually stored by a DataAbstract object.
61 The specific instance of DataAbstract used may vary over the lifetime
62 of the Data object.
63 Some methods on this class return references (eg getShape()).
64 These references should not be used after an operation which changes the underlying DataAbstract object.
65 Doing so will lead to invalid memory access.
66 This should not affect any methods exposed via boost::python.
67 */
68 class Data {
69
70 public:
71
72 /**
73 Constructors.
74 */
75
76 /**
77 \brief
78 Default constructor.
79 Creates a DataEmpty object.
80 */
81 ESCRIPT_DLL_API
82 Data();
83
84 /**
85 \brief
86 Copy constructor.
87 WARNING: Only performs a shallow copy.
88 */
89 ESCRIPT_DLL_API
90 Data(const Data& inData);
91
92 /**
93 \brief
94 Constructor from another Data object. If "what" is different from the
95 function space of inData the inData are tried to be interpolated to what,
96 otherwise a shallow copy of inData is returned.
97 */
98 ESCRIPT_DLL_API
99 Data(const Data& inData,
100 const FunctionSpace& what);
101
102 /**
103 \brief Copy Data from an existing vector
104 */
105
106 ESCRIPT_DLL_API
107 Data(const DataTypes::RealVectorType& value,
108 const DataTypes::ShapeType& shape,
109 const FunctionSpace& what=FunctionSpace(),
110 bool expanded=false);
111
112 /**
113 \brief
114 Constructor which creates a Data with points having the specified shape.
115
116 \param value - Input - Single value applied to all Data.
117 \param dataPointShape - Input - The shape of each data point.
118 \param what - Input - A description of what this data represents.
119 \param expanded - Input - Flag, if true fill the entire container with
120 the given value. Otherwise a more efficient storage
121 mechanism will be used.
122 */
123 ESCRIPT_DLL_API
124 Data(DataTypes::real_t value,
125 const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
126 const FunctionSpace& what=FunctionSpace(),
127 bool expanded=false);
128
129 /**
130 \brief
131 Constructor which performs a deep copy of a region from another Data object.
132
133 \param inData - Input - Input Data object.
134 \param region - Input - Region to copy.
135 */
136 ESCRIPT_DLL_API
137 Data(const Data& inData,
138 const DataTypes::RegionType& region);
139
140 /**
141 \brief
142 Constructor which copies data from any object that can be treated like a python array/sequence.
143
144 \param value - Input - Input data.
145 \param what - Input - A description of what this data represents.
146 \param expanded - Input - Flag, if true fill the entire container with
147 the value. Otherwise a more efficient storage
148 mechanism will be used.
149 */
150 ESCRIPT_DLL_API
151 Data(const boost::python::object& value,
152 const FunctionSpace& what=FunctionSpace(),
153 bool expanded=false);
154
155 /**
156 \brief
157 Constructor which copies data from a wrapped array.
158
159 \param w - Input - Input data.
160 \param what - Input - A description of what this data represents.
161 \param expanded - Input - Flag, if true fill the entire container with
162 the value. Otherwise a more efficient storage
163 mechanism will be used.
164 */
165 ESCRIPT_DLL_API
166 Data(const WrappedArray& w, const FunctionSpace& what,
167 bool expanded=false);
168
169
170 /**
171 \brief
172 Constructor which creates a DataConstant.
173 Copies data from any object that can be treated like a python array/sequence.
174 All other parameters are copied from other.
175
176 \param value - Input - Input data.
177 \param other - Input - contains all other parameters.
178 */
179 ESCRIPT_DLL_API
180 Data(const boost::python::object& value,
181 const Data& other);
182
183 /**
184 \brief
185 Constructor which creates a DataConstant of "shape" with constant value.
186 */
187 ESCRIPT_DLL_API
188 Data(DataTypes::real_t value,
189 const boost::python::tuple& shape=boost::python::make_tuple(),
190 const FunctionSpace& what=FunctionSpace(),
191 bool expanded=false);
192
193 /**
194 \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
195 Once you have passed the pointer, do not delete it.
196 */
197 ESCRIPT_DLL_API
198 explicit Data(DataAbstract* underlyingdata);
199
200 /**
201 \brief Create a Data based on the supplied DataAbstract
202 */
203 ESCRIPT_DLL_API
204 explicit Data(DataAbstract_ptr underlyingdata);
205
206 /**
207 \brief
208 Destructor
209 */
210 ESCRIPT_DLL_API
211 ~Data();
212
213 /**
214 \brief Make this object a deep copy of "other".
215 */
216 ESCRIPT_DLL_API
217 void
218 copy(const Data& other);
219
220 /**
221 \brief Return a pointer to a deep copy of this object.
222 */
223 ESCRIPT_DLL_API
224 Data
225 copySelf() const;
226
227
228 /**
229 \brief produce a delayed evaluation version of this Data.
230 */
231 ESCRIPT_DLL_API
232 Data
233 delay();
234
235 /**
236 \brief convert the current data into lazy data.
237 */
238 ESCRIPT_DLL_API
239 void
240 delaySelf();
241
242
243 /**
244 Member access methods.
245 */
246
247 /**
248 \brief
249 switches on update protection
250
251 */
252 ESCRIPT_DLL_API
253 void
254 setProtection();
255
256 /**
257 \brief
258 Returns true, if the data object is protected against update
259
260 */
261 ESCRIPT_DLL_API
262 bool
263 isProtected() const;
264
265
266 /**
267 \brief
268 Return the value of a data point as a python tuple.
269 */
270 ESCRIPT_DLL_API
271 const boost::python::object
272 getValueOfDataPointAsTuple(int dataPointNo);
273
274 /**
275 \brief
276 sets the values of a data-point from a python object on this process
277 */
278 ESCRIPT_DLL_API
279 void
280 setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
281
282 /**
283 \brief
284 sets the values of a data-point from a array-like object on this process
285 */
286 ESCRIPT_DLL_API
287 void
288 setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
289
290 /**
291 \brief
292 sets the values of a data-point on this process
293 */
294 ESCRIPT_DLL_API
295 void
296 setValueOfDataPoint(int dataPointNo, const DataTypes::real_t);
297
298 ESCRIPT_DLL_API
299 void
300 setValueOfDataPointC(int dataPointNo, const DataTypes::cplx_t);
301
302
303 /**
304 \brief Return a data point across all processors as a python tuple.
305 */
306 ESCRIPT_DLL_API
307 const boost::python::object
308 getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
309
310
311 /**
312 \brief Set the value of a global data point
313 */
314 ESCRIPT_DLL_API
315 void
316 setTupleForGlobalDataPoint(int id, int proc, boost::python::object);
317
318 /**
319 \brief
320 Return the tag number associated with the given data-point.
321
322 */
323 ESCRIPT_DLL_API
324 int
325 getTagNumber(int dpno);
326
327
328 /**
329 \brief
330 Write the data as a string. For large amounts of data, a summary is printed.
331 */
332 ESCRIPT_DLL_API
333 std::string
334 toString() const;
335
336 /**
337 \brief
338 Whatever the current Data type make this into a DataExpanded.
339 */
340 ESCRIPT_DLL_API
341 void
342 expand();
343
344 /**
345 \brief
346 If possible convert this Data to DataTagged. This will only allow
347 Constant data to be converted to tagged. An attempt to convert
348 Expanded data to tagged will throw an exception.
349 */
350 ESCRIPT_DLL_API
351 void
352 tag();
353
354 /**
355 \brief If this data is lazy, then convert it to ready data.
356 What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
357 */
358 ESCRIPT_DLL_API
359 void
360 resolve();
361
362 /**
363 \brief returns return true if data contains NaN.
364 \warning This is dependent on the ability to reliably detect NaNs on your compiler.
365 See the nancheck function in LocalOps for details.
366 */
367 ESCRIPT_DLL_API
368 bool
369 hasNaN();
370
371 /**
372 \brief replaces all NaN values with value
373 */
374 ESCRIPT_DLL_API
375 void
376 replaceNaN(DataTypes::real_t value);
377
378 /**
379 \brief Ensures data is ready for write access.
380 This means that the data will be resolved if lazy and will be copied if shared with another Data object.
381 \warning This method should only be called in single threaded sections of code. (It modifies m_data).
382 Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
383 Doing so might introduce additional sharing.
384 */
385 ESCRIPT_DLL_API
386 void
387 requireWrite();
388
389 /**
390 \brief
391 Return true if this Data is expanded.
392 \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
393 */
394 ESCRIPT_DLL_API
395 bool
396 isExpanded() const;
397
398 /**
399 \brief
400 Return true if this Data is expanded or resolves to expanded.
401 That is, if it has a separate value for each datapoint in the sample.
402 */
403 ESCRIPT_DLL_API
404 bool
405 actsExpanded() const;
406
407
408 /**
409 \brief
410 Return true if this Data is tagged.
411 */
412 ESCRIPT_DLL_API
413 bool
414 isTagged() const;
415
416 /**
417 \brief
418 Return true if this Data is constant.
419 */
420 ESCRIPT_DLL_API
421 bool
422 isConstant() const;
423
424 /**
425 \brief Return true if this Data is lazy.
426 */
427 ESCRIPT_DLL_API
428 bool
429 isLazy() const;
430
431 /**
432 \brief Return true if this data is ready.
433 */
434 ESCRIPT_DLL_API
435 bool
436 isReady() const;
437
438 /**
439 \brief
440 Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
441 contains datapoints.
442 */
443 ESCRIPT_DLL_API
444 bool
445 isEmpty() const;
446
447 /**
448 \brief
449 True if components of this data are stored as complex
450 */
451 ESCRIPT_DLL_API
452 bool
453 isComplex() const;
454
455 /**
456 \brief
457 Return the function space.
458 */
459 ESCRIPT_DLL_API
460 inline
461 const FunctionSpace&
462 getFunctionSpace() const
463 {
464 return m_data->getFunctionSpace();
465 }
466
467 /**
468 \brief
469 Return the domain.
470 */
471 ESCRIPT_DLL_API
472 inline
473 // const AbstractDomain&
474 const_Domain_ptr
475 getDomain() const
476 {
477 return getFunctionSpace().getDomain();
478 }
479
480
481 /**
482 \brief
483 Return the domain.
484 TODO: For internal use only. This should be removed.
485 */
486 ESCRIPT_DLL_API
487 inline
488 // const AbstractDomain&
489 Domain_ptr
490 getDomainPython() const
491 {
492 return getFunctionSpace().getDomainPython();
493 }
494
495 /**
496 \brief
497 Return the rank of the point data.
498 */
499 ESCRIPT_DLL_API
500 inline
501 unsigned int
502 getDataPointRank() const
503 {
504 return m_data->getRank();
505 }
506
507 /**
508 \brief
509 Return the number of data points
510 */
511 ESCRIPT_DLL_API
512 inline
513 int
514 getNumDataPoints() const
515 {
516 return getNumSamples() * getNumDataPointsPerSample();
517 }
518 /**
519 \brief
520 Return the number of samples.
521 */
522 ESCRIPT_DLL_API
523 inline
524 int
525 getNumSamples() const
526 {
527 return m_data->getNumSamples();
528 }
529
530 /**
531 \brief
532 Return the number of data points per sample.
533 */
534 ESCRIPT_DLL_API
535 inline
536 int
537 getNumDataPointsPerSample() const
538 {
539 return m_data->getNumDPPSample();
540 }
541
542 /**
543 \brief
544 Returns true if the number of data points per sample and the number of
545 samples match the respective argument. DataEmpty always returns true.
546 */
547 ESCRIPT_DLL_API
548 inline
549 bool numSamplesEqual(int numDataPointsPerSample, int numSamples) const
550 {
551 return (isEmpty() ||
552 (numDataPointsPerSample==getNumDataPointsPerSample() && numSamples==getNumSamples()));
553 }
554
555 /**
556 \brief
557 Returns true if the shape matches the vector (dimensions[0],...,
558 dimensions[rank-1]). DataEmpty always returns true.
559 */
560 inline
561 bool isDataPointShapeEqual(int rank, const int* dimensions) const
562 {
563 if (isEmpty())
564 return true;
565 const DataTypes::ShapeType givenShape(&dimensions[0],&dimensions[rank]);
566 return (getDataPointShape()==givenShape);
567 }
568
569 /**
570 \brief
571 Return the number of values in the shape for this object.
572 */
573 ESCRIPT_DLL_API
574 int
575 getNoValues() const
576 {
577 return m_data->getNoValues();
578 }
579
580
581 /**
582 \brief
583 dumps the object into a netCDF file
584 */
585 ESCRIPT_DLL_API
586 void
587 dump(const std::string fileName) const;
588
589 /**
590 \brief returns the values of the object as a list of tuples (one for each datapoint).
591
592 \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
593 If false, the result is a list of scalars [1, 2, ...]
594 */
595 ESCRIPT_DLL_API
596 const boost::python::object
597 toListOfTuples(bool scalarastuple=true);
598
599
600 /**
601 \brief
602 Return the sample data for the given sample no.
603 Please do not use this unless you NEED to access samples individually
604 \param sampleNo - Input - the given sample no.
605 \return pointer to the sample data.
606 */
607 ESCRIPT_DLL_API
608 inline
609 const DataTypes::real_t*
610 getSampleDataRO(DataTypes::RealVectorType::size_type sampleNo) const;
611
612
613 /**
614 \brief
615 Return the sample data for the given sample no.
616 Please do not use this unless you NEED to access samples individually
617 \param sampleNo - Input - the given sample no.
618 \return pointer to the sample data.
619 */
620 ESCRIPT_DLL_API
621 inline
622 DataTypes::real_t*
623 getSampleDataRW(DataTypes::RealVectorType::size_type sampleNo);
624
625
626 /**
627 \brief
628 Return a pointer to the beginning of the underlying data
629 \warning please avoid using this method since it by-passes possible lazy improvements. May be removed without notice.
630 \return pointer to the data.
631 */
632 ESCRIPT_DLL_API
633 const DataTypes::real_t*
634 getDataRO(DataTypes::real_t dummy=0) const;
635
636 ESCRIPT_DLL_API
637 const DataTypes::cplx_t*
638 getDataRO(DataTypes::cplx_t dummy) const;
639
640
641
642 /**
643 \brief
644 Return the sample data for the given tag. If an attempt is made to
645 access data that isn't tagged an exception will be thrown.
646 \param tag - Input - the tag key.
647 */
648 ESCRIPT_DLL_API
649 inline
650 DataTypes::real_t*
651 getSampleDataByTag(int tag)
652 {
653 return m_data->getSampleDataByTag(tag);
654 }
655
656 /**
657 \brief
658 Return a reference into the DataVector which points to the specified data point.
659 \param sampleNo - Input -
660 \param dataPointNo - Input -
661 */
662 ESCRIPT_DLL_API
663 DataTypes::RealVectorType::const_reference
664 getDataPointRO(int sampleNo, int dataPointNo);
665
666 /**
667 \brief
668 Return a reference into the DataVector which points to the specified data point.
669 \param sampleNo - Input -
670 \param dataPointNo - Input -
671 */
672 ESCRIPT_DLL_API
673 DataTypes::RealVectorType::reference
674 getDataPointRW(int sampleNo, int dataPointNo);
675
676
677
678 /**
679 \brief
680 Return the offset for the given sample and point within the sample
681 */
682 ESCRIPT_DLL_API
683 inline
684 DataTypes::RealVectorType::size_type
685 getDataOffset(int sampleNo,
686 int dataPointNo)
687 {
688 return m_data->getPointOffset(sampleNo,dataPointNo);
689 }
690
691 /**
692 \brief
693 Return a reference to the data point shape.
694 */
695 ESCRIPT_DLL_API
696 inline
697 const DataTypes::ShapeType&
698 getDataPointShape() const
699 {
700 return m_data->getShape();
701 }
702
703 /**
704 \brief
705 Return the data point shape as a tuple of integers.
706 */
707 ESCRIPT_DLL_API
708 const boost::python::tuple
709 getShapeTuple() const;
710
711 /**
712 \brief
713 Return the size of the data point. It is the product of the
714 data point shape dimensions.
715 */
716 ESCRIPT_DLL_API
717 int
718 getDataPointSize() const;
719
720 /**
721 \brief
722 Return the number of doubles stored for this Data.
723 */
724 ESCRIPT_DLL_API
725 DataTypes::RealVectorType::size_type
726 getLength() const;
727
728 /**
729 \brief Return true if this object contains no samples.
730 This is not the same as isEmpty()
731 */
732 ESCRIPT_DLL_API
733 bool
734 hasNoSamples() const
735 {
736 return getLength()==0;
737 }
738
739 /**
740 \brief
741 Assign the given value to the tag assocciated with name. Implicitly converts this
742 object to type DataTagged. Throws an exception if this object
743 cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
744 \param name - Input - name of tag.
745 \param value - Input - Value to associate with given key.
746 */
747 ESCRIPT_DLL_API
748 void
749 setTaggedValueByName(std::string name,
750 const boost::python::object& value);
751
752 /**
753 \brief
754 Assign the given value to the tag. Implicitly converts this
755 object to type DataTagged if it is constant.
756
757 \param tagKey - Input - Integer key.
758 \param value - Input - Value to associate with given key.
759 ==>*
760 */
761 ESCRIPT_DLL_API
762 void
763 setTaggedValue(int tagKey,
764 const boost::python::object& value);
765
766 /**
767 \brief
768 Assign the given value to the tag. Implicitly converts this
769 object to type DataTagged if it is constant.
770
771 \param tagKey - Input - Integer key.
772 \param pointshape - Input - The shape of the value parameter
773 \param value - Input - Value to associate with given key.
774 \param dataOffset - Input - Offset of the begining of the point within the value parameter
775 */
776 ESCRIPT_DLL_API
777 void
778 setTaggedValueFromCPP(int tagKey,
779 const DataTypes::ShapeType& pointshape,
780 const DataTypes::RealVectorType& value,
781 int dataOffset=0);
782
783
784 ESCRIPT_DLL_API
785 void
786 setTaggedValueFromCPP(int tagKey,
787 const DataTypes::ShapeType& pointshape,
788 const DataTypes::CplxVectorType& value,
789 int dataOffset=0);
790
791 /**
792 \brief
793 Copy other Data object into this Data object where mask is positive.
794 */
795 ESCRIPT_DLL_API
796 void
797 copyWithMask(const Data& other,
798 const Data& mask);
799
800 /**
801 Data object operation methods and operators.
802 */
803
804 /**
805 \brief
806 set all values to zero
807 *
808 */
809 ESCRIPT_DLL_API
810 void
811 setToZero();
812
813 /**
814 \brief
815 Interpolates this onto the given functionspace and returns
816 the result as a Data object.
817 *
818 */
819 ESCRIPT_DLL_API
820 Data
821 interpolate(const FunctionSpace& functionspace) const;
822
823 ESCRIPT_DLL_API
824 Data
825 interpolateFromTable3D(const WrappedArray& table, DataTypes::real_t Amin, DataTypes::real_t Astep,
826 DataTypes::real_t undef, Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep, Data& C,
827 DataTypes::real_t Cmin, DataTypes::real_t Cstep, bool check_boundaries);
828
829 ESCRIPT_DLL_API
830 Data
831 interpolateFromTable2D(const WrappedArray& table, DataTypes::real_t Amin, DataTypes::real_t Astep,
832 DataTypes::real_t undef, Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep,bool check_boundaries);
833
834 ESCRIPT_DLL_API
835 Data
836 interpolateFromTable1D(const WrappedArray& table, DataTypes::real_t Amin, DataTypes::real_t Astep,
837 DataTypes::real_t undef,bool check_boundaries);
838
839
840 ESCRIPT_DLL_API
841 Data
842 interpolateFromTable3DP(boost::python::object table, DataTypes::real_t Amin, DataTypes::real_t Astep,
843 Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep, Data& C, DataTypes::real_t Cmin, DataTypes::real_t Cstep, DataTypes::real_t undef,bool check_boundaries);
844
845
846 ESCRIPT_DLL_API
847 Data
848 interpolateFromTable2DP(boost::python::object table, DataTypes::real_t Amin, DataTypes::real_t Astep,
849 Data& B, DataTypes::real_t Bmin, DataTypes::real_t Bstep, DataTypes::real_t undef,bool check_boundaries);
850
851 ESCRIPT_DLL_API
852 Data
853 interpolateFromTable1DP(boost::python::object table, DataTypes::real_t Amin, DataTypes::real_t Astep,
854 DataTypes::real_t undef,bool check_boundaries);
855
856 ESCRIPT_DLL_API
857 Data
858 nonuniforminterp(boost::python::object in, boost::python::object out, bool check_boundaries);
859
860 ESCRIPT_DLL_API
861 Data
862 nonuniformslope(boost::python::object in, boost::python::object out, bool check_boundaries);
863
864 /**
865 \brief
866 Calculates the gradient of the data at the data points of functionspace.
867 If functionspace is not present the function space of Function(getDomain()) is used.
868 *
869 */
870 ESCRIPT_DLL_API
871 Data
872 gradOn(const FunctionSpace& functionspace) const;
873
874 ESCRIPT_DLL_API
875 Data
876 grad() const;
877
878 /**
879 \brief
880 Calculate the integral over the function space domain as a python tuple.
881 */
882 ESCRIPT_DLL_API
883 boost::python::object
884 integrateToTuple_const() const;
885
886
887 /**
888 \brief
889 Calculate the integral over the function space domain as a python tuple.
890 */
891 ESCRIPT_DLL_API
892 boost::python::object
893 integrateToTuple();
894
895
896
897 /**
898 \brief
899 Returns 1./ Data object
900 *
901 */
902 ESCRIPT_DLL_API
903 Data
904 oneOver() const;
905 /**
906 \brief
907 Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.
908 *
909 */
910 ESCRIPT_DLL_API
911 Data
912 wherePositive() const;
913
914 /**
915 \brief
916 Return a Data with a 1 for -ive values and a 0 for +ive or 0 values.
917 *
918 */
919 ESCRIPT_DLL_API
920 Data
921 whereNegative() const;
922
923 /**
924 \brief
925 Return a Data with a 1 for +ive or 0 values and a 0 for -ive values.
926 *
927 */
928 ESCRIPT_DLL_API
929 Data
930 whereNonNegative() const;
931
932 /**
933 \brief
934 Return a Data with a 1 for -ive or 0 values and a 0 for +ive values.
935 *
936 */
937 ESCRIPT_DLL_API
938 Data
939 whereNonPositive() const;
940
941 /**
942 \brief
943 Return a Data with a 1 for 0 values and a 0 for +ive or -ive values.
944 *
945 */
946 ESCRIPT_DLL_API
947 Data
948 whereZero(DataTypes::real_t tol=0.0) const;
949
950 /**
951 \brief
952 Return a Data with a 0 for 0 values and a 1 for +ive or -ive values.
953 *
954 */
955 ESCRIPT_DLL_API
956 Data
957 whereNonZero(DataTypes::real_t tol=0.0) const;
958
959 /**
960 \brief
961 Return the maximum absolute value of this Data object.
962
963 The method is not const because lazy data needs to be expanded before Lsup can be computed.
964 The _const form can be used when the Data object is const, however this will only work for
965 Data which is not Lazy.
966
967 For Data which contain no samples (or tagged Data for which no tags in use have a value)
968 zero is returned.
969 */
970 ESCRIPT_DLL_API
971 DataTypes::real_t
972 Lsup();
973
974 ESCRIPT_DLL_API
975 DataTypes::real_t
976 Lsup_const() const;
977
978
979 /**
980 \brief
981 Return the maximum value of this Data object.
982
983 The method is not const because lazy data needs to be expanded before sup can be computed.
984 The _const form can be used when the Data object is const, however this will only work for
985 Data which is not Lazy.
986
987 For Data which contain no samples (or tagged Data for which no tags in use have a value)
988 a large negative value is returned.
989 */
990 ESCRIPT_DLL_API
991 DataTypes::real_t
992 sup();
993
994 ESCRIPT_DLL_API
995 DataTypes::real_t
996 sup_const() const;
997
998
999 /**
1000 \brief
1001 Return the minimum value of this Data object.
1002
1003 The method is not const because lazy data needs to be expanded before inf can be computed.
1004 The _const form can be used when the Data object is const, however this will only work for
1005 Data which is not Lazy.
1006
1007 For Data which contain no samples (or tagged Data for which no tags in use have a value)
1008 a large positive value is returned.
1009 */
1010 ESCRIPT_DLL_API
1011 DataTypes::real_t
1012 inf();
1013
1014 ESCRIPT_DLL_API
1015 DataTypes::real_t
1016 inf_const() const;
1017
1018
1019
1020 /**
1021 \brief
1022 Return the absolute value of each data point of this Data object.
1023 *
1024 */
1025 ESCRIPT_DLL_API
1026 Data
1027 abs() const;
1028
1029 /**
1030 \brief
1031 Return the maximum value of each data point of this Data object.
1032 *
1033 */
1034 ESCRIPT_DLL_API
1035 Data
1036 maxval() const;
1037
1038 /**
1039 \brief
1040 Return the minimum value of each data point of this Data object.
1041 *
1042 */
1043 ESCRIPT_DLL_API
1044 Data
1045 minval() const;
1046
1047 /**
1048 \brief
1049 Return the (sample number, data-point number) of the data point with
1050 the minimum component value in this Data object.
1051 \note If you are working in python, please consider using Locator
1052 instead of manually manipulating process and point IDs.
1053 */
1054 ESCRIPT_DLL_API
1055 const boost::python::tuple
1056 minGlobalDataPoint() const;
1057
1058 /**
1059 \brief
1060 Return the (sample number, data-point number) of the data point with
1061 the minimum component value in this Data object.
1062 \note If you are working in python, please consider using Locator
1063 instead of manually manipulating process and point IDs.
1064 */
1065 ESCRIPT_DLL_API
1066 const boost::python::tuple
1067 maxGlobalDataPoint() const;
1068
1069
1070
1071 /**
1072 \brief
1073 Return the sign of each data point of this Data object.
1074 -1 for negative values, zero for zero values, 1 for positive values.
1075 *
1076 */
1077 ESCRIPT_DLL_API
1078 Data
1079 sign() const;
1080
1081 /**
1082 \brief
1083 Return the symmetric part of a matrix which is half the matrix plus its transpose.
1084 *
1085 */
1086 ESCRIPT_DLL_API
1087 Data
1088 symmetric() const;
1089
1090 /**
1091 \brief
1092 Return the nonsymmetric part of a matrix which is half the matrix minus its transpose.
1093 *
1094 */
1095 ESCRIPT_DLL_API
1096 Data
1097 nonsymmetric() const;
1098
1099 /**
1100 \brief
1101 Return the trace of a matrix
1102 *
1103 */
1104 ESCRIPT_DLL_API
1105 Data
1106 trace(int axis_offset) const;
1107
1108 /**
1109 \brief
1110 Transpose each data point of this Data object around the given axis.
1111 *
1112 */
1113 ESCRIPT_DLL_API
1114 Data
1115 transpose(int axis_offset) const;
1116
1117 /**
1118 \brief
1119 Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
1120 Currently this function is restricted to rank 2, square shape, and dimension 3.
1121 *
1122 */
1123 ESCRIPT_DLL_API
1124 Data
1125 eigenvalues() const;
1126
1127 /**
1128 \brief
1129 Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1130 the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1131 tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1132 first non-zero entry is positive.
1133 Currently this function is restricted to rank 2, square shape, and dimension 3
1134 *
1135 */
1136 ESCRIPT_DLL_API
1137 const boost::python::tuple
1138 eigenvalues_and_eigenvectors(const DataTypes::real_t tol=1.e-12) const;
1139
1140 /**
1141 \brief
1142 swaps the components axis0 and axis1
1143 *
1144 */
1145 ESCRIPT_DLL_API
1146 Data
1147 swapaxes(const int axis0, const int axis1) const;
1148
1149 /**
1150 \brief
1151 Return the error function erf of each data point of this Data object.
1152 *
1153 */
1154 ESCRIPT_DLL_API
1155 Data
1156 erf() const;
1157
1158
1159 /**
1160 \brief
1161 For complex values return the conjugate values.
1162 For non-complex data return a copy
1163 */
1164 ESCRIPT_DLL_API
1165 Data
1166 conjugate() const;
1167
1168 ESCRIPT_DLL_API
1169 Data
1170 real() const;
1171
1172 ESCRIPT_DLL_API
1173 Data
1174 imag() const;
1175
1176 /**
1177 \brief
1178 Return the sin of each data point of this Data object.
1179 *
1180 */
1181 ESCRIPT_DLL_API
1182 Data
1183 sin() const;
1184
1185 /**
1186 \brief
1187 Return the cos of each data point of this Data object.
1188 *
1189 */
1190 ESCRIPT_DLL_API
1191 Data
1192 cos() const;
1193
1194 /**
1195 \brief
1196 Bessel worker function.
1197 *
1198 */
1199 ESCRIPT_DLL_API
1200 Data
1201 bessel(int order, DataTypes::real_t (*besselfunc) (int,DataTypes::real_t) );
1202
1203
1204 /**
1205 \brief
1206 Return the Bessel function of the first kind for each data point of this Data object.
1207 *
1208 */
1209 ESCRIPT_DLL_API
1210 Data
1211 besselFirstKind(int order);
1212
1213 /**
1214 \brief
1215 Return the Bessel function of the second kind for each data point of this Data object.
1216 *
1217 */
1218 ESCRIPT_DLL_API
1219 Data
1220 besselSecondKind(int order);
1221
1222
1223 /**
1224 \brief
1225 Return the tan of each data point of this Data object.
1226 *
1227 */
1228 ESCRIPT_DLL_API
1229 Data
1230 tan() const;
1231
1232 /**
1233 \brief
1234 Return the asin of each data point of this Data object.
1235 *
1236 */
1237 ESCRIPT_DLL_API
1238 Data
1239 asin() const;
1240
1241 /**
1242 \brief
1243 Return the acos of each data point of this Data object.
1244 *
1245 */
1246 ESCRIPT_DLL_API
1247 Data
1248 acos() const;
1249
1250 /**
1251 \brief
1252 Return the atan of each data point of this Data object.
1253 *
1254 */
1255 ESCRIPT_DLL_API
1256 Data
1257 atan() const;
1258
1259 /**
1260 \brief
1261 Return the sinh of each data point of this Data object.
1262 *
1263 */
1264 ESCRIPT_DLL_API
1265 Data
1266 sinh() const;
1267
1268 /**
1269 \brief
1270 Return the cosh of each data point of this Data object.
1271 *
1272 */
1273 ESCRIPT_DLL_API
1274 Data
1275 cosh() const;
1276
1277 /**
1278 \brief
1279 Return the tanh of each data point of this Data object.
1280 *
1281 */
1282 ESCRIPT_DLL_API
1283 Data
1284 tanh() const;
1285
1286 /**
1287 \brief
1288 Return the asinh of each data point of this Data object.
1289 *
1290 */
1291 ESCRIPT_DLL_API
1292 Data
1293 asinh() const;
1294
1295 /**
1296 \brief
1297 Return the acosh of each data point of this Data object.
1298 *
1299 */
1300 ESCRIPT_DLL_API
1301 Data
1302 acosh() const;
1303
1304 /**
1305 \brief
1306 Return the atanh of each data point of this Data object.
1307 *
1308 */
1309 ESCRIPT_DLL_API
1310 Data
1311 atanh() const;
1312
1313 /**
1314 \brief
1315 Return the log to base 10 of each data point of this Data object.
1316 *
1317 */
1318 ESCRIPT_DLL_API
1319 Data
1320 log10() const;
1321
1322 /**
1323 \brief
1324 Return the natural log of each data point of this Data object.
1325 *
1326 */
1327 ESCRIPT_DLL_API
1328 Data
1329 log() const;
1330
1331 /**
1332 \brief
1333 Return the exponential function of each data point of this Data object.
1334 *
1335 */
1336 ESCRIPT_DLL_API
1337 Data
1338 exp() const;
1339
1340 /**
1341 \brief
1342 Return the square root of each data point of this Data object.
1343 *
1344 */
1345 ESCRIPT_DLL_API
1346 Data
1347 sqrt() const;
1348
1349 /**
1350 \brief
1351 Return the negation of each data point of this Data object.
1352 *
1353 */
1354 ESCRIPT_DLL_API
1355 Data
1356 neg() const;
1357
1358 /**
1359 \brief
1360 Return the identity of each data point of this Data object.
1361 Simply returns this object unmodified.
1362 *
1363 */
1364 ESCRIPT_DLL_API
1365 Data
1366 pos() const;
1367
1368 /**
1369 \brief
1370 Return the given power of each data point of this Data object.
1371
1372 \param right Input - the power to raise the object to.
1373 *
1374 */
1375 ESCRIPT_DLL_API
1376 Data
1377 powD(const Data& right) const;
1378
1379 /**
1380 \brief
1381 Return the given power of each data point of this boost python object.
1382
1383 \param right Input - the power to raise the object to.
1384 *
1385 */
1386 ESCRIPT_DLL_API
1387 Data
1388 powO(const boost::python::object& right) const;
1389
1390 /**
1391 \brief
1392 Return the given power of each data point of this boost python object.
1393
1394 \param left Input - the bases
1395 *
1396 */
1397
1398 ESCRIPT_DLL_API
1399 Data
1400 rpowO(const boost::python::object& left) const;
1401
1402 /**
1403 \brief
1404 Overloaded operator +=
1405 \param right - Input - The right hand side.
1406 *
1407 */
1408 ESCRIPT_DLL_API
1409 Data& operator+=(const Data& right);
1410 ESCRIPT_DLL_API
1411 Data& operator+=(const boost::python::object& right);
1412
1413 ESCRIPT_DLL_API
1414 Data& operator=(const Data& other);
1415
1416 /**
1417 \brief
1418 Overloaded operator -=
1419 \param right - Input - The right hand side.
1420 *
1421 */
1422 ESCRIPT_DLL_API
1423 Data& operator-=(const Data& right);
1424 ESCRIPT_DLL_API
1425 Data& operator-=(const boost::python::object& right);
1426
1427 /**
1428 \brief
1429 Overloaded operator *=
1430 \param right - Input - The right hand side.
1431 *
1432 */
1433 ESCRIPT_DLL_API
1434 Data& operator*=(const Data& right);
1435 ESCRIPT_DLL_API
1436 Data& operator*=(const boost::python::object& right);
1437
1438 /**
1439 \brief
1440 Overloaded operator /=
1441 \param right - Input - The right hand side.
1442 *
1443 */
1444 ESCRIPT_DLL_API
1445 Data& operator/=(const Data& right);
1446 ESCRIPT_DLL_API
1447 Data& operator/=(const boost::python::object& right);
1448
1449 /**
1450 \brief
1451 Newer style division operator for python
1452 */
1453 ESCRIPT_DLL_API
1454 Data truedivD(const Data& right);
1455
1456 /**
1457 \brief
1458 Newer style division operator for python
1459 */
1460 ESCRIPT_DLL_API
1461 Data truedivO(const boost::python::object& right);
1462
1463 /**
1464 \brief
1465 Newer style division operator for python
1466 */
1467 ESCRIPT_DLL_API
1468 Data rtruedivO(const boost::python::object& left);
1469
1470 /**
1471 \brief
1472 wrapper for python add operation
1473 */
1474 ESCRIPT_DLL_API
1475 boost::python::object __add__(const boost::python::object& right);
1476
1477
1478 /**
1479 \brief
1480 wrapper for python subtract operation
1481 */
1482 ESCRIPT_DLL_API
1483 boost::python::object __sub__(const boost::python::object& right);
1484
1485 /**
1486 \brief
1487 wrapper for python reverse subtract operation
1488 */
1489 ESCRIPT_DLL_API
1490 boost::python::object __rsub__(const boost::python::object& right);
1491
1492 /**
1493 \brief
1494 wrapper for python multiply operation
1495 */
1496 ESCRIPT_DLL_API
1497 boost::python::object __mul__(const boost::python::object& right);
1498
1499 /**
1500 \brief
1501 wrapper for python divide operation
1502 */
1503 ESCRIPT_DLL_API
1504 boost::python::object __div__(const boost::python::object& right);
1505
1506 /**
1507 \brief
1508 wrapper for python reverse divide operation
1509 */
1510 ESCRIPT_DLL_API
1511 boost::python::object __rdiv__(const boost::python::object& right);
1512
1513 /**
1514 \brief return inverse of matricies.
1515 */
1516 ESCRIPT_DLL_API
1517 Data
1518 matrixInverse() const;
1519
1520 /**
1521 \brief
1522 Returns true if this can be interpolated to functionspace.
1523 */
1524 ESCRIPT_DLL_API
1525 bool
1526 probeInterpolation(const FunctionSpace& functionspace) const;
1527
1528 /**
1529 Data object slicing methods.
1530 */
1531
1532 /**
1533 \brief
1534 Returns a slice from this Data object.
1535
1536 /description
1537 Implements the [] get operator in python.
1538 Calls getSlice.
1539
1540 \param key - Input - python slice tuple specifying
1541 slice to return.
1542 */
1543 ESCRIPT_DLL_API
1544 Data
1545 getItem(const boost::python::object& key) const;
1546
1547 /**
1548 \brief
1549 Copies slice from value into this Data object.
1550
1551 Implements the [] set operator in python.
1552 Calls setSlice.
1553
1554 \param key - Input - python slice tuple specifying
1555 slice to copy from value.
1556 \param value - Input - Data object to copy from.
1557 */
1558 ESCRIPT_DLL_API
1559 void
1560 setItemD(const boost::python::object& key,
1561 const Data& value);
1562
1563 ESCRIPT_DLL_API
1564 void
1565 setItemO(const boost::python::object& key,
1566 const boost::python::object& value);
1567
1568 // These following public methods should be treated as private.
1569
1570 /**
1571 \brief
1572 Perform the given unary operation on every element of every data point in
1573 this Data object.
1574 */
1575 template <class UnaryFunction>
1576 ESCRIPT_DLL_API
1577 inline
1578 void
1579 unaryOp2(UnaryFunction operation);
1580
1581 /**
1582 \brief
1583 Return a Data object containing the specified slice of
1584 this Data object.
1585 \param region - Input - Region to copy.
1586 *
1587 */
1588 ESCRIPT_DLL_API
1589 Data
1590 getSlice(const DataTypes::RegionType& region) const;
1591
1592 /**
1593 \brief
1594 Copy the specified slice from the given value into this
1595 Data object.
1596 \param value - Input - Data to copy from.
1597 \param region - Input - Region to copy.
1598 *
1599 */
1600 ESCRIPT_DLL_API
1601 void
1602 setSlice(const Data& value,
1603 const DataTypes::RegionType& region);
1604
1605 /**
1606 \brief
1607 print the data values to stdout. Used for debugging
1608 */
1609 ESCRIPT_DLL_API
1610 void
1611 print(void);
1612
1613 /**
1614 \brief
1615 return the MPI rank number of the local data
1616 MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1617 is returned
1618 */
1619 ESCRIPT_DLL_API
1620 int
1621 get_MPIRank(void) const;
1622
1623 /**
1624 \brief
1625 return the MPI rank number of the local data
1626 MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1627 is returned
1628 */
1629 ESCRIPT_DLL_API
1630 int
1631 get_MPISize(void) const;
1632
1633 /**
1634 \brief
1635 return the MPI rank number of the local data
1636 MPI_COMM_WORLD is assumed and returned.
1637 */
1638 ESCRIPT_DLL_API
1639 MPI_Comm
1640 get_MPIComm(void) const;
1641
1642 /**
1643 \brief
1644 return the object produced by the factory, which is a DataConstant or DataExpanded
1645 TODO Ownership of this object should be explained in doco.
1646 */
1647 ESCRIPT_DLL_API
1648 DataAbstract*
1649 borrowData(void) const;
1650
1651 ESCRIPT_DLL_API
1652 DataAbstract_ptr
1653 borrowDataPtr(void) const;
1654
1655 ESCRIPT_DLL_API
1656 DataReady_ptr
1657 borrowReadyPtr(void) const;
1658
1659
1660
1661 /**
1662 \brief
1663 Return a pointer to the beginning of the datapoint at the specified offset.
1664 TODO Eventually these should be inlined.
1665 \param i - position(offset) in the underlying datastructure
1666 */
1667
1668 ESCRIPT_DLL_API
1669 DataTypes::RealVectorType::const_reference
1670 getDataAtOffsetRO(DataTypes::RealVectorType::size_type i, DataTypes::real_t dummy=0);
1671
1672
1673 ESCRIPT_DLL_API
1674 DataTypes::RealVectorType::reference
1675 getDataAtOffsetRW(DataTypes::RealVectorType::size_type i, DataTypes::real_t dummy=0);
1676
1677 ESCRIPT_DLL_API
1678 DataTypes::CplxVectorType::const_reference
1679 getDataAtOffsetRO(DataTypes::CplxVectorType::size_type i, DataTypes::cplx_t dummy);
1680
1681
1682 ESCRIPT_DLL_API
1683 DataTypes::CplxVectorType::reference
1684 getDataAtOffsetRW(DataTypes::CplxVectorType::size_type i, DataTypes::cplx_t dummy);
1685
1686
1687 /**
1688 \brief Ensures that the Data is expanded and returns its underlying vector
1689 Does not check for exclusive write so do that before calling if sharing
1690 Is a posibility.
1691 \warning For domain implementors only. Using this function will
1692 avoid using optimisations like lazy evaluation. It is intended
1693 to allow quick initialisation of Data by domain; not as a bypass around
1694 escript's other mechanisms.
1695 */
1696 ESCRIPT_DLL_API
1697 DataTypes::RealVectorType&
1698 getExpandedVectorReference(DataTypes::real_t dummy=0);
1699
1700 ESCRIPT_DLL_API
1701 DataTypes::CplxVectorType&
1702 getExpandedVectorReference(DataTypes::cplx_t dummy);
1703
1704
1705 /**
1706 * \brief For tagged Data returns the number of tags with values.
1707 * For non-tagged data will return 0 (even Data which has been expanded from tagged).
1708 */
1709 ESCRIPT_DLL_API
1710 size_t
1711 getNumberOfTaggedValues() const;
1712
1713 /*
1714 * \brief make the data complex
1715 */
1716 ESCRIPT_DLL_API
1717 void complicate();
1718
1719 protected:
1720
1721 private:
1722
1723 template <class BinaryOp>
1724 DataTypes::real_t
1725 #ifdef ESYS_MPI
1726 lazyAlgWorker(DataTypes::real_t init, MPI_Op mpiop_type);
1727 #else
1728 lazyAlgWorker(DataTypes::real_t init);
1729 #endif
1730
1731 DataTypes::real_t
1732 LsupWorker() const;
1733
1734 DataTypes::real_t
1735 supWorker() const;
1736
1737 DataTypes::real_t
1738 infWorker() const;
1739
1740 boost::python::object
1741 integrateWorker() const;
1742
1743 void
1744 calc_minGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1745
1746 void
1747 calc_maxGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1748
1749 // For internal use in Data.cpp only!
1750 // other uses should call the main entry points and allow laziness
1751 Data
1752 minval_nonlazy() const;
1753
1754 // For internal use in Data.cpp only!
1755 Data
1756 maxval_nonlazy() const;
1757
1758
1759 /**
1760 \brief
1761 Check *this and the right operand are compatible. Throws
1762 an exception if they aren't.
1763 \param right - Input - The right hand side.
1764 */
1765 inline
1766 void
1767 operandCheck(const Data& right) const
1768 {
1769 return m_data->operandCheck(*(right.m_data.get()));
1770 }
1771
1772 /**
1773 \brief
1774 Perform the specified reduction algorithm on every element of every data point in
1775 this Data object according to the given function and return the single value result.
1776 */
1777 template <class BinaryFunction>
1778 inline
1779 DataTypes::real_t
1780 algorithm(BinaryFunction operation,
1781 DataTypes::real_t initial_value) const;
1782
1783 /**
1784 \brief
1785 Reduce each data-point in this Data object using the given operation. Return a Data
1786 object with the same number of data-points, but with each data-point containing only
1787 one value - the result of the reduction operation on the corresponding data-point in
1788 this Data object
1789 */
1790 template <class BinaryFunction>
1791 inline
1792 Data
1793 dp_algorithm(BinaryFunction operation,
1794 DataTypes::real_t initial_value) const;
1795
1796 /**
1797 \brief
1798 Perform the given binary operation on all of the data's elements.
1799 The underlying type of the right hand side (right) determines the final
1800 type of *this after the operation. For example if the right hand side
1801 is expanded *this will be expanded if necessary.
1802 RHS is a Data object.
1803 */
1804 template <class BinaryFunction>
1805 inline
1806 void
1807 binaryOp(const Data& right,
1808 BinaryFunction operation);
1809
1810 void
1811 binaryDataOp(const Data& right,
1812 escript::ESFunction operation);
1813
1814
1815 /**
1816 \brief
1817 Convert the data type of the RHS to match this.
1818 \param right - Input - data type to match.
1819 */
1820 void
1821 typeMatchLeft(Data& right) const;
1822
1823 /**
1824 \brief
1825 Convert the data type of this to match the RHS.
1826 \param right - Input - data type to match.
1827 */
1828 void
1829 typeMatchRight(const Data& right);
1830
1831 /**
1832 \brief
1833 Construct a Data object of the appropriate type.
1834 */
1835
1836 void
1837 initialise(const DataTypes::RealVectorType& value,
1838 const DataTypes::ShapeType& shape,
1839 const FunctionSpace& what,
1840 bool expanded);
1841
1842 void
1843 initialise(const WrappedArray& value,
1844 const FunctionSpace& what,
1845 bool expanded);
1846
1847 void
1848 initialise(const DataTypes::real_t value,
1849 const DataTypes::ShapeType& shape,
1850 const FunctionSpace& what,
1851 bool expanded);
1852
1853 //
1854 // flag to protect the data object against any update
1855 bool m_protected;
1856 mutable bool m_shared;
1857 bool m_lazy;
1858
1859 //
1860 // pointer to the actual data object
1861 // boost::shared_ptr<DataAbstract> m_data;
1862 DataAbstract_ptr m_data;
1863
1864 // If possible please use getReadyPtr instead.
1865 // But see warning below.
1866 const DataReady*
1867 getReady() const;
1868
1869 DataReady*
1870 getReady();
1871
1872
1873 // Be wary of using this for local operations since it (temporarily) increases reference count.
1874 // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1875 // getReady() instead
1876 DataReady_ptr
1877 getReadyPtr();
1878
1879 const_DataReady_ptr
1880 getReadyPtr() const;
1881
1882
1883 /**
1884 \brief Update the Data's shared flag
1885 This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1886 For internal use only.
1887 */
1888 void updateShareStatus(bool nowshared) const
1889 {
1890 m_shared=nowshared; // m_shared is mutable
1891 }
1892
1893 // In the isShared() method below:
1894 // A problem would occur if m_data (the address pointed to) were being modified
1895 // while the call m_data->is_shared is being executed.
1896 //
1897 // Q: So why do I think this code can be thread safe/correct?
1898 // A: We need to make some assumptions.
1899 // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1900 // 2. We assume that no constructions or assignments which will share previously unshared
1901 // will occur while this call is executing. This is consistent with the way Data:: and C are written.
1902 //
1903 // This means that the only transition we need to consider, is when a previously shared object is
1904 // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1905 // In those cases the m_shared flag changes to false after m_data has completed changing.
1906 // For any threads executing before the flag switches they will assume the object is still shared.
1907 bool isShared() const
1908 {
1909 return m_shared;
1910 /* if (m_shared) return true;
1911 if (m_data->isShared())
1912 {
1913 updateShareStatus(true);
1914 return true;
1915 }
1916 return false;*/
1917 }
1918
1919 void forceResolve()
1920 {
1921 if (isLazy())
1922 {
1923 #ifdef _OPENMP
1924 if (omp_in_parallel())
1925 { // Yes this is throwing an exception out of an omp thread which is forbidden.
1926 throw DataException("Please do not call forceResolve() in a parallel region.");
1927 }
1928 #endif
1929 resolve();
1930 }
1931 }
1932
1933 /**
1934 \brief if another object is sharing out member data make a copy to work with instead.
1935 This code should only be called from single threaded sections of code.
1936 */
1937 void exclusiveWrite()
1938 {
1939 #ifdef _OPENMP
1940 if (omp_in_parallel())
1941 {
1942 // *((int*)0)=17;
1943 throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1944 }
1945 #endif
1946 forceResolve();
1947 if (isShared())
1948 {
1949 DataAbstract* t=m_data->deepCopy();
1950 set_m_data(DataAbstract_ptr(t));
1951 }
1952 #ifdef EXWRITECHK
1953 m_data->exclusivewritecalled=true;
1954 #endif
1955 }
1956
1957 /**
1958 \brief checks if caller can have exclusive write to the object
1959 */
1960 void checkExclusiveWrite()
1961 {
1962 if (isLazy() || isShared())
1963 {
1964 throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1965 }
1966 }
1967
1968 /**
1969 \brief Modify the data abstract hosted by this Data object
1970 For internal use only.
1971 Passing a pointer to null is permitted (do this in the destructor)
1972 \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1973 */
1974 void set_m_data(DataAbstract_ptr p);
1975
1976 friend class DataAbstract; // To allow calls to updateShareStatus
1977 friend class TestDomain; // so its getX will work quickly
1978 #ifdef IKNOWWHATIMDOING
1979 friend ESCRIPT_DLL_API Data applyBinaryCFunction(boost::python::object cfunc, boost::python::tuple shape, escript::Data& d, escript::Data& e);
1980 #endif
1981 friend ESCRIPT_DLL_API Data condEval(escript::Data& mask, escript::Data& trueval, escript::Data& falseval);
1982 friend ESCRIPT_DLL_API Data randomData(const boost::python::tuple& shape, const FunctionSpace& what, long seed, const boost::python::tuple& filter);
1983
1984 };
1985
1986
1987 #ifdef IKNOWWHATIMDOING
1988 ESCRIPT_DLL_API
1989 Data
1990 applyBinaryCFunction(boost::python::object func, boost::python::tuple shape, escript::Data& d, escript::Data& e);
1991 #endif
1992
1993
1994 ESCRIPT_DLL_API
1995 Data
1996 condEval(escript::Data& mask, escript::Data& trueval, escript::Data& falseval);
1997
1998
1999
2000 /**
2001 \brief Create a new Expanded Data object filled with pseudo-random data.
2002 */
2003 ESCRIPT_DLL_API
2004 Data randomData(const boost::python::tuple& shape,
2005 const FunctionSpace& what,
2006 long seed, const boost::python::tuple& filter);
2007
2008
2009 } // end namespace escript
2010
2011
2012 // No, this is not supposed to be at the top of the file
2013 // DataAbstact needs to be declared first, then DataReady needs to be fully declared
2014 // so that I can dynamic cast between them below.
2015 #include "DataReady.h"
2016 #include "DataLazy.h"
2017
2018 namespace escript
2019 {
2020
2021 inline
2022 const DataReady*
2023 Data::getReady() const
2024 {
2025 const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
2026 EsysAssert((dr!=0), "Error - casting to DataReady.");
2027 return dr;
2028 }
2029
2030 inline
2031 DataReady*
2032 Data::getReady()
2033 {
2034 DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
2035 EsysAssert((dr!=0), "Error - casting to DataReady.");
2036 return dr;
2037 }
2038
2039 // Be wary of using this for local operations since it (temporarily) increases reference count.
2040 // If you are just using this to call a method on DataReady instead of DataAbstract consider using
2041 // getReady() instead
2042 inline
2043 DataReady_ptr
2044 Data::getReadyPtr()
2045 {
2046 DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
2047 EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
2048 return dr;
2049 }
2050
2051
2052 inline
2053 const_DataReady_ptr
2054 Data::getReadyPtr() const
2055 {
2056 const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
2057 EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
2058 return dr;
2059 }
2060
2061 inline
2062 DataTypes::real_t*
2063 Data::getSampleDataRW(DataTypes::RealVectorType::size_type sampleNo)
2064 {
2065 if (isLazy())
2066 {
2067 throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
2068 }
2069 #ifdef EXWRITECHK
2070 if (!getReady()->exclusivewritecalled)
2071 {
2072 throw DataException("Error, call to Data::getSampleDataRW without a preceeding call to requireWrite/exclusiveWrite.");
2073 }
2074 #endif
2075 return getReady()->getSampleDataRW(sampleNo);
2076 }
2077
2078 inline
2079 const DataTypes::real_t*
2080 Data::getSampleDataRO(DataTypes::RealVectorType::size_type sampleNo) const
2081 {
2082 DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
2083 if (l!=0)
2084 {
2085 size_t offset=0;
2086 const DataTypes::RealVectorType* res=l->resolveSample(sampleNo,offset);
2087 return &((*res)[offset]);
2088 }
2089 return getReady()->getSampleDataRO(sampleNo);
2090 }
2091
2092 inline
2093 const DataTypes::real_t*
2094 Data::getDataRO(DataTypes::real_t dummy) const
2095 {
2096 if (isLazy())
2097 {
2098 throw DataException("Programmer error - getDataRO must not be called on Lazy Data.");
2099 }
2100 if (getNumSamples()==0)
2101 {
2102 return 0;
2103 }
2104 else
2105 {
2106 return &(getReady()->getTypedVectorRO(0)[0]);
2107 }
2108 }
2109
2110 inline
2111 const DataTypes::cplx_t*
2112 Data::getDataRO(DataTypes::cplx_t dummy) const
2113 {
2114 if (isLazy())
2115 {
2116 throw DataException("Programmer error - getDataRO must not be called on Lazy Data.");
2117 }
2118 if (getNumSamples()==0)
2119 {
2120 return 0;
2121 }
2122 else
2123 {
2124 return &(getReady()->getTypedVectorRO(dummy)[0]);
2125 }
2126 }
2127
2128
2129 /**
2130 Binary Data object operators.
2131 */
2132 inline DataTypes::real_t rpow(DataTypes::real_t x,DataTypes::real_t y)
2133 {
2134 return pow(y,x);
2135 }
2136
2137 /**
2138 \brief
2139 Operator+
2140 Takes two Data objects.
2141 */
2142 ESCRIPT_DLL_API Data operator+(const Data& left, const Data& right);
2143
2144 /**
2145 \brief
2146 Operator-
2147 Takes two Data objects.
2148 */
2149 ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);
2150
2151 /**
2152 \brief
2153 Operator*
2154 Takes two Data objects.
2155 */
2156 ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);
2157
2158 /**
2159 \brief
2160 Operator/
2161 Takes two Data objects.
2162 */
2163 ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);
2164
2165 /**
2166 \brief
2167 Operator+
2168 Takes LHS Data object and RHS python::object.
2169 python::object must be convertable to Data type.
2170 */
2171 ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);
2172
2173 /**
2174 \brief
2175 Operator-
2176 Takes LHS Data object and RHS python::object.
2177 python::object must be convertable to Data type.
2178 */
2179 ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);
2180
2181 /**
2182 \brief
2183 Operator*
2184 Takes LHS Data object and RHS python::object.
2185 python::object must be convertable to Data type.
2186 */
2187 ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);
2188
2189 /**
2190 \brief
2191 Operator/
2192 Takes LHS Data object and RHS python::object.
2193 python::object must be convertable to Data type.
2194 */
2195 ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);
2196
2197 /**
2198 \brief
2199 Operator+
2200 Takes LHS python::object and RHS Data object.
2201 python::object must be convertable to Data type.
2202 */
2203 ESCRIPT_DLL_API Data operator+(const boost::python::object& left, const Data& right);
2204
2205 /**
2206 \brief
2207 Operator-
2208 Takes LHS python::object and RHS Data object.
2209 python::object must be convertable to Data type.
2210 */
2211 ESCRIPT_DLL_API Data operator-(const boost::python::object& left, const Data& right);
2212
2213 /**
2214 \brief
2215 Operator*
2216 Takes LHS python::object and RHS Data object.
2217 python::object must be convertable to Data type.
2218 */
2219 ESCRIPT_DLL_API Data operator*(const boost::python::object& left, const Data& right);
2220
2221 /**
2222 \brief
2223 Operator/
2224 Takes LHS python::object and RHS Data object.
2225 python::object must be convertable to Data type.
2226 */
2227 ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
2228
2229
2230
2231 /**
2232 \brief
2233 Output operator
2234 */
2235 ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);
2236
2237 /**
2238 \brief
2239 Compute a tensor product of two Data objects
2240 \param arg_0 - Input - Data object
2241 \param arg_1 - Input - Data object
2242 \param axis_offset - Input - axis offset
2243 \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
2244 */
2245 ESCRIPT_DLL_API
2246 Data
2247 C_GeneralTensorProduct(Data& arg_0,
2248 Data& arg_1,
2249 int axis_offset=0,
2250 int transpose=0);
2251
2252 /**
2253 \brief
2254 Operator/
2255 Takes RHS Data object.
2256 */
2257 inline
2258 Data
2259 Data::truedivD(const Data& right)
2260 {
2261 return *this / right;
2262 }
2263
2264 /**
2265 \brief
2266 Operator/
2267 Takes RHS python::object.
2268 */
2269 inline
2270 Data
2271 Data::truedivO(const boost::python::object& right)
2272 {
2273 Data tmp(right, getFunctionSpace(), false);
2274 return truedivD(tmp);
2275 }
2276
2277 /**
2278 \brief
2279 Operator/
2280 Takes LHS python::object.
2281 */
2282 inline
2283 Data
2284 Data::rtruedivO(const boost::python::object& left)
2285 {
2286 Data tmp(left, getFunctionSpace(), false);
2287 return tmp.truedivD(*this);
2288 }
2289
2290
2291 inline
2292 void
2293 Data::binaryDataOp(const Data& right,
2294 escript::ESFunction operation)
2295 {
2296 //
2297 // if this has a rank of zero promote it to the rank of the RHS
2298 if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
2299 throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
2300 }
2301
2302 if (isLazy() || right.isLazy())
2303 {
2304 throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
2305 }
2306 //
2307 // initially make the temporary a shallow copy
2308 Data tempRight(right);
2309 FunctionSpace fsl=getFunctionSpace();
2310 FunctionSpace fsr=right.getFunctionSpace();
2311 if (fsl!=fsr) {
2312 signed char intres=fsl.getDomain()->preferredInterpolationOnDomain(fsr.getTypeCode(), fsl.getTypeCode());
2313 if (intres==0)
2314 {
2315 std::string msg="Error - attempt to combine incompatible FunctionSpaces.";
2316 msg+=fsl.toString();
2317 msg+=" ";
2318 msg+=fsr.toString();
2319 throw DataException(msg.c_str());
2320 }
2321 else if (intres==1)
2322 {
2323 // an interpolation is required so create a new Data
2324 tempRight=Data(right,fsl);
2325 }
2326 else // reverse interpolation preferred
2327 {
2328 // interpolate onto the RHS function space
2329 Data tempLeft(*this,fsr);
2330 set_m_data(tempLeft.m_data);
2331 }
2332 }
2333 operandCheck(tempRight);
2334 //
2335 // ensure this has the right type for the RHS
2336 typeMatchRight(tempRight);
2337 //
2338 // Need to cast to the concrete types so that the correct binaryOp
2339 // is called.
2340 if (isExpanded()) {
2341 //
2342 // Expanded data will be done in parallel, the right hand side can be
2343 // of any data type
2344 DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2345 EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2346 escript::binaryOpDataReady(*leftC,*(tempRight.getReady()),operation);
2347 } else if (isTagged()) {
2348 //
2349 // Tagged data is operated on serially, the right hand side can be
2350 // either DataConstant or DataTagged
2351 DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2352 EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2353 if (right.isTagged()) {
2354 DataTagged* rightC=dynamic_cast<DataTagged*>(tempRight.m_data.get());
2355 EsysAssert((rightC!=0), "Programming error - casting to DataTagged.");
2356 escript::binaryOpDataReady(*leftC,*rightC,operation);
2357 } else {
2358 DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2359 EsysAssert((rightC!=0), "Programming error - casting to DataConstant.");
2360 escript::binaryOpDataReady(*leftC,*rightC,operation);
2361 }
2362 } else if (isConstant()) {
2363 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2364 DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2365 EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2366 escript::binaryOpDataReady(*leftC,*rightC,operation);
2367 }
2368 }
2369
2370 /**
2371 \brief
2372 Perform the given binary operation with this and right as operands.
2373 Right is a Data object.
2374 */
2375 template <class BinaryFunction>
2376 inline
2377 void
2378 Data::binaryOp(const Data& right,
2379 BinaryFunction operation)
2380 {
2381 //
2382 // if this has a rank of zero promote it to the rank of the RHS
2383 if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
2384 throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
2385 }
2386
2387 if (isLazy() || right.isLazy())
2388 {
2389 throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
2390 }
2391 //
2392 // initially make the temporary a shallow copy
2393 //Data tempRight(right); // unfortunately, this adds an owner and breaks d+=d
2394
2395 DataAbstract_ptr rp=right.m_data;
2396
2397 FunctionSpace fsl=getFunctionSpace();
2398 FunctionSpace fsr=right.getFunctionSpace();
2399 if (fsl!=fsr) {
2400 signed char intres=fsl.getDomain()->preferredInterpolationOnDomain(fsr.getTypeCode(), fsl.getTypeCode());
2401 if (intres==0)
2402 {
2403 std::string msg="Error - attempt to combine incompatible FunctionSpaces.";
2404 msg+=fsl.toString();
2405 msg+=" ";
2406 msg+=fsr.toString();
2407 throw DataException(msg.c_str());
2408 }
2409 else if (intres==1)
2410 {
2411 // an interpolation is required so create a new Data
2412 rp=Data(right,fsl).m_data;
2413 }
2414 else // reverse interpolation preferred
2415 {
2416 // interpolate onto the RHS function space
2417 Data tempLeft(*this,fsr);
2418 set_m_data(tempLeft.m_data);
2419 }
2420 }
2421 m_data->operandCheck(*rp.get());
2422 //operandCheck(tempRight);
2423 //
2424 // ensure this has the right type for the RHS
2425 typeMatchRight(right); // yes we are actually processing rp
2426 // (which may be pointing at a different object) not right.
2427 // but rp and right will be referring to the same type of data
2428 //
2429 // Need to cast to the concrete types so that the correct binaryOp
2430 // is called.
2431 if (isExpanded()) {
2432 //
2433 // Expanded data will be done in parallel, the right hand side can be
2434 // of any data type
2435 DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2436 EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2437 DataReady* dr=dynamic_cast<DataReady*>(rp.get());
2438 escript::binaryOp(*leftC,*dr,operation);
2439 } else if (isTagged()) {
2440 //
2441 // Tagged data is operated on serially, the right hand side can be
2442 // either DataConstant or DataTagged
2443 DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2444 EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2445 if (right.isTagged()) {
2446 DataTagged* rightC=dynamic_cast<DataTagged*>(rp.get());
2447 EsysAssert((rightC!=0), "Programming error - casting to DataTagged.");
2448 escript::binaryOp(*leftC,*rightC,operation);
2449 } else {
2450 DataConstant* rightC=dynamic_cast<DataConstant*>(rp.get());
2451 EsysAssert((rightC!=0), "Programming error - casting to DataConstant.");
2452 escript::binaryOp(*leftC,*rightC,operation);
2453 }
2454 } else if (isConstant()) {
2455 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2456 DataConstant* rightC=dynamic_cast<DataConstant*>(rp.get());
2457 EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2458 escript::binaryOp(*leftC,*rightC,operation);
2459 }
2460 }
2461
2462 /**
2463 \brief
2464 Perform the given Data object reduction algorithm on this and return the result.
2465 Given operation combines each element of each data point, thus argument
2466 object (*this) is a rank n Data object, and returned object is a scalar.
2467 Calls escript::algorithm.
2468 */
2469 template <class BinaryFunction>
2470 inline
2471 DataTypes::real_t
2472 Data::algorithm(BinaryFunction operation, DataTypes::real_t initial_value) const
2473 {
2474 if (isExpanded()) {
2475 DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2476 EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2477 return escript::algorithm(*leftC,operation,initial_value);
2478 } else if (isTagged()) {
2479 DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2480 EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2481 return escript::algorithm(*leftC,operation,initial_value);
2482 } else if (isConstant()) {
2483 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2484 EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2485 return escript::algorithm(*leftC,operation,initial_value);
2486 } else if (isEmpty()) {
2487 throw DataException("Error - Operations (algorithm) not permitted on instances of DataEmpty.");
2488 } else if (isLazy()) {
2489 throw DataException("Error - Operations not permitted on instances of DataLazy.");
2490 } else {
2491 throw DataException("Error - Data encapsulates an unknown type.");
2492 }
2493 }
2494
2495 /**
2496 \brief
2497 Perform the given data point reduction algorithm on data and return the result.
2498 Given operation combines each element within each data point into a scalar,
2499 thus argument object is a rank n Data object, and returned object is a
2500 rank 0 Data object.
2501 Calls escript::dp_algorithm.
2502 */
2503 template <class BinaryFunction>
2504 inline
2505 Data
2506 Data::dp_algorithm(BinaryFunction operation, DataTypes::real_t initial_value) const
2507 {
2508 if (isEmpty()) {
2509 throw DataException("Error - Operations (dp_algorithm) not permitted on instances of DataEmpty.");
2510 }
2511 else if (isExpanded()) {
2512 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2513 DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2514 DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2515 EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2516 EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2517 escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2518 return result;
2519 }
2520 else if (isTagged()) {
2521 DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2522 EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2523 DataTypes::RealVectorType defval(1);
2524 defval[0]=0;
2525 DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2526 escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2527 return Data(resultT); // note: the Data object now owns the resultT pointer
2528 }
2529 else if (isConstant()) {
2530 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2531 DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2532 DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2533 EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2534 EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2535 escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2536 return result;
2537 } else if (isLazy()) {
2538 throw DataException("Error - Operations not permitted on instances of DataLazy.");
2539 } else {
2540 throw DataException("Error - Data encapsulates an unknown type.");
2541 }
2542 }
2543
2544 /**
2545 \brief
2546 Compute a tensor operation with two Data objects
2547 \param arg_0 - Input - Data object
2548 \param arg_1 - Input - Data object
2549 \param operation - Input - Binary op functor
2550 */
2551 template <typename BinaryFunction>
2552 inline
2553 Data
2554 C_TensorBinaryOperation(Data const &arg_0,
2555 Data const &arg_1,
2556 BinaryFunction operation)
2557 {
2558 if (arg_0.isEmpty() || arg_1.isEmpty())
2559 {
2560 throw DataException("Error - Operations (C_TensorBinaryOperation) not permitted on instances of DataEmpty.");
2561 }
2562 if (arg_0.isLazy() || arg_1.isLazy())
2563 {
2564 throw DataException("Error - Operations not permitted on lazy data.");
2565 }
2566
2567 // for disambiguation of some methods later
2568 typename decltype(operation)::first_argument_type dummy0=0;
2569 typename decltype(operation)::second_argument_type dummy1=0;
2570 typename decltype(operation)::result_type dummyr=0;
2571
2572 // Interpolate if necessary and find an appropriate function space
2573 Data arg_0_Z, arg_1_Z;
2574 FunctionSpace fsl=arg_0.getFunctionSpace();
2575 FunctionSpace fsr=arg_1.getFunctionSpace();
2576 if (fsl!=fsr) {
2577 signed char intres=fsl.getDomain()->preferredInterpolationOnDomain(fsr.getTypeCode(), fsl.getTypeCode());
2578 if (intres==0)
2579 {
2580 std::string msg="Error - C_TensorBinaryOperation: arguments have incompatible function spaces.";
2581 msg+=fsl.toString();
2582 msg+=" ";
2583 msg+=fsr.toString();
2584 throw DataException(msg.c_str());
2585 }
2586 else if (intres==1)
2587 {
2588 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2589 arg_0_Z =Data(arg_0);
2590 }
2591 else // reverse interpolation preferred
2592 {
2593 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2594 arg_1_Z = Data(arg_1);
2595 }
2596 } else {
2597 arg_0_Z = Data(arg_0);
2598 arg_1_Z = Data(arg_1);
2599 }
2600 // Get rank and shape of inputs
2601 int rank0 = arg_0_Z.getDataPointRank();
2602 int rank1 = arg_1_Z.getDataPointRank();
2603 DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2604 DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2605 int size0 = arg_0_Z.getDataPointSize();
2606 int size1 = arg_1_Z.getDataPointSize();
2607
2608
2609
2610 // Declare output Data object
2611 Data res;
2612
2613 if (shape0 == shape1) {
2614 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2615 typename decltype(operation)::first_argument_type dummy=0;
2616 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
2617 if (arg_0_Z.isComplex() || arg_1_Z.isComplex())
2618 {
2619 res.complicate(); // It would be much better to create the Data object as complex to start with
2620 } // But that would require more work so let's just get this case working first
2621 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0, dummy));
2622 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0, dummy1));
2623 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(0, dummyr));
2624
2625 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2626 }
2627 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2628
2629 // Prepare the DataConstant input
2630 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2631
2632 // Borrow DataTagged input from Data object
2633 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2634
2635 // Prepare a DataTagged output 2
2636 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
2637 res.tag();
2638 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2639
2640
2641
2642 // Prepare offset into DataConstant
2643 int offset_0 = tmp_0->getPointOffset(0,0);
2644 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2645
2646 // Get the pointers to the actual data
2647 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDefaultValueRO(0, dummy1));
2648 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
2649
2650 // Compute a result for the default
2651 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2652 // Compute a result for each tag
2653 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2654 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2655 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2656 tmp_2->addTag(i->first);
2657 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0, dummy1));
2658 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
2659
2660 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2661 }
2662
2663 }
2664 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2665 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2666 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2667 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2668 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2669
2670 int sampleNo_1,dataPointNo_1;
2671 int numSamples_1 = arg_1_Z.getNumSamples();
2672 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2673 int offset_0 = tmp_0->getPointOffset(0,0);
2674 res.requireWrite();
2675 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2676 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2677 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2678 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2679 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2680 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2681 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2682 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
2683 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2684 }
2685 }
2686
2687 }
2688 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2689 // Borrow DataTagged input from Data object
2690 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2691
2692 // Prepare the DataConstant input
2693 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2694
2695 // Prepare a DataTagged output 2
2696 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2697 res.tag();
2698 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2699
2700 // Prepare offset into DataConstant
2701 int offset_1 = tmp_1->getPointOffset(0,0);
2702
2703 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2704 // Get the pointers to the actual data
2705 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDefaultValueRO(0, dummy0));
2706 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
2707 // Compute a result for the default
2708 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2709 // Compute a result for each tag
2710 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2711 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2712 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2713 tmp_2->addTag(i->first);
2714 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0, dummy0));
2715 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0,dummyr));
2716 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2717 }
2718
2719 }
2720 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2721 // Borrow DataTagged input from Data object
2722 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2723
2724 // Borrow DataTagged input from Data object
2725 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2726
2727 // Prepare a DataTagged output 2
2728 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2729 res.tag(); // DataTagged output
2730 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2731
2732 // Get the pointers to the actual data
2733 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDefaultValueRO(0, dummy0));
2734 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDefaultValueRO(0, dummy1));
2735 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
2736
2737 // Compute a result for the default
2738 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2739 // Merge the tags
2740 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2741 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2742 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2743 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2744 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2745 }
2746 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2747 tmp_2->addTag(i->first);
2748 }
2749 // Compute a result for each tag
2750 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2751 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2752
2753 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0, dummy0));
2754 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0, dummy1));
2755 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
2756
2757 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2758 }
2759
2760 }
2761 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2762 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2763 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2764 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2765 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2766 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2767
2768 int sampleNo_0,dataPointNo_0;
2769 int numSamples_0 = arg_0_Z.getNumSamples();
2770 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2771 res.requireWrite();
2772 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2773 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2774 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2775 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2776 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2777 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2778 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2779 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2780 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
2781 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2782 }
2783 }
2784
2785 }
2786 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2787 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2788 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2789 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2790 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2791
2792 int sampleNo_0,dataPointNo_0;
2793 int numSamples_0 = arg_0_Z.getNumSamples();
2794 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2795 int offset_1 = tmp_1->getPointOffset(0,0);
2796 res.requireWrite();
2797 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2798 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2799 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2800 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2801 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2802
2803 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2804 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2805 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
2806
2807
2808 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2809 }
2810 }
2811
2812 }
2813 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2814 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2815 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2816 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2817 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2818 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2819
2820 int sampleNo_0,dataPointNo_0;
2821 int numSamples_0 = arg_0_Z.getNumSamples();
2822 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2823 res.requireWrite();
2824 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2825 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2826 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2827 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2828 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2829 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2830 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2831 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2832 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
2833 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2834 }
2835 }
2836
2837 }
2838 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2839 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2840 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2841 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2842 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2843 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2844
2845 int sampleNo_0,dataPointNo_0;
2846 int numSamples_0 = arg_0_Z.getNumSamples();
2847 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2848 res.requireWrite();
2849 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2850 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2851 dataPointNo_0=0;
2852 // for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2853 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2854 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2855 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2856 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2857 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2858 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
2859 tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2860 // }
2861 }
2862
2863 }
2864 else {
2865 throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2866 }
2867
2868 } else if (0 == rank0) {
2869 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2870 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataConstant output
2871 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0, dummy0));
2872 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0, dummy1));
2873 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(0, dummyr));
2874 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2875 }
2876 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2877
2878 // Prepare the DataConstant input
2879 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2880
2881 // Borrow DataTagged input from Data object
2882 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2883
2884 // Prepare a DataTagged output 2
2885 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataTagged output
2886 res.tag();
2887 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2888
2889 // Prepare offset into DataConstant
2890 int offset_0 = tmp_0->getPointOffset(0,0);
2891 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2892
2893 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDefaultValueRO(0, dummy1));
2894 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
2895
2896 // Compute a result for the default
2897 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2898 // Compute a result for each tag
2899 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2900 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2901 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2902 tmp_2->addTag(i->first);
2903 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0, dummy1));
2904 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
2905 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2906 }
2907
2908 }
2909 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2910
2911 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2912 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2913 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2914 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2915
2916 int sampleNo_1;
2917 int numSamples_1 = arg_1_Z.getNumSamples();
2918 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2919 int offset_0 = tmp_0->getPointOffset(0,0);
2920 const typename decltype(operation)::first_argument_type *ptr_src = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
2921 typename decltype(operation)::first_argument_type ptr_0 = ptr_src[0];
2922 int size = size1*numDataPointsPerSample_1;
2923 res.requireWrite();
2924 #pragma omp parallel for private(sampleNo_1) schedule(static)
2925 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2926 int offset_1 = tmp_1->getPointOffset(sampleNo_1,0);
2927 int offset_2 = tmp_2->getPointOffset(sampleNo_1,0);
2928 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2929 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
2930 tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);
2931 }
2932
2933 }
2934 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2935
2936 // Borrow DataTagged input from Data object
2937 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2938
2939 // Prepare the DataConstant input
2940 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2941
2942 // Prepare a DataTagged output 2
2943 res = Data(0.0, shape1, arg_0_Z.getFunctionSpace()); // DataTagged output
2944 res.tag();
2945 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2946
2947 // Prepare offset into DataConstant
2948 int offset_1 = tmp_1->getPointOffset(0,0);
2949 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
2950
2951 // Get the pointers to the actual data
2952 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDefaultValueRO(0, dummy0));
2953 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
2954
2955
2956 // Compute a result for the default
2957 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2958 // Compute a result for each tag
2959 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2960 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2961 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2962 tmp_2->addTag(i->first);
2963 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0, dummy0));
2964 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
2965
2966 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2967 }
2968
2969 }
2970 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2971
2972 // Borrow DataTagged input from Data object
2973 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2974
2975 // Borrow DataTagged input from Data object
2976 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2977
2978 // Prepare a DataTagged output 2
2979 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2980 res.tag(); // DataTagged output
2981 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2982
2983 // Get the pointers to the actual data
2984 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDefaultValueRO(0, dummy0));
2985 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDefaultValueRO(0, dummy1));
2986 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
2987
2988 // Compute a result for the default
2989 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2990 // Merge the tags
2991 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2992 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2993 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2994 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2995 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2996 }
2997 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2998 tmp_2->addTag(i->first);
2999 }
3000 // Compute a result for each tag
3001 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
3002 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
3003 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0, dummy0));
3004 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0, dummy1));
3005 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
3006
3007 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
3008 }
3009
3010 }
3011 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
3012
3013 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3014 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3015 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3016 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3017 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3018
3019 int sampleNo_0,dataPointNo_0;
3020 int numSamples_0 = arg_0_Z.getNumSamples();
3021 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3022 res.requireWrite();
3023 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3024 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3025 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
3026 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3027 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3028 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3029 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3030 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3031 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3032 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
3033 }
3034 }
3035
3036 }
3037 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
3038 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3039 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3040 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
3041 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3042
3043 int sampleNo_0,dataPointNo_0;
3044 int numSamples_0 = arg_0_Z.getNumSamples();
3045 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3046 int offset_1 = tmp_1->getPointOffset(0,0);
3047 res.requireWrite();
3048 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3049 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3050 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3051 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3052 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3053 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3054 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3055 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3056 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
3057 }
3058 }
3059
3060
3061 }
3062 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
3063
3064 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3065 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3066 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3067 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3068 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3069
3070 int sampleNo_0,dataPointNo_0;
3071 int numSamples_0 = arg_0_Z.getNumSamples();
3072 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3073 res.requireWrite();
3074 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3075 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3076 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
3077 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3078 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3079 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3080 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3081 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3082 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3083 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
3084 }
3085 }
3086
3087 }
3088 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
3089
3090 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3091 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3092 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3093 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3094 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3095
3096 int sampleNo_0,dataPointNo_0;
3097 int numSamples_0 = arg_0_Z.getNumSamples();
3098 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3099 res.requireWrite();
3100 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3101 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3102 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3103 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3104 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3105 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3106 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3107 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3108 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3109 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
3110 }
3111 }
3112
3113 }
3114 else {
3115 throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
3116 }
3117
3118 } else if (0 == rank1) {
3119 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
3120 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
3121 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0, dummy0));
3122 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0, dummy1));
3123 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(0, dummyr));
3124 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3125 }
3126 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
3127
3128 // Prepare the DataConstant input
3129 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
3130
3131 // Borrow DataTagged input from Data object
3132 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3133
3134 // Prepare a DataTagged output 2
3135 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
3136 res.tag();
3137 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3138
3139 // Prepare offset into DataConstant
3140 int offset_0 = tmp_0->getPointOffset(0,0);
3141 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3142
3143 //Get the pointers to the actual data
3144 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDefaultValueRO(0, dummy1));
3145 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
3146
3147 // Compute a result for the default
3148 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3149 // Compute a result for each tag
3150 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
3151 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3152 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
3153 tmp_2->addTag(i->first);
3154 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0, dummy1));
3155 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
3156 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3157 }
3158 }
3159 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
3160
3161 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3162 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
3163 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3164 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3165
3166 int sampleNo_1,dataPointNo_1;
3167 int numSamples_1 = arg_1_Z.getNumSamples();
3168 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
3169 int offset_0 = tmp_0->getPointOffset(0,0);
3170 res.requireWrite();
3171 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
3172 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
3173 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
3174 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
3175 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
3176 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3177 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3178 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3179 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3180 }
3181 }
3182
3183 }
3184 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
3185
3186 // Borrow DataTagged input from Data object
3187 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3188
3189 // Prepare the DataConstant input
3190 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
3191
3192 // Prepare a DataTagged output 2
3193 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
3194 res.tag();
3195 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3196
3197 // Prepare offset into DataConstant
3198 int offset_1 = tmp_1->getPointOffset(0,0);
3199 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3200 // Get the pointers to the actual data
3201 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDefaultValueRO(0, dummy0));
3202 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
3203 // Compute a result for the default
3204 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3205 // Compute a result for each tag
3206 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3207 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3208 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3209 tmp_2->addTag(i->first);
3210 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0, dummy0));
3211 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
3212 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3213 }
3214
3215 }
3216 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
3217
3218 // Borrow DataTagged input from Data object
3219 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3220
3221 // Borrow DataTagged input from Data object
3222 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3223
3224 // Prepare a DataTagged output 2
3225 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
3226 res.tag(); // DataTagged output
3227 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3228
3229 // Get the pointers to the actual data
3230 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDefaultValueRO(0, dummy0));
3231 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDefaultValueRO(0, dummy1));
3232 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0, dummyr));
3233
3234 // Compute a result for the default
3235 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3236 // Merge the tags
3237 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3238 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3239 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
3240 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3241 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
3242 }
3243 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
3244 tmp_2->addTag(i->first);
3245 }
3246 // Compute a result for each tag
3247 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
3248 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
3249 const typename decltype(operation)::first_argument_type *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0, dummy0));
3250 const typename decltype(operation)::second_argument_type *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0, dummy1));
3251 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0, dummyr));
3252 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3253 }
3254
3255 }
3256 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
3257
3258 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3259 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3260 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3261 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3262 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3263
3264 int sampleNo_0,dataPointNo_0;
3265 int numSamples_0 = arg_0_Z.getNumSamples();
3266 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3267 res.requireWrite();
3268 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3269 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3270 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
3271 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3272 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3273 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3274 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3275 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3276 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3277 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3278 }
3279 }
3280
3281 }
3282 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
3283 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3284 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3285 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
3286 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3287
3288 int sampleNo_0;
3289 int numSamples_0 = arg_0_Z.getNumSamples();
3290 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3291 int offset_1 = tmp_1->getPointOffset(0,0);
3292 const typename decltype(operation)::second_argument_type *ptr_src = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3293 typename decltype(operation)::second_argument_type ptr_1 = ptr_src[0];
3294 int size = size0 * numDataPointsPerSample_0;
3295 res.requireWrite();
3296 #pragma omp parallel for private(sampleNo_0) schedule(static)
3297 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3298 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0);
3299 int offset_2 = tmp_2->getPointOffset(sampleNo_0,0);
3300 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3301 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3302 tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);
3303 }
3304
3305
3306 }
3307 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
3308
3309 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3310 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3311 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3312 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3313 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3314
3315 int sampleNo_0,dataPointNo_0;
3316 int numSamples_0 = arg_0_Z.getNumSamples();
3317 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3318 res.requireWrite();
3319 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3320 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3321 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
3322 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3323 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3324 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3325 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3326 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3327 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3328 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3329 }
3330 }
3331
3332 }
3333 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
3334
3335 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3336 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3337 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3338 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3339 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3340
3341 int sampleNo_0,dataPointNo_0;
3342 int numSamples_0 = arg_0_Z.getNumSamples();
3343 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3344 res.requireWrite();
3345 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3346 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3347 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3348 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3349 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3350 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3351 const typename decltype(operation)::first_argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, dummy0));
3352 const typename decltype(operation)::second_argument_type *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1, dummy1));
3353 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2, dummyr));
3354 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
3355 }
3356 }
3357
3358 }
3359 else {
3360 throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
3361 }
3362
3363 } else {
3364 throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
3365 }
3366
3367 return res;
3368 }
3369
3370 template <typename UnaryFunction>
3371 Data
3372 C_TensorUnaryOperation(Data const &arg_0,
3373 UnaryFunction operation)
3374 {
3375 if (arg_0.isEmpty()) // do this before we attempt to interpolate
3376 {
3377 throw DataException("Error - Operations (C_TensorUnaryOperation) not permitted on instances of DataEmpty.");
3378 }
3379 if (arg_0.isLazy())
3380 {
3381 throw DataException("Error - Operations not permitted on lazy data.");
3382 }
3383 // Interpolate if necessary and find an appropriate function space
3384 Data arg_0_Z = Data(arg_0);
3385
3386 // Get rank and shape of inputs
3387 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
3388 int size0 = arg_0_Z.getDataPointSize();
3389
3390 // Sanity check on the types of the function
3391 if (typeid(typename UnaryFunction::argument_type)!=typeid(typename UnaryFunction::result_type))
3392 {
3393 throw DataException("Error - Types for C_TensorUnaryOperation must be consistent");
3394 }
3395
3396 // Declare output Data object
3397 Data res;
3398
3399 if (arg_0_Z.isConstant()) {
3400 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataConstant output
3401 const typename UnaryFunction::argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
3402 typename UnaryFunction::result_type *ptr_2 = &(res.getDataAtOffsetRW(0));
3403 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
3404 }
3405 else if (arg_0_Z.isTagged()) {
3406
3407 // Borrow DataTagged input from Data object
3408 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3409
3410 // Prepare a DataTagged output 2
3411 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
3412 res.tag();
3413 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3414
3415 // Get the pointers to the actual data
3416 const typename decltype(operation)::argument_type *ptr_0 = &(tmp_0->getDefaultValueRO(0));
3417 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDefaultValueRW(0));
3418 // Compute a result for the default
3419 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
3420 // Compute a result for each tag
3421 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3422 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3423 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3424 tmp_2->addTag(i->first);
3425 const typename decltype(operation)::argument_type *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
3426 typename decltype(operation)::result_type *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
3427 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
3428 }
3429
3430 }
3431 else if (arg_0_Z.isExpanded()) {
3432
3433 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
3434 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3435 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3436
3437 int sampleNo_0,dataPointNo_0;
3438 int numSamples_0 = arg_0_Z.getNumSamples();
3439 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3440 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3441 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3442 dataPointNo_0=0;
3443 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3444 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3445 const typename decltype(operation)::argument_type *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3446 typename decltype(operation)::result_type *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3447 tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3448 }
3449 }
3450 else {
3451 throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3452 }
3453
3454 return res;
3455 }
3456
3457
3458 Data
3459 C_TensorUnaryOperation(Data const &arg_0,
3460 escript::ESFunction operation,
3461 DataTypes::real_t tol=0
3462 );
3463
3464
3465
3466
3467
3468 }
3469 #endif

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26