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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 854 - (show annotations)
Thu Sep 21 05:29:42 2006 UTC (13 years, 2 months ago) by gross
File size: 82904 byte(s)
Some modifications to the binary operations +,-,*/, pow. 
The code is a bit simpler now and more efficient has there is
no reseising required now. the resizing method has been removed as
it is very, very inefficient. Even serial code should be faster now.
It is now forbidden to do an inplace update of scalar data object with an object 
of rank >0 as this is very slow (and does not make much sense). 


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::sin() const
1036 {
1037 #if defined DOPROF
1038 profData->unary++;
1039 #endif
1040 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
1041 }
1042
1043 Data
1044 Data::cos() const
1045 {
1046 #if defined DOPROF
1047 profData->unary++;
1048 #endif
1049 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
1050 }
1051
1052 Data
1053 Data::tan() const
1054 {
1055 #if defined DOPROF
1056 profData->unary++;
1057 #endif
1058 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
1059 }
1060
1061 Data
1062 Data::asin() const
1063 {
1064 #if defined DOPROF
1065 profData->unary++;
1066 #endif
1067 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
1068 }
1069
1070 Data
1071 Data::acos() const
1072 {
1073 #if defined DOPROF
1074 profData->unary++;
1075 #endif
1076 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1077 }
1078
1079 Data
1080 Data::atan() const
1081 {
1082 #if defined DOPROF
1083 profData->unary++;
1084 #endif
1085 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
1086 }
1087
1088 Data
1089 Data::sinh() const
1090 {
1091 #if defined DOPROF
1092 profData->unary++;
1093 #endif
1094 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
1095 }
1096
1097 Data
1098 Data::cosh() const
1099 {
1100 #if defined DOPROF
1101 profData->unary++;
1102 #endif
1103 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
1104 }
1105
1106 Data
1107 Data::tanh() const
1108 {
1109 #if defined DOPROF
1110 profData->unary++;
1111 #endif
1112 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1113 }
1114
1115 Data
1116 Data::asinh() const
1117 {
1118 #if defined DOPROF
1119 profData->unary++;
1120 #endif
1121 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1122 }
1123
1124 Data
1125 Data::acosh() const
1126 {
1127 #if defined DOPROF
1128 profData->unary++;
1129 #endif
1130 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1131 }
1132
1133 Data
1134 Data::atanh() const
1135 {
1136 #if defined DOPROF
1137 profData->unary++;
1138 #endif
1139 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1140 }
1141
1142 Data
1143 Data::log10() const
1144 {
1145 #if defined DOPROF
1146 profData->unary++;
1147 #endif
1148 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1149 }
1150
1151 Data
1152 Data::log() const
1153 {
1154 #if defined DOPROF
1155 profData->unary++;
1156 #endif
1157 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1158 }
1159
1160 Data
1161 Data::sign() const
1162 {
1163 #if defined DOPROF
1164 profData->unary++;
1165 #endif
1166 return escript::unaryOp(*this,escript::fsign);
1167 }
1168
1169 Data
1170 Data::abs() const
1171 {
1172 #if defined DOPROF
1173 profData->unary++;
1174 #endif
1175 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1176 }
1177
1178 Data
1179 Data::neg() const
1180 {
1181 #if defined DOPROF
1182 profData->unary++;
1183 #endif
1184 return escript::unaryOp(*this,negate<double>());
1185 }
1186
1187 Data
1188 Data::pos() const
1189 {
1190 #if defined DOPROF
1191 profData->unary++;
1192 #endif
1193 Data result;
1194 // perform a deep copy
1195 result.copy(*this);
1196 return result;
1197 }
1198
1199 Data
1200 Data::exp() const
1201 {
1202 #if defined DOPROF
1203 profData->unary++;
1204 #endif
1205 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1206 }
1207
1208 Data
1209 Data::sqrt() const
1210 {
1211 #if defined DOPROF
1212 profData->unary++;
1213 #endif
1214 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1215 }
1216
1217 double
1218 Data::Lsup() const
1219 {
1220 double localValue, globalValue;
1221 #if defined DOPROF
1222 profData->reduction1++;
1223 #endif
1224 //
1225 // set the initial absolute maximum value to zero
1226
1227 AbsMax abs_max_func;
1228 localValue = algorithm(abs_max_func,0);
1229 #ifdef PASO_MPI
1230 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1231 return globalValue;
1232 #else
1233 return localValue;
1234 #endif
1235 }
1236
1237 double
1238 Data::Linf() const
1239 {
1240 double localValue, globalValue;
1241 #if defined DOPROF
1242 profData->reduction1++;
1243 #endif
1244 //
1245 // set the initial absolute minimum value to max double
1246 AbsMin abs_min_func;
1247 localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1248
1249 #ifdef PASO_MPI
1250 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1251 return globalValue;
1252 #else
1253 return localValue;
1254 #endif
1255 }
1256
1257 double
1258 Data::sup() const
1259 {
1260 double localValue, globalValue;
1261 #if defined DOPROF
1262 profData->reduction1++;
1263 #endif
1264 //
1265 // set the initial maximum value to min possible double
1266 FMax fmax_func;
1267 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1268 #ifdef PASO_MPI
1269 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1270 return globalValue;
1271 #else
1272 return localValue;
1273 #endif
1274 }
1275
1276 double
1277 Data::inf() const
1278 {
1279 double localValue, globalValue;
1280 #if defined DOPROF
1281 profData->reduction1++;
1282 #endif
1283 //
1284 // set the initial minimum value to max possible double
1285 FMin fmin_func;
1286 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1287 #ifdef PASO_MPI
1288 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1289 return globalValue;
1290 #else
1291 return localValue;
1292 #endif
1293 }
1294
1295 /* TODO */
1296 /* global reduction */
1297 Data
1298 Data::maxval() const
1299 {
1300 #if defined DOPROF
1301 profData->reduction2++;
1302 #endif
1303 //
1304 // set the initial maximum value to min possible double
1305 FMax fmax_func;
1306 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1307 }
1308
1309 Data
1310 Data::minval() const
1311 {
1312 #if defined DOPROF
1313 profData->reduction2++;
1314 #endif
1315 //
1316 // set the initial minimum value to max possible double
1317 FMin fmin_func;
1318 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1319 }
1320
1321 Data
1322 Data::swapaxes(const int axis0, const int axis1) const
1323 {
1324 int axis0_tmp,axis1_tmp;
1325 #if defined DOPROF
1326 profData->unary++;
1327 #endif
1328 DataArrayView::ShapeType s=getDataPointShape();
1329 DataArrayView::ShapeType ev_shape;
1330 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1331 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1332 int rank=getDataPointRank();
1333 if (rank<2) {
1334 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1335 }
1336 if (axis0<0 || axis0>rank-1) {
1337 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1338 }
1339 if (axis1<0 || axis1>rank-1) {
1340 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1341 }
1342 if (axis0 == axis1) {
1343 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1344 }
1345 if (axis0 > axis1) {
1346 axis0_tmp=axis1;
1347 axis1_tmp=axis0;
1348 } else {
1349 axis0_tmp=axis0;
1350 axis1_tmp=axis1;
1351 }
1352 for (int i=0; i<rank; i++) {
1353 if (i == axis0_tmp) {
1354 ev_shape.push_back(s[axis1_tmp]);
1355 } else if (i == axis1_tmp) {
1356 ev_shape.push_back(s[axis0_tmp]);
1357 } else {
1358 ev_shape.push_back(s[i]);
1359 }
1360 }
1361 Data ev(0.,ev_shape,getFunctionSpace());
1362 ev.typeMatchRight(*this);
1363 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1364 return ev;
1365
1366 }
1367
1368 Data
1369 Data::symmetric() const
1370 {
1371 #if defined DOPROF
1372 profData->unary++;
1373 #endif
1374 // check input
1375 DataArrayView::ShapeType s=getDataPointShape();
1376 if (getDataPointRank()==2) {
1377 if(s[0] != s[1])
1378 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1379 }
1380 else if (getDataPointRank()==4) {
1381 if(!(s[0] == s[2] && s[1] == s[3]))
1382 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1383 }
1384 else {
1385 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1386 }
1387 Data ev(0.,getDataPointShape(),getFunctionSpace());
1388 ev.typeMatchRight(*this);
1389 m_data->symmetric(ev.m_data.get());
1390 return ev;
1391 }
1392
1393 Data
1394 Data::nonsymmetric() const
1395 {
1396 #if defined DOPROF
1397 profData->unary++;
1398 #endif
1399 // check input
1400 DataArrayView::ShapeType s=getDataPointShape();
1401 if (getDataPointRank()==2) {
1402 if(s[0] != s[1])
1403 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1404 DataArrayView::ShapeType ev_shape;
1405 ev_shape.push_back(s[0]);
1406 ev_shape.push_back(s[1]);
1407 Data ev(0.,ev_shape,getFunctionSpace());
1408 ev.typeMatchRight(*this);
1409 m_data->nonsymmetric(ev.m_data.get());
1410 return ev;
1411 }
1412 else if (getDataPointRank()==4) {
1413 if(!(s[0] == s[2] && s[1] == s[3]))
1414 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1415 DataArrayView::ShapeType ev_shape;
1416 ev_shape.push_back(s[0]);
1417 ev_shape.push_back(s[1]);
1418 ev_shape.push_back(s[2]);
1419 ev_shape.push_back(s[3]);
1420 Data ev(0.,ev_shape,getFunctionSpace());
1421 ev.typeMatchRight(*this);
1422 m_data->nonsymmetric(ev.m_data.get());
1423 return ev;
1424 }
1425 else {
1426 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1427 }
1428 }
1429
1430 Data
1431 Data::trace(int axis_offset) const
1432 {
1433 #if defined DOPROF
1434 profData->unary++;
1435 #endif
1436 DataArrayView::ShapeType s=getDataPointShape();
1437 if (getDataPointRank()==2) {
1438 DataArrayView::ShapeType ev_shape;
1439 Data ev(0.,ev_shape,getFunctionSpace());
1440 ev.typeMatchRight(*this);
1441 m_data->trace(ev.m_data.get(), axis_offset);
1442 return ev;
1443 }
1444 if (getDataPointRank()==3) {
1445 DataArrayView::ShapeType ev_shape;
1446 if (axis_offset==0) {
1447 int s2=s[2];
1448 ev_shape.push_back(s2);
1449 }
1450 else if (axis_offset==1) {
1451 int s0=s[0];
1452 ev_shape.push_back(s0);
1453 }
1454 Data ev(0.,ev_shape,getFunctionSpace());
1455 ev.typeMatchRight(*this);
1456 m_data->trace(ev.m_data.get(), axis_offset);
1457 return ev;
1458 }
1459 if (getDataPointRank()==4) {
1460 DataArrayView::ShapeType ev_shape;
1461 if (axis_offset==0) {
1462 ev_shape.push_back(s[2]);
1463 ev_shape.push_back(s[3]);
1464 }
1465 else if (axis_offset==1) {
1466 ev_shape.push_back(s[0]);
1467 ev_shape.push_back(s[3]);
1468 }
1469 else if (axis_offset==2) {
1470 ev_shape.push_back(s[0]);
1471 ev_shape.push_back(s[1]);
1472 }
1473 Data ev(0.,ev_shape,getFunctionSpace());
1474 ev.typeMatchRight(*this);
1475 m_data->trace(ev.m_data.get(), axis_offset);
1476 return ev;
1477 }
1478 else {
1479 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1480 }
1481 }
1482
1483 Data
1484 Data::transpose(int axis_offset) const
1485 {
1486 #if defined DOPROF
1487 profData->unary++;
1488 #endif
1489 DataArrayView::ShapeType s=getDataPointShape();
1490 DataArrayView::ShapeType ev_shape;
1491 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1492 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1493 int rank=getDataPointRank();
1494 if (axis_offset<0 || axis_offset>rank) {
1495 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1496 }
1497 for (int i=0; i<rank; i++) {
1498 int index = (axis_offset+i)%rank;
1499 ev_shape.push_back(s[index]); // Append to new shape
1500 }
1501 Data ev(0.,ev_shape,getFunctionSpace());
1502 ev.typeMatchRight(*this);
1503 m_data->transpose(ev.m_data.get(), axis_offset);
1504 return ev;
1505 }
1506
1507 Data
1508 Data::eigenvalues() const
1509 {
1510 #if defined DOPROF
1511 profData->unary++;
1512 #endif
1513 // check input
1514 DataArrayView::ShapeType s=getDataPointShape();
1515 if (getDataPointRank()!=2)
1516 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1517 if(s[0] != s[1])
1518 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1519 // create return
1520 DataArrayView::ShapeType ev_shape(1,s[0]);
1521 Data ev(0.,ev_shape,getFunctionSpace());
1522 ev.typeMatchRight(*this);
1523 m_data->eigenvalues(ev.m_data.get());
1524 return ev;
1525 }
1526
1527 const boost::python::tuple
1528 Data::eigenvalues_and_eigenvectors(const double tol) const
1529 {
1530 #if defined DOPROF
1531 profData->unary++;
1532 #endif
1533 DataArrayView::ShapeType s=getDataPointShape();
1534 if (getDataPointRank()!=2)
1535 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1536 if(s[0] != s[1])
1537 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1538 // create return
1539 DataArrayView::ShapeType ev_shape(1,s[0]);
1540 Data ev(0.,ev_shape,getFunctionSpace());
1541 ev.typeMatchRight(*this);
1542 DataArrayView::ShapeType V_shape(2,s[0]);
1543 Data V(0.,V_shape,getFunctionSpace());
1544 V.typeMatchRight(*this);
1545 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1546 return make_tuple(boost::python::object(ev),boost::python::object(V));
1547 }
1548
1549 const boost::python::tuple
1550 Data::mindp() const
1551 {
1552 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1553 // abort (for unknown reasons) if there are openmp directives with it in the
1554 // surrounding function
1555
1556 int SampleNo;
1557 int DataPointNo;
1558 int ProcNo;
1559 calc_mindp(ProcNo,SampleNo,DataPointNo);
1560 return make_tuple(ProcNo,SampleNo,DataPointNo);
1561 }
1562
1563 void
1564 Data::calc_mindp( int& ProcNo,
1565 int& SampleNo,
1566 int& DataPointNo) const
1567 {
1568 int i,j;
1569 int lowi=0,lowj=0;
1570 double min=numeric_limits<double>::max();
1571
1572 Data temp=minval();
1573
1574 int numSamples=temp.getNumSamples();
1575 int numDPPSample=temp.getNumDataPointsPerSample();
1576
1577 double next,local_min;
1578 int local_lowi,local_lowj;
1579
1580 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1581 {
1582 local_min=min;
1583 #pragma omp for private(i,j) schedule(static)
1584 for (i=0; i<numSamples; i++) {
1585 for (j=0; j<numDPPSample; j++) {
1586 next=temp.getDataPoint(i,j)();
1587 if (next<local_min) {
1588 local_min=next;
1589 local_lowi=i;
1590 local_lowj=j;
1591 }
1592 }
1593 }
1594 #pragma omp critical
1595 if (local_min<min) {
1596 min=local_min;
1597 lowi=local_lowi;
1598 lowj=local_lowj;
1599 }
1600 }
1601
1602 #ifdef PASO_MPI
1603 // determine the processor on which the minimum occurs
1604 next = temp.getDataPoint(lowi,lowj)();
1605 int lowProc = 0;
1606 double *globalMins = new double[get_MPISize()+1];
1607 int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1608
1609 if( get_MPIRank()==0 ){
1610 next = globalMins[lowProc];
1611 for( i=1; i<get_MPISize(); i++ )
1612 if( next>globalMins[i] ){
1613 lowProc = i;
1614 next = globalMins[i];
1615 }
1616 }
1617 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1618
1619 delete [] globalMins;
1620 ProcNo = lowProc;
1621 #else
1622 ProcNo = 0;
1623 #endif
1624 SampleNo = lowi;
1625 DataPointNo = lowj;
1626 }
1627
1628 void
1629 Data::saveDX(std::string fileName) const
1630 {
1631 boost::python::dict args;
1632 args["data"]=boost::python::object(this);
1633 getDomain().saveDX(fileName,args);
1634 return;
1635 }
1636
1637 void
1638 Data::saveVTK(std::string fileName) const
1639 {
1640 boost::python::dict args;
1641 args["data"]=boost::python::object(this);
1642 getDomain().saveVTK(fileName,args);
1643 return;
1644 }
1645
1646 Data&
1647 Data::operator+=(const Data& right)
1648 {
1649 if (isProtected()) {
1650 throw DataException("Error - attempt to update protected Data object.");
1651 }
1652 binaryOp(right,plus<double>());
1653 return (*this);
1654 }
1655
1656 Data&
1657 Data::operator+=(const boost::python::object& right)
1658 {
1659 Data tmp(right,getFunctionSpace(),false);
1660 binaryOp(tmp,plus<double>());
1661 return (*this);
1662 }
1663
1664 Data&
1665 Data::operator-=(const Data& right)
1666 {
1667 if (isProtected()) {
1668 throw DataException("Error - attempt to update protected Data object.");
1669 }
1670 binaryOp(right,minus<double>());
1671 return (*this);
1672 }
1673
1674 Data&
1675 Data::operator-=(const boost::python::object& right)
1676 {
1677 Data tmp(right,getFunctionSpace(),false);
1678 binaryOp(tmp,minus<double>());
1679 return (*this);
1680 }
1681
1682 Data&
1683 Data::operator*=(const Data& right)
1684 {
1685 if (isProtected()) {
1686 throw DataException("Error - attempt to update protected Data object.");
1687 }
1688 binaryOp(right,multiplies<double>());
1689 return (*this);
1690 }
1691
1692 Data&
1693 Data::operator*=(const boost::python::object& right)
1694 {
1695 Data tmp(right,getFunctionSpace(),false);
1696 binaryOp(tmp,multiplies<double>());
1697 return (*this);
1698 }
1699
1700 Data&
1701 Data::operator/=(const Data& right)
1702 {
1703 if (isProtected()) {
1704 throw DataException("Error - attempt to update protected Data object.");
1705 }
1706 binaryOp(right,divides<double>());
1707 return (*this);
1708 }
1709
1710 Data&
1711 Data::operator/=(const boost::python::object& right)
1712 {
1713 Data tmp(right,getFunctionSpace(),false);
1714 binaryOp(tmp,divides<double>());
1715 return (*this);
1716 }
1717
1718 Data
1719 Data::rpowO(const boost::python::object& left) const
1720 {
1721 Data left_d(left,*this);
1722 return left_d.powD(*this);
1723 }
1724
1725 Data
1726 Data::powO(const boost::python::object& right) const
1727 {
1728 Data tmp(right,getFunctionSpace(),false);
1729 return powD(tmp);
1730 }
1731
1732 Data
1733 Data::powD(const Data& right) const
1734 {
1735 Data result;
1736 if (getDataPointRank()<right.getDataPointRank()) {
1737 result.copy(right);
1738 result.binaryOp(*this,escript::rpow);
1739 } else {
1740 result.copy(*this);
1741 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1742 }
1743 return result;
1744 }
1745
1746
1747 //
1748 // NOTE: It is essential to specify the namespace this operator belongs to
1749 Data
1750 escript::operator+(const Data& left, const Data& right)
1751 {
1752 Data result;
1753 //
1754 // perform a deep copy
1755 if (left.getDataPointRank()<right.getDataPointRank()) {
1756 result.copy(right);
1757 result+=left;
1758 } else {
1759 result.copy(left);
1760 result+=right;
1761 }
1762 return result;
1763 }
1764
1765 //
1766 // NOTE: It is essential to specify the namespace this operator belongs to
1767 Data
1768 escript::operator-(const Data& left, const Data& right)
1769 {
1770 Data result;
1771 //
1772 // perform a deep copy
1773 if (left.getDataPointRank()<right.getDataPointRank()) {
1774 result=right.neg();
1775 result+=left;
1776 } else {
1777 result.copy(left);
1778 result-=right;
1779 }
1780 return result;
1781 }
1782
1783 //
1784 // NOTE: It is essential to specify the namespace this operator belongs to
1785 Data
1786 escript::operator*(const Data& left, const Data& right)
1787 {
1788 Data result;
1789 //
1790 // perform a deep copy
1791 if (left.getDataPointRank()<right.getDataPointRank()) {
1792 result.copy(right);
1793 result*=left;
1794 } else {
1795 result.copy(left);
1796 result*=right;
1797 }
1798 return result;
1799 }
1800
1801 //
1802 // NOTE: It is essential to specify the namespace this operator belongs to
1803 Data
1804 escript::operator/(const Data& left, const Data& right)
1805 {
1806 Data result;
1807 //
1808 // perform a deep copy
1809 if (left.getDataPointRank()<right.getDataPointRank()) {
1810 result=right.oneOver();
1811 result*=left;
1812 } else {
1813 result.copy(left);
1814 result/=right;
1815 }
1816 return result;
1817 }
1818
1819 //
1820 // NOTE: It is essential to specify the namespace this operator belongs to
1821 Data
1822 escript::operator+(const Data& left, const boost::python::object& right)
1823 {
1824 return left+Data(right,left.getFunctionSpace(),false);
1825 }
1826
1827 //
1828 // NOTE: It is essential to specify the namespace this operator belongs to
1829 Data
1830 escript::operator-(const Data& left, const boost::python::object& right)
1831 {
1832 return left-Data(right,left.getFunctionSpace(),false);
1833 }
1834
1835 //
1836 // NOTE: It is essential to specify the namespace this operator belongs to
1837 Data
1838 escript::operator*(const Data& left, const boost::python::object& right)
1839 {
1840 return left*Data(right,left.getFunctionSpace(),false);
1841 }
1842
1843 //
1844 // NOTE: It is essential to specify the namespace this operator belongs to
1845 Data
1846 escript::operator/(const Data& left, const boost::python::object& right)
1847 {
1848 return left/Data(right,left.getFunctionSpace(),false);
1849 }
1850
1851 //
1852 // NOTE: It is essential to specify the namespace this operator belongs to
1853 Data
1854 escript::operator+(const boost::python::object& left, const Data& right)
1855 {
1856 return Data(left,right.getFunctionSpace(),false)+right;
1857 }
1858
1859 //
1860 // NOTE: It is essential to specify the namespace this operator belongs to
1861 Data
1862 escript::operator-(const boost::python::object& left, const Data& right)
1863 {
1864 return Data(left,right.getFunctionSpace(),false)-right;
1865 }
1866
1867 //
1868 // NOTE: It is essential to specify the namespace this operator belongs to
1869 Data
1870 escript::operator*(const boost::python::object& left, const Data& right)
1871 {
1872 return Data(left,right.getFunctionSpace(),false)*right;
1873 }
1874
1875 //
1876 // NOTE: It is essential to specify the namespace this operator belongs to
1877 Data
1878 escript::operator/(const boost::python::object& left, const Data& right)
1879 {
1880 return Data(left,right.getFunctionSpace(),false)/right;
1881 }
1882
1883 //
1884 //bool escript::operator==(const Data& left, const Data& right)
1885 //{
1886 // /*
1887 // NB: this operator does very little at this point, and isn't to
1888 // be relied on. Requires further implementation.
1889 // */
1890 //
1891 // bool ret;
1892 //
1893 // if (left.isEmpty()) {
1894 // if(!right.isEmpty()) {
1895 // ret = false;
1896 // } else {
1897 // ret = true;
1898 // }
1899 // }
1900 //
1901 // if (left.isConstant()) {
1902 // if(!right.isConstant()) {
1903 // ret = false;
1904 // } else {
1905 // ret = true;
1906 // }
1907 // }
1908 //
1909 // if (left.isTagged()) {
1910 // if(!right.isTagged()) {
1911 // ret = false;
1912 // } else {
1913 // ret = true;
1914 // }
1915 // }
1916 //
1917 // if (left.isExpanded()) {
1918 // if(!right.isExpanded()) {
1919 // ret = false;
1920 // } else {
1921 // ret = true;
1922 // }
1923 // }
1924 //
1925 // return ret;
1926 //}
1927
1928 /* TODO */
1929 /* global reduction */
1930 Data
1931 Data::getItem(const boost::python::object& key) const
1932 {
1933 const DataArrayView& view=getPointDataView();
1934
1935 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1936
1937 if (slice_region.size()!=view.getRank()) {
1938 throw DataException("Error - slice size does not match Data rank.");
1939 }
1940
1941 return getSlice(slice_region);
1942 }
1943
1944 /* TODO */
1945 /* global reduction */
1946 Data
1947 Data::getSlice(const DataArrayView::RegionType& region) const
1948 {
1949 #if defined DOPROF
1950 profData->slicing++;
1951 #endif
1952 return Data(*this,region);
1953 }
1954
1955 /* TODO */
1956 /* global reduction */
1957 void
1958 Data::setItemO(const boost::python::object& key,
1959 const boost::python::object& value)
1960 {
1961 Data tempData(value,getFunctionSpace());
1962 setItemD(key,tempData);
1963 }
1964
1965 /* TODO */
1966 /* global reduction */
1967 void
1968 Data::setItemD(const boost::python::object& key,
1969 const Data& value)
1970 {
1971 const DataArrayView& view=getPointDataView();
1972
1973 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1974 if (slice_region.size()!=view.getRank()) {
1975 throw DataException("Error - slice size does not match Data rank.");
1976 }
1977 if (getFunctionSpace()!=value.getFunctionSpace()) {
1978 setSlice(Data(value,getFunctionSpace()),slice_region);
1979 } else {
1980 setSlice(value,slice_region);
1981 }
1982 }
1983
1984 /* TODO */
1985 /* global reduction */
1986 void
1987 Data::setSlice(const Data& value,
1988 const DataArrayView::RegionType& region)
1989 {
1990 if (isProtected()) {
1991 throw DataException("Error - attempt to update protected Data object.");
1992 }
1993 #if defined DOPROF
1994 profData->slicing++;
1995 #endif
1996 Data tempValue(value);
1997 typeMatchLeft(tempValue);
1998 typeMatchRight(tempValue);
1999 m_data->setSlice(tempValue.m_data.get(),region);
2000 }
2001
2002 void
2003 Data::typeMatchLeft(Data& right) const
2004 {
2005 if (isExpanded()){
2006 right.expand();
2007 } else if (isTagged()) {
2008 if (right.isConstant()) {
2009 right.tag();
2010 }
2011 }
2012 }
2013
2014 void
2015 Data::typeMatchRight(const Data& right)
2016 {
2017 if (isTagged()) {
2018 if (right.isExpanded()) {
2019 expand();
2020 }
2021 } else if (isConstant()) {
2022 if (right.isExpanded()) {
2023 expand();
2024 } else if (right.isTagged()) {
2025 tag();
2026 }
2027 }
2028 }
2029
2030 /* TODO */
2031 /* global reduction */
2032 void
2033 Data::setTaggedValue(int tagKey,
2034 const boost::python::object& value)
2035 {
2036 if (isProtected()) {
2037 throw DataException("Error - attempt to update protected Data object.");
2038 }
2039 //
2040 // Ensure underlying data object is of type DataTagged
2041 tag();
2042
2043 if (!isTagged()) {
2044 throw DataException("Error - DataTagged conversion failed!!");
2045 }
2046
2047 //
2048 // Construct DataArray from boost::python::object input value
2049 DataArray valueDataArray(value);
2050
2051 //
2052 // Call DataAbstract::setTaggedValue
2053 m_data->setTaggedValue(tagKey,valueDataArray.getView());
2054 }
2055
2056 /* TODO */
2057 /* global reduction */
2058 void
2059 Data::setTaggedValueFromCPP(int tagKey,
2060 const DataArrayView& value)
2061 {
2062 if (isProtected()) {
2063 throw DataException("Error - attempt to update protected Data object.");
2064 }
2065 //
2066 // Ensure underlying data object is of type DataTagged
2067 tag();
2068
2069 if (!isTagged()) {
2070 throw DataException("Error - DataTagged conversion failed!!");
2071 }
2072
2073 //
2074 // Call DataAbstract::setTaggedValue
2075 m_data->setTaggedValue(tagKey,value);
2076 }
2077
2078 /* TODO */
2079 /* global reduction */
2080 int
2081 Data::getTagNumber(int dpno)
2082 {
2083 return m_data->getTagNumber(dpno);
2084 }
2085
2086 /* TODO */
2087 /* global reduction */
2088 void
2089 Data::setRefValue(int ref,
2090 const boost::python::numeric::array& value)
2091 {
2092 if (isProtected()) {
2093 throw DataException("Error - attempt to update protected Data object.");
2094 }
2095 //
2096 // Construct DataArray from boost::python::object input value
2097 DataArray valueDataArray(value);
2098
2099 //
2100 // Call DataAbstract::setRefValue
2101 m_data->setRefValue(ref,valueDataArray);
2102 }
2103
2104 /* TODO */
2105 /* global reduction */
2106 void
2107 Data::getRefValue(int ref,
2108 boost::python::numeric::array& value)
2109 {
2110 //
2111 // Construct DataArray for boost::python::object return value
2112 DataArray valueDataArray(value);
2113
2114 //
2115 // Load DataArray with values from data-points specified by ref
2116 m_data->getRefValue(ref,valueDataArray);
2117
2118 //
2119 // Load values from valueDataArray into return numarray
2120
2121 // extract the shape of the numarray
2122 int rank = value.getrank();
2123 DataArrayView::ShapeType shape;
2124 for (int i=0; i < rank; i++) {
2125 shape.push_back(extract<int>(value.getshape()[i]));
2126 }
2127
2128 // and load the numarray with the data from the DataArray
2129 DataArrayView valueView = valueDataArray.getView();
2130
2131 if (rank==0) {
2132 boost::python::numeric::array temp_numArray(valueView());
2133 value = temp_numArray;
2134 }
2135 if (rank==1) {
2136 for (int i=0; i < shape[0]; i++) {
2137 value[i] = valueView(i);
2138 }
2139 }
2140 if (rank==2) {
2141 for (int i=0; i < shape[0]; i++) {
2142 for (int j=0; j < shape[1]; j++) {
2143 value[i][j] = valueView(i,j);
2144 }
2145 }
2146 }
2147 if (rank==3) {
2148 for (int i=0; i < shape[0]; i++) {
2149 for (int j=0; j < shape[1]; j++) {
2150 for (int k=0; k < shape[2]; k++) {
2151 value[i][j][k] = valueView(i,j,k);
2152 }
2153 }
2154 }
2155 }
2156 if (rank==4) {
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 for (int l=0; l < shape[3]; l++) {
2161 value[i][j][k][l] = valueView(i,j,k,l);
2162 }
2163 }
2164 }
2165 }
2166 }
2167
2168 }
2169
2170 void
2171 Data::archiveData(const std::string fileName)
2172 {
2173 cout << "Archiving Data object to: " << fileName << endl;
2174
2175 //
2176 // Determine type of this Data object
2177 int dataType = -1;
2178
2179 if (isEmpty()) {
2180 dataType = 0;
2181 cout << "\tdataType: DataEmpty" << endl;
2182 }
2183 if (isConstant()) {
2184 dataType = 1;
2185 cout << "\tdataType: DataConstant" << endl;
2186 }
2187 if (isTagged()) {
2188 dataType = 2;
2189 cout << "\tdataType: DataTagged" << endl;
2190 }
2191 if (isExpanded()) {
2192 dataType = 3;
2193 cout << "\tdataType: DataExpanded" << endl;
2194 }
2195
2196 if (dataType == -1) {
2197 throw DataException("archiveData Error: undefined dataType");
2198 }
2199
2200 //
2201 // Collect data items common to all Data types
2202 int noSamples = getNumSamples();
2203 int noDPPSample = getNumDataPointsPerSample();
2204 int functionSpaceType = getFunctionSpace().getTypeCode();
2205 int dataPointRank = getDataPointRank();
2206 int dataPointSize = getDataPointSize();
2207 int dataLength = getLength();
2208 DataArrayView::ShapeType dataPointShape = getDataPointShape();
2209 vector<int> referenceNumbers(noSamples);
2210 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2211 referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2212 }
2213 vector<int> tagNumbers(noSamples);
2214 if (isTagged()) {
2215 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2216 tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2217 }
2218 }
2219
2220 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2221 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2222 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2223
2224 //
2225 // Flatten Shape to an array of integers suitable for writing to file
2226 int flatShape[4] = {0,0,0,0};
2227 cout << "\tshape: < ";
2228 for (int dim=0; dim<dataPointRank; dim++) {
2229 flatShape[dim] = dataPointShape[dim];
2230 cout << dataPointShape[dim] << " ";
2231 }
2232 cout << ">" << endl;
2233
2234 //
2235 // Open archive file
2236 ofstream archiveFile;
2237 archiveFile.open(fileName.data(), ios::out);
2238
2239 if (!archiveFile.good()) {
2240 throw DataException("archiveData Error: problem opening archive file");
2241 }
2242
2243 //
2244 // Write common data items to archive file
2245 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2246 archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2247 archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2248 archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2249 archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2250 archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2251 archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2252 for (int dim = 0; dim < 4; dim++) {
2253 archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2254 }
2255 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2256 archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2257 }
2258 if (isTagged()) {
2259 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2260 archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2261 }
2262 }
2263
2264 if (!archiveFile.good()) {
2265 throw DataException("archiveData Error: problem writing to archive file");
2266 }
2267
2268 //
2269 // Archive underlying data values for each Data type
2270 int noValues;
2271 switch (dataType) {
2272 case 0:
2273 // DataEmpty
2274 noValues = 0;
2275 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2276 cout << "\tnoValues: " << noValues << endl;
2277 break;
2278 case 1:
2279 // DataConstant
2280 noValues = m_data->getLength();
2281 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2282 cout << "\tnoValues: " << noValues << endl;
2283 if (m_data->archiveData(archiveFile,noValues)) {
2284 throw DataException("archiveData Error: problem writing data to archive file");
2285 }
2286 break;
2287 case 2:
2288 // DataTagged
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 3:
2297 // DataExpanded
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 }
2306
2307 if (!archiveFile.good()) {
2308 throw DataException("archiveData Error: problem writing data to archive file");
2309 }
2310
2311 //
2312 // Close archive file
2313 archiveFile.close();
2314
2315 if (!archiveFile.good()) {
2316 throw DataException("archiveData Error: problem closing archive file");
2317 }
2318
2319 }
2320
2321 void
2322 Data::extractData(const std::string fileName,
2323 const FunctionSpace& fspace)
2324 {
2325 //
2326 // Can only extract Data to an object which is initially DataEmpty
2327 if (!isEmpty()) {
2328 throw DataException("extractData Error: can only extract to DataEmpty object");
2329 }
2330
2331 cout << "Extracting Data object from: " << fileName << endl;
2332
2333 int dataType;
2334 int noSamples;
2335 int noDPPSample;
2336 int functionSpaceType;
2337 int dataPointRank;
2338 int dataPointSize;
2339 int dataLength;
2340 DataArrayView::ShapeType dataPointShape;
2341 int flatShape[4];
2342
2343 //
2344 // Open the archive file
2345 ifstream archiveFile;
2346 archiveFile.open(fileName.data(), ios::in);
2347
2348 if (!archiveFile.good()) {
2349 throw DataException("extractData Error: problem opening archive file");
2350 }
2351
2352 //
2353 // Read common data items from archive file
2354 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2355 archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2356 archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2357 archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2358 archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2359 archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2360 archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2361 for (int dim = 0; dim < 4; dim++) {
2362 archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2363 if (flatShape[dim]>0) {
2364 dataPointShape.push_back(flatShape[dim]);
2365 }
2366 }
2367 vector<int> referenceNumbers(noSamples);
2368 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2369 archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2370 }
2371 vector<int> tagNumbers(noSamples);
2372 if (dataType==2) {
2373 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2374 archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2375 }
2376 }
2377
2378 if (!archiveFile.good()) {
2379 throw DataException("extractData Error: problem reading from archive file");
2380 }
2381
2382 //
2383 // Verify the values just read from the archive file
2384 switch (dataType) {
2385 case 0:
2386 cout << "\tdataType: DataEmpty" << endl;
2387 break;
2388 case 1:
2389 cout << "\tdataType: DataConstant" << endl;
2390 break;
2391 case 2:
2392 cout << "\tdataType: DataTagged" << endl;
2393 break;
2394 case 3:
2395 cout << "\tdataType: DataExpanded" << endl;
2396 break;
2397 default:
2398 throw DataException("extractData Error: undefined dataType read from archive file");
2399 break;
2400 }
2401
2402 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2403 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2404 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2405 cout << "\tshape: < ";
2406 for (int dim = 0; dim < dataPointRank; dim++) {
2407 cout << dataPointShape[dim] << " ";
2408 }
2409 cout << ">" << endl;
2410
2411 //
2412 // Verify that supplied FunctionSpace object is compatible with this Data object.
2413 if ( (fspace.getTypeCode()!=functionSpaceType) ||
2414 (fspace.getNumSamples()!=noSamples) ||
2415 (fspace.getNumDPPSample()!=noDPPSample)
2416 ) {
2417 throw DataException("extractData Error: incompatible FunctionSpace");
2418 }
2419 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2420 if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2421 throw DataException("extractData Error: incompatible FunctionSpace");
2422 }
2423 }
2424 if (dataType==2) {
2425 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2426 if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2427 throw DataException("extractData Error: incompatible FunctionSpace");
2428 }
2429 }
2430 }
2431
2432 //
2433 // Construct a DataVector to hold underlying data values
2434 DataVector dataVec(dataLength);
2435
2436 //
2437 // Load this DataVector with the appropriate values
2438 int noValues;
2439 archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2440 cout << "\tnoValues: " << noValues << endl;
2441 switch (dataType) {
2442 case 0:
2443 // DataEmpty
2444 if (noValues != 0) {
2445 throw DataException("extractData Error: problem reading data from archive file");
2446 }
2447 break;
2448 case 1:
2449 // DataConstant
2450 if (dataVec.extractData(archiveFile,noValues)) {
2451 throw DataException("extractData Error: problem reading data from archive file");
2452 }
2453 break;
2454 case 2:
2455 // DataTagged
2456 if (dataVec.extractData(archiveFile,noValues)) {
2457 throw DataException("extractData Error: problem reading data from archive file");
2458 }
2459 break;
2460 case 3:
2461 // DataExpanded
2462 if (dataVec.extractData(archiveFile,noValues)) {
2463 throw DataException("extractData Error: problem reading data from archive file");
2464 }
2465 break;
2466 }
2467
2468 if (!archiveFile.good()) {
2469 throw DataException("extractData Error: problem reading from archive file");
2470 }
2471
2472 //
2473 // Close archive file
2474 archiveFile.close();
2475
2476 if (!archiveFile.good()) {
2477 throw DataException("extractData Error: problem closing archive file");
2478 }
2479
2480 //
2481 // Construct an appropriate Data object
2482 DataAbstract* tempData;
2483 switch (dataType) {
2484 case 0:
2485 // DataEmpty
2486 tempData=new DataEmpty();
2487 break;
2488 case 1:
2489 // DataConstant
2490 tempData=new DataConstant(fspace,dataPointShape,dataVec);
2491 break;
2492 case 2:
2493 // DataTagged
2494 tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2495 break;
2496 case 3:
2497 // DataExpanded
2498 tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2499 break;
2500 }
2501 shared_ptr<DataAbstract> temp_data(tempData);
2502 m_data=temp_data;
2503 }
2504
2505 ostream& escript::operator<<(ostream& o, const Data& data)
2506 {
2507 o << data.toString();
2508 return o;
2509 }
2510
2511 Data
2512 escript::C_GeneralTensorProduct(Data& arg_0,
2513 Data& arg_1,
2514 int axis_offset,
2515 int transpose)
2516 {
2517 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2518 // SM is the product of the last axis_offset entries in arg_0.getShape().
2519
2520 #if defined DOPROF
2521 // profData->binary++;
2522 #endif
2523
2524 // Interpolate if necessary and find an appropriate function space
2525 Data arg_0_Z, arg_1_Z;
2526 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2527 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2528 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2529 arg_1_Z = Data(arg_1);
2530 }
2531 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2532 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2533 arg_0_Z =Data(arg_0);
2534 }
2535 else {
2536 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2537 }
2538 } else {
2539 arg_0_Z = Data(arg_0);
2540 arg_1_Z = Data(arg_1);
2541 }
2542 // Get rank and shape of inputs
2543 int rank0 = arg_0_Z.getDataPointRank();
2544 int rank1 = arg_1_Z.getDataPointRank();
2545 DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2546 DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2547
2548 // Prepare for the loops of the product and verify compatibility of shapes
2549 int start0=0, start1=0;
2550 if (transpose == 0) {}
2551 else if (transpose == 1) { start0 = axis_offset; }
2552 else if (transpose == 2) { start1 = rank1-axis_offset; }
2553 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2554
2555 // Adjust the shapes for transpose
2556 DataArrayView::ShapeType tmpShape0;
2557 DataArrayView::ShapeType tmpShape1;
2558 for (int i=0; i<rank0; i++) { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2559 for (int i=0; i<rank1; i++) { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2560
2561 #if 0
2562 // For debugging: show shape after transpose
2563 char tmp[100];
2564 std::string shapeStr;
2565 shapeStr = "(";
2566 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2567 shapeStr += ")";
2568 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2569 shapeStr = "(";
2570 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2571 shapeStr += ")";
2572 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2573 #endif
2574
2575 // Prepare for the loops of the product
2576 int SL=1, SM=1, SR=1;
2577 for (int i=0; i<rank0-axis_offset; i++) {
2578 SL *= tmpShape0[i];
2579 }
2580 for (int i=rank0-axis_offset; i<rank0; i++) {
2581 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2582 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2583 }
2584 SM *= tmpShape0[i];
2585 }
2586 for (int i=axis_offset; i<rank1; i++) {
2587 SR *= tmpShape1[i];
2588 }
2589
2590 // Define the shape of the output
2591 DataArrayView::ShapeType shape2;
2592 for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2593 for (int i=axis_offset; i<rank1; i++) { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2594
2595 // Declare output Data object
2596 Data res;
2597
2598 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2599 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2600 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2601 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2602 double *ptr_2 = &((res.getPointDataView().getData())[0]);
2603 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2604 }
2605 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2606
2607 // Prepare the DataConstant input
2608 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2609 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2610
2611 // Borrow DataTagged input from Data object
2612 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2613 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2614
2615 // Prepare a DataTagged output 2
2616 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2617 res.tag();
2618 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2619 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2620
2621 // Prepare offset into DataConstant
2622 int offset_0 = tmp_0->getPointOffset(0,0);
2623 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2624 // Get the views
2625 DataArrayView view_1 = tmp_1->getDefaultValue();
2626 DataArrayView view_2 = tmp_2->getDefaultValue();
2627 // Get the pointers to the actual data
2628 double *ptr_1 = &((view_1.getData())[0]);
2629 double *ptr_2 = &((view_2.getData())[0]);
2630 // Compute an MVP for the default
2631 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2632 // Compute an MVP for each tag
2633 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2634 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2635 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2636 tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2637 DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2638 DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2639 double *ptr_1 = &view_1.getData(0);
2640 double *ptr_2 = &view_2.getData(0);
2641 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2642 }
2643
2644 }
2645 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2646
2647 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2648 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2649 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2650 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2651 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2652 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2653 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2654 int sampleNo_1,dataPointNo_1;
2655 int numSamples_1 = arg_1_Z.getNumSamples();
2656 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2657 int offset_0 = tmp_0->getPointOffset(0,0);
2658 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2659 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2660 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2661 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2662 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2663 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2664 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2665 double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2666 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2667 }
2668 }
2669
2670 }
2671 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2672
2673 // Borrow DataTagged input from Data object
2674 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2675 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2676
2677 // Prepare the DataConstant input
2678 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2679 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2680
2681 // Prepare a DataTagged output 2
2682 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2683 res.tag();
2684 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2685 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2686
2687 // Prepare offset into DataConstant
2688 int offset_1 = tmp_1->getPointOffset(0,0);
2689 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2690 // Get the views
2691 DataArrayView view_0 = tmp_0->getDefaultValue();
2692 DataArrayView view_2 = tmp_2->getDefaultValue();
2693 // Get the pointers to the actual data
2694 double *ptr_0 = &((view_0.getData())[0]);
2695 double *ptr_2 = &((view_2.getData())[0]);
2696 // Compute an MVP for the default
2697 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2698 // Compute an MVP for each tag
2699 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2700 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2701 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2702 tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2703 DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2704 DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2705 double *ptr_0 = &view_0.getData(0);
2706 double *ptr_2 = &view_2.getData(0);
2707 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2708 }
2709
2710 }
2711 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2712
2713 // Borrow DataTagged input from Data object
2714 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2715 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2716
2717 // Borrow DataTagged input from Data object
2718 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2719 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2720
2721 // Prepare a DataTagged output 2
2722 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2723 res.tag(); // DataTagged output
2724 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2725 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2726
2727 // Get the views
2728 DataArrayView view_0 = tmp_0->getDefaultValue();
2729 DataArrayView view_1 = tmp_1->getDefaultValue();
2730 DataArrayView view_2 = tmp_2->getDefaultValue();
2731 // Get the pointers to the actual data
2732 double *ptr_0 = &((view_0.getData())[0]);
2733 double *ptr_1 = &((view_1.getData())[0]);
2734 double *ptr_2 = &((view_2.getData())[0]);
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->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2743 }
2744 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2745 tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
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 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2757 }
2758
2759 }
2760 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2761
2762 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2763 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2764 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2765 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2766 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2767 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2768 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2769 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2770 int sampleNo_0,dataPointNo_0;
2771 int numSamples_0 = arg_0_Z.getNumSamples();
2772 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2773 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2774 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2775 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2776 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2777 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2778 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2779 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2780 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2781 double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2782 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2783 }
2784 }
2785
2786 }
2787 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2788
2789 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2790 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2791 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2792 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2793 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2794 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2795 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2796 int sampleNo_0,dataPointNo_0;
2797 int numSamples_0 = arg_0_Z.getNumSamples();
2798 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2799 int offset_1 = tmp_1->getPointOffset(0,0);
2800 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2801 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2802 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2803 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2804 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2805 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2806 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2807 double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2808 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2809 }
2810 }
2811
2812
2813 }
2814 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2815
2816 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2817 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2818 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2819 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2820 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2821 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2822 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2823 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2824 int sampleNo_0,dataPointNo_0;
2825 int numSamples_0 = arg_0_Z.getNumSamples();
2826 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2827 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2828 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2829 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2830 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2831 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2832 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2833 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2834 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2835 double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2836 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2837 }
2838 }
2839
2840 }
2841 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2842
2843 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2844 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2845 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2846 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2847 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2848 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2849 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2850 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2851 int sampleNo_0,dataPointNo_0;
2852 int numSamples_0 = arg_0_Z.getNumSamples();
2853 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2854 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2855 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2856 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2857 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2858 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2859 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2860 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2861 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2862 double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2863 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2864 }
2865 }
2866
2867 }
2868 else {
2869 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2870 }
2871
2872 return res;
2873 }
2874
2875 DataAbstract*
2876 Data::borrowData() const
2877 {
2878 return m_data.get();
2879 }
2880
2881 /* Member functions specific to the MPI implementation */
2882
2883 void
2884 Data::print()
2885 {
2886 int i,j;
2887
2888 printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2889 for( i=0; i<getNumSamples(); i++ )
2890 {
2891 printf( "[%6d]", i );
2892 for( j=0; j<getNumDataPointsPerSample(); j++ )
2893 printf( "\t%10.7g", (getSampleData(i))[j] );
2894 printf( "\n" );
2895 }
2896 }
2897
2898 int
2899 Data::get_MPISize() const
2900 {
2901 int error, size;
2902 #ifdef PASO_MPI
2903 error = MPI_Comm_size( get_MPIComm(), &size );
2904 #else
2905 size = 1;
2906 #endif
2907 return size;
2908 }
2909
2910 int
2911 Data::get_MPIRank() const
2912 {
2913 int error, rank;
2914 #ifdef PASO_MPI
2915 error = MPI_Comm_rank( get_MPIComm(), &rank );
2916 #else
2917 rank = 0;
2918 #endif
2919 return rank;
2920 }
2921
2922 MPI_Comm
2923 Data::get_MPIComm() const
2924 {
2925 #ifdef PASO_MPI
2926 return MPI_COMM_WORLD;
2927 #else
2928 return -1;
2929 #endif
2930 }
2931

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26