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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1949 - (show annotations)
Thu Oct 30 00:10:35 2008 UTC (11 years, 5 months ago) by jfenwick
File size: 85483 byte(s)
Branch commit.
Fixed some compile errors.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26