/[escript]/trunk/escript/src/Data.h
ViewVC logotype

Contents of /trunk/escript/src/Data.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3390 - (show annotations)
Thu Dec 2 00:34:37 2010 UTC (10 years ago) by jfenwick
File MIME type: text/plain
File size: 93400 byte(s)
RandomData added

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26