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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3390 - (show annotations)
Thu Dec 2 00:34:37 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 28567 byte(s)
RandomData added

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 void DataExpanded::randomFill(double seed)
820 {
821 CHECK_FOR_EX_WRITE
822 if (seed==0)
823 {
824 time_t s=time(0);
825 seed=s;
826 }
827 srand(seed);
828 DataVector& dv=getVectorRW();
829 for (long i=0;i<dv.size();++i)
830 {
831 dv[i]=(double)rand()/RAND_MAX;
832 }
833 }
834
835 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26