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