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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 8 months ago) by jfenwick
File size: 28315 byte(s)
Don't panic.
Updating copyright stamps

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 jfenwick 2766 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 jgs 102 string
286     DataExpanded::toString() const
287     {
288     stringstream temp;
289 gross 856 FunctionSpace fs=getFunctionSpace();
290 jfenwick 1796
291     int offset=0;
292 jgs 117 for (int i=0;i<m_data.getNumRows();i++) {
293     for (int j=0;j<m_data.getNumCols();j++) {
294 jfenwick 1796 offset=getPointOffset(i,j);
295 jgs 102 stringstream suffix;
296 gross 964 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
297 jfenwick 2271 temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
298 jgs 102 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
299     temp << endl;
300 jgs 82 }
301     }
302     }
303 jfenwick 1808 string result=temp.str();
304     if (result.empty())
305     {
306     return "(data contains no samples)\n";
307     }
308 ksteube 1312 return temp.str();
309 jgs 102 }
310 jgs 82
311 jfenwick 1796 DataTypes::ValueType::size_type
312 jgs 102 DataExpanded::getPointOffset(int sampleNo,
313     int dataPointNo) const
314     {
315     return m_data.index(sampleNo,dataPointNo);
316     }
317 jgs 82
318 jfenwick 1796 DataTypes::ValueType::size_type
319 jfenwick 2005 DataExpanded::getPointOffset(int sampleNo,
320     int dataPointNo)
321     {
322     return m_data.index(sampleNo,dataPointNo);
323     }
324    
325     DataTypes::ValueType::size_type
326 jgs 102 DataExpanded::getLength() const
327     {
328 jgs 117 return m_data.size();
329 jgs 102 }
330 jgs 82
331 jgs 123
332    
333 matt 1319 void
334 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
335 jfenwick 2271 CHECK_FOR_EX_WRITE
336 gross 922 //
337     // Get the number of samples and data-points per sample.
338     int numSamples = getNumSamples();
339     int numDataPointsPerSample = getNumDPPSample();
340 jfenwick 1796 int dataPointRank = getRank();
341     ShapeType dataPointShape = getShape();
342 gross 922 if (numSamples*numDataPointsPerSample > 0) {
343     //TODO: global error handling
344 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
345 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
346     }
347 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
348 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
349     }
350 jfenwick 1796 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
351 jfenwick 2271 ValueType& vec=getVectorRW();
352 gross 922 if (dataPointRank==0) {
353 gross 1846 vec[offset]=value;
354 gross 922 } else if (dataPointRank==1) {
355     for (int i=0; i<dataPointShape[0]; i++) {
356 jfenwick 1796 vec[offset+i]=value;
357 gross 922 }
358     } else if (dataPointRank==2) {
359     for (int i=0; i<dataPointShape[0]; i++) {
360     for (int j=0; j<dataPointShape[1]; j++) {
361 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
362 gross 922 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
369 gross 922 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
378 gross 922 }
379     }
380     }
381     }
382 matt 1319 }
383 gross 922 }
384     }
385 jfenwick 2271
386 matt 1319 void
387 jfenwick 2271 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
388     CHECK_FOR_EX_WRITE
389 gross 921 //
390     // Get the number of samples and data-points per sample.
391     int numSamples = getNumSamples();
392     int numDataPointsPerSample = getNumDPPSample();
393     //
394     // check rank:
395 jfenwick 2271 if (value.getRank()!=getRank())
396 jfenwick 2458 throw DataException("Rank of value does not match Data object rank");
397 gross 921 if (numSamples*numDataPointsPerSample > 0) {
398     //TODO: global error handling
399 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
400 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
401     }
402 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
403 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
404     }
405 gross 1962 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
406 jfenwick 2271 ValueType& vec=getVectorRW();
407     vec.copyFromArrayToOffset(value,offset,1);
408 gross 921 }
409     }
410 jfenwick 2271
411 jgs 126 void
412 ksteube 775 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 jfenwick 2271 const ValueType& vec=getVectorRO();
422 jfenwick 1796 const ShapeType& shape=getShape();
423 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
424 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
425 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
426     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
427     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
428 jfenwick 1796 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
429     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
430 ksteube 775 }
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 jfenwick 2271 const ValueType& vec=getVectorRO();
444 jfenwick 1796 const ShapeType& shape=getShape();
445 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
446 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
447 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
448     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
449     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
450 jfenwick 1796 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
451     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
452 ksteube 775 }
453     }
454     }
455     void
456 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
457 ksteube 775 {
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 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
464 ksteube 775 }
465 jfenwick 2271 const ValueType& vec=getVectorRO();
466 jfenwick 1796 const ShapeType& shape=getShape();
467 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
468 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
469 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
470     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
471     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
472 jfenwick 1796 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
473     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
474 ksteube 775 }
475     }
476     }
477 gross 800
478 ksteube 775 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 jfenwick 2271 const ValueType& vec=getVectorRO();
489 jfenwick 1796 const ShapeType& shape=getShape();
490 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
491 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
492 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
493     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
494     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
495 jfenwick 1796 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
496     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
497 ksteube 775 }
498     }
499     }
500 gross 800
501 ksteube 775 void
502 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
503 gross 800 {
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 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
510 gross 800 }
511 jfenwick 2271 const ValueType& vec=getVectorRO();
512 jfenwick 1796 const ShapeType& shape=getShape();
513 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
514 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
515 gross 800 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
516     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
517     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
518 jfenwick 1796 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
519     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
520 gross 800 }
521     }
522     }
523     void
524 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
525     {
526 gross 584 int sampleNo,dataPointNo;
527 gross 580 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 jfenwick 2271 const ValueType& vec=getVectorRO();
534 jfenwick 1796 const ShapeType& shape=getShape();
535 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
536 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
537 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
538 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
539     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
540 jfenwick 1796 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
541     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
542 gross 580 }
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 gross 584 int sampleNo,dataPointNo;
551 gross 580 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 jfenwick 2271 const ValueType& vec=getVectorRO();
560 jfenwick 1796 const ShapeType& shape=getShape();
561 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
562 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
563 jfenwick 2271 ValueType& VVec=temp_V->getVectorRW();
564 jfenwick 1796 const ShapeType& VShape=temp_V->getShape();
565 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
566 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
567     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
568 jfenwick 1796 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 gross 580 }
572     }
573     }
574    
575 jfenwick 2742
576 jfenwick 2792 int
577 jfenwick 2742 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 jfenwick 2792 return errcode;
616 jfenwick 2742 if (errcode)
617     {
618     DataMaths::matrixInverseError(errcode); // throws exceptions
619     }
620     }
621    
622     void
623 gross 1118 DataExpanded::setToZero(){
624 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
625     // Why is there no memset here? Parallel issues?
626 jfenwick 2742 // A: This ensures that memory is touched by the correct thread.
627 jfenwick 2271 CHECK_FOR_EX_WRITE
628 gross 1118 int numSamples = getNumSamples();
629     int numDataPointsPerSample = getNumDPPSample();
630 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
631 gross 1118 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 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
638 gross 1118 }
639     }
640     }
641    
642 ksteube 1748 /* 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 gross 1118
651     void
652 gross 950 DataExpanded::dump(const std::string fileName) const
653     {
654 gross 1023 #ifdef USE_NETCDF
655 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
656 gross 967 const NcDim* ncdims[ldims];
657     NcVar *var, *ids;
658 jfenwick 1796 int rank = getRank();
659 gross 967 int type= getFunctionSpace().getTypeCode();
660     int ndims =0;
661     long dims[ldims];
662 gross 1141 const double* d_ptr=&(m_data[0]);
663 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
664 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
665     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
666 ksteube 1800 #ifdef PASO_MPI
667     MPI_Status status;
668     #endif
669 gross 965
670 ksteube 1827 #ifdef PASO_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 gross 967 // netCDF error handler
676     NcError err(NcError::verbose_nonfatal);
677     // Create the file.
678 ksteube 1748 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
679     NcFile dataFile(newFileName, NcFile::Replace);
680 gross 967 // check if writing was successful
681     if (!dataFile.is_valid())
682     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
683 gross 1141 if (!dataFile.add_att("type_id",2) )
684 gross 967 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 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
694 gross 967 }
695     if ( rank >1 ) {
696     dims[1]=shape[1];
697     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
698 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
699 gross 967 }
700     if ( rank >2 ) {
701     dims[2]=shape[2];
702     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
703 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
704 gross 967 }
705     if ( rank >3 ) {
706     dims[3]=shape[3];
707     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
708 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
709 gross 967 }
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 jfenwick 1808 if (getFunctionSpace().getNumSamples()>0)
718     {
719    
720     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
721 gross 967 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
722 jfenwick 1808 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
723     if (! (ids->put(ids_p,dims[rank+1])) )
724 gross 967 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
725 jfenwick 1808 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
726 gross 967 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
727 jfenwick 1808 if (! (var->put(d_ptr,dims)) )
728 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
729 jfenwick 1808 }
730 ksteube 1800 #ifdef PASO_MPI
731 ksteube 1801 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
732 ksteube 1800 #endif
733 gross 1023 #else
734     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
735     #endif
736 gross 950 }
737 gross 1358
738 jfenwick 1796 void
739 gross 1358 DataExpanded::setTaggedValue(int tagKey,
740 jfenwick 1796 const DataTypes::ShapeType& pointshape,
741     const DataTypes::ValueType& value,
742     int dataOffset)
743 gross 1358 {
744 jfenwick 2271 CHECK_FOR_EX_WRITE
745 gross 1358 int numSamples = getNumSamples();
746     int numDataPointsPerSample = getNumDPPSample();
747     int sampleNo,dataPointNo, i;
748 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
749     double* p;
750     const double* in=&value[0+dataOffset];
751 gross 1358
752 jfenwick 1796 if (value.size() != n) {
753 gross 1358 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 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
762 gross 1358 }
763     }
764     }
765     }
766    
767 jfenwick 1796
768 gross 1487 void
769     DataExpanded::reorderByReferenceIDs(int *reference_ids)
770     {
771 jfenwick 2271 CHECK_FOR_EX_WRITE
772 gross 1487 int numSamples = getNumSamples();
773 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
774 gross 1487 int sampleNo, sampleNo2,i;
775     double* p,*p2;
776     register double rtmp;
777     FunctionSpace fs=getFunctionSpace();
778 gross 1358
779 gross 1487 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 gross 1513 for (i=0; i<n ;i++) {
789 gross 1487 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 phornby 1628 if (! matched) {
800 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
801     }
802     }
803     }
804     }
805    
806 jfenwick 1796 DataTypes::ValueType&
807 jfenwick 2271 DataExpanded::getVectorRW()
808 jfenwick 1796 {
809 jfenwick 2271 CHECK_FOR_EX_WRITE
810 jfenwick 1796 return m_data.getData();
811     }
812 gross 1487
813 jfenwick 1796 const DataTypes::ValueType&
814 jfenwick 2271 DataExpanded::getVectorRO() const
815 jfenwick 1796 {
816     return m_data.getData();
817     }
818    
819 jfenwick 2271
820 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