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