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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2396 - (show annotations)
Thu Apr 23 23:58:29 2009 UTC (10 years ago) by jfenwick
File size: 26748 byte(s)
Branching to port escript to use numpy rather than numarray

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
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 string
271 DataExpanded::toString() const
272 {
273 stringstream temp;
274 FunctionSpace fs=getFunctionSpace();
275
276 int offset=0;
277 for (int i=0;i<m_data.getNumRows();i++) {
278 for (int j=0;j<m_data.getNumCols();j++) {
279 offset=getPointOffset(i,j);
280 stringstream suffix;
281 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
282 temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
283 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
284 temp << endl;
285 }
286 }
287 }
288 string result=temp.str();
289 if (result.empty())
290 {
291 return "(data contains no samples)\n";
292 }
293 return temp.str();
294 }
295
296 DataTypes::ValueType::size_type
297 DataExpanded::getPointOffset(int sampleNo,
298 int dataPointNo) const
299 {
300 return m_data.index(sampleNo,dataPointNo);
301 }
302
303 DataTypes::ValueType::size_type
304 DataExpanded::getPointOffset(int sampleNo,
305 int dataPointNo)
306 {
307 return m_data.index(sampleNo,dataPointNo);
308 }
309
310 DataTypes::ValueType::size_type
311 DataExpanded::getLength() const
312 {
313 return m_data.size();
314 }
315
316
317
318 void
319 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
320 CHECK_FOR_EX_WRITE
321 //
322 // Get the number of samples and data-points per sample.
323 int numSamples = getNumSamples();
324 int numDataPointsPerSample = getNumDPPSample();
325 int dataPointRank = getRank();
326 ShapeType dataPointShape = getShape();
327 if (numSamples*numDataPointsPerSample > 0) {
328 //TODO: global error handling
329 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
330 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
331 }
332 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
333 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
334 }
335 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
336 ValueType& vec=getVectorRW();
337 if (dataPointRank==0) {
338 vec[offset]=value;
339 } else if (dataPointRank==1) {
340 for (int i=0; i<dataPointShape[0]; i++) {
341 vec[offset+i]=value;
342 }
343 } else if (dataPointRank==2) {
344 for (int i=0; i<dataPointShape[0]; i++) {
345 for (int j=0; j<dataPointShape[1]; j++) {
346 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
347 }
348 }
349 } else if (dataPointRank==3) {
350 for (int i=0; i<dataPointShape[0]; i++) {
351 for (int j=0; j<dataPointShape[1]; j++) {
352 for (int k=0; k<dataPointShape[2]; k++) {
353 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
354 }
355 }
356 }
357 } else if (dataPointRank==4) {
358 for (int i=0; i<dataPointShape[0]; i++) {
359 for (int j=0; j<dataPointShape[1]; j++) {
360 for (int k=0; k<dataPointShape[2]; k++) {
361 for (int l=0; l<dataPointShape[3]; l++) {
362 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
363 }
364 }
365 }
366 }
367 }
368 }
369 }
370
371 void
372 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
373 CHECK_FOR_EX_WRITE
374 //
375 // Get the number of samples and data-points per sample.
376 int numSamples = getNumSamples();
377 int numDataPointsPerSample = getNumDPPSample();
378 //
379 // check rank:
380 if (value.getRank()!=getRank())
381 throw DataException("Rank of numarray does not match Data object rank");
382 if (numSamples*numDataPointsPerSample > 0) {
383 //TODO: global error handling
384 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
385 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
386 }
387 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
388 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
389 }
390 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
391 ValueType& vec=getVectorRW();
392 vec.copyFromArrayToOffset(value,offset,1);
393 }
394 }
395
396 void
397 DataExpanded::symmetric(DataAbstract* ev)
398 {
399 int sampleNo,dataPointNo;
400 int numSamples = getNumSamples();
401 int numDataPointsPerSample = getNumDPPSample();
402 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
403 if (temp_ev==0) {
404 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
405 }
406 const ValueType& vec=getVectorRO();
407 const ShapeType& shape=getShape();
408 ValueType& evVec=temp_ev->getVectorRW();
409 const ShapeType& evShape=temp_ev->getShape();
410 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
411 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
412 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
413 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
414 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
415 }
416 }
417 }
418 void
419 DataExpanded::nonsymmetric(DataAbstract* ev)
420 {
421 int sampleNo,dataPointNo;
422 int numSamples = getNumSamples();
423 int numDataPointsPerSample = getNumDPPSample();
424 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
425 if (temp_ev==0) {
426 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
427 }
428 const ValueType& vec=getVectorRO();
429 const ShapeType& shape=getShape();
430 ValueType& evVec=temp_ev->getVectorRW();
431 const ShapeType& evShape=temp_ev->getShape();
432 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
433 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
434 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
435 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
436 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
437 }
438 }
439 }
440 void
441 DataExpanded::trace(DataAbstract* ev, int axis_offset)
442 {
443 int sampleNo,dataPointNo;
444 int numSamples = getNumSamples();
445 int numDataPointsPerSample = getNumDPPSample();
446 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
447 if (temp_ev==0) {
448 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
449 }
450 const ValueType& vec=getVectorRO();
451 const ShapeType& shape=getShape();
452 ValueType& evVec=temp_ev->getVectorRW();
453 const ShapeType& evShape=temp_ev->getShape();
454 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
455 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
456 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
457 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
458 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
459 }
460 }
461 }
462
463 void
464 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
465 {
466 int sampleNo,dataPointNo;
467 int numSamples = getNumSamples();
468 int numDataPointsPerSample = getNumDPPSample();
469 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
470 if (temp_ev==0) {
471 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
472 }
473 const ValueType& vec=getVectorRO();
474 const ShapeType& shape=getShape();
475 ValueType& evVec=temp_ev->getVectorRW();
476 const ShapeType& evShape=temp_ev->getShape();
477 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
478 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
479 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
480 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
481 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
482 }
483 }
484 }
485
486 void
487 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
488 {
489 int sampleNo,dataPointNo;
490 int numSamples = getNumSamples();
491 int numDataPointsPerSample = getNumDPPSample();
492 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
493 if (temp_ev==0) {
494 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
495 }
496 const ValueType& vec=getVectorRO();
497 const ShapeType& shape=getShape();
498 ValueType& evVec=temp_ev->getVectorRW();
499 const ShapeType& evShape=temp_ev->getShape();
500 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
501 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
504 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
505 }
506 }
507 }
508 void
509 DataExpanded::eigenvalues(DataAbstract* ev)
510 {
511 int sampleNo,dataPointNo;
512 int numSamples = getNumSamples();
513 int numDataPointsPerSample = getNumDPPSample();
514 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
515 if (temp_ev==0) {
516 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
517 }
518 const ValueType& vec=getVectorRO();
519 const ShapeType& shape=getShape();
520 ValueType& evVec=temp_ev->getVectorRW();
521 const ShapeType& evShape=temp_ev->getShape();
522 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
523 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
524 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
525 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
526 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
527 }
528 }
529 }
530 void
531 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
532 {
533 int numSamples = getNumSamples();
534 int numDataPointsPerSample = getNumDPPSample();
535 int sampleNo,dataPointNo;
536 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
537 if (temp_ev==0) {
538 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
539 }
540 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
541 if (temp_V==0) {
542 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
543 }
544 const ValueType& vec=getVectorRO();
545 const ShapeType& shape=getShape();
546 ValueType& evVec=temp_ev->getVectorRW();
547 const ShapeType& evShape=temp_ev->getShape();
548 ValueType& VVec=temp_V->getVectorRW();
549 const ShapeType& VShape=temp_V->getShape();
550 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
551 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
552 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
553 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
554 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
555 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
556 }
557 }
558 }
559
560 void
561 DataExpanded::setToZero(){
562 // TODO: Surely there is a more efficient way to do this????
563 // Why is there no memset here? Parallel issues?
564 CHECK_FOR_EX_WRITE
565 int numSamples = getNumSamples();
566 int numDataPointsPerSample = getNumDPPSample();
567 DataTypes::ValueType::size_type n = getNoValues();
568 double* p;
569 int sampleNo,dataPointNo, i;
570 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
571 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
572 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
573 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
574 for (i=0; i<n ;++i) p[i]=0.;
575 }
576 }
577 }
578
579 /* Append MPI rank to file name if multiple MPI processes */
580 char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
581 /* Make plenty of room for the mpi_rank number and terminating '\0' */
582 char *newFileName = (char *)malloc(strlen(fileName)+20);
583 strncpy(newFileName, fileName, strlen(fileName)+1);
584 if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
585 return(newFileName);
586 }
587
588 void
589 DataExpanded::dump(const std::string fileName) const
590 {
591 #ifdef USE_NETCDF
592 const int ldims=2+DataTypes::maxRank;
593 const NcDim* ncdims[ldims];
594 NcVar *var, *ids;
595 int rank = getRank();
596 int type= getFunctionSpace().getTypeCode();
597 int ndims =0;
598 long dims[ldims];
599 const double* d_ptr=&(m_data[0]);
600 const DataTypes::ShapeType& shape = getShape();
601 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
602 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
603 #ifdef PASO_MPI
604 MPI_Status status;
605 #endif
606
607 #ifdef PASO_MPI
608 /* Serialize NetCDF I/O */
609 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
610 #endif
611
612 // netCDF error handler
613 NcError err(NcError::verbose_nonfatal);
614 // Create the file.
615 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
616 NcFile dataFile(newFileName, NcFile::Replace);
617 // check if writing was successful
618 if (!dataFile.is_valid())
619 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
620 if (!dataFile.add_att("type_id",2) )
621 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
622 if (!dataFile.add_att("rank",rank) )
623 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
624 if (!dataFile.add_att("function_space_type",type))
625 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
626 ndims=rank+2;
627 if ( rank >0 ) {
628 dims[0]=shape[0];
629 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
630 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
631 }
632 if ( rank >1 ) {
633 dims[1]=shape[1];
634 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
635 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
636 }
637 if ( rank >2 ) {
638 dims[2]=shape[2];
639 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
640 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
641 }
642 if ( rank >3 ) {
643 dims[3]=shape[3];
644 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
645 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
646 }
647 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
648 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
649 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
650 dims[rank+1]=getFunctionSpace().getNumSamples();
651 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
652 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
653
654 if (getFunctionSpace().getNumSamples()>0)
655 {
656
657 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
658 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
659 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
660 if (! (ids->put(ids_p,dims[rank+1])) )
661 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
662 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
663 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
664 if (! (var->put(d_ptr,dims)) )
665 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
666 }
667 #ifdef PASO_MPI
668 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
669 #endif
670 #else
671 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
672 #endif
673 }
674
675 void
676 DataExpanded::setTaggedValue(int tagKey,
677 const DataTypes::ShapeType& pointshape,
678 const DataTypes::ValueType& value,
679 int dataOffset)
680 {
681 CHECK_FOR_EX_WRITE
682 int numSamples = getNumSamples();
683 int numDataPointsPerSample = getNumDPPSample();
684 int sampleNo,dataPointNo, i;
685 DataTypes::ValueType::size_type n = getNoValues();
686 double* p;
687 const double* in=&value[0+dataOffset];
688
689 if (value.size() != n) {
690 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
691 }
692
693 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
694 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
695 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
696 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
697 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
698 for (i=0; i<n ;++i) p[i]=in[i];
699 }
700 }
701 }
702 }
703
704
705 void
706 DataExpanded::reorderByReferenceIDs(int *reference_ids)
707 {
708 CHECK_FOR_EX_WRITE
709 int numSamples = getNumSamples();
710 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
711 int sampleNo, sampleNo2,i;
712 double* p,*p2;
713 register double rtmp;
714 FunctionSpace fs=getFunctionSpace();
715
716 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
717 const int id_in=reference_ids[sampleNo];
718 const int id=fs.getReferenceIDOfSample(sampleNo);
719 if (id!=id_in) {
720 bool matched=false;
721 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
722 if (id == reference_ids[sampleNo2]) {
723 p=&(m_data[getPointOffset(sampleNo,0)]);
724 p2=&(m_data[getPointOffset(sampleNo2,0)]);
725 for (i=0; i<n ;i++) {
726 rtmp=p[i];
727 p[i]=p2[i];
728 p2[i]=rtmp;
729 }
730 reference_ids[sampleNo]=id;
731 reference_ids[sampleNo2]=id_in;
732 matched=true;
733 break;
734 }
735 }
736 if (! matched) {
737 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
738 }
739 }
740 }
741 }
742
743 DataTypes::ValueType&
744 DataExpanded::getVectorRW()
745 {
746 CHECK_FOR_EX_WRITE
747 return m_data.getData();
748 }
749
750 const DataTypes::ValueType&
751 DataExpanded::getVectorRO() const
752 {
753 return m_data.getData();
754 }
755
756
757 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26