/[escript]/branches/arrexp_trunk2098/escript/src/Data.cpp
ViewVC logotype

Contents of /branches/arrexp_trunk2098/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2119 - (show annotations)
Tue Dec 2 06:06:04 2008 UTC (11 years, 10 months ago) by jfenwick
File size: 92926 byte(s)
Branch commit.
Threading Wrapped Array through the code.
<
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "Data.h"
16
17 #include "DataExpanded.h"
18 #include "DataConstant.h"
19 #include "DataTagged.h"
20 #include "DataEmpty.h"
21 #include "DataLazy.h"
22 #include "FunctionSpaceFactory.h"
23 #include "AbstractContinuousDomain.h"
24 #include "UnaryFuncs.h"
25 #include "FunctionSpaceException.h"
26 #include "EscriptParams.h"
27
28 extern "C" {
29 #include "esysUtils/blocktimer.h"
30 }
31
32 #include <fstream>
33 #include <algorithm>
34 #include <vector>
35 #include <functional>
36 #include <sstream> // so we can throw messages about ranks
37
38 #include <boost/python/dict.hpp>
39 #include <boost/python/extract.hpp>
40 #include <boost/python/long.hpp>
41
42 #include "WrappedArray.h"
43
44 using namespace std;
45 using namespace boost::python;
46 using namespace boost;
47 using namespace escript;
48
49 // ensure the current object is not a DataLazy
50 // The idea was that we could add an optional warning whenever a resolve is forced
51 #define FORCERESOLVE if (isLazy()) {resolve();}
52
53
54 void escript::jtest(boost::python::object& obj)
55 {
56 cout << "Testing1\n";
57 try
58 {
59 double v=extract<double>(obj);
60 cout << "Found a double " << v << endl;
61 return;
62 }
63 catch (...)
64 {
65 PyErr_Clear();
66 // There isn't actually an error here
67 }
68 cout << "Testing2\n";
69 try
70 {
71 double v=extract<double>(obj[make_tuple()]);
72 cout << "Found a () double " << v << endl;
73 return;
74 }
75 catch (...)
76 {
77 PyErr_Clear();
78 }
79 cout << "Testing 3\n";
80 WrappedArray w(obj);
81 DataTypes::ShapeType shape=w.getShape();
82 size_t rank=shape.size();
83 cout << "Rank=" << rank << endl;
84 for (unsigned int i=0;i<rank;++i)
85 {
86 cout << shape[i] << ",";
87 }
88 DataVector vec(DataTypes::noValues(shape));
89 cout << endl;
90 if (rank==1)
91 {
92 int i=0;
93 try
94 {
95 for (i=0;i<shape[0];++i)
96 {
97 vec[i]=w.getElt(i);
98 }
99 }
100 catch (...)
101 {
102 PyErr_Clear();
103 ostringstream oss;
104 oss << "Array filter - Error reading element " << i;
105 throw DataException(oss.str());
106 }
107 }
108 else if (rank==2)
109 {
110 int i=0;
111 int j=0;
112 try
113 {
114 for (i=0;i<shape[0];++i) // element
115 {
116 for (j=0;j<shape[1];++j) // row
117 {
118 vec[DataTypes::getRelIndex(shape,i,j)]=w.getElt(i,j);
119 }
120 }
121 }
122 catch (...)
123 {
124 PyErr_Clear();
125 ostringstream oss;
126 oss << "Array filter - Error reading element " << i << ',' << j;
127 throw DataException(oss.str());
128 }
129 }
130 else if (rank==3)
131 {
132 int i=0;
133 int j=0;
134 int k=0;
135 try
136 {
137 for (i=0;i<shape[0];++i)
138 {
139 for (j=0;j<shape[1];++j)
140 {
141 for (k=0;k<shape[2];++k)
142 {
143 vec[DataTypes::getRelIndex(shape,i,j,k)]=w.getElt(i,j,k);
144 }
145 }
146 }
147 }
148 catch (...)
149 {
150 PyErr_Clear();
151 ostringstream oss;
152 oss << "Array filter - Error reading element " << i << ',' << j << ',' << k;
153 throw DataException(oss.str());
154 }
155 }
156 else if (rank==4)
157 {
158 int i=0;
159 int j=0;
160 int k=0;
161 int m=0;
162 try
163 {
164 for (i=0;i<shape[0];++i)
165 {
166 for (j=0;j<shape[1];++j)
167 {
168 for (k=0;k<shape[2];++k)
169 {
170 for (m=0;m<shape[3];++m)
171 {
172 vec[DataTypes::getRelIndex(shape,i,j,k,m)]=w.getElt(i,j,k,m);
173 }
174 }
175 }
176 }
177 }
178 catch (...)
179 {
180 PyErr_Clear();
181 ostringstream oss;
182 oss << "Array filter - Error reading element " << i << ',' << j << ',' << k << ',' << m;
183 throw DataException(oss.str());
184 }
185 }
186 for (int i=0;i<vec.size();++i)
187 {
188 cout << vec[i] << " ";
189 }
190 cout << endl;
191 }
192
193
194
195
196
197
198 /*
199 ArrayTools::getObjShape(obj,shape);
200 size_t rank=shape.size();
201 cout << "Rank=" << rank << endl;
202 for (unsigned int i=0;i<rank;++i)
203 {
204 cout << shape[i] << ",";
205 }
206 DataVector vec(DataTypes::noValues(shape));
207 cout << endl;
208 if (rank==1)
209 {
210 int i=0;
211 try
212 {
213 for (i=0;i<shape[0];++i)
214 {
215 vec[i]=extract<double>(obj[i]);
216 }
217 }
218 catch (...)
219 {
220 PyErr_Clear();
221 ostringstream oss;
222 oss << "Array filter - Error reading element " << i;
223 throw DataException(oss.str());
224 }
225 }
226 else if (rank==2)
227 {
228 int i=0;
229 int j=0;
230 try
231 {
232 for (i=0;i<shape[0];++i) // element
233 {
234 for (j=0;j<shape[1];++j) // row
235 {
236 vec[DataTypes::getRelIndex(shape,i,j)]=extract<double>(obj[i][j]);
237 }
238 }
239 }
240 catch (...)
241 {
242 PyErr_Clear();
243 ostringstream oss;
244 oss << "Array filter - Error reading element " << i << ',' << j;
245 throw DataException(oss.str());
246 }
247 }
248 else if (rank==3)
249 {
250 int i=0;
251 int j=0;
252 int k=0;
253 try
254 {
255 for (i=0;i<shape[0];++i)
256 {
257 for (j=0;j<shape[1];++j)
258 {
259 for (k=0;k<shape[2];++k)
260 {
261 vec[DataTypes::getRelIndex(shape,i,j,k)]=extract<double>(obj[i][j][k]);
262 }
263 }
264 }
265 }
266 catch (...)
267 {
268 PyErr_Clear();
269 ostringstream oss;
270 oss << "Array filter - Error reading element " << i << ',' << j << ',' << k;
271 throw DataException(oss.str());
272 }
273 }
274 else if (rank==4)
275 {
276 int i=0;
277 int j=0;
278 int k=0;
279 int m=0;
280 try
281 {
282 for (i=0;i<shape[0];++i)
283 {
284 for (j=0;j<shape[1];++j)
285 {
286 for (k=0;k<shape[2];++k)
287 {
288 for (m=0;m<shape[3];++m)
289 {
290 vec[DataTypes::getRelIndex(shape,i,j,k,m)]=extract<double>(obj[i][j][k][m]);
291 }
292 }
293 }
294 }
295 }
296 catch (...)
297 {
298 PyErr_Clear();
299 ostringstream oss;
300 oss << "Array filter - Error reading element " << i << ',' << j << ',' << k << ',' << m;
301 throw DataException(oss.str());
302 }
303 }
304 */
305
306
307
308
309
310
311
312
313
314
315
316
317 // void escript::jtest()
318 // {
319 // // string s=extract<std::string>(obj.attr("__str__")());
320 // // cout << "[" << s << "]" << endl;
321 // // int len=extract<int>(obj.attr("__len__")());
322 // // cerr << "Length=" << len << endl;
323 // // for (int j=0;j<len;++j)
324 // // {
325 // // object it=obj.attr("__getitem__")(j);
326 // // string s=extract<std::string>(it.attr("__str__")());
327 // // cerr << "Item: " << j << s;
328 // // }
329 //
330 // // try
331 // // {
332 // // object z=obj.attr("nothere");
333 // // }
334 // // catch(...)
335 // // {
336 // // cout << "Caught exception\n";
337 // // }
338 //
339 // // 1 test capabilities -- check for required methods - write this later
340 // // 2 compute rank and shape
341 // // 3 get values
342 //
343 // DataTypes::ShapeType shape;
344 // object t=obj;
345 // // int len=0;
346 // // while (len=extract<int>(obj.attr("__len__")()),len>=1)
347 // // {
348 // // shape.push_back(len);
349 // //
350 // // }
351 //
352 // try
353 // {
354 // object z=t.attr("__getitem__");
355 // if (z==0)
356 // {
357 // cerr << "Found null" << endl;
358 // }
359 // }
360 // catch (...)
361 // {
362 // cerr << "Type does not support __getitem__" << endl;
363 // boost::python::handle_exception();
364 // return;
365 // }
366 // try
367 // {
368 // object z=t.attr("__len__")();
369 // }
370 // catch (...)
371 // {
372 // cerr << "Object does not support __len__" << endl;
373 // boost::python::handle_exception();
374 // return;
375 // }
376 //
377 // cerr << "Loop entry" << endl;
378 // bool done=false;
379 // do
380 // {
381 // string s=extract<std::string>(t.attr("__str__")());
382 // cerr << "Processing: " << s<< endl;
383 // int len=extract<int>(t.attr("__len__")());
384 // cerr << "Length is " << len << endl;
385 // if (len<1)
386 // {
387 // throw DataException("Array filter - No empty elements please.");
388 // }
389 // shape.push_back(len);
390 // if (shape.size()>ESCRIPT_MAX_DATA_RANK)
391 // {
392 // throw DataException("ArrayFilter - Maximum permitted rank is 4");
393 // }
394 // cerr << "About to call [0]" << endl;
395 // // t=t.attr("__getitem__")(0);
396 // t=t[0];
397 // try
398 // {
399 // cerr << "About to test __getitem__ " << endl;
400 // // object o=t.attr("__getitem__"); // test if the object has the method
401 // object o=t[0];
402 // cerr << "Call to test succeeded" << endl;
403 // }
404 // catch (...)
405 // {
406 // done=true;
407 // cerr << "In catch" << endl;
408 // boost::python::handle_exception();
409 // }
410 // cerr << "Bottom of the loop" << endl;
411 // } while (!done);
412 // cout << "Rank=" << shape.size() << endl;
413 // for (int i=0;i<shape.size(); ++i)
414 // {
415 // cout << shape[i] << endl;
416 // }
417 // }
418
419 Data::Data()
420 {
421 //
422 // Default data is type DataEmpty
423 DataAbstract* temp=new DataEmpty();
424 m_data=temp->getPtr();
425 m_protected=false;
426 }
427
428 Data::Data(double value,
429 const tuple& shape,
430 const FunctionSpace& what,
431 bool expanded)
432 {
433 DataTypes::ShapeType dataPointShape;
434 for (int i = 0; i < shape.attr("__len__")(); ++i) {
435 dataPointShape.push_back(extract<const int>(shape[i]));
436 }
437
438 int len = DataTypes::noValues(dataPointShape);
439 DataVector temp_data(len,value,len);
440 initialise(temp_data, dataPointShape, what, expanded);
441 m_protected=false;
442 }
443
444 Data::Data(double value,
445 const DataTypes::ShapeType& dataPointShape,
446 const FunctionSpace& what,
447 bool expanded)
448 {
449 int len = DataTypes::noValues(dataPointShape);
450
451 DataVector temp_data(len,value,len);
452 // DataArrayView temp_dataView(temp_data, dataPointShape);
453
454 // initialise(temp_dataView, what, expanded);
455 initialise(temp_data, dataPointShape, what, expanded);
456
457 m_protected=false;
458 }
459
460 Data::Data(const Data& inData)
461 {
462 m_data=inData.m_data;
463 m_protected=inData.isProtected();
464 }
465
466
467 Data::Data(const Data& inData,
468 const DataTypes::RegionType& region)
469 {
470 DataAbstract_ptr dat=inData.m_data;
471 if (inData.isLazy())
472 {
473 dat=inData.m_data->resolve();
474 }
475 else
476 {
477 dat=inData.m_data;
478 }
479 //
480 // Create Data which is a slice of another Data
481 DataAbstract* tmp = dat->getSlice(region);
482 m_data=DataAbstract_ptr(tmp);
483 m_protected=false;
484 }
485
486 Data::Data(const Data& inData,
487 const FunctionSpace& functionspace)
488 {
489 if (inData.isEmpty())
490 {
491 throw DataException("Error - will not interpolate for instances of DataEmpty.");
492 }
493 if (inData.getFunctionSpace()==functionspace) {
494 m_data=inData.m_data;
495 }
496 else
497 {
498
499 if (inData.isConstant()) { // for a constant function, we just need to use the new function space
500 if (!inData.probeInterpolation(functionspace))
501 { // Even though this is constant, we still need to check whether interpolation is allowed
502 throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
503 }
504 // if the data is not lazy, this will just be a cast to DataReady
505 DataReady_ptr dr=inData.m_data->resolve();
506 DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVector());
507 m_data=DataAbstract_ptr(dc);
508 } else {
509 Data tmp(0,inData.getDataPointShape(),functionspace,true);
510 // Note: Must use a reference or pointer to a derived object
511 // in order to get polymorphic behaviour. Shouldn't really
512 // be able to create an instance of AbstractDomain but that was done
513 // as a boost:python work around which may no longer be required.
514 /*const AbstractDomain& inDataDomain=inData.getDomain();*/
515 const_Domain_ptr inDataDomain=inData.getDomain();
516 if (inDataDomain==functionspace.getDomain()) {
517 inDataDomain->interpolateOnDomain(tmp,inData);
518 } else {
519 inDataDomain->interpolateACross(tmp,inData);
520 }
521 m_data=tmp.m_data;
522 }
523 }
524 m_protected=false;
525 }
526
527 Data::Data(DataAbstract* underlyingdata)
528 {
529 // m_data=shared_ptr<DataAbstract>(underlyingdata);
530 m_data=underlyingdata->getPtr();
531 m_protected=false;
532 }
533
534 Data::Data(DataAbstract_ptr underlyingdata)
535 {
536 m_data=underlyingdata;
537 m_protected=false;
538 }
539
540
541 // Data::Data(const numeric::array& value,
542 // const FunctionSpace& what,
543 // bool expanded)
544 // {
545 // initialise(value,what,expanded);
546 // m_protected=false;
547 // }
548
549 Data::Data(const DataTypes::ValueType& value,
550 const DataTypes::ShapeType& shape,
551 const FunctionSpace& what,
552 bool expanded)
553 {
554 initialise(value,shape,what,expanded);
555 m_protected=false;
556 }
557
558
559 Data::Data(const object& value,
560 const FunctionSpace& what,
561 bool expanded)
562 {
563 //numeric::array asNumArray(value);
564 WrappedArray w(value);
565 initialise(w,what,expanded);
566 m_protected=false;
567 }
568
569
570 Data::Data(const object& value,
571 const Data& other)
572 {
573 numeric::array asNumArray(value);
574
575 // extract the shape of the numarray
576 DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);
577 if (DataTypes::getRank(tempShape)==0) {
578
579
580 // get the space for the data vector
581 int len1 = DataTypes::noValues(tempShape);
582 DataVector temp_data(len1, 0.0, len1);
583 temp_data.copyFromNumArray(asNumArray);
584
585 int len = DataTypes::noValues(other.getDataPointShape());
586
587 DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);
588 DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
589 m_data=DataAbstract_ptr(t);
590
591 } else {
592 //
593 // Create a DataConstant with the same sample shape as other
594 DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());
595 // boost::shared_ptr<DataAbstract> sp(t);
596 // m_data=sp;
597 m_data=DataAbstract_ptr(t);
598 }
599 m_protected=false;
600 }
601
602 Data::~Data()
603 {
604
605 }
606
607
608 void Data::initialise(const WrappedArray& value,
609 const FunctionSpace& what,
610 bool expanded)
611 {
612 //
613 // Construct a Data object of the appropriate type.
614 // Construct the object first as there seems to be a bug which causes
615 // undefined behaviour if an exception is thrown during construction
616 // within the shared_ptr constructor.
617 if (expanded) {
618 DataAbstract* temp=new DataExpanded(value, what);
619 m_data=temp->getPtr();
620 } else {
621 DataAbstract* temp=new DataConstant(value, what);
622 m_data=temp->getPtr();
623 }
624 }
625
626
627 void
628 Data::initialise(const boost::python::numeric::array& value,
629 const FunctionSpace& what,
630 bool expanded)
631 {
632 //
633 // Construct a Data object of the appropriate type.
634 // Construct the object first as there seems to be a bug which causes
635 // undefined behaviour if an exception is thrown during construction
636 // within the shared_ptr constructor.
637 if (expanded) {
638 DataAbstract* temp=new DataExpanded(value, what);
639 // boost::shared_ptr<DataAbstract> temp_data(temp);
640 // m_data=temp_data;
641 m_data=temp->getPtr();
642 } else {
643 DataAbstract* temp=new DataConstant(value, what);
644 // boost::shared_ptr<DataAbstract> temp_data(temp);
645 // m_data=temp_data;
646 m_data=temp->getPtr();
647 }
648 }
649
650
651 void
652 Data::initialise(const DataTypes::ValueType& value,
653 const DataTypes::ShapeType& shape,
654 const FunctionSpace& what,
655 bool expanded)
656 {
657 //
658 // Construct a Data object of the appropriate type.
659 // Construct the object first as there seems to be a bug which causes
660 // undefined behaviour if an exception is thrown during construction
661 // within the shared_ptr constructor.
662 if (expanded) {
663 DataAbstract* temp=new DataExpanded(what, shape, value);
664 // boost::shared_ptr<DataAbstract> temp_data(temp);
665 // m_data=temp_data;
666 m_data=temp->getPtr();
667 } else {
668 DataAbstract* temp=new DataConstant(what, shape, value);
669 // boost::shared_ptr<DataAbstract> temp_data(temp);
670 // m_data=temp_data;
671 m_data=temp->getPtr();
672 }
673 }
674
675
676 // void
677 // Data::CompareDebug(const Data& rd)
678 // {
679 // using namespace std;
680 // bool mismatch=false;
681 // std::cout << "Comparing left and right" << endl;
682 // const DataTagged* left=dynamic_cast<DataTagged*>(m_data.get());
683 // const DataTagged* right=dynamic_cast<DataTagged*>(rd.m_data.get());
684 //
685 // if (left==0)
686 // {
687 // cout << "left arg is not a DataTagged\n";
688 // return;
689 // }
690 //
691 // if (right==0)
692 // {
693 // cout << "right arg is not a DataTagged\n";
694 // return;
695 // }
696 // cout << "Num elements=" << left->getVector().size() << ":" << right->getVector().size() << std::endl;
697 // cout << "Shapes ";
698 // if (left->getShape()==right->getShape())
699 // {
700 // cout << "ok\n";
701 // }
702 // else
703 // {
704 // cout << "Problem: shapes do not match\n";
705 // mismatch=true;
706 // }
707 // int lim=left->getVector().size();
708 // if (right->getVector().size()) lim=right->getVector().size();
709 // for (int i=0;i<lim;++i)
710 // {
711 // if (left->getVector()[i]!=right->getVector()[i])
712 // {
713 // cout << "[" << i << "] value mismatch " << left->getVector()[i] << ":" << right->getVector()[i] << endl;
714 // mismatch=true;
715 // }
716 // }
717 //
718 // // still need to check the tag map
719 // // also need to watch what is happening to function spaces, are they copied or what?
720 //
721 // const DataTagged::DataMapType& mapleft=left->getTagLookup();
722 // const DataTagged::DataMapType& mapright=right->getTagLookup();
723 //
724 // if (mapleft.size()!=mapright.size())
725 // {
726 // cout << "Maps are different sizes " << mapleft.size() << ":" << mapright.size() << endl;
727 // mismatch=true;
728 // cout << "Left map\n";
729 // DataTagged::DataMapType::const_iterator i,j;
730 // for (i=mapleft.begin();i!=mapleft.end();++i) {
731 // cout << "(" << i->first << "=>" << i->second << ")\n";
732 // }
733 // cout << "Right map\n";
734 // for (i=mapright.begin();i!=mapright.end();++i) {
735 // cout << "(" << i->first << "=>" << i->second << ")\n";
736 // }
737 // cout << "End map\n";
738 //
739 // }
740 //
741 // DataTagged::DataMapType::const_iterator i,j;
742 // for (i=mapleft.begin(),j=mapright.begin();i!=mapleft.end() && j!=mapright.end();++i,++j) {
743 // if ((i->first!=j->first) || (i->second!=j->second))
744 // {
745 // cout << "(" << i->first << "=>" << i->second << ")";
746 // cout << ":(" << j->first << "=>" << j->second << ") ";
747 // mismatch=true;
748 // }
749 // }
750 // if (mismatch)
751 // {
752 // cout << "#Mismatch\n";
753 // }
754 // }
755
756 escriptDataC
757 Data::getDataC()
758 {
759 escriptDataC temp;
760 temp.m_dataPtr=(void*)this;
761 return temp;
762 }
763
764 escriptDataC
765 Data::getDataC() const
766 {
767 escriptDataC temp;
768 temp.m_dataPtr=(void*)this;
769 return temp;
770 }
771
772 const boost::python::tuple
773 Data::getShapeTuple() const
774 {
775 const DataTypes::ShapeType& shape=getDataPointShape();
776 switch(getDataPointRank()) {
777 case 0:
778 return make_tuple();
779 case 1:
780 return make_tuple(long_(shape[0]));
781 case 2:
782 return make_tuple(long_(shape[0]),long_(shape[1]));
783 case 3:
784 return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
785 case 4:
786 return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
787 default:
788 throw DataException("Error - illegal Data rank.");
789 }
790 }
791
792
793 // The different name is needed because boost has trouble with overloaded functions.
794 // It can't work out what type the function is based soley on its name.
795 // There are ways to fix this involving creating function pointer variables for each form
796 // but there doesn't seem to be a need given that the methods have the same name from the python point of view
797 Data*
798 Data::copySelf()
799 {
800 DataAbstract* temp=m_data->deepCopy();
801 return new Data(temp);
802 }
803
804 void
805 Data::copy(const Data& other)
806 {
807 DataAbstract* temp=other.m_data->deepCopy();
808 DataAbstract_ptr p=temp->getPtr();
809 m_data=p;
810 }
811
812
813 Data
814 Data::delay()
815 {
816 DataLazy* dl=new DataLazy(m_data);
817 return Data(dl);
818 }
819
820 void
821 Data::delaySelf()
822 {
823 if (!isLazy())
824 {
825 m_data=(new DataLazy(m_data))->getPtr();
826 }
827 }
828
829 void
830 Data::setToZero()
831 {
832 if (isEmpty())
833 {
834 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
835 }
836 m_data->setToZero();
837 }
838
839 void
840 Data::copyWithMask(const Data& other,
841 const Data& mask)
842 {
843 // 1. Interpolate if required so all Datas use the same FS as this
844 // 2. Tag or Expand so that all Data's are the same type
845 // 3. Iterate over the data vectors copying values where mask is >0
846 if (other.isEmpty() || mask.isEmpty())
847 {
848 throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
849 }
850 Data other2(other);
851 Data mask2(mask);
852 other2.resolve();
853 mask2.resolve();
854 this->resolve();
855 FunctionSpace myFS=getFunctionSpace();
856 FunctionSpace oFS=other2.getFunctionSpace();
857 FunctionSpace mFS=mask2.getFunctionSpace();
858 if (oFS!=myFS)
859 {
860 if (other2.probeInterpolation(myFS))
861 {
862 other2=other2.interpolate(myFS);
863 }
864 else
865 {
866 throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
867 }
868 }
869 if (mFS!=myFS)
870 {
871 if (mask2.probeInterpolation(myFS))
872 {
873 mask2=mask2.interpolate(myFS);
874 }
875 else
876 {
877 throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
878 }
879 }
880 // Ensure that all args have the same type
881 if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
882 {
883 this->expand();
884 other2.expand();
885 mask2.expand();
886 }
887 else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
888 {
889 this->tag();
890 other2.tag();
891 mask2.tag();
892 }
893 else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
894 {
895 }
896 else
897 {
898 throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
899 }
900 // Now we iterate over the elements
901 DataVector& self=getReadyPtr()->getVector();
902 const DataVector& ovec=other2.getReadyPtr()->getVector();
903 const DataVector& mvec=mask2.getReadyPtr()->getVector();
904 if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))
905 {
906 throw DataException("Error - size mismatch in arguments to copyWithMask.");
907 }
908 size_t num_points=self.size();
909
910 // OPENMP 3.0 allows unsigned loop vars.
911 #if defined(_OPENMP) && (_OPENMP < 200805)
912 long i;
913 #else
914 size_t i;
915 #endif
916 #pragma omp parallel for private(i) schedule(static)
917 for (i=0;i<num_points;++i)
918 {
919 if (mvec[i]>0)
920 {
921 self[i]=ovec[i];
922 }
923 }
924 }
925
926
927
928 bool
929 Data::isExpanded() const
930 {
931 DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
932 return (temp!=0);
933 }
934
935 bool
936 Data::isTagged() const
937 {
938 DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
939 return (temp!=0);
940 }
941
942 bool
943 Data::isEmpty() const
944 {
945 DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
946 return (temp!=0);
947 }
948
949 bool
950 Data::isConstant() const
951 {
952 DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
953 return (temp!=0);
954 }
955
956 bool
957 Data::isLazy() const
958 {
959 return m_data->isLazy();
960 }
961
962 // at the moment this is synonymous with !isLazy() but that could change
963 bool
964 Data::isReady() const
965 {
966 return (dynamic_cast<DataReady*>(m_data.get())!=0);
967 }
968
969
970 void
971 Data::setProtection()
972 {
973 m_protected=true;
974 }
975
976 bool
977 Data::isProtected() const
978 {
979 return m_protected;
980 }
981
982
983
984 void
985 Data::expand()
986 {
987 if (isConstant()) {
988 DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
989 DataAbstract* temp=new DataExpanded(*tempDataConst);
990 // shared_ptr<DataAbstract> temp_data(temp);
991 // m_data=temp_data;
992 m_data=temp->getPtr();
993 } else if (isTagged()) {
994 DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
995 DataAbstract* temp=new DataExpanded(*tempDataTag);
996 // shared_ptr<DataAbstract> temp_data(temp);
997 // m_data=temp_data;
998 m_data=temp->getPtr();
999 } else if (isExpanded()) {
1000 //
1001 // do nothing
1002 } else if (isEmpty()) {
1003 throw DataException("Error - Expansion of DataEmpty not possible.");
1004 } else if (isLazy()) {
1005 resolve();
1006 expand(); // resolve might not give us expanded data
1007 } else {
1008 throw DataException("Error - Expansion not implemented for this Data type.");
1009 }
1010 }
1011
1012 void
1013 Data::tag()
1014 {
1015 if (isConstant()) {
1016 DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
1017 DataAbstract* temp=new DataTagged(*tempDataConst);
1018 // shared_ptr<DataAbstract> temp_data(temp);
1019 // m_data=temp_data;
1020 m_data=temp->getPtr();
1021 } else if (isTagged()) {
1022 // do nothing
1023 } else if (isExpanded()) {
1024 throw DataException("Error - Creating tag data from DataExpanded not possible.");
1025 } else if (isEmpty()) {
1026 throw DataException("Error - Creating tag data from DataEmpty not possible.");
1027 } else if (isLazy()) {
1028 DataAbstract_ptr res=m_data->resolve();
1029 if (m_data->isExpanded())
1030 {
1031 throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
1032 }
1033 m_data=res;
1034 tag();
1035 } else {
1036 throw DataException("Error - Tagging not implemented for this Data type.");
1037 }
1038 }
1039
1040 void
1041 Data::resolve()
1042 {
1043 if (isLazy())
1044 {
1045 m_data=m_data->resolve();
1046 }
1047 }
1048
1049
1050 Data
1051 Data::oneOver() const
1052 {
1053 if (isLazy())
1054 {
1055 DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);
1056 return Data(c);
1057 }
1058 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
1059 }
1060
1061 Data
1062 Data::wherePositive() const
1063 {
1064 if (isLazy())
1065 {
1066 DataLazy* c=new DataLazy(borrowDataPtr(),GZ);
1067 return Data(c);
1068 }
1069 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
1070 }
1071
1072 Data
1073 Data::whereNegative() const
1074 {
1075 if (isLazy())
1076 {
1077 DataLazy* c=new DataLazy(borrowDataPtr(),LZ);
1078 return Data(c);
1079 }
1080 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
1081 }
1082
1083 Data
1084 Data::whereNonNegative() const
1085 {
1086 if (isLazy())
1087 {
1088 DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);
1089 return Data(c);
1090 }
1091 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
1092 }
1093
1094 Data
1095 Data::whereNonPositive() const
1096 {
1097 if (isLazy())
1098 {
1099 DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);
1100 return Data(c);
1101 }
1102 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
1103 }
1104
1105 Data
1106 Data::whereZero(double tol) const
1107 {
1108 Data dataAbs=abs();
1109 return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1110 }
1111
1112 Data
1113 Data::whereNonZero(double tol) const
1114 {
1115 Data dataAbs=abs();
1116 return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1117 }
1118
1119 Data
1120 Data::interpolate(const FunctionSpace& functionspace) const
1121 {
1122 return Data(*this,functionspace);
1123 }
1124
1125 bool
1126 Data::probeInterpolation(const FunctionSpace& functionspace) const
1127 {
1128 return getFunctionSpace().probeInterpolation(functionspace);
1129 // if (getFunctionSpace()==functionspace) {
1130 // return true;
1131 // } else {
1132 // const_Domain_ptr domain=getDomain();
1133 // if (*domain==*functionspace.getDomain()) {
1134 // return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
1135 // } else {
1136 // return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());
1137 // }
1138 // }
1139 }
1140
1141 Data
1142 Data::gradOn(const FunctionSpace& functionspace) const
1143 {
1144 if (isEmpty())
1145 {
1146 throw DataException("Error - operation not permitted on instances of DataEmpty.");
1147 }
1148 double blocktimer_start = blocktimer_time();
1149 if (functionspace.getDomain()!=getDomain())
1150 throw DataException("Error - gradient cannot be calculated on different domains.");
1151 DataTypes::ShapeType grad_shape=getDataPointShape();
1152 grad_shape.push_back(functionspace.getDim());
1153 Data out(0.0,grad_shape,functionspace,true);
1154 getDomain()->setToGradient(out,*this);
1155 blocktimer_increment("grad()", blocktimer_start);
1156 return out;
1157 }
1158
1159 Data
1160 Data::grad() const
1161 {
1162 if (isEmpty())
1163 {
1164 throw DataException("Error - operation not permitted on instances of DataEmpty.");
1165 }
1166 return gradOn(escript::function(*getDomain()));
1167 }
1168
1169 int
1170 Data::getDataPointSize() const
1171 {
1172 return m_data->getNoValues();
1173 }
1174
1175 DataTypes::ValueType::size_type
1176 Data::getLength() const
1177 {
1178 return m_data->getLength();
1179 }
1180
1181 const
1182 boost::python::numeric::array
1183 Data:: getValueOfDataPoint(int dataPointNo)
1184 {
1185 int i, j, k, l;
1186
1187 FORCERESOLVE;
1188
1189 //
1190 // determine the rank and shape of each data point
1191 int dataPointRank = getDataPointRank();
1192 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1193
1194 //
1195 // create the numeric array to be returned
1196 boost::python::numeric::array numArray(0.0);
1197
1198 //
1199 // the shape of the returned numeric array will be the same
1200 // as that of the data point
1201 int arrayRank = dataPointRank;
1202 const DataTypes::ShapeType& arrayShape = dataPointShape;
1203
1204 //
1205 // resize the numeric array to the shape just calculated
1206 if (arrayRank==0) {
1207 numArray.resize(1);
1208 }
1209 if (arrayRank==1) {
1210 numArray.resize(arrayShape[0]);
1211 }
1212 if (arrayRank==2) {
1213 numArray.resize(arrayShape[0],arrayShape[1]);
1214 }
1215 if (arrayRank==3) {
1216 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
1217 }
1218 if (arrayRank==4) {
1219 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
1220 }
1221
1222 if (getNumDataPointsPerSample()>0) {
1223 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1224 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1225 //
1226 // Check a valid sample number has been supplied
1227 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1228 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
1229 }
1230
1231 //
1232 // Check a valid data point number has been supplied
1233 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1234 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
1235 }
1236 // TODO: global error handling
1237 // create a view of the data if it is stored locally
1238 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
1239 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1240
1241
1242 switch( dataPointRank ){
1243 case 0 :
1244 numArray[0] = getDataAtOffset(offset);
1245 break;
1246 case 1 :
1247 for( i=0; i<dataPointShape[0]; i++ )
1248 numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
1249 break;
1250 case 2 :
1251 for( i=0; i<dataPointShape[0]; i++ )
1252 for( j=0; j<dataPointShape[1]; j++)
1253 numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
1254 break;
1255 case 3 :
1256 for( i=0; i<dataPointShape[0]; i++ )
1257 for( j=0; j<dataPointShape[1]; j++ )
1258 for( k=0; k<dataPointShape[2]; k++)
1259 numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
1260 break;
1261 case 4 :
1262 for( i=0; i<dataPointShape[0]; i++ )
1263 for( j=0; j<dataPointShape[1]; j++ )
1264 for( k=0; k<dataPointShape[2]; k++ )
1265 for( l=0; l<dataPointShape[3]; l++)
1266 numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
1267 break;
1268 }
1269 }
1270 //
1271 // return the array
1272 return numArray;
1273
1274 }
1275
1276 void
1277 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1278 {
1279 // this will throw if the value cannot be represented
1280 boost::python::numeric::array num_array(py_object);
1281 setValueOfDataPointToArray(dataPointNo,num_array);
1282 }
1283
1284 void
1285 Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
1286 {
1287 if (isProtected()) {
1288 throw DataException("Error - attempt to update protected Data object.");
1289 }
1290 FORCERESOLVE;
1291 //
1292 // check rank
1293 if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())
1294 throw DataException("Rank of numarray does not match Data object rank");
1295
1296 //
1297 // check shape of num_array
1298 for (unsigned int i=0; i<getDataPointRank(); i++) {
1299 if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
1300 throw DataException("Shape of numarray does not match Data object rank");
1301 }
1302 //
1303 // make sure data is expanded:
1304 //
1305 if (!isExpanded()) {
1306 expand();
1307 }
1308 if (getNumDataPointsPerSample()>0) {
1309 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1310 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1311 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
1312 } else {
1313 m_data->copyToDataPoint(-1, 0,num_array);
1314 }
1315 }
1316
1317 void
1318 Data::setValueOfDataPoint(int dataPointNo, const double value)
1319 {
1320 if (isProtected()) {
1321 throw DataException("Error - attempt to update protected Data object.");
1322 }
1323 //
1324 // make sure data is expanded:
1325 FORCERESOLVE;
1326 if (!isExpanded()) {
1327 expand();
1328 }
1329 if (getNumDataPointsPerSample()>0) {
1330 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1331 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1332 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
1333 } else {
1334 m_data->copyToDataPoint(-1, 0,value);
1335 }
1336 }
1337
1338 const
1339 boost::python::numeric::array
1340 Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
1341 {
1342 size_t length=0;
1343 int i, j, k, l, pos;
1344 FORCERESOLVE;
1345 //
1346 // determine the rank and shape of each data point
1347 int dataPointRank = getDataPointRank();
1348 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1349
1350 //
1351 // create the numeric array to be returned
1352 boost::python::numeric::array numArray(0.0);
1353
1354 //
1355 // the shape of the returned numeric array will be the same
1356 // as that of the data point
1357 int arrayRank = dataPointRank;
1358 const DataTypes::ShapeType& arrayShape = dataPointShape;
1359
1360 //
1361 // resize the numeric array to the shape just calculated
1362 if (arrayRank==0) {
1363 numArray.resize(1);
1364 }
1365 if (arrayRank==1) {
1366 numArray.resize(arrayShape[0]);
1367 }
1368 if (arrayRank==2) {
1369 numArray.resize(arrayShape[0],arrayShape[1]);
1370 }
1371 if (arrayRank==3) {
1372 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
1373 }
1374 if (arrayRank==4) {
1375 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
1376 }
1377
1378 // added for the MPI communication
1379 length=1;
1380 for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
1381 double *tmpData = new double[length];
1382
1383 //
1384 // load the values for the data point into the numeric array.
1385
1386 // updated for the MPI case
1387 if( get_MPIRank()==procNo ){
1388 if (getNumDataPointsPerSample()>0) {
1389 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1390 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1391 //
1392 // Check a valid sample number has been supplied
1393 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1394 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
1395 }
1396
1397 //
1398 // Check a valid data point number has been supplied
1399 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1400 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
1401 }
1402 // TODO: global error handling
1403 // create a view of the data if it is stored locally
1404 //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
1405 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1406
1407 // pack the data from the view into tmpData for MPI communication
1408 pos=0;
1409 switch( dataPointRank ){
1410 case 0 :
1411 tmpData[0] = getDataAtOffset(offset);
1412 break;
1413 case 1 :
1414 for( i=0; i<dataPointShape[0]; i++ )
1415 tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
1416 break;
1417 case 2 :
1418 for( i=0; i<dataPointShape[0]; i++ )
1419 for( j=0; j<dataPointShape[1]; j++, pos++ )
1420 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
1421 break;
1422 case 3 :
1423 for( i=0; i<dataPointShape[0]; i++ )
1424 for( j=0; j<dataPointShape[1]; j++ )
1425 for( k=0; k<dataPointShape[2]; k++, pos++ )
1426 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
1427 break;
1428 case 4 :
1429 for( i=0; i<dataPointShape[0]; i++ )
1430 for( j=0; j<dataPointShape[1]; j++ )
1431 for( k=0; k<dataPointShape[2]; k++ )
1432 for( l=0; l<dataPointShape[3]; l++, pos++ )
1433 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
1434 break;
1435 }
1436 }
1437 }
1438 #ifdef PASO_MPI
1439 // broadcast the data to all other processes
1440 MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1441 #endif
1442
1443 // unpack the data
1444 switch( dataPointRank ){
1445 case 0 :
1446 numArray[0]=tmpData[0];
1447 break;
1448 case 1 :
1449 for( i=0; i<dataPointShape[0]; i++ )
1450 numArray[i]=tmpData[i];
1451 break;
1452 case 2 :
1453 for( i=0; i<dataPointShape[0]; i++ )
1454 for( j=0; j<dataPointShape[1]; j++ )
1455 numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
1456 break;
1457 case 3 :
1458 for( i=0; i<dataPointShape[0]; i++ )
1459 for( j=0; j<dataPointShape[1]; j++ )
1460 for( k=0; k<dataPointShape[2]; k++ )
1461 numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
1462 break;
1463 case 4 :
1464 for( i=0; i<dataPointShape[0]; i++ )
1465 for( j=0; j<dataPointShape[1]; j++ )
1466 for( k=0; k<dataPointShape[2]; k++ )
1467 for( l=0; l<dataPointShape[3]; l++ )
1468 numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
1469 break;
1470 }
1471
1472 delete [] tmpData;
1473 //
1474 // return the loaded array
1475 return numArray;
1476 }
1477
1478
1479 boost::python::numeric::array
1480 Data::integrate_const() const
1481 {
1482 if (isLazy())
1483 {
1484 throw DataException("Error - cannot integrate for constant lazy data.");
1485 }
1486 return integrateWorker();
1487 }
1488
1489 boost::python::numeric::array
1490 Data::integrate()
1491 {
1492 if (isLazy())
1493 {
1494 expand();
1495 }
1496 return integrateWorker();
1497 }
1498
1499
1500
1501 boost::python::numeric::array
1502 Data::integrateWorker() const
1503 {
1504 int index;
1505 int rank = getDataPointRank();
1506 DataTypes::ShapeType shape = getDataPointShape();
1507 int dataPointSize = getDataPointSize();
1508
1509 //
1510 // calculate the integral values
1511 vector<double> integrals(dataPointSize);
1512 vector<double> integrals_local(dataPointSize);
1513 const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1514 if (dom==0)
1515 {
1516 throw DataException("Can not integrate over non-continuous domains.");
1517 }
1518 #ifdef PASO_MPI
1519 dom->setToIntegrals(integrals_local,*this);
1520 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1521 double *tmp = new double[dataPointSize];
1522 double *tmp_local = new double[dataPointSize];
1523 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1524 MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1525 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1526 delete[] tmp;
1527 delete[] tmp_local;
1528 #else
1529 dom->setToIntegrals(integrals,*this);
1530 #endif
1531
1532 //
1533 // create the numeric array to be returned
1534 // and load the array with the integral values
1535 boost::python::numeric::array bp_array(1.0);
1536 if (rank==0) {
1537 bp_array.resize(1);
1538 index = 0;
1539 bp_array[0] = integrals[index];
1540 }
1541 if (rank==1) {
1542 bp_array.resize(shape[0]);
1543 for (int i=0; i<shape[0]; i++) {
1544 index = i;
1545 bp_array[i] = integrals[index];
1546 }
1547 }
1548 if (rank==2) {
1549 bp_array.resize(shape[0],shape[1]);
1550 for (int i=0; i<shape[0]; i++) {
1551 for (int j=0; j<shape[1]; j++) {
1552 index = i + shape[0] * j;
1553 bp_array[make_tuple(i,j)] = integrals[index];
1554 }
1555 }
1556 }
1557 if (rank==3) {
1558 bp_array.resize(shape[0],shape[1],shape[2]);
1559 for (int i=0; i<shape[0]; i++) {
1560 for (int j=0; j<shape[1]; j++) {
1561 for (int k=0; k<shape[2]; k++) {
1562 index = i + shape[0] * ( j + shape[1] * k );
1563 bp_array[make_tuple(i,j,k)] = integrals[index];
1564 }
1565 }
1566 }
1567 }
1568 if (rank==4) {
1569 bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1570 for (int i=0; i<shape[0]; i++) {
1571 for (int j=0; j<shape[1]; j++) {
1572 for (int k=0; k<shape[2]; k++) {
1573 for (int l=0; l<shape[3]; l++) {
1574 index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1575 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1576 }
1577 }
1578 }
1579 }
1580 }
1581
1582 //
1583 // return the loaded array
1584 return bp_array;
1585 }
1586
1587 Data
1588 Data::sin() const
1589 {
1590 if (isLazy())
1591 {
1592 DataLazy* c=new DataLazy(borrowDataPtr(),SIN);
1593 return Data(c);
1594 }
1595 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1596 }
1597
1598 Data
1599 Data::cos() const
1600 {
1601 if (isLazy())
1602 {
1603 DataLazy* c=new DataLazy(borrowDataPtr(),COS);
1604 return Data(c);
1605 }
1606 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1607 }
1608
1609 Data
1610 Data::tan() const
1611 {
1612 if (isLazy())
1613 {
1614 DataLazy* c=new DataLazy(borrowDataPtr(),TAN);
1615 return Data(c);
1616 }
1617 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1618 }
1619
1620 Data
1621 Data::asin() const
1622 {
1623 if (isLazy())
1624 {
1625 DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);
1626 return Data(c);
1627 }
1628 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1629 }
1630
1631 Data
1632 Data::acos() const
1633 {
1634 if (isLazy())
1635 {
1636 DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);
1637 return Data(c);
1638 }
1639 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1640 }
1641
1642
1643 Data
1644 Data::atan() const
1645 {
1646 if (isLazy())
1647 {
1648 DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);
1649 return Data(c);
1650 }
1651 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1652 }
1653
1654 Data
1655 Data::sinh() const
1656 {
1657 if (isLazy())
1658 {
1659 DataLazy* c=new DataLazy(borrowDataPtr(),SINH);
1660 return Data(c);
1661 }
1662 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1663 }
1664
1665 Data
1666 Data::cosh() const
1667 {
1668 if (isLazy())
1669 {
1670 DataLazy* c=new DataLazy(borrowDataPtr(),COSH);
1671 return Data(c);
1672 }
1673 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1674 }
1675
1676 Data
1677 Data::tanh() const
1678 {
1679 if (isLazy())
1680 {
1681 DataLazy* c=new DataLazy(borrowDataPtr(),TANH);
1682 return Data(c);
1683 }
1684 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1685 }
1686
1687
1688 Data
1689 Data::erf() const
1690 {
1691 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1692 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1693 #else
1694 if (isLazy())
1695 {
1696 DataLazy* c=new DataLazy(borrowDataPtr(),ERF);
1697 return Data(c);
1698 }
1699 return C_TensorUnaryOperation(*this, ::erf);
1700 #endif
1701 }
1702
1703 Data
1704 Data::asinh() const
1705 {
1706 if (isLazy())
1707 {
1708 DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);
1709 return Data(c);
1710 }
1711 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1712 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1713 #else
1714 return C_TensorUnaryOperation(*this, ::asinh);
1715 #endif
1716 }
1717
1718 Data
1719 Data::acosh() const
1720 {
1721 if (isLazy())
1722 {
1723 DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);
1724 return Data(c);
1725 }
1726 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1727 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1728 #else
1729 return C_TensorUnaryOperation(*this, ::acosh);
1730 #endif
1731 }
1732
1733 Data
1734 Data::atanh() const
1735 {
1736 if (isLazy())
1737 {
1738 DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);
1739 return Data(c);
1740 }
1741 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1742 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1743 #else
1744 return C_TensorUnaryOperation(*this, ::atanh);
1745 #endif
1746 }
1747
1748 Data
1749 Data::log10() const
1750 { if (isLazy())
1751 {
1752 DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);
1753 return Data(c);
1754 }
1755 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1756 }
1757
1758 Data
1759 Data::log() const
1760 {
1761 if (isLazy())
1762 {
1763 DataLazy* c=new DataLazy(borrowDataPtr(),LOG);
1764 return Data(c);
1765 }
1766 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1767 }
1768
1769 Data
1770 Data::sign() const
1771 {
1772 if (isLazy())
1773 {
1774 DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);
1775 return Data(c);
1776 }
1777 return C_TensorUnaryOperation(*this, escript::fsign);
1778 }
1779
1780 Data
1781 Data::abs() const
1782 {
1783 if (isLazy())
1784 {
1785 DataLazy* c=new DataLazy(borrowDataPtr(),ABS);
1786 return Data(c);
1787 }
1788 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1789 }
1790
1791 Data
1792 Data::neg() const
1793 {
1794 if (isLazy())
1795 {
1796 DataLazy* c=new DataLazy(borrowDataPtr(),NEG);
1797 return Data(c);
1798 }
1799 return C_TensorUnaryOperation(*this, negate<double>());
1800 }
1801
1802 Data
1803 Data::pos() const
1804 {
1805 // not doing lazy check here is deliberate.
1806 // since a deep copy of lazy data should be cheap, I'll just let it happen now
1807 Data result;
1808 // perform a deep copy
1809 result.copy(*this);
1810 return result;
1811 }
1812
1813 Data
1814 Data::exp() const
1815 {
1816 if (isLazy())
1817 {
1818 DataLazy* c=new DataLazy(borrowDataPtr(),EXP);
1819 return Data(c);
1820 }
1821 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1822 }
1823
1824 Data
1825 Data::sqrt() const
1826 {
1827 if (isLazy())
1828 {
1829 DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);
1830 return Data(c);
1831 }
1832 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1833 }
1834
1835 double
1836 Data::Lsup_const() const
1837 {
1838 if (isLazy())
1839 {
1840 throw DataException("Error - cannot compute Lsup for constant lazy data.");
1841 }
1842 return LsupWorker();
1843 }
1844
1845 double
1846 Data::Lsup()
1847 {
1848 if (isLazy())
1849 {
1850 resolve();
1851 }
1852 return LsupWorker();
1853 }
1854
1855 double
1856 Data::sup_const() const
1857 {
1858 if (isLazy())
1859 {
1860 throw DataException("Error - cannot compute sup for constant lazy data.");
1861 }
1862 return supWorker();
1863 }
1864
1865 double
1866 Data::sup()
1867 {
1868 if (isLazy())
1869 {
1870 resolve();
1871 }
1872 return supWorker();
1873 }
1874
1875 double
1876 Data::inf_const() const
1877 {
1878 if (isLazy())
1879 {
1880 throw DataException("Error - cannot compute inf for constant lazy data.");
1881 }
1882 return infWorker();
1883 }
1884
1885 double
1886 Data::inf()
1887 {
1888 if (isLazy())
1889 {
1890 resolve();
1891 }
1892 return infWorker();
1893 }
1894
1895 double
1896 Data::LsupWorker() const
1897 {
1898 double localValue;
1899 //
1900 // set the initial absolute maximum value to zero
1901
1902 AbsMax abs_max_func;
1903 localValue = algorithm(abs_max_func,0);
1904 #ifdef PASO_MPI
1905 double globalValue;
1906 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1907 return globalValue;
1908 #else
1909 return localValue;
1910 #endif
1911 }
1912
1913 double
1914 Data::supWorker() const
1915 {
1916 double localValue;
1917 //
1918 // set the initial maximum value to min possible double
1919 FMax fmax_func;
1920 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1921 #ifdef PASO_MPI
1922 double globalValue;
1923 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1924 return globalValue;
1925 #else
1926 return localValue;
1927 #endif
1928 }
1929
1930 double
1931 Data::infWorker() const
1932 {
1933 double localValue;
1934 //
1935 // set the initial minimum value to max possible double
1936 FMin fmin_func;
1937 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1938 #ifdef PASO_MPI
1939 double globalValue;
1940 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1941 return globalValue;
1942 #else
1943 return localValue;
1944 #endif
1945 }
1946
1947 /* TODO */
1948 /* global reduction */
1949 Data
1950 Data::maxval() const
1951 {
1952 if (isLazy())
1953 {
1954 Data temp(*this); // to get around the fact that you can't resolve a const Data
1955 temp.resolve();
1956 return temp.maxval();
1957 }
1958 //
1959 // set the initial maximum value to min possible double
1960 FMax fmax_func;
1961 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1962 }
1963
1964 Data
1965 Data::minval() const
1966 {
1967 if (isLazy())
1968 {
1969 Data temp(*this); // to get around the fact that you can't resolve a const Data
1970 temp.resolve();
1971 return temp.minval();
1972 }
1973 //
1974 // set the initial minimum value to max possible double
1975 FMin fmin_func;
1976 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1977 }
1978
1979 Data
1980 Data::swapaxes(const int axis0, const int axis1) const
1981 {
1982 if (isLazy())
1983 {
1984 Data temp(*this);
1985 temp.resolve();
1986 return temp.swapaxes(axis0,axis1);
1987 }
1988 int axis0_tmp,axis1_tmp;
1989 DataTypes::ShapeType s=getDataPointShape();
1990 DataTypes::ShapeType ev_shape;
1991 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1992 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1993 int rank=getDataPointRank();
1994 if (rank<2) {
1995 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1996 }
1997 if (axis0<0 || axis0>rank-1) {
1998 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1999 }
2000 if (axis1<0 || axis1>rank-1) {
2001 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
2002 }
2003 if (axis0 == axis1) {
2004 throw DataException("Error - Data::swapaxes: axis indices must be different.");
2005 }
2006 if (axis0 > axis1) {
2007 axis0_tmp=axis1;
2008 axis1_tmp=axis0;
2009 } else {
2010 axis0_tmp=axis0;
2011 axis1_tmp=axis1;
2012 }
2013 for (int i=0; i<rank; i++) {
2014 if (i == axis0_tmp) {
2015 ev_shape.push_back(s[axis1_tmp]);
2016 } else if (i == axis1_tmp) {
2017 ev_shape.push_back(s[axis0_tmp]);
2018 } else {
2019 ev_shape.push_back(s[i]);
2020 }
2021 }
2022 Data ev(0.,ev_shape,getFunctionSpace());
2023 ev.typeMatchRight(*this);
2024 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
2025 return ev;
2026
2027 }
2028
2029 Data
2030 Data::symmetric() const
2031 {
2032 // check input
2033 DataTypes::ShapeType s=getDataPointShape();
2034 if (getDataPointRank()==2) {
2035 if(s[0] != s[1])
2036 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
2037 }
2038 else if (getDataPointRank()==4) {
2039 if(!(s[0] == s[2] && s[1] == s[3]))
2040 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
2041 }
2042 else {
2043 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
2044 }
2045 if (isLazy())
2046 {
2047 DataLazy* c=new DataLazy(borrowDataPtr(),SYM);
2048 return Data(c);
2049 }
2050 Data ev(0.,getDataPointShape(),getFunctionSpace());
2051 ev.typeMatchRight(*this);
2052 m_data->symmetric(ev.m_data.get());
2053 return ev;
2054 }
2055
2056 Data
2057 Data::nonsymmetric() const
2058 {
2059 if (isLazy())
2060 {
2061 DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);
2062 return Data(c);
2063 }
2064 // check input
2065 DataTypes::ShapeType s=getDataPointShape();
2066 if (getDataPointRank()==2) {
2067 if(s[0] != s[1])
2068 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
2069 DataTypes::ShapeType ev_shape;
2070 ev_shape.push_back(s[0]);
2071 ev_shape.push_back(s[1]);
2072 Data ev(0.,ev_shape,getFunctionSpace());
2073 ev.typeMatchRight(*this);
2074 m_data->nonsymmetric(ev.m_data.get());
2075 return ev;
2076 }
2077 else if (getDataPointRank()==4) {
2078 if(!(s[0] == s[2] && s[1] == s[3]))
2079 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
2080 DataTypes::ShapeType ev_shape;
2081 ev_shape.push_back(s[0]);
2082 ev_shape.push_back(s[1]);
2083 ev_shape.push_back(s[2]);
2084 ev_shape.push_back(s[3]);
2085 Data ev(0.,ev_shape,getFunctionSpace());
2086 ev.typeMatchRight(*this);
2087 m_data->nonsymmetric(ev.m_data.get());
2088 return ev;
2089 }
2090 else {
2091 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
2092 }
2093 }
2094
2095 Data
2096 Data::trace(int axis_offset) const
2097 {
2098 if (isLazy())
2099 {
2100 DataLazy* c=new DataLazy(borrowDataPtr(),TRACE,axis_offset);
2101 return Data(c);
2102 }
2103 DataTypes::ShapeType s=getDataPointShape();
2104 if (getDataPointRank()==2) {
2105 DataTypes::ShapeType ev_shape;
2106 Data ev(0.,ev_shape,getFunctionSpace());
2107 ev.typeMatchRight(*this);
2108 m_data->trace(ev.m_data.get(), axis_offset);
2109 return ev;
2110 }
2111 if (getDataPointRank()==3) {
2112 DataTypes::ShapeType ev_shape;
2113 if (axis_offset==0) {
2114 int s2=s[2];
2115 ev_shape.push_back(s2);
2116 }
2117 else if (axis_offset==1) {
2118 int s0=s[0];
2119 ev_shape.push_back(s0);
2120 }
2121 Data ev(0.,ev_shape,getFunctionSpace());
2122 ev.typeMatchRight(*this);
2123 m_data->trace(ev.m_data.get(), axis_offset);
2124 return ev;
2125 }
2126 if (getDataPointRank()==4) {
2127 DataTypes::ShapeType ev_shape;
2128 if (axis_offset==0) {
2129 ev_shape.push_back(s[2]);
2130 ev_shape.push_back(s[3]);
2131 }
2132 else if (axis_offset==1) {
2133 ev_shape.push_back(s[0]);
2134 ev_shape.push_back(s[3]);
2135 }
2136 else if (axis_offset==2) {
2137 ev_shape.push_back(s[0]);
2138 ev_shape.push_back(s[1]);
2139 }
2140 Data ev(0.,ev_shape,getFunctionSpace());
2141 ev.typeMatchRight(*this);
2142 m_data->trace(ev.m_data.get(), axis_offset);
2143 return ev;
2144 }
2145 else {
2146 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
2147 }
2148 }
2149
2150 Data
2151 Data::transpose(int axis_offset) const
2152 {
2153 if (isLazy())
2154 {
2155 DataLazy* c=new DataLazy(borrowDataPtr(),TRANS,axis_offset);
2156 return Data(c);
2157 }
2158 DataTypes::ShapeType s=getDataPointShape();
2159 DataTypes::ShapeType ev_shape;
2160 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
2161 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
2162 int rank=getDataPointRank();
2163 if (axis_offset<0 || axis_offset>rank) {
2164 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
2165 }
2166 for (int i=0; i<rank; i++) {
2167 int index = (axis_offset+i)%rank;
2168 ev_shape.push_back(s[index]); // Append to new shape
2169 }
2170 Data ev(0.,ev_shape,getFunctionSpace());
2171 ev.typeMatchRight(*this);
2172 m_data->transpose(ev.m_data.get(), axis_offset);
2173 return ev;
2174 }
2175
2176 Data
2177 Data::eigenvalues() const
2178 {
2179 if (isLazy())
2180 {
2181 Data temp(*this); // to get around the fact that you can't resolve a const Data
2182 temp.resolve();
2183 return temp.eigenvalues();
2184 }
2185 // check input
2186 DataTypes::ShapeType s=getDataPointShape();
2187 if (getDataPointRank()!=2)
2188 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
2189 if(s[0] != s[1])
2190 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
2191 // create return
2192 DataTypes::ShapeType ev_shape(1,s[0]);
2193 Data ev(0.,ev_shape,getFunctionSpace());
2194 ev.typeMatchRight(*this);
2195 m_data->eigenvalues(ev.m_data.get());
2196 return ev;
2197 }
2198
2199 const boost::python::tuple
2200 Data::eigenvalues_and_eigenvectors(const double tol) const
2201 {
2202 if (isLazy())
2203 {
2204 Data temp(*this); // to get around the fact that you can't resolve a const Data
2205 temp.resolve();
2206 return temp.eigenvalues_and_eigenvectors(tol);
2207 }
2208 DataTypes::ShapeType s=getDataPointShape();
2209 if (getDataPointRank()!=2)
2210 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
2211 if(s[0] != s[1])
2212 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
2213 // create return
2214 DataTypes::ShapeType ev_shape(1,s[0]);
2215 Data ev(0.,ev_shape,getFunctionSpace());
2216 ev.typeMatchRight(*this);
2217 DataTypes::ShapeType V_shape(2,s[0]);
2218 Data V(0.,V_shape,getFunctionSpace());
2219 V.typeMatchRight(*this);
2220 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
2221 return make_tuple(boost::python::object(ev),boost::python::object(V));
2222 }
2223
2224 const boost::python::tuple
2225 Data::minGlobalDataPoint() const
2226 {
2227 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
2228 // abort (for unknown reasons) if there are openmp directives with it in the
2229 // surrounding function
2230
2231 int DataPointNo;
2232 int ProcNo;
2233 calc_minGlobalDataPoint(ProcNo,DataPointNo);
2234 return make_tuple(ProcNo,DataPointNo);
2235 }
2236
2237 void
2238 Data::calc_minGlobalDataPoint(int& ProcNo,
2239 int& DataPointNo) const
2240 {
2241 if (isLazy())
2242 {
2243 Data temp(*this); // to get around the fact that you can't resolve a const Data
2244 temp.resolve();
2245 return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
2246 }
2247 int i,j;
2248 int lowi=0,lowj=0;
2249 double min=numeric_limits<double>::max();
2250
2251 Data temp=minval();
2252
2253 int numSamples=temp.getNumSamples();
2254 int numDPPSample=temp.getNumDataPointsPerSample();
2255
2256 double next,local_min;
2257 int local_lowi=0,local_lowj=0;
2258
2259 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
2260 {
2261 local_min=min;
2262 #pragma omp for private(i,j) schedule(static)
2263 for (i=0; i<numSamples; i++) {
2264 for (j=0; j<numDPPSample; j++) {
2265 next=temp.getDataAtOffset(temp.getDataOffset(i,j));
2266 if (next<local_min) {
2267 local_min=next;
2268 local_lowi=i;
2269 local_lowj=j;
2270 }
2271 }
2272 }
2273 #pragma omp critical
2274 if (local_min<min) {
2275 min=local_min;
2276 lowi=local_lowi;
2277 lowj=local_lowj;
2278 }
2279 }
2280
2281 #ifdef PASO_MPI
2282 // determine the processor on which the minimum occurs
2283 next = temp.getDataPoint(lowi,lowj);
2284 int lowProc = 0;
2285 double *globalMins = new double[get_MPISize()+1];
2286 int error;
2287 error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
2288
2289 if( get_MPIRank()==0 ){
2290 next = globalMins[lowProc];
2291 for( i=1; i<get_MPISize(); i++ )
2292 if( next>globalMins[i] ){
2293 lowProc = i;
2294 next = globalMins[i];
2295 }
2296 }
2297 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
2298
2299 delete [] globalMins;
2300 ProcNo = lowProc;
2301 #else
2302 ProcNo = 0;
2303 #endif
2304 DataPointNo = lowj + lowi * numDPPSample;
2305 }
2306
2307 void
2308 Data::saveDX(std::string fileName) const
2309 {
2310 if (isEmpty())
2311 {
2312 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2313 }
2314 if (isLazy())
2315 {
2316 Data temp(*this); // to get around the fact that you can't resolve a const Data
2317 temp.resolve();
2318 temp.saveDX(fileName);
2319 return;
2320 }
2321 boost::python::dict args;
2322 args["data"]=boost::python::object(this);
2323 getDomain()->saveDX(fileName,args);
2324 return;
2325 }
2326
2327 void
2328 Data::saveVTK(std::string fileName) const
2329 {
2330 if (isEmpty())
2331 {
2332 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2333 }
2334 if (isLazy())
2335 {
2336 Data temp(*this); // to get around the fact that you can't resolve a const Data
2337 temp.resolve();
2338 temp.saveVTK(fileName);
2339 return;
2340 }
2341 boost::python::dict args;
2342 args["data"]=boost::python::object(this);
2343 getDomain()->saveVTK(fileName,args);
2344 return;
2345 }
2346
2347 Data&
2348 Data::operator+=(const Data& right)
2349 {
2350 if (isProtected()) {
2351 throw DataException("Error - attempt to update protected Data object.");
2352 }
2353 if (isLazy() || right.isLazy())
2354 {
2355 DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=
2356 m_data=c->getPtr();
2357 return (*this);
2358 }
2359 else
2360 {
2361 binaryOp(right,plus<double>());
2362 return (*this);
2363 }
2364 }
2365
2366 Data&
2367 Data::operator+=(const boost::python::object& right)
2368 {
2369 if (isProtected()) {
2370 throw DataException("Error - attempt to update protected Data object.");
2371 }
2372 Data tmp(right,getFunctionSpace(),false);
2373 if (isLazy())
2374 {
2375 DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),ADD); // for lazy + is equivalent to +=
2376 m_data=c->getPtr();
2377 return (*this);
2378 }
2379 else
2380 {
2381 binaryOp(tmp,plus<double>());
2382 return (*this);
2383 }
2384 }
2385
2386 // Hmmm, operator= makes a deep copy but the copy constructor does not?
2387 Data&
2388 Data::operator=(const Data& other)
2389 {
2390 copy(other);
2391 return (*this);
2392 }
2393
2394 Data&
2395 Data::operator-=(const Data& right)
2396 {
2397 if (isProtected()) {
2398 throw DataException("Error - attempt to update protected Data object.");
2399 }
2400 if (isLazy() || right.isLazy())
2401 {
2402 DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=
2403 m_data=c->getPtr();
2404 return (*this);
2405 }
2406 else
2407 {
2408 binaryOp(right,minus<double>());
2409 return (*this);
2410 }
2411 }
2412
2413 Data&
2414 Data::operator-=(const boost::python::object& right)
2415 {
2416 if (isProtected()) {
2417 throw DataException("Error - attempt to update protected Data object.");
2418 }
2419 Data tmp(right,getFunctionSpace(),false);
2420 if (isLazy())
2421 {
2422 DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),SUB); // for lazy - is equivalent to -=
2423 m_data=c->getPtr();
2424 return (*this);
2425 }
2426 else
2427 {
2428 binaryOp(tmp,minus<double>());
2429 return (*this);
2430 }
2431 }
2432
2433 Data&
2434 Data::operator*=(const Data& right)
2435 {
2436 if (isProtected()) {
2437 throw DataException("Error - attempt to update protected Data object.");
2438 }
2439 if (isLazy() || right.isLazy())
2440 {
2441 DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=
2442 m_data=c->getPtr();
2443 return (*this);
2444 }
2445 else
2446 {
2447 binaryOp(right,multiplies<double>());
2448 return (*this);
2449 }
2450 }
2451
2452 Data&
2453 Data::operator*=(const boost::python::object& right)
2454 {
2455 if (isProtected()) {
2456 throw DataException("Error - attempt to update protected Data object.");
2457 }
2458 Data tmp(right,getFunctionSpace(),false);
2459 if (isLazy())
2460 {
2461 DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),MUL); // for lazy * is equivalent to *=
2462 m_data=c->getPtr();
2463 return (*this);
2464 }
2465 else
2466 {
2467 binaryOp(tmp,multiplies<double>());
2468 return (*this);
2469 }
2470 }
2471
2472 Data&
2473 Data::operator/=(const Data& right)
2474 {
2475 if (isProtected()) {
2476 throw DataException("Error - attempt to update protected Data object.");
2477 }
2478 if (isLazy() || right.isLazy())
2479 {
2480 DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=
2481 m_data=c->getPtr();
2482 return (*this);
2483 }
2484 else
2485 {
2486 binaryOp(right,divides<double>());
2487 return (*this);
2488 }
2489 }
2490
2491 Data&
2492 Data::operator/=(const boost::python::object& right)
2493 {
2494 if (isProtected()) {
2495 throw DataException("Error - attempt to update protected Data object.");
2496 }
2497 Data tmp(right,getFunctionSpace(),false);
2498 if (isLazy())
2499 {
2500 DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),DIV); // for lazy / is equivalent to /=
2501 m_data=c->getPtr();
2502 return (*this);
2503 }
2504 else
2505 {
2506 binaryOp(tmp,divides<double>());
2507 return (*this);
2508 }
2509 }
2510
2511 Data
2512 Data::rpowO(const boost::python::object& left) const
2513 {
2514 Data left_d(left,*this);
2515 return left_d.powD(*this);
2516 }
2517
2518 Data
2519 Data::powO(const boost::python::object& right) const
2520 {
2521 Data tmp(right,getFunctionSpace(),false);
2522 return powD(tmp);
2523 }
2524
2525 Data
2526 Data::powD(const Data& right) const
2527 {
2528 if (isLazy() || right.isLazy())
2529 {
2530 DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);
2531 return Data(c);
2532 }
2533 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2534 }
2535
2536 //
2537 // NOTE: It is essential to specify the namespace this operator belongs to
2538 Data
2539 escript::operator+(const Data& left, const Data& right)
2540 {
2541 if (left.isLazy() || right.isLazy())
2542 {
2543 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);
2544 return Data(c);
2545 }
2546 return C_TensorBinaryOperation(left, right, plus<double>());
2547 }
2548
2549 //
2550 // NOTE: It is essential to specify the namespace this operator belongs to
2551 Data
2552 escript::operator-(const Data& left, const Data& right)
2553 {
2554 if (left.isLazy() || right.isLazy())
2555 {
2556 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);
2557 return Data(c);
2558 }
2559 return C_TensorBinaryOperation(left, right, minus<double>());
2560 }
2561
2562 //
2563 // NOTE: It is essential to specify the namespace this operator belongs to
2564 Data
2565 escript::operator*(const Data& left, const Data& right)
2566 {
2567 if (left.isLazy() || right.isLazy())
2568 {
2569 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);
2570 return Data(c);
2571 }
2572 return C_TensorBinaryOperation(left, right, multiplies<double>());
2573 }
2574
2575 //
2576 // NOTE: It is essential to specify the namespace this operator belongs to
2577 Data
2578 escript::operator/(const Data& left, const Data& right)
2579 {
2580 if (left.isLazy() || right.isLazy())
2581 {
2582 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);
2583 return Data(c);
2584 }
2585 return C_TensorBinaryOperation(left, right, divides<double>());
2586 }
2587
2588 //
2589 // NOTE: It is essential to specify the namespace this operator belongs to
2590 Data
2591 escript::operator+(const Data& left, const boost::python::object& right)
2592 {
2593 if (left.isLazy())
2594 {
2595 DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),ADD);
2596 return Data(c);
2597 }
2598 return left+Data(right,left.getFunctionSpace(),false);
2599 }
2600
2601 //
2602 // NOTE: It is essential to specify the namespace this operator belongs to
2603 Data
2604 escript::operator-(const Data& left, const boost::python::object& right)
2605 {
2606 if (left.isLazy())
2607 {
2608 DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),SUB);
2609 return Data(c);
2610 }
2611 return left-Data(right,left.getFunctionSpace(),false);
2612 }
2613
2614 //
2615 // NOTE: It is essential to specify the namespace this operator belongs to
2616 Data
2617 escript::operator*(const Data& left, const boost::python::object& right)
2618 {
2619 if (left.isLazy())
2620 {
2621 DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),MUL);
2622 return Data(c);
2623 }
2624 return left*Data(right,left.getFunctionSpace(),false);
2625 }
2626
2627 //
2628 // NOTE: It is essential to specify the namespace this operator belongs to
2629 Data
2630 escript::operator/(const Data& left, const boost::python::object& right)
2631 {
2632 if (left.isLazy())
2633 {
2634 DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),DIV);
2635 return Data(c);
2636 }
2637 return left/Data(right,left.getFunctionSpace(),false);
2638 }
2639
2640 //
2641 // NOTE: It is essential to specify the namespace this operator belongs to
2642 Data
2643 escript::operator+(const boost::python::object& left, const Data& right)
2644 {
2645 if (right.isLazy())
2646 {
2647 DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),ADD);
2648 return Data(c);
2649 }
2650 return Data(left,right.getFunctionSpace(),false)+right;
2651 }
2652
2653 //
2654 // NOTE: It is essential to specify the namespace this operator belongs to
2655 Data
2656 escript::operator-(const boost::python::object& left, const Data& right)
2657 {
2658 if (right.isLazy())
2659 {
2660 DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),SUB);
2661 return Data(c);
2662 }
2663 return Data(left,right.getFunctionSpace(),false)-right;
2664 }
2665
2666 //
2667 // NOTE: It is essential to specify the namespace this operator belongs to
2668 Data
2669 escript::operator*(const boost::python::object& left, const Data& right)
2670 {
2671 if (right.isLazy())
2672 {
2673 DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),MUL);
2674 return Data(c);
2675 }
2676 return Data(left,right.getFunctionSpace(),false)*right;
2677 }
2678
2679 //
2680 // NOTE: It is essential to specify the namespace this operator belongs to
2681 Data
2682 escript::operator/(const boost::python::object& left, const Data& right)
2683 {
2684 if (right.isLazy())
2685 {
2686 DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),DIV);
2687 return Data(c);
2688 }
2689 return Data(left,right.getFunctionSpace(),false)/right;
2690 }
2691
2692
2693 /* TODO */
2694 /* global reduction */
2695 Data
2696 Data::getItem(const boost::python::object& key) const
2697 {
2698
2699 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2700
2701 if (slice_region.size()!=getDataPointRank()) {
2702 throw DataException("Error - slice size does not match Data rank.");
2703 }
2704
2705 return getSlice(slice_region);
2706 }
2707
2708 /* TODO */
2709 /* global reduction */
2710 Data
2711 Data::getSlice(const DataTypes::RegionType& region) const
2712 {
2713 return Data(*this,region);
2714 }
2715
2716 /* TODO */
2717 /* global reduction */
2718 void
2719 Data::setItemO(const boost::python::object& key,
2720 const boost::python::object& value)
2721 {
2722 Data tempData(value,getFunctionSpace());
2723 setItemD(key,tempData);
2724 }
2725
2726 void
2727 Data::setItemD(const boost::python::object& key,
2728 const Data& value)
2729 {
2730 // const DataArrayView& view=getPointDataView();
2731
2732 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2733 if (slice_region.size()!=getDataPointRank()) {
2734 throw DataException("Error - slice size does not match Data rank.");
2735 }
2736 if (getFunctionSpace()!=value.getFunctionSpace()) {
2737 setSlice(Data(value,getFunctionSpace()),slice_region);
2738 } else {
2739 setSlice(value,slice_region);
2740 }
2741 }
2742
2743 void
2744 Data::setSlice(const Data& value,
2745 const DataTypes::RegionType& region)
2746 {
2747 if (isProtected()) {
2748 throw DataException("Error - attempt to update protected Data object.");
2749 }
2750 FORCERESOLVE;
2751 /* if (isLazy())
2752 {
2753 throw DataException("Error - setSlice not permitted on lazy data.");
2754 }*/
2755 Data tempValue(value);
2756 typeMatchLeft(tempValue);
2757 typeMatchRight(tempValue);
2758 getReady()->setSlice(tempValue.m_data.get(),region);
2759 }
2760
2761 void
2762 Data::typeMatchLeft(Data& right) const
2763 {
2764 if (right.isLazy() && !isLazy())
2765 {
2766 right.resolve();
2767 }
2768 if (isExpanded()){
2769 right.expand();
2770 } else if (isTagged()) {
2771 if (right.isConstant()) {
2772 right.tag();
2773 }
2774 }
2775 }
2776
2777 void
2778 Data::typeMatchRight(const Data& right)
2779 {
2780 if (isLazy() && !right.isLazy())
2781 {
2782 resolve();
2783 }
2784 if (isTagged()) {
2785 if (right.isExpanded()) {
2786 expand();
2787 }
2788 } else if (isConstant()) {
2789 if (right.isExpanded()) {
2790 expand();
2791 } else if (right.isTagged()) {
2792 tag();
2793 }
2794 }
2795 }
2796
2797 void
2798 Data::setTaggedValueByName(std::string name,
2799 const boost::python::object& value)
2800 {
2801 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2802 FORCERESOLVE;
2803 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2804 setTaggedValue(tagKey,value);
2805 }
2806 }
2807 void
2808 Data::setTaggedValue(int tagKey,
2809 const boost::python::object& value)
2810 {
2811 if (isProtected()) {
2812 throw DataException("Error - attempt to update protected Data object.");
2813 }
2814 //
2815 // Ensure underlying data object is of type DataTagged
2816 FORCERESOLVE;
2817 if (isConstant()) tag();
2818 numeric::array asNumArray(value);
2819
2820 // extract the shape of the numarray
2821 DataTypes::ShapeType tempShape;
2822 for (int i=0; i < asNumArray.getrank(); i++) {
2823 tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2824 }
2825
2826 DataVector temp_data2;
2827 temp_data2.copyFromNumArray(asNumArray);
2828
2829 m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2830 }
2831
2832
2833 void
2834 Data::setTaggedValueFromCPP(int tagKey,
2835 const DataTypes::ShapeType& pointshape,
2836 const DataTypes::ValueType& value,
2837 int dataOffset)
2838 {
2839 if (isProtected()) {
2840 throw DataException("Error - attempt to update protected Data object.");
2841 }
2842 //
2843 // Ensure underlying data object is of type DataTagged
2844 FORCERESOLVE;
2845 if (isConstant()) tag();
2846 //
2847 // Call DataAbstract::setTaggedValue
2848 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2849 }
2850
2851 int
2852 Data::getTagNumber(int dpno)
2853 {
2854 if (isEmpty())
2855 {
2856 throw DataException("Error - operation not permitted on instances of DataEmpty.");
2857 }
2858 return getFunctionSpace().getTagFromDataPointNo(dpno);
2859 }
2860
2861
2862 ostream& escript::operator<<(ostream& o, const Data& data)
2863 {
2864 o << data.toString();
2865 return o;
2866 }
2867
2868 Data
2869 escript::C_GeneralTensorProduct(Data& arg_0,
2870 Data& arg_1,
2871 int axis_offset,
2872 int transpose)
2873 {
2874 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2875 // SM is the product of the last axis_offset entries in arg_0.getShape().
2876
2877 // deal with any lazy data
2878 // if (arg_0.isLazy()) {arg_0.resolve();}
2879 // if (arg_1.isLazy()) {arg_1.resolve();}
2880 if (arg_0.isLazy() || arg_1.isLazy())
2881 {
2882 DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2883 return Data(c);
2884 }
2885
2886 // Interpolate if necessary and find an appropriate function space
2887 Data arg_0_Z, arg_1_Z;
2888 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2889 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2890 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2891 arg_1_Z = Data(arg_1);
2892 }
2893 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2894 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2895 arg_0_Z =Data(arg_0);
2896 }
2897 else {
2898 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2899 }
2900 } else {
2901 arg_0_Z = Data(arg_0);
2902 arg_1_Z = Data(arg_1);
2903 }
2904 // Get rank and shape of inputs
2905 int rank0 = arg_0_Z.getDataPointRank();
2906 int rank1 = arg_1_Z.getDataPointRank();
2907 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2908 const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2909
2910 // Prepare for the loops of the product and verify compatibility of shapes
2911 int start0=0, start1=0;
2912 if (transpose == 0) {}
2913 else if (transpose == 1) { start0 = axis_offset; }
2914 else if (transpose == 2) { start1 = rank1-axis_offset; }
2915 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2916
2917
2918 // Adjust the shapes for transpose
2919 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2920 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2921 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2922 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2923
2924 #if 0
2925 // For debugging: show shape after transpose
2926 char tmp[100];
2927 std::string shapeStr;
2928 shapeStr = "(";
2929 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2930 shapeStr += ")";
2931 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2932 shapeStr = "(";
2933 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2934 shapeStr += ")";
2935 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2936 #endif
2937
2938 // Prepare for the loops of the product
2939 int SL=1, SM=1, SR=1;
2940 for (int i=0; i<rank0-axis_offset; i++) {
2941 SL *= tmpShape0[i];
2942 }
2943 for (int i=rank0-axis_offset; i<rank0; i++) {
2944 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2945 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2946 }
2947 SM *= tmpShape0[i];
2948 }
2949 for (int i=axis_offset; i<rank1; i++) {
2950 SR *= tmpShape1[i];
2951 }
2952
2953 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2954 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2955 { // block to limit the scope of out_index
2956 int out_index=0;
2957 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2958 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2959 }
2960
2961 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2962 {
2963 ostringstream os;
2964 os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2965 throw DataException(os.str());
2966 }
2967
2968 // Declare output Data object
2969 Data res;
2970
2971 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2972 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2973 double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2974 double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2975 double *ptr_2 = &(res.getDataAtOffset(0));
2976 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2977 }
2978 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2979
2980 // Prepare the DataConstant input
2981 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2982 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2983
2984 // Borrow DataTagged input from Data object
2985 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2986 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2987
2988 // Prepare a DataTagged output 2
2989 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2990 res.tag();
2991 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2992 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2993
2994 // Prepare offset into DataConstant
2995 int offset_0 = tmp_0->getPointOffset(0,0);
2996 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2997 // Get the views
2998 // DataArrayView view_1 = tmp_1->getDefaultValue();
2999 // DataArrayView view_2 = tmp_2->getDefaultValue();
3000 // // Get the pointers to the actual data
3001 // double *ptr_1 = &((view_1.getData())[0]);
3002 // double *ptr_2 = &((view_2.getData())[0]);
3003
3004 double *ptr_1 = &(tmp_1->getDefaultValue(0));
3005 double *ptr_2 = &(tmp_2->getDefaultValue(0));
3006
3007
3008 // Compute an MVP for the default
3009 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3010 // Compute an MVP for each tag
3011 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
3012 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3013 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
3014 tmp_2->addTag(i->first);
3015 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
3016 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
3017 // double *ptr_1 = &view_1.getData(0);
3018 // double *ptr_2 = &view_2.getData(0);
3019
3020 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
3021 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
3022
3023 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3024 }
3025
3026 }
3027 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
3028
3029 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3030 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
3031 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3032 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3033 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
3034 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3035 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3036 int sampleNo_1,dataPointNo_1;
3037 int numSamples_1 = arg_1_Z.getNumSamples();
3038 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
3039 int offset_0 = tmp_0->getPointOffset(0,0);
3040 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
3041 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
3042 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
3043 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
3044 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
3045 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
3046 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
3047 double *ptr_2 = &(res.getDataAtOffset(offset_2));
3048 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3049 }
3050 }
3051
3052 }
3053 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
3054
3055 // Borrow DataTagged input from Data object
3056 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3057 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
3058
3059 // Prepare the DataConstant input
3060 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
3061 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
3062
3063 // Prepare a DataTagged output 2
3064 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
3065 res.tag();
3066 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3067 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3068
3069 // Prepare offset into DataConstant
3070 int offset_1 = tmp_1->getPointOffset(0,0);
3071 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
3072 // Get the views
3073 // DataArrayView view_0 = tmp_0->getDefaultValue();
3074 // DataArrayView view_2 = tmp_2->getDefaultValue();
3075 // // Get the pointers to the actual data
3076 // double *ptr_0 = &((view_0.getData())[0]);
3077 // double *ptr_2 = &((view_2.getData())[0]);
3078
3079 double *ptr_0 = &(tmp_0->getDefaultValue(0));
3080 double *ptr_2 = &(tmp_2->getDefaultValue(0));
3081
3082 // Compute an MVP for the default
3083 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3084 // Compute an MVP for each tag
3085 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3086 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3087 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3088 // tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
3089 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
3090 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
3091 // double *ptr_0 = &view_0.getData(0);
3092 // double *ptr_2 = &view_2.getData(0);
3093
3094 tmp_2->addTag(i->first);
3095 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
3096 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
3097 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3098 }
3099
3100 }
3101 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
3102
3103 // Borrow DataTagged input from Data object
3104 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3105 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3106
3107 // Borrow DataTagged input from Data object
3108 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3109 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3110
3111 // Prepare a DataTagged output 2
3112 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
3113 res.tag(); // DataTagged output
3114 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3115 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3116
3117 // // Get the views
3118 // DataArrayView view_0 = tmp_0->getDefaultValue();
3119 // DataArrayView view_1 = tmp_1->getDefaultValue();
3120 // DataArrayView view_2 = tmp_2->getDefaultValue();
3121 // // Get the pointers to the actual data
3122 // double *ptr_0 = &((view_0.getData())[0]);
3123 // double *ptr_1 = &((view_1.getData())[0]);
3124 // double *ptr_2 = &((view_2.getData())[0]);
3125
3126 double *ptr_0 = &(tmp_0->getDefaultValue(0));
3127 double *ptr_1 = &(tmp_1->getDefaultValue(0));
3128 double *ptr_2 = &(tmp_2->getDefaultValue(0));
3129
3130
3131 // Compute an MVP for the default
3132 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3133 // Merge the tags
3134 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3135 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3136 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
3137 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3138 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
3139 }
3140 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
3141 tmp_2->addTag(i->first);
3142 }
3143 // Compute an MVP for each tag
3144 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
3145 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
3146 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
3147 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
3148 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
3149 // double *ptr_0 = &view_0.getData(0);
3150 // double *ptr_1 = &view_1.getData(0);
3151 // double *ptr_2 = &view_2.getData(0);
3152
3153 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
3154 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
3155 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
3156
3157 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3158 }
3159
3160 }
3161 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
3162
3163 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3164 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3165 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3166 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3167 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3168 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3169 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3170 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3171 int sampleNo_0,dataPointNo_0;
3172 int numSamples_0 = arg_0_Z.getNumSamples();
3173 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3174 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3175 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3176 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
3177 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
3178 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3179 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3180 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3181 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
3182 double *ptr_2 = &(res.getDataAtOffset(offset_2));
3183 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3184 }
3185 }
3186
3187 }
3188 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
3189
3190 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3191 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3192 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
3193 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3194 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3195 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
3196 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3197 int sampleNo_0,dataPointNo_0;
3198 int numSamples_0 = arg_0_Z.getNumSamples();
3199 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3200 int offset_1 = tmp_1->getPointOffset(0,0);
3201 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3202 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3203 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3204 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3205 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3206 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
3207 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
3208 double *ptr_2 = &(res.getDataAtOffset(offset_2));
3209 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3210 }
3211 }
3212
3213
3214 }
3215 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
3216
3217 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3218 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3219 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3220 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3221 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3222 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3223 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3224 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3225 int sampleNo_0,dataPointNo_0;
3226 int numSamples_0 = arg_0_Z.getNumSamples();
3227 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3228 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3229 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3230 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
3231 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
3232 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3233 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3234 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3235 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
3236 double *ptr_2 = &(res.getDataAtOffset(offset_2));
3237 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3238 }
3239 }
3240
3241 }
3242 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
3243
3244 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3245 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3246 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3247 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3248 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3249 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3250 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3251 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3252 int sampleNo_0,dataPointNo_0;
3253 int numSamples_0 = arg_0_Z.getNumSamples();
3254 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3255 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3256 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3257 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3258 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3259 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3260 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3261 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
3262 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
3263 double *ptr_2 = &(res.getDataAtOffset(offset_2));
3264 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3265 }
3266 }
3267
3268 }
3269 else {
3270 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
3271 }
3272
3273 return res;
3274 }
3275
3276 DataAbstract*
3277 Data::borrowData() const
3278 {
3279 return m_data.get();
3280 }
3281
3282 // Not all that happy about returning a non-const from a const
3283 DataAbstract_ptr
3284 Data::borrowDataPtr() const
3285 {
3286 return m_data;
3287 }
3288
3289 // Not all that happy about returning a non-const from a const
3290 DataReady_ptr
3291 Data::borrowReadyPtr() const
3292 {
3293 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(m_data);
3294 EsysAssert((dr!=0), "Error - casting to DataReady.");
3295 return dr;
3296 }
3297
3298 std::string
3299 Data::toString() const
3300 {
3301 if (!m_data->isEmpty() &&
3302 !m_data->isLazy() &&
3303 getLength()>escriptParams.getInt("TOO_MANY_LINES"))
3304 {
3305 stringstream temp;
3306 temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
3307 return temp.str();
3308 }
3309 return m_data->toString();
3310 }
3311
3312
3313
3314 DataTypes::ValueType::const_reference
3315 Data::getDataAtOffset(DataTypes::ValueType::size_type i) const
3316 {
3317 if (isLazy())
3318 {
3319 throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");
3320 }
3321 return getReady()->getDataAtOffset(i);
3322 }
3323
3324
3325 DataTypes::ValueType::reference
3326 Data::getDataAtOffset(DataTypes::ValueType::size_type i)
3327 {
3328 // if (isLazy())
3329 // {
3330 // throw DataException("getDataAtOffset not permitted on lazy data.");
3331 // }
3332 FORCERESOLVE;
3333 return getReady()->getDataAtOffset(i);
3334 }
3335
3336 DataTypes::ValueType::const_reference
3337 Data::getDataPoint(int sampleNo, int dataPointNo) const
3338 {
3339