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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3317 - (show annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 11 months ago) by caltinay
File size: 28309 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26