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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 921 - (show annotations)
Fri Jan 5 00:54:37 2007 UTC (12 years, 9 months ago) by gross
File size: 83284 byte(s)
I have done some clarification on functions that allow to access individual data point values in a Data object. 
The term "data point number" is always local on a MPI process and referes to the value (data_point_in_sample, sample)
as a single identifyer (data_point_in_sample + sample * number_data_points_per_sample). a "global data point number"
referes to a tuple of a processour id and local data point number.

The function convertToNumArrayFromSampleNo has been removed now and convertToNumArrayFromDPNo renamed to getValueOfDataPoint.
There are two new functions:

   getNumberOfDataPoints
   setValueOfDataPoint

This allows you to do things like:

  in=Data(..)
  out=Data(..)
   for i in xrange(in.getNumberOfDataPoints())
       in_loc=in.getValueOfDataPoint(i)
       out_loc=< some operations on in_loc>
       out.setValueOfDataPoint(i,out_loc)


Also mindp  is renamed to  minGlobalDataPoint and there is a new function getValueOfGlobalDataPoint. While in MPI the functions getNumberOfDataPoints and getValueOfDataPoint are working locally on each process (so the code above is executed in parallel).
the latter allows getting a single value across all processors. 


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26