/[escript]/branches/arrexp_trunk2098/escript/src/DataExpanded.cpp
ViewVC logotype

Contents of /branches/arrexp_trunk2098/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2119 - (show annotations)
Tue Dec 2 06:06:04 2008 UTC (11 years, 10 months ago) by jfenwick
File size: 31949 byte(s)
Branch commit.
Threading Wrapped Array through the code.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
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
15 #include "Data.h"
16 #include "DataExpanded.h"
17 #include "DataException.h"
18 #include "DataConstant.h"
19 #include "DataTagged.h"
20 #ifdef USE_NETCDF
21 #include <netcdfcpp.h>
22 #endif
23 #ifdef PASO_MPI
24 #include <mpi.h>
25 #endif
26
27 #include <boost/python/extract.hpp>
28 #include "DataMaths.h"
29
30 using namespace std;
31 using namespace boost::python;
32 using namespace boost;
33 using namespace escript::DataTypes;
34
35 namespace escript {
36
37 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
38 const FunctionSpace& what)
39 : parent(what,DataTypes::shapeFromNumArray(value))
40 {
41 //
42 // initialise the data array for this object
43 initialise(what.getNumSamples(),what.getNumDPPSample());
44 //
45 // copy the given value to every data point
46 copy(value);
47 }
48
49 DataExpanded::DataExpanded(const WrappedArray& value,
50 const FunctionSpace& what)
51 : parent(what,value.getShape())
52 {
53 //
54 // initialise the data array for this object
55 initialise(what.getNumSamples(),what.getNumDPPSample());
56 //
57 // copy the given value to every data point
58 copy(value);
59 }
60
61 DataExpanded::DataExpanded(const DataExpanded& other)
62 : parent(other.getFunctionSpace(), other.getShape()),
63 m_data(other.m_data)
64 {
65 }
66
67 DataExpanded::DataExpanded(const DataConstant& other)
68 : parent(other.getFunctionSpace(), other.getShape())
69 {
70 //
71 // initialise the data array for this object
72 initialise(other.getNumSamples(),other.getNumDPPSample());
73 //
74 // DataConstant only has one value, copy this to every data point
75 copy(other);
76 }
77
78 DataExpanded::DataExpanded(const DataTagged& other)
79 : parent(other.getFunctionSpace(), other.getShape())
80 {
81 //
82 // initialise the data array for this object
83 initialise(other.getNumSamples(),other.getNumDPPSample());
84 //
85 // for each data point in this object, extract and copy the corresponding data
86 // value from the given DataTagged object
87 int i,j;
88 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
89 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
90 #pragma omp parallel for private(i,j) schedule(static)
91 for (i=0;i<numRows;i++) {
92 for (j=0;j<numCols;j++) {
93 try {
94 DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
95 other.getVector(),
96 other.getPointOffset(i,j));
97 }
98 catch (std::exception& e) {
99 cout << e.what() << endl;
100 }
101 }
102 }
103 }
104
105 DataExpanded::DataExpanded(const DataExpanded& other,
106 const DataTypes::RegionType& region)
107 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
108 {
109 //
110 // get the shape of the slice
111 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
112 //
113 // initialise this Data object to the shape of the slice
114 initialise(other.getNumSamples(),other.getNumDPPSample());
115 //
116 // copy the data
117 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
118 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
119 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
120 int i,j;
121 #pragma omp parallel for private(i,j) schedule(static)
122 for (i=0;i<numRows;i++) {
123 for (j=0;j<numCols;j++) {
124 try {
125 // getPointDataView().copySlice(getPointOffset(i,j),
126 // other.getPointDataView(),
127 // other.getPointOffset(i,j),
128 // region_loop_range);
129 DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
130 other.getVector(),
131 other.getShape(),
132 other.getPointOffset(i,j),
133 region_loop_range);
134 }
135 catch (std::exception& e) {
136 cout << e.what() << endl;
137 }
138 }
139 }
140 }
141
142 // DataExpanded::DataExpanded(const DataArrayView& value,
143 // const FunctionSpace& what)
144 // : DataAbstract(what)
145 // {
146 // //
147 // // get the shape of the given data value
148 // DataTypes::ShapeType tempShape=value.getShape();
149 // //
150 // // initialise this Data object to the shape of the given data value
151 // initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
152 // //
153 // // copy the given value to every data point
154 // copy(value);
155 // }
156
157 DataExpanded::DataExpanded(const FunctionSpace& what,
158 const DataTypes::ShapeType &shape,
159 const DataTypes::ValueType &data)
160 : parent(what,shape)
161 {
162 EsysAssert(data.size()%getNoValues()==0,
163 "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
164
165 if (data.size()==getNoValues())
166 {
167 ValueType& vec=m_data.getData();
168 //
169 // create the view of the data
170 initialise(what.getNumSamples(),what.getNumDPPSample());
171 // now we copy this value to all elements
172 for (int i=0;i<getLength();)
173 {
174 for (unsigned int j=0;j<getNoValues();++j,++i)
175 {
176 vec[i]=data[j];
177 }
178 }
179 }
180 else
181 {
182 //
183 // copy the data in the correct format
184 m_data.getData()=data;
185 }
186
187
188 }
189
190 DataExpanded::~DataExpanded()
191 {
192 }
193
194 DataAbstract*
195 DataExpanded::deepCopy()
196 {
197 return new DataExpanded(*this);
198 }
199
200
201 DataAbstract*
202 DataExpanded::getSlice(const DataTypes::RegionType& region) const
203 {
204 return new DataExpanded(*this,region);
205 }
206
207 void
208 DataExpanded::setSlice(const DataAbstract* value,
209 const DataTypes::RegionType& region)
210 {
211 const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
212 if (tempDataExp==0) {
213 throw DataException("Programming error - casting to DataExpanded.");
214 }
215 //
216 // get shape of slice
217 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
218 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
219 //
220 // check shape
221 if (getRank()!=region.size()) {
222 throw DataException("Error - Invalid slice region.");
223 }
224 if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
225 throw DataException (DataTypes::createShapeErrorMessage(
226 "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
227 }
228 //
229 // copy the data from the slice into this object
230 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
231 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
232 int i, j;
233 ValueType& vec=getVector();
234 const ShapeType& mshape=getShape();
235 const ValueType& tVec=tempDataExp->getVector();
236 const ShapeType& tShape=tempDataExp->getShape();
237 #pragma omp parallel for private(i,j) schedule(static)
238 for (i=0;i<numRows;i++) {
239 for (j=0;j<numCols;j++) {
240 /* getPointDataView().copySliceFrom(getPointOffset(i,j),
241 tempDataExp->getPointDataView(),
242 tempDataExp->getPointOffset(i,j),
243 region_loop_range);*/
244 DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
245 tVec,
246 tShape,
247 tempDataExp->getPointOffset(i,j),
248 region_loop_range);
249
250 }
251 }
252 }
253
254 void
255 DataExpanded::copy(const DataConstant& value)
256 {
257 EsysAssert((checkShape(getShape(), value.getShape())),
258 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
259
260 //
261 // copy a single value to every data point in this object
262 int nRows=m_data.getNumRows();
263 int nCols=m_data.getNumCols();
264 int i,j;
265 #pragma omp parallel for private(i,j) schedule(static)
266 for (i=0;i<nRows;i++) {
267 for (j=0;j<nCols;j++) {
268 // NOTE: An exception may be thown from this call if
269 // DOASSERT is on which of course will play
270 // havoc with the omp threads. Run single threaded
271 // if using DOASSERT.
272 //getPointDataView().copy(getPointOffset(i,j),value);
273 DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
274 }
275 }
276 }
277
278
279
280 void
281 DataExpanded::copy(const boost::python::numeric::array& value)
282 {
283
284 // extract the shape of the numarray
285 DataTypes::ShapeType tempShape;
286 for (int i=0; i < value.getrank(); i++) {
287 tempShape.push_back(extract<int>(value.getshape()[i]));
288 }
289
290 // get the space for the data vector
291 // int len = DataTypes::noValues(tempShape);
292 // DataVector temp_data(len, 0.0, len);
293 // DataArrayView temp_dataView(temp_data, tempShape);
294 // temp_dataView.copy(value);
295
296 //
297 // check the input shape matches this shape
298 if (!DataTypes::checkShape(getShape(),tempShape)) {
299 throw DataException(DataTypes::createShapeErrorMessage(
300 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
301 tempShape,getShape()));
302 }
303 //
304 // now copy over the data
305 //copy(temp_dataView);
306 getVector().copyFromNumArray(value);
307 }
308
309 void
310 DataExpanded::copy(const WrappedArray& value)
311 {
312 // check the input shape matches this shape
313 if (!DataTypes::checkShape(getShape(),value.getShape())) {
314 throw DataException(DataTypes::createShapeErrorMessage(
315 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
316 value.getShape(),getShape()));
317 }
318 getVector().copyFromArray(value);
319 }
320
321
322 void
323 DataExpanded::initialise(int noSamples,
324 int noDataPointsPerSample)
325 {
326 if (noSamples==0) //retain the default empty object
327 {
328 return;
329 }
330 //
331 // resize data array to the required size
332 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
333 }
334
335 string
336 DataExpanded::toString() const
337 {
338 stringstream temp;
339 FunctionSpace fs=getFunctionSpace();
340
341 int offset=0;
342 for (int i=0;i<m_data.getNumRows();i++) {
343 for (int j=0;j<m_data.getNumCols();j++) {
344 offset=getPointOffset(i,j);
345 stringstream suffix;
346 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
347 temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
348 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
349 temp << endl;
350 }
351 }
352 }
353 string result=temp.str();
354 if (result.empty())
355 {
356 return "(data contains no samples)\n";
357 }
358 return temp.str();
359 }
360
361 DataTypes::ValueType::size_type
362 DataExpanded::getPointOffset(int sampleNo,
363 int dataPointNo) const
364 {
365 return m_data.index(sampleNo,dataPointNo);
366 }
367
368 DataTypes::ValueType::size_type
369 DataExpanded::getPointOffset(int sampleNo,
370 int dataPointNo)
371 {
372 return m_data.index(sampleNo,dataPointNo);
373 }
374
375 DataTypes::ValueType::size_type
376 DataExpanded::getLength() const
377 {
378 return m_data.size();
379 }
380
381
382
383 void
384 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
385 //
386 // Get the number of samples and data-points per sample.
387 int numSamples = getNumSamples();
388 int numDataPointsPerSample = getNumDPPSample();
389 int dataPointRank = getRank();
390 ShapeType dataPointShape = getShape();
391 if (numSamples*numDataPointsPerSample > 0) {
392 //TODO: global error handling
393 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
394 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
395 }
396 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
397 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
398 }
399 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
400 ValueType& vec=getVector();
401 if (dataPointRank==0) {
402 vec[offset]=value;
403 } else if (dataPointRank==1) {
404 for (int i=0; i<dataPointShape[0]; i++) {
405 vec[offset+i]=value;
406 }
407 } else if (dataPointRank==2) {
408 for (int i=0; i<dataPointShape[0]; i++) {
409 for (int j=0; j<dataPointShape[1]; j++) {
410 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
411 }
412 }
413 } else if (dataPointRank==3) {
414 for (int i=0; i<dataPointShape[0]; i++) {
415 for (int j=0; j<dataPointShape[1]; j++) {
416 for (int k=0; k<dataPointShape[2]; k++) {
417 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
418 }
419 }
420 }
421 } else if (dataPointRank==4) {
422 for (int i=0; i<dataPointShape[0]; i++) {
423 for (int j=0; j<dataPointShape[1]; j++) {
424 for (int k=0; k<dataPointShape[2]; k++) {
425 for (int l=0; l<dataPointShape[3]; l++) {
426 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
427 }
428 }
429 }
430 }
431 }
432 }
433 }
434 void
435 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
436 //
437 // Get the number of samples and data-points per sample.
438 int numSamples = getNumSamples();
439 int numDataPointsPerSample = getNumDPPSample();
440 int dataPointRank = getRank();
441 const ShapeType& shape = getShape();
442 //
443 // check rank:
444 if (value.getrank()!=dataPointRank)
445 throw DataException("Rank of numarray does not match Data object rank");
446 if (numSamples*numDataPointsPerSample > 0) {
447 //TODO: global error handling
448 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
449 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
450 }
451 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
452 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
453 }
454 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
455 ValueType& vec=getVector();
456 if (dataPointRank==0) {
457 vec[offset]=extract<double>(value[0]);
458 } else if (dataPointRank==1) {
459 for (int i=0; i<shape[0]; i++) {
460 vec[offset+i]=extract<double>(value[i]);
461 }
462 } else if (dataPointRank==2) {
463 for (int i=0; i<shape[0]; i++) {
464 for (int j=0; j<shape[1]; j++) {
465 vec[offset+getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
466 }
467 }
468 } else if (dataPointRank==3) {
469 for (int i=0; i<shape[0]; i++) {
470 for (int j=0; j<shape[1]; j++) {
471 for (int k=0; k<shape[2]; k++) {
472 vec[offset+getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
473 }
474 }
475 }
476 } else if (dataPointRank==4) {
477 for (int i=0; i<shape[0]; i++) {
478 for (int j=0; j<shape[1]; j++) {
479 for (int k=0; k<shape[2]; k++) {
480 for (int l=0; l<shape[3]; l++) {
481 vec[offset+getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
482 }
483 }
484 }
485 }
486 }
487 }
488 }
489 void
490 DataExpanded::copyAll(const boost::python::numeric::array& value) {
491 //
492 // Get the number of samples and data-points per sample.
493 int numSamples = getNumSamples();
494 int numDataPointsPerSample = getNumDPPSample();
495 int dataPointRank = getRank();
496 const ShapeType& dataPointShape = getShape();
497 //
498 // check rank:
499 if (value.getrank()!=dataPointRank+1)
500 throw DataException("Rank of numarray does not match Data object rank");
501 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
502 throw DataException("leading dimension of numarray is too small");
503 //
504 ValueType& vec=getVector();
505 int dataPoint = 0;
506 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
507 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
508 ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
509 if (dataPointRank==0) {
510 vec[offset]=extract<double>(value[dataPoint]);
511 } else if (dataPointRank==1) {
512 for (int i=0; i<dataPointShape[0]; i++) {
513 vec[offset+i]=extract<double>(value[dataPoint][i]);
514 }
515 } else if (dataPointRank==2) {
516 for (int i=0; i<dataPointShape[0]; i++) {
517 for (int j=0; j<dataPointShape[1]; j++) {
518 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
519 }
520 }
521 } else if (dataPointRank==3) {
522 for (int i=0; i<dataPointShape[0]; i++) {
523 for (int j=0; j<dataPointShape[1]; j++) {
524 for (int k=0; k<dataPointShape[2]; k++) {
525 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
526 }
527 }
528 }
529 } else if (dataPointRank==4) {
530 for (int i=0; i<dataPointShape[0]; i++) {
531 for (int j=0; j<dataPointShape[1]; j++) {
532 for (int k=0; k<dataPointShape[2]; k++) {
533 for (int l=0; l<dataPointShape[3]; l++) {
534 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
535 }
536 }
537 }
538 }
539 }
540 dataPoint++;
541 }
542 }
543 }
544 void
545 DataExpanded::symmetric(DataAbstract* ev)
546 {
547 int sampleNo,dataPointNo;
548 int numSamples = getNumSamples();
549 int numDataPointsPerSample = getNumDPPSample();
550 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
551 if (temp_ev==0) {
552 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
553 }
554 ValueType& vec=getVector();
555 const ShapeType& shape=getShape();
556 ValueType& evVec=temp_ev->getVector();
557 const ShapeType& evShape=temp_ev->getShape();
558 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
559 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
560 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
561 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
562 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
563 }
564 }
565 }
566 void
567 DataExpanded::nonsymmetric(DataAbstract* ev)
568 {
569 int sampleNo,dataPointNo;
570 int numSamples = getNumSamples();
571 int numDataPointsPerSample = getNumDPPSample();
572 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
573 if (temp_ev==0) {
574 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
575 }
576 ValueType& vec=getVector();
577 const ShapeType& shape=getShape();
578 ValueType& evVec=temp_ev->getVector();
579 const ShapeType& evShape=temp_ev->getShape();
580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
581 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
582 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
583 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
584 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
585 }
586 }
587 }
588 void
589 DataExpanded::trace(DataAbstract* ev, int axis_offset)
590 {
591 int sampleNo,dataPointNo;
592 int numSamples = getNumSamples();
593 int numDataPointsPerSample = getNumDPPSample();
594 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
595 if (temp_ev==0) {
596 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
597 }
598 ValueType& vec=getVector();
599 const ShapeType& shape=getShape();
600 ValueType& evVec=temp_ev->getVector();
601 const ShapeType& evShape=temp_ev->getShape();
602 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
603 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
604 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
605 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
606 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
607 }
608 }
609 }
610
611 void
612 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
613 {
614 int sampleNo,dataPointNo;
615 int numSamples = getNumSamples();
616 int numDataPointsPerSample = getNumDPPSample();
617 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
618 if (temp_ev==0) {
619 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
620 }
621 ValueType& vec=getVector();
622 const ShapeType& shape=getShape();
623 ValueType& evVec=temp_ev->getVector();
624 const ShapeType& evShape=temp_ev->getShape();
625 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
626 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
627 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
628 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
629 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
630 }
631 }
632 }
633
634 void
635 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
636 {
637 int sampleNo,dataPointNo;
638 int numSamples = getNumSamples();
639 int numDataPointsPerSample = getNumDPPSample();
640 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
641 if (temp_ev==0) {
642 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
643 }
644 ValueType& vec=getVector();
645 const ShapeType& shape=getShape();
646 ValueType& evVec=temp_ev->getVector();
647 const ShapeType& evShape=temp_ev->getShape();
648 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
649 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
650 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
651 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
652 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
653 }
654 }
655 }
656 void
657 DataExpanded::eigenvalues(DataAbstract* ev)
658 {
659 int sampleNo,dataPointNo;
660 int numSamples = getNumSamples();
661 int numDataPointsPerSample = getNumDPPSample();
662 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
663 if (temp_ev==0) {
664 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
665 }
666 ValueType& vec=getVector();
667 const ShapeType& shape=getShape();
668 ValueType& evVec=temp_ev->getVector();
669 const ShapeType& evShape=temp_ev->getShape();
670 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
671 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
672 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
673 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
674 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
675 }
676 }
677 }
678 void
679 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
680 {
681 int numSamples = getNumSamples();
682 int numDataPointsPerSample = getNumDPPSample();
683 int sampleNo,dataPointNo;
684 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
685 if (temp_ev==0) {
686 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
687 }
688 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
689 if (temp_V==0) {
690 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
691 }
692 ValueType& vec=getVector();
693 const ShapeType& shape=getShape();
694 ValueType& evVec=temp_ev->getVector();
695 const ShapeType& evShape=temp_ev->getShape();
696 ValueType& VVec=temp_V->getVector();
697 const ShapeType& VShape=temp_V->getShape();
698 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
699 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
700 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
701 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
702 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
703 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
704 }
705 }
706 }
707
708 void
709 DataExpanded::setToZero(){
710 // TODO: Surely there is a more efficient way to do this????
711 // Why is there no memset here? Parallel issues?
712 int numSamples = getNumSamples();
713 int numDataPointsPerSample = getNumDPPSample();
714 DataTypes::ValueType::size_type n = getNoValues();
715 double* p;
716 int sampleNo,dataPointNo, i;
717 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
718 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
719 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
720 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
721 for (i=0; i<n ;++i) p[i]=0.;
722 }
723 }
724 }
725
726 /* Append MPI rank to file name if multiple MPI processes */
727 char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
728 /* Make plenty of room for the mpi_rank number and terminating '\0' */
729 char *newFileName = (char *)malloc(strlen(fileName)+20);
730 strncpy(newFileName, fileName, strlen(fileName)+1);
731 if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
732 return(newFileName);
733 }
734
735 void
736 DataExpanded::dump(const std::string fileName) const
737 {
738 #ifdef USE_NETCDF
739 const int ldims=2+DataTypes::maxRank;
740 const NcDim* ncdims[ldims];
741 NcVar *var, *ids;
742 int rank = getRank();
743 int type= getFunctionSpace().getTypeCode();
744 int ndims =0;
745 long dims[ldims];
746 const double* d_ptr=&(m_data[0]);
747 const DataTypes::ShapeType& shape = getShape();
748 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
749 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
750 #ifdef PASO_MPI
751 MPI_Status status;
752 #endif
753
754 #ifdef PASO_MPI
755 /* Serialize NetCDF I/O */
756 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
757 #endif
758
759 // netCDF error handler
760 NcError err(NcError::verbose_nonfatal);
761 // Create the file.
762 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
763 NcFile dataFile(newFileName, NcFile::Replace);
764 // check if writing was successful
765 if (!dataFile.is_valid())
766 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
767 if (!dataFile.add_att("type_id",2) )
768 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
769 if (!dataFile.add_att("rank",rank) )
770 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
771 if (!dataFile.add_att("function_space_type",type))
772 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
773 ndims=rank+2;
774 if ( rank >0 ) {
775 dims[0]=shape[0];
776 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
777 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
778 }
779 if ( rank >1 ) {
780 dims[1]=shape[1];
781 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
782 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
783 }
784 if ( rank >2 ) {
785 dims[2]=shape[2];
786 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
787 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
788 }
789 if ( rank >3 ) {
790 dims[3]=shape[3];
791 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
792 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
793 }
794 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
795 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
796 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
797 dims[rank+1]=getFunctionSpace().getNumSamples();
798 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
799 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
800
801 if (getFunctionSpace().getNumSamples()>0)
802 {
803
804 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
805 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
806 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
807 if (! (ids->put(ids_p,dims[rank+1])) )
808 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
809 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
810 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
811 if (! (var->put(d_ptr,dims)) )
812 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
813 }
814 #ifdef PASO_MPI
815 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
816 #endif
817 #else
818 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
819 #endif
820 }
821
822 void
823 DataExpanded::setTaggedValue(int tagKey,
824 const DataTypes::ShapeType& pointshape,
825 const DataTypes::ValueType& value,
826 int dataOffset)
827 {
828 int numSamples = getNumSamples();
829 int numDataPointsPerSample = getNumDPPSample();
830 int sampleNo,dataPointNo, i;
831 DataTypes::ValueType::size_type n = getNoValues();
832 double* p;
833 const double* in=&value[0+dataOffset];
834
835 if (value.size() != n) {
836 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
837 }
838
839 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
840 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
841 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
842 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
843 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
844 for (i=0; i<n ;++i) p[i]=in[i];
845 }
846 }
847 }
848 }
849
850
851 void
852 DataExpanded::reorderByReferenceIDs(int *reference_ids)
853 {
854 int numSamples = getNumSamples();
855 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
856 int sampleNo, sampleNo2,i;
857 double* p,*p2;
858 register double rtmp;
859 FunctionSpace fs=getFunctionSpace();
860
861 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
862 const int id_in=reference_ids[sampleNo];
863 const int id=fs.getReferenceIDOfSample(sampleNo);
864 if (id!=id_in) {
865 bool matched=false;
866 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
867 if (id == reference_ids[sampleNo2]) {
868 p=&(m_data[getPointOffset(sampleNo,0)]);
869 p2=&(m_data[getPointOffset(sampleNo2,0)]);
870 for (i=0; i<n ;i++) {
871 rtmp=p[i];
872 p[i]=p2[i];
873 p2[i]=rtmp;
874 }
875 reference_ids[sampleNo]=id;
876 reference_ids[sampleNo2]=id_in;
877 matched=true;
878 break;
879 }
880 }
881 if (! matched) {
882 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
883 }
884 }
885 }
886 }
887
888 DataTypes::ValueType&
889 DataExpanded::getVector()
890 {
891 return m_data.getData();
892 }
893
894 const DataTypes::ValueType&
895 DataExpanded::getVector() const
896 {
897 return m_data.getData();
898 }
899
900 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26