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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26