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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 876 - (show annotations)
Thu Oct 19 03:50:23 2006 UTC (12 years, 6 months ago) by ksteube
File size: 83039 byte(s)
Added erf (error function) implementation

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26