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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 813 - (show annotations)
Mon Aug 21 02:08:47 2006 UTC (12 years, 10 months ago) by ksteube
File size: 83823 byte(s)
Tensor products for Data objects are now computed by a C++ method
C_GeneralTensorProduct, which calls C function matrix_matrix_product
to do the actual calculation.

Can perform product with either input transposed in place, meaning
without first computing the transpose in a separate step.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26