/[escript]/trunk-mpi-branch/escript/src/Data.cpp
ViewVC logotype

Contents of /trunk-mpi-branch/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1087 - (show annotations)
Thu Apr 12 10:01:47 2007 UTC (12 years, 2 months ago) by gross
File size: 77910 byte(s)
the MPI version of PASO.PCG is running. 
There is a bug in the rectangular mesh generators but they need to be
revised in any case to clean up the code.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26