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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1196 - (hide annotations)
Fri Jun 15 03:45:48 2007 UTC (11 years, 11 months ago) by ksteube
File size: 76223 byte(s)
Use of PAPI on solver is now enabled with papi_instrument_solver=1 in scons/ess_options.py.
Can instrument other blocks of code with blockpapi.c.
Added interval timers to grad, integrate and Assemble_PDE.

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