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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1031 - (show annotations)
Wed Mar 14 06:03:21 2007 UTC (12 years, 8 months ago) by phornby
File size: 81959 byte(s)
Finer control of the inverse hyp. functions.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26