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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2271 - (hide annotations)
Mon Feb 16 05:08:29 2009 UTC (10 years, 8 months ago) by jfenwick
File size: 26748 byte(s)
Merging version 2269 to trunk

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
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 ksteube 1312
14 ksteube 1811
15 ksteube 1800 #include "Data.h"
16 jgs 480 #include "DataExpanded.h"
17 jgs 474 #include "DataException.h"
18     #include "DataConstant.h"
19     #include "DataTagged.h"
20 gross 1023 #ifdef USE_NETCDF
21 ksteube 1312 #include <netcdfcpp.h>
22 gross 1023 #endif
23 ksteube 1748 #ifdef PASO_MPI
24     #include <mpi.h>
25     #endif
26 jgs 82
27     #include <boost/python/extract.hpp>
28 jfenwick 1796 #include "DataMaths.h"
29 jgs 82
30     using namespace std;
31     using namespace boost::python;
32     using namespace boost;
33 jfenwick 1796 using namespace escript::DataTypes;
34 jgs 82
35 jfenwick 2271
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 jgs 82 namespace escript {
41    
42 jfenwick 2271 DataExpanded::DataExpanded(const WrappedArray& value,
43 jgs 102 const FunctionSpace& what)
44 jfenwick 2271 : parent(what,value.getShape())
45 jgs 102 {
46     //
47 jgs 117 // initialise the data array for this object
48 jfenwick 1796 initialise(what.getNumSamples(),what.getNumDPPSample());
49 jgs 102 //
50     // copy the given value to every data point
51     copy(value);
52     }
53 jgs 82
54 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other)
55 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape()),
56 jgs 102 m_data(other.m_data)
57 jgs 117 {
58 jgs 102 }
59 jgs 82
60 jgs 102 DataExpanded::DataExpanded(const DataConstant& other)
61 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape())
62 jgs 117 {
63     //
64     // initialise the data array for this object
65 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
66 jgs 102 //
67 jgs 117 // DataConstant only has one value, copy this to every data point
68 jfenwick 1796 copy(other);
69 jgs 102 }
70 jgs 82
71 jgs 102 DataExpanded::DataExpanded(const DataTagged& other)
72 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape())
73 jgs 117 {
74     //
75 jgs 119 // initialise the data array for this object
76 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
77 jgs 117 //
78     // for each data point in this object, extract and copy the corresponding data
79     // value from the given DataTagged object
80 jgs 102 int i,j;
81 jfenwick 1796 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
82     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
83 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
84 jgs 117 for (i=0;i<numRows;i++) {
85     for (j=0;j<numCols;j++) {
86 jgs 102 try {
87 jfenwick 2271 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
88     other.getVectorRO(),
89 jgs 122 other.getPointOffset(i,j));
90 jgs 82 }
91 jgs 102 catch (std::exception& e) {
92     cout << e.what() << endl;
93     }
94 jgs 82 }
95     }
96 jgs 102 }
97 jgs 82
98 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other,
99 jfenwick 1796 const DataTypes::RegionType& region)
100 jfenwick 2005 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
101 jgs 102 {
102     //
103     // get the shape of the slice
104 jfenwick 1796 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
105 jgs 102 //
106     // initialise this Data object to the shape of the slice
107 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
108 jgs 102 //
109 jgs 117 // copy the data
110 jfenwick 1796 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 jgs 102 int i,j;
114 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
115 jgs 117 for (i=0;i<numRows;i++) {
116     for (j=0;j<numCols;j++) {
117 jgs 102 try {
118 jfenwick 2271 DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
119     other.getVectorRO(),
120 jfenwick 1796 other.getShape(),
121 jgs 122 other.getPointOffset(i,j),
122     region_loop_range);
123 jgs 82 }
124 jgs 102 catch (std::exception& e) {
125     cout << e.what() << endl;
126     }
127 jgs 82 }
128     }
129 jgs 102 }
130 jgs 82
131 jgs 119 DataExpanded::DataExpanded(const FunctionSpace& what,
132 jfenwick 1796 const DataTypes::ShapeType &shape,
133     const DataTypes::ValueType &data)
134 jfenwick 2005 : parent(what,shape)
135 jgs 119 {
136 jfenwick 1796 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 phornby 2008 for (unsigned int j=0;j<getNoValues();++j,++i)
149 jfenwick 1796 {
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 jgs 119 }
163    
164 jgs 102 DataExpanded::~DataExpanded()
165     {
166     }
167 jgs 82
168 jgs 102 DataAbstract*
169 jfenwick 1799 DataExpanded::deepCopy()
170     {
171     return new DataExpanded(*this);
172     }
173    
174    
175     DataAbstract*
176 jfenwick 1796 DataExpanded::getSlice(const DataTypes::RegionType& region) const
177 jgs 102 {
178     return new DataExpanded(*this,region);
179     }
180    
181     void
182     DataExpanded::setSlice(const DataAbstract* value,
183 jfenwick 1796 const DataTypes::RegionType& region)
184 jgs 102 {
185     const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
186     if (tempDataExp==0) {
187     throw DataException("Programming error - casting to DataExpanded.");
188     }
189 jfenwick 2271 CHECK_FOR_EX_WRITE
190 jgs 108 //
191 jgs 117 // get shape of slice
192 jfenwick 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
193     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
194 jgs 108 //
195 jgs 117 // check shape
196 jfenwick 1796 if (getRank()!=region.size()) {
197 jgs 102 throw DataException("Error - Invalid slice region.");
198     }
199 jfenwick 1796 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 jgs 102 }
203     //
204 jgs 117 // copy the data from the slice into this object
205 jfenwick 1796 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
206     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
207 jgs 102 int i, j;
208 jfenwick 2271 ValueType& vec=getVectorRW();
209 jfenwick 1796 const ShapeType& mshape=getShape();
210 jfenwick 2271 const ValueType& tVec=tempDataExp->getVectorRO();
211 jfenwick 1796 const ShapeType& tShape=tempDataExp->getShape();
212 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
213 jgs 102 for (i=0;i<numRows;i++) {
214     for (j=0;j<numCols;j++) {
215 jfenwick 1796 DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
216     tVec,
217     tShape,
218     tempDataExp->getPointOffset(i,j),
219 jgs 122 region_loop_range);
220 jfenwick 1796
221 jgs 82 }
222     }
223 jgs 102 }
224 jgs 82
225 jgs 102 void
226 jfenwick 1796 DataExpanded::copy(const DataConstant& value)
227 jgs 102 {
228 jfenwick 1796 EsysAssert((checkShape(getShape(), value.getShape())),
229     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
230    
231 jgs 102 //
232 jgs 117 // copy a single value to every data point in this object
233 jgs 102 int nRows=m_data.getNumRows();
234     int nCols=m_data.getNumCols();
235 jgs 117 int i,j;
236 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
237 jgs 117 for (i=0;i<nRows;i++) {
238     for (j=0;j<nCols;j++) {
239 jfenwick 2271 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(), value.getVectorRO(), 0);
240 jgs 82 }
241     }
242 jgs 102 }
243 jgs 82
244 jfenwick 2271 void
245     DataExpanded::copy(const WrappedArray& value)
246 jgs 102 {
247 jgs 117 // check the input shape matches this shape
248 jfenwick 2271 if (!DataTypes::checkShape(getShape(),value.getShape())) {
249 jfenwick 1796 throw DataException(DataTypes::createShapeErrorMessage(
250 jgs 102 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
251 jfenwick 2271 value.getShape(),getShape()));
252 jgs 82 }
253 jfenwick 2271 getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
254 jgs 102 }
255 jgs 82
256 jfenwick 1796
257 jgs 102 void
258 jfenwick 1796 DataExpanded::initialise(int noSamples,
259 jgs 102 int noDataPointsPerSample)
260     {
261 jfenwick 1808 if (noSamples==0) //retain the default empty object
262     {
263     return;
264     }
265 jgs 102 //
266 jgs 117 // resize data array to the required size
267 jfenwick 1796 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
268 jgs 102 }
269 jgs 82
270 jgs 102 string
271     DataExpanded::toString() const
272     {
273     stringstream temp;
274 gross 856 FunctionSpace fs=getFunctionSpace();
275 jfenwick 1796
276     int offset=0;
277 jgs 117 for (int i=0;i<m_data.getNumRows();i++) {
278     for (int j=0;j<m_data.getNumCols();j++) {
279 jfenwick 1796 offset=getPointOffset(i,j);
280 jgs 102 stringstream suffix;
281 gross 964 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
282 jfenwick 2271 temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
283 jgs 102 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
284     temp << endl;
285 jgs 82 }
286     }
287     }
288 jfenwick 1808 string result=temp.str();
289     if (result.empty())
290     {
291     return "(data contains no samples)\n";
292     }
293 ksteube 1312 return temp.str();
294 jgs 102 }
295 jgs 82
296 jfenwick 1796 DataTypes::ValueType::size_type
297 jgs 102 DataExpanded::getPointOffset(int sampleNo,
298     int dataPointNo) const
299     {
300     return m_data.index(sampleNo,dataPointNo);
301     }
302 jgs 82
303 jfenwick 1796 DataTypes::ValueType::size_type
304 jfenwick 2005 DataExpanded::getPointOffset(int sampleNo,
305     int dataPointNo)
306     {
307     return m_data.index(sampleNo,dataPointNo);
308     }
309    
310     DataTypes::ValueType::size_type
311 jgs 102 DataExpanded::getLength() const
312     {
313 jgs 117 return m_data.size();
314 jgs 102 }
315 jgs 82
316 jgs 123
317    
318 matt 1319 void
319 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
320 jfenwick 2271 CHECK_FOR_EX_WRITE
321 gross 922 //
322     // Get the number of samples and data-points per sample.
323     int numSamples = getNumSamples();
324     int numDataPointsPerSample = getNumDPPSample();
325 jfenwick 1796 int dataPointRank = getRank();
326     ShapeType dataPointShape = getShape();
327 gross 922 if (numSamples*numDataPointsPerSample > 0) {
328     //TODO: global error handling
329 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
330 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
331     }
332 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
333 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
334     }
335 jfenwick 1796 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
336 jfenwick 2271 ValueType& vec=getVectorRW();
337 gross 922 if (dataPointRank==0) {
338 gross 1846 vec[offset]=value;
339 gross 922 } else if (dataPointRank==1) {
340     for (int i=0; i<dataPointShape[0]; i++) {
341 jfenwick 1796 vec[offset+i]=value;
342 gross 922 }
343     } else if (dataPointRank==2) {
344     for (int i=0; i<dataPointShape[0]; i++) {
345     for (int j=0; j<dataPointShape[1]; j++) {
346 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
347 gross 922 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
354 gross 922 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
363 gross 922 }
364     }
365     }
366     }
367 matt 1319 }
368 gross 922 }
369     }
370 jfenwick 2271
371 matt 1319 void
372 jfenwick 2271 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
373     CHECK_FOR_EX_WRITE
374 gross 921 //
375     // Get the number of samples and data-points per sample.
376     int numSamples = getNumSamples();
377     int numDataPointsPerSample = getNumDPPSample();
378     //
379     // check rank:
380 jfenwick 2271 if (value.getRank()!=getRank())
381 gross 921 throw DataException("Rank of numarray does not match Data object rank");
382     if (numSamples*numDataPointsPerSample > 0) {
383     //TODO: global error handling
384 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
385 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
386     }
387 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
388 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
389     }
390 gross 1962 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
391 jfenwick 2271 ValueType& vec=getVectorRW();
392     vec.copyFromArrayToOffset(value,offset,1);
393 gross 921 }
394     }
395 jfenwick 2271
396 jgs 126 void
397 ksteube 775 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 jfenwick 2271 const ValueType& vec=getVectorRO();
407 jfenwick 1796 const ShapeType& shape=getShape();
408 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
409 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
410 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
411     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
412     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
413 jfenwick 1796 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
414     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
415 ksteube 775 }
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 jfenwick 2271 const ValueType& vec=getVectorRO();
429 jfenwick 1796 const ShapeType& shape=getShape();
430 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
431 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
432 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
433     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
434     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
435 jfenwick 1796 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
436     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
437 ksteube 775 }
438     }
439     }
440     void
441 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
442 ksteube 775 {
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 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
449 ksteube 775 }
450 jfenwick 2271 const ValueType& vec=getVectorRO();
451 jfenwick 1796 const ShapeType& shape=getShape();
452 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
453 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
454 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
455     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
456     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
457 jfenwick 1796 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
458     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
459 ksteube 775 }
460     }
461     }
462 gross 800
463 ksteube 775 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 jfenwick 2271 const ValueType& vec=getVectorRO();
474 jfenwick 1796 const ShapeType& shape=getShape();
475 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
476 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
477 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
478     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
479     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
480 jfenwick 1796 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
481     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
482 ksteube 775 }
483     }
484     }
485 gross 800
486 ksteube 775 void
487 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
488 gross 800 {
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 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
495 gross 800 }
496 jfenwick 2271 const ValueType& vec=getVectorRO();
497 jfenwick 1796 const ShapeType& shape=getShape();
498 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
499 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
500 gross 800 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
501     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503 jfenwick 1796 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
504     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
505 gross 800 }
506     }
507     }
508     void
509 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
510     {
511 gross 584 int sampleNo,dataPointNo;
512 gross 580 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 jfenwick 2271 const ValueType& vec=getVectorRO();
519 jfenwick 1796 const ShapeType& shape=getShape();
520 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
521 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
522 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
523 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
524     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
525 jfenwick 1796 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
526     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
527 gross 580 }
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 gross 584 int sampleNo,dataPointNo;
536 gross 580 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 jfenwick 2271 const ValueType& vec=getVectorRO();
545 jfenwick 1796 const ShapeType& shape=getShape();
546 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
547 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
548 jfenwick 2271 ValueType& VVec=temp_V->getVectorRW();
549 jfenwick 1796 const ShapeType& VShape=temp_V->getShape();
550 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
551 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
552     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
553 jfenwick 1796 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 gross 580 }
557     }
558     }
559    
560 gross 950 void
561 gross 1118 DataExpanded::setToZero(){
562 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
563     // Why is there no memset here? Parallel issues?
564 jfenwick 2271 CHECK_FOR_EX_WRITE
565 gross 1118 int numSamples = getNumSamples();
566     int numDataPointsPerSample = getNumDPPSample();
567 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
568 gross 1118 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 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
575 gross 1118 }
576     }
577     }
578    
579 ksteube 1748 /* 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 gross 1118
588     void
589 gross 950 DataExpanded::dump(const std::string fileName) const
590     {
591 gross 1023 #ifdef USE_NETCDF
592 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
593 gross 967 const NcDim* ncdims[ldims];
594     NcVar *var, *ids;
595 jfenwick 1796 int rank = getRank();
596 gross 967 int type= getFunctionSpace().getTypeCode();
597     int ndims =0;
598     long dims[ldims];
599 gross 1141 const double* d_ptr=&(m_data[0]);
600 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
601 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
602     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
603 ksteube 1800 #ifdef PASO_MPI
604     MPI_Status status;
605     #endif
606 gross 965
607 ksteube 1827 #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 gross 967 // netCDF error handler
613     NcError err(NcError::verbose_nonfatal);
614     // Create the file.
615 ksteube 1748 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
616     NcFile dataFile(newFileName, NcFile::Replace);
617 gross 967 // check if writing was successful
618     if (!dataFile.is_valid())
619     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
620 gross 1141 if (!dataFile.add_att("type_id",2) )
621 gross 967 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 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
631 gross 967 }
632     if ( rank >1 ) {
633     dims[1]=shape[1];
634     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
635 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
636 gross 967 }
637     if ( rank >2 ) {
638     dims[2]=shape[2];
639     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
640 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
641 gross 967 }
642     if ( rank >3 ) {
643     dims[3]=shape[3];
644     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
645 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
646 gross 967 }
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 jfenwick 1808 if (getFunctionSpace().getNumSamples()>0)
655     {
656    
657     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
658 gross 967 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
659 jfenwick 1808 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
660     if (! (ids->put(ids_p,dims[rank+1])) )
661 gross 967 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
662 jfenwick 1808 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
663 gross 967 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
664 jfenwick 1808 if (! (var->put(d_ptr,dims)) )
665 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
666 jfenwick 1808 }
667 ksteube 1800 #ifdef PASO_MPI
668 ksteube 1801 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
669 ksteube 1800 #endif
670 gross 1023 #else
671     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
672     #endif
673 gross 950 }
674 gross 1358
675 jfenwick 1796 void
676 gross 1358 DataExpanded::setTaggedValue(int tagKey,
677 jfenwick 1796 const DataTypes::ShapeType& pointshape,
678     const DataTypes::ValueType& value,
679     int dataOffset)
680 gross 1358 {
681 jfenwick 2271 CHECK_FOR_EX_WRITE
682 gross 1358 int numSamples = getNumSamples();
683     int numDataPointsPerSample = getNumDPPSample();
684     int sampleNo,dataPointNo, i;
685 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
686     double* p;
687     const double* in=&value[0+dataOffset];
688 gross 1358
689 jfenwick 1796 if (value.size() != n) {
690 gross 1358 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 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
699 gross 1358 }
700     }
701     }
702     }
703    
704 jfenwick 1796
705 gross 1487 void
706     DataExpanded::reorderByReferenceIDs(int *reference_ids)
707     {
708 jfenwick 2271 CHECK_FOR_EX_WRITE
709 gross 1487 int numSamples = getNumSamples();
710 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
711 gross 1487 int sampleNo, sampleNo2,i;
712     double* p,*p2;
713     register double rtmp;
714     FunctionSpace fs=getFunctionSpace();
715 gross 1358
716 gross 1487 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 gross 1513 for (i=0; i<n ;i++) {
726 gross 1487 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 phornby 1628 if (! matched) {
737 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
738     }
739     }
740     }
741     }
742    
743 jfenwick 1796 DataTypes::ValueType&
744 jfenwick 2271 DataExpanded::getVectorRW()
745 jfenwick 1796 {
746 jfenwick 2271 CHECK_FOR_EX_WRITE
747 jfenwick 1796 return m_data.getData();
748     }
749 gross 1487
750 jfenwick 1796 const DataTypes::ValueType&
751 jfenwick 2271 DataExpanded::getVectorRO() const
752 jfenwick 1796 {
753     return m_data.getData();
754     }
755    
756 jfenwick 2271
757 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26