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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2742 - (hide annotations)
Thu Nov 12 06:03:37 2009 UTC (10 years ago) by jfenwick
File size: 28105 byte(s)
Merging changes from the lapack branch.

The inverse() operation has been moved into c++. [No lazy support for this operation yet.]
Optional Lapack support has been added for matrices larger than 3x3. 
service0 is set to use mkl_lapack.



1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 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 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 jfenwick 2458 throw DataException("Rank of value does not match Data object rank");
382 gross 921 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 jfenwick 2742
561 gross 950 void
562 jfenwick 2742 DataExpanded::matrixInverse(DataAbstract* out) const
563     {
564     DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
565     if (temp==0)
566     {
567     throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
568     }
569    
570     if (getRank()!=2)
571     {
572     throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
573    
574     }
575     int sampleNo;
576     const int numdpps=getNumDPPSample();
577     const int numSamples = getNumSamples();
578     const ValueType& vec=m_data.getData();
579     int errcode=0;
580     #pragma omp parallel private(sampleNo)
581     {
582     int errorcode=0;
583     LapackInverseHelper h(getShape()[0]);
584     #pragma omp for schedule(static)
585     for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
586     {
587     // not sure I like all those virtual calls to getPointOffset
588     DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
589     int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
590     if (res>errorcode)
591     {
592     errorcode=res;
593     #pragma omp critical
594     {
595     errcode=errorcode; // I'm not especially concerned which error gets reported as long as one is
596     }
597     }
598     }
599     }
600     if (errcode)
601     {
602     DataMaths::matrixInverseError(errcode); // throws exceptions
603     }
604     }
605    
606     void
607 gross 1118 DataExpanded::setToZero(){
608 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
609     // Why is there no memset here? Parallel issues?
610 jfenwick 2742 // A: This ensures that memory is touched by the correct thread.
611 jfenwick 2271 CHECK_FOR_EX_WRITE
612 gross 1118 int numSamples = getNumSamples();
613     int numDataPointsPerSample = getNumDPPSample();
614 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
615 gross 1118 double* p;
616     int sampleNo,dataPointNo, i;
617     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
618     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
619     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
620     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
621 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
622 gross 1118 }
623     }
624     }
625    
626 ksteube 1748 /* Append MPI rank to file name if multiple MPI processes */
627     char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
628     /* Make plenty of room for the mpi_rank number and terminating '\0' */
629     char *newFileName = (char *)malloc(strlen(fileName)+20);
630     strncpy(newFileName, fileName, strlen(fileName)+1);
631     if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
632     return(newFileName);
633     }
634 gross 1118
635     void
636 gross 950 DataExpanded::dump(const std::string fileName) const
637     {
638 gross 1023 #ifdef USE_NETCDF
639 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
640 gross 967 const NcDim* ncdims[ldims];
641     NcVar *var, *ids;
642 jfenwick 1796 int rank = getRank();
643 gross 967 int type= getFunctionSpace().getTypeCode();
644     int ndims =0;
645     long dims[ldims];
646 gross 1141 const double* d_ptr=&(m_data[0]);
647 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
648 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
649     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
650 ksteube 1800 #ifdef PASO_MPI
651     MPI_Status status;
652     #endif
653 gross 965
654 ksteube 1827 #ifdef PASO_MPI
655     /* Serialize NetCDF I/O */
656     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
657     #endif
658    
659 gross 967 // netCDF error handler
660     NcError err(NcError::verbose_nonfatal);
661     // Create the file.
662 ksteube 1748 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
663     NcFile dataFile(newFileName, NcFile::Replace);
664 gross 967 // check if writing was successful
665     if (!dataFile.is_valid())
666     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
667 gross 1141 if (!dataFile.add_att("type_id",2) )
668 gross 967 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
669     if (!dataFile.add_att("rank",rank) )
670     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
671     if (!dataFile.add_att("function_space_type",type))
672     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
673     ndims=rank+2;
674     if ( rank >0 ) {
675     dims[0]=shape[0];
676     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
677 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
678 gross 967 }
679     if ( rank >1 ) {
680     dims[1]=shape[1];
681     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
682 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
683 gross 967 }
684     if ( rank >2 ) {
685     dims[2]=shape[2];
686     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
687 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
688 gross 967 }
689     if ( rank >3 ) {
690     dims[3]=shape[3];
691     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
692 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
693 gross 967 }
694     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
695     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
696     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
697     dims[rank+1]=getFunctionSpace().getNumSamples();
698     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
699     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
700    
701 jfenwick 1808 if (getFunctionSpace().getNumSamples()>0)
702     {
703    
704     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
705 gross 967 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
706 jfenwick 1808 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
707     if (! (ids->put(ids_p,dims[rank+1])) )
708 gross 967 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
709 jfenwick 1808 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
710 gross 967 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
711 jfenwick 1808 if (! (var->put(d_ptr,dims)) )
712 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
713 jfenwick 1808 }
714 ksteube 1800 #ifdef PASO_MPI
715 ksteube 1801 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
716 ksteube 1800 #endif
717 gross 1023 #else
718     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
719     #endif
720 gross 950 }
721 gross 1358
722 jfenwick 1796 void
723 gross 1358 DataExpanded::setTaggedValue(int tagKey,
724 jfenwick 1796 const DataTypes::ShapeType& pointshape,
725     const DataTypes::ValueType& value,
726     int dataOffset)
727 gross 1358 {
728 jfenwick 2271 CHECK_FOR_EX_WRITE
729 gross 1358 int numSamples = getNumSamples();
730     int numDataPointsPerSample = getNumDPPSample();
731     int sampleNo,dataPointNo, i;
732 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
733     double* p;
734     const double* in=&value[0+dataOffset];
735 gross 1358
736 jfenwick 1796 if (value.size() != n) {
737 gross 1358 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
738     }
739    
740     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
741     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
742     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
743     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
744     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
745 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
746 gross 1358 }
747     }
748     }
749     }
750    
751 jfenwick 1796
752 gross 1487 void
753     DataExpanded::reorderByReferenceIDs(int *reference_ids)
754     {
755 jfenwick 2271 CHECK_FOR_EX_WRITE
756 gross 1487 int numSamples = getNumSamples();
757 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
758 gross 1487 int sampleNo, sampleNo2,i;
759     double* p,*p2;
760     register double rtmp;
761     FunctionSpace fs=getFunctionSpace();
762 gross 1358
763 gross 1487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
764     const int id_in=reference_ids[sampleNo];
765     const int id=fs.getReferenceIDOfSample(sampleNo);
766     if (id!=id_in) {
767     bool matched=false;
768     for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
769     if (id == reference_ids[sampleNo2]) {
770     p=&(m_data[getPointOffset(sampleNo,0)]);
771     p2=&(m_data[getPointOffset(sampleNo2,0)]);
772 gross 1513 for (i=0; i<n ;i++) {
773 gross 1487 rtmp=p[i];
774     p[i]=p2[i];
775     p2[i]=rtmp;
776     }
777     reference_ids[sampleNo]=id;
778     reference_ids[sampleNo2]=id_in;
779     matched=true;
780     break;
781     }
782     }
783 phornby 1628 if (! matched) {
784 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
785     }
786     }
787     }
788     }
789    
790 jfenwick 1796 DataTypes::ValueType&
791 jfenwick 2271 DataExpanded::getVectorRW()
792 jfenwick 1796 {
793 jfenwick 2271 CHECK_FOR_EX_WRITE
794 jfenwick 1796 return m_data.getData();
795     }
796 gross 1487
797 jfenwick 1796 const DataTypes::ValueType&
798 jfenwick 2271 DataExpanded::getVectorRO() const
799 jfenwick 1796 {
800     return m_data.getData();
801     }
802    
803 jfenwick 2271
804 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