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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1962 - (show annotations)
Tue Nov 4 07:17:31 2008 UTC (10 years, 9 months ago) by gross
File size: 31705 byte(s)
there is an offet missing when copy data to a data point.


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "Data.h"
16 #include "DataExpanded.h"
17 #include "DataException.h"
18 #include "DataConstant.h"
19 #include "DataTagged.h"
20 #ifdef USE_NETCDF
21 #include <netcdfcpp.h>
22 #endif
23 #ifdef PASO_MPI
24 #include <mpi.h>
25 #endif
26
27 #include <boost/python/extract.hpp>
28 #include "DataMaths.h"
29
30 using namespace std;
31 using namespace boost::python;
32 using namespace boost;
33 using namespace escript::DataTypes;
34
35 namespace escript {
36
37 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
38 const FunctionSpace& what)
39 : DataAbstract(what,DataTypes::shapeFromNumArray(value))
40 {
41 //
42 // initialise the data array for this object
43 initialise(what.getNumSamples(),what.getNumDPPSample());
44 //
45 // copy the given value to every data point
46 copy(value);
47 }
48
49 DataExpanded::DataExpanded(const DataExpanded& other)
50 : DataAbstract(other.getFunctionSpace(), other.getShape()),
51 m_data(other.m_data)
52 {
53 }
54
55 DataExpanded::DataExpanded(const DataConstant& other)
56 : DataAbstract(other.getFunctionSpace(), other.getShape())
57 {
58 //
59 // initialise the data array for this object
60 initialise(other.getNumSamples(),other.getNumDPPSample());
61 //
62 // DataConstant only has one value, copy this to every data point
63 copy(other);
64 }
65
66 DataExpanded::DataExpanded(const DataTagged& other)
67 : DataAbstract(other.getFunctionSpace(), other.getShape())
68 {
69 //
70 // initialise the data array for this object
71 initialise(other.getNumSamples(),other.getNumDPPSample());
72 //
73 // for each data point in this object, extract and copy the corresponding data
74 // value from the given DataTagged object
75 int i,j;
76 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
77 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
78 #pragma omp parallel for private(i,j) schedule(static)
79 for (i=0;i<numRows;i++) {
80 for (j=0;j<numCols;j++) {
81 try {
82 DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
83 other.getVector(),
84 other.getPointOffset(i,j));
85 }
86 catch (std::exception& e) {
87 cout << e.what() << endl;
88 }
89 }
90 }
91 }
92
93 DataExpanded::DataExpanded(const DataExpanded& other,
94 const DataTypes::RegionType& region)
95 : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
96 {
97 //
98 // get the shape of the slice
99 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
100 //
101 // initialise this Data object to the shape of the slice
102 initialise(other.getNumSamples(),other.getNumDPPSample());
103 //
104 // copy the data
105 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
106 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
107 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
108 int i,j;
109 #pragma omp parallel for private(i,j) schedule(static)
110 for (i=0;i<numRows;i++) {
111 for (j=0;j<numCols;j++) {
112 try {
113 // getPointDataView().copySlice(getPointOffset(i,j),
114 // other.getPointDataView(),
115 // other.getPointOffset(i,j),
116 // region_loop_range);
117 DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
118 other.getVector(),
119 other.getShape(),
120 other.getPointOffset(i,j),
121 region_loop_range);
122 }
123 catch (std::exception& e) {
124 cout << e.what() << endl;
125 }
126 }
127 }
128 }
129
130 // DataExpanded::DataExpanded(const DataArrayView& value,
131 // const FunctionSpace& what)
132 // : DataAbstract(what)
133 // {
134 // //
135 // // get the shape of the given data value
136 // DataTypes::ShapeType tempShape=value.getShape();
137 // //
138 // // initialise this Data object to the shape of the given data value
139 // initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
140 // //
141 // // copy the given value to every data point
142 // copy(value);
143 // }
144
145 DataExpanded::DataExpanded(const FunctionSpace& what,
146 const DataTypes::ShapeType &shape,
147 const DataTypes::ValueType &data)
148 : DataAbstract(what,shape)
149 {
150 EsysAssert(data.size()%getNoValues()==0,
151 "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
152
153 if (data.size()==getNoValues())
154 {
155 ValueType& vec=m_data.getData();
156 //
157 // create the view of the data
158 initialise(what.getNumSamples(),what.getNumDPPSample());
159 // now we copy this value to all elements
160 for (int i=0;i<getLength();)
161 {
162 for (int j=0;j<getNoValues();++j,++i)
163 {
164 vec[i]=data[j];
165 }
166 }
167 }
168 else
169 {
170 //
171 // copy the data in the correct format
172 m_data.getData()=data;
173 }
174
175
176 }
177
178 DataExpanded::~DataExpanded()
179 {
180 }
181
182 DataAbstract*
183 DataExpanded::deepCopy()
184 {
185 return new DataExpanded(*this);
186 }
187
188
189 DataAbstract*
190 DataExpanded::getSlice(const DataTypes::RegionType& region) const
191 {
192 return new DataExpanded(*this,region);
193 }
194
195 void
196 DataExpanded::setSlice(const DataAbstract* value,
197 const DataTypes::RegionType& region)
198 {
199 const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
200 if (tempDataExp==0) {
201 throw DataException("Programming error - casting to DataExpanded.");
202 }
203 //
204 // get shape of slice
205 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
206 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
207 //
208 // check shape
209 if (getRank()!=region.size()) {
210 throw DataException("Error - Invalid slice region.");
211 }
212 if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
213 throw DataException (DataTypes::createShapeErrorMessage(
214 "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
215 }
216 //
217 // copy the data from the slice into this object
218 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
219 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
220 int i, j;
221 ValueType& vec=getVector();
222 const ShapeType& mshape=getShape();
223 const ValueType& tVec=tempDataExp->getVector();
224 const ShapeType& tShape=tempDataExp->getShape();
225 #pragma omp parallel for private(i,j) schedule(static)
226 for (i=0;i<numRows;i++) {
227 for (j=0;j<numCols;j++) {
228 /* getPointDataView().copySliceFrom(getPointOffset(i,j),
229 tempDataExp->getPointDataView(),
230 tempDataExp->getPointOffset(i,j),
231 region_loop_range);*/
232 DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
233 tVec,
234 tShape,
235 tempDataExp->getPointOffset(i,j),
236 region_loop_range);
237
238 }
239 }
240 }
241
242 void
243 DataExpanded::copy(const DataConstant& value)
244 {
245 EsysAssert((checkShape(getShape(), value.getShape())),
246 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
247
248 //
249 // copy a single value to every data point in this object
250 int nRows=m_data.getNumRows();
251 int nCols=m_data.getNumCols();
252 int i,j;
253 #pragma omp parallel for private(i,j) schedule(static)
254 for (i=0;i<nRows;i++) {
255 for (j=0;j<nCols;j++) {
256 // NOTE: An exception may be thown from this call if
257 // DOASSERT is on which of course will play
258 // havoc with the omp threads. Run single threaded
259 // if using DOASSERT.
260 //getPointDataView().copy(getPointOffset(i,j),value);
261 DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
262 }
263 }
264 }
265
266
267 // void
268 // DataExpanded::copy(const DataArrayView& value)
269 // {
270 // //
271 // // copy a single value to every data point in this object
272 // int nRows=m_data.getNumRows();
273 // int nCols=m_data.getNumCols();
274 // int i,j;
275 // #pragma omp parallel for private(i,j) schedule(static)
276 // for (i=0;i<nRows;i++) {
277 // for (j=0;j<nCols;j++) {
278 // // NOTE: An exception may be thown from this call if
279 // // DOASSERT is on which of course will play
280 // // havoc with the omp threads. Run single threaded
281 // // if using DOASSERT.
282 // getPointDataView().copy(getPointOffset(i,j),value);
283 // }
284 // }
285 // }
286
287 void
288 DataExpanded::copy(const boost::python::numeric::array& value)
289 {
290
291 // extract the shape of the numarray
292 DataTypes::ShapeType tempShape;
293 for (int i=0; i < value.getrank(); i++) {
294 tempShape.push_back(extract<int>(value.getshape()[i]));
295 }
296
297 // get the space for the data vector
298 // int len = DataTypes::noValues(tempShape);
299 // DataVector temp_data(len, 0.0, len);
300 // DataArrayView temp_dataView(temp_data, tempShape);
301 // temp_dataView.copy(value);
302
303 //
304 // check the input shape matches this shape
305 if (!DataTypes::checkShape(getShape(),tempShape)) {
306 throw DataException(DataTypes::createShapeErrorMessage(
307 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
308 tempShape,getShape()));
309 }
310 //
311 // now copy over the data
312 //copy(temp_dataView);
313 getVector().copyFromNumArray(value);
314 }
315
316
317 void
318 DataExpanded::initialise(int noSamples,
319 int noDataPointsPerSample)
320 {
321 if (noSamples==0) //retain the default empty object
322 {
323 return;
324 }
325 //
326 // resize data array to the required size
327 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
328 }
329
330 string
331 DataExpanded::toString() const
332 {
333 stringstream temp;
334 FunctionSpace fs=getFunctionSpace();
335
336 int offset=0;
337 for (int i=0;i<m_data.getNumRows();i++) {
338 for (int j=0;j<m_data.getNumCols();j++) {
339 offset=getPointOffset(i,j);
340 stringstream suffix;
341 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
342 temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
343 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
344 temp << endl;
345 }
346 }
347 }
348 string result=temp.str();
349 if (result.empty())
350 {
351 return "(data contains no samples)\n";
352 }
353 return temp.str();
354 }
355
356 DataTypes::ValueType::size_type
357 DataExpanded::getPointOffset(int sampleNo,
358 int dataPointNo) const
359 {
360 return m_data.index(sampleNo,dataPointNo);
361 }
362
363 DataTypes::ValueType::size_type
364 DataExpanded::getLength() const
365 {
366 return m_data.size();
367 }
368
369
370
371 void
372 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
373 //
374 // Get the number of samples and data-points per sample.
375 int numSamples = getNumSamples();
376 int numDataPointsPerSample = getNumDPPSample();
377 int dataPointRank = getRank();
378 ShapeType dataPointShape = getShape();
379 if (numSamples*numDataPointsPerSample > 0) {
380 //TODO: global error handling
381 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
382 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
383 }
384 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
385 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
386 }
387 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
388 ValueType& vec=getVector();
389 if (dataPointRank==0) {
390 vec[offset]=value;
391 } else if (dataPointRank==1) {
392 for (int i=0; i<dataPointShape[0]; i++) {
393 vec[offset+i]=value;
394 }
395 } else if (dataPointRank==2) {
396 for (int i=0; i<dataPointShape[0]; i++) {
397 for (int j=0; j<dataPointShape[1]; j++) {
398 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
399 }
400 }
401 } else if (dataPointRank==3) {
402 for (int i=0; i<dataPointShape[0]; i++) {
403 for (int j=0; j<dataPointShape[1]; j++) {
404 for (int k=0; k<dataPointShape[2]; k++) {
405 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
406 }
407 }
408 }
409 } else if (dataPointRank==4) {
410 for (int i=0; i<dataPointShape[0]; i++) {
411 for (int j=0; j<dataPointShape[1]; j++) {
412 for (int k=0; k<dataPointShape[2]; k++) {
413 for (int l=0; l<dataPointShape[3]; l++) {
414 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
415 }
416 }
417 }
418 }
419 }
420 }
421 }
422 void
423 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
424 //
425 // Get the number of samples and data-points per sample.
426 int numSamples = getNumSamples();
427 int numDataPointsPerSample = getNumDPPSample();
428 int dataPointRank = getRank();
429 const ShapeType& shape = getShape();
430 //
431 // check rank:
432 if (value.getrank()!=dataPointRank)
433 throw DataException("Rank of numarray does not match Data object rank");
434 if (numSamples*numDataPointsPerSample > 0) {
435 //TODO: global error handling
436 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
437 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
438 }
439 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
440 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
441 }
442 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
443 ValueType& vec=getVector();
444 if (dataPointRank==0) {
445 vec[offset]=extract<double>(value[0]);
446 } else if (dataPointRank==1) {
447 for (int i=0; i<shape[0]; i++) {
448 vec[offset+i]=extract<double>(value[i]);
449 }
450 } else if (dataPointRank==2) {
451 for (int i=0; i<shape[0]; i++) {
452 for (int j=0; j<shape[1]; j++) {
453 vec[offset+getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
454 }
455 }
456 } else if (dataPointRank==3) {
457 for (int i=0; i<shape[0]; i++) {
458 for (int j=0; j<shape[1]; j++) {
459 for (int k=0; k<shape[2]; k++) {
460 vec[offset+getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
461 }
462 }
463 }
464 } else if (dataPointRank==4) {
465 for (int i=0; i<shape[0]; i++) {
466 for (int j=0; j<shape[1]; j++) {
467 for (int k=0; k<shape[2]; k++) {
468 for (int l=0; l<shape[3]; l++) {
469 vec[offset+getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
470 }
471 }
472 }
473 }
474 }
475 }
476 }
477 void
478 DataExpanded::copyAll(const boost::python::numeric::array& value) {
479 //
480 // Get the number of samples and data-points per sample.
481 int numSamples = getNumSamples();
482 int numDataPointsPerSample = getNumDPPSample();
483 int dataPointRank = getRank();
484 const ShapeType& dataPointShape = getShape();
485 //
486 // check rank:
487 if (value.getrank()!=dataPointRank+1)
488 throw DataException("Rank of numarray does not match Data object rank");
489 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
490 throw DataException("leading dimension of numarray is too small");
491 //
492 ValueType& vec=getVector();
493 int dataPoint = 0;
494 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
495 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
496 ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
497 if (dataPointRank==0) {
498 vec[offset]=extract<double>(value[dataPoint]);
499 } else if (dataPointRank==1) {
500 for (int i=0; i<dataPointShape[0]; i++) {
501 vec[offset+i]=extract<double>(value[dataPoint][i]);
502 }
503 } else if (dataPointRank==2) {
504 for (int i=0; i<dataPointShape[0]; i++) {
505 for (int j=0; j<dataPointShape[1]; j++) {
506 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
507 }
508 }
509 } else if (dataPointRank==3) {
510 for (int i=0; i<dataPointShape[0]; i++) {
511 for (int j=0; j<dataPointShape[1]; j++) {
512 for (int k=0; k<dataPointShape[2]; k++) {
513 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
514 }
515 }
516 }
517 } else if (dataPointRank==4) {
518 for (int i=0; i<dataPointShape[0]; i++) {
519 for (int j=0; j<dataPointShape[1]; j++) {
520 for (int k=0; k<dataPointShape[2]; k++) {
521 for (int l=0; l<dataPointShape[3]; l++) {
522 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
523 }
524 }
525 }
526 }
527 }
528 dataPoint++;
529 }
530 }
531 }
532 void
533 DataExpanded::symmetric(DataAbstract* ev)
534 {
535 int sampleNo,dataPointNo;
536 int numSamples = getNumSamples();
537 int numDataPointsPerSample = getNumDPPSample();
538 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
539 if (temp_ev==0) {
540 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
541 }
542 ValueType& vec=getVector();
543 const ShapeType& shape=getShape();
544 ValueType& evVec=temp_ev->getVector();
545 const ShapeType& evShape=temp_ev->getShape();
546 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
547 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
548 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
549 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
550 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
551 }
552 }
553 }
554 void
555 DataExpanded::nonsymmetric(DataAbstract* ev)
556 {
557 int sampleNo,dataPointNo;
558 int numSamples = getNumSamples();
559 int numDataPointsPerSample = getNumDPPSample();
560 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
561 if (temp_ev==0) {
562 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
563 }
564 ValueType& vec=getVector();
565 const ShapeType& shape=getShape();
566 ValueType& evVec=temp_ev->getVector();
567 const ShapeType& evShape=temp_ev->getShape();
568 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
569 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
570 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
571 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
572 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
573 }
574 }
575 }
576 void
577 DataExpanded::trace(DataAbstract* ev, int axis_offset)
578 {
579 int sampleNo,dataPointNo;
580 int numSamples = getNumSamples();
581 int numDataPointsPerSample = getNumDPPSample();
582 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
583 if (temp_ev==0) {
584 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
585 }
586 ValueType& vec=getVector();
587 const ShapeType& shape=getShape();
588 ValueType& evVec=temp_ev->getVector();
589 const ShapeType& evShape=temp_ev->getShape();
590 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
591 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
592 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
593 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
594 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
595 }
596 }
597 }
598
599 void
600 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
601 {
602 int sampleNo,dataPointNo;
603 int numSamples = getNumSamples();
604 int numDataPointsPerSample = getNumDPPSample();
605 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
606 if (temp_ev==0) {
607 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
608 }
609 ValueType& vec=getVector();
610 const ShapeType& shape=getShape();
611 ValueType& evVec=temp_ev->getVector();
612 const ShapeType& evShape=temp_ev->getShape();
613 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
614 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
615 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
616 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
617 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
618 }
619 }
620 }
621
622 void
623 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
624 {
625 int sampleNo,dataPointNo;
626 int numSamples = getNumSamples();
627 int numDataPointsPerSample = getNumDPPSample();
628 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
629 if (temp_ev==0) {
630 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
631 }
632 ValueType& vec=getVector();
633 const ShapeType& shape=getShape();
634 ValueType& evVec=temp_ev->getVector();
635 const ShapeType& evShape=temp_ev->getShape();
636 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
637 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
638 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
639 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
640 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
641 }
642 }
643 }
644 void
645 DataExpanded::eigenvalues(DataAbstract* ev)
646 {
647 int sampleNo,dataPointNo;
648 int numSamples = getNumSamples();
649 int numDataPointsPerSample = getNumDPPSample();
650 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
651 if (temp_ev==0) {
652 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
653 }
654 ValueType& vec=getVector();
655 const ShapeType& shape=getShape();
656 ValueType& evVec=temp_ev->getVector();
657 const ShapeType& evShape=temp_ev->getShape();
658 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
659 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
660 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
661 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
662 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
663 }
664 }
665 }
666 void
667 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
668 {
669 int numSamples = getNumSamples();
670 int numDataPointsPerSample = getNumDPPSample();
671 int sampleNo,dataPointNo;
672 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
673 if (temp_ev==0) {
674 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
675 }
676 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
677 if (temp_V==0) {
678 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
679 }
680 ValueType& vec=getVector();
681 const ShapeType& shape=getShape();
682 ValueType& evVec=temp_ev->getVector();
683 const ShapeType& evShape=temp_ev->getShape();
684 ValueType& VVec=temp_V->getVector();
685 const ShapeType& VShape=temp_V->getShape();
686 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
687 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
688 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
689 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
690 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
691 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
692 }
693 }
694 }
695
696 void
697 DataExpanded::setToZero(){
698 // TODO: Surely there is a more efficient way to do this????
699 // Why is there no memset here? Parallel issues?
700 int numSamples = getNumSamples();
701 int numDataPointsPerSample = getNumDPPSample();
702 DataTypes::ValueType::size_type n = getNoValues();
703 double* p;
704 int sampleNo,dataPointNo, i;
705 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
706 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
707 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
708 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
709 for (i=0; i<n ;++i) p[i]=0.;
710 }
711 }
712 }
713
714 /* Append MPI rank to file name if multiple MPI processes */
715 char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
716 /* Make plenty of room for the mpi_rank number and terminating '\0' */
717 char *newFileName = (char *)malloc(strlen(fileName)+20);
718 strncpy(newFileName, fileName, strlen(fileName)+1);
719 if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
720 return(newFileName);
721 }
722
723 void
724 DataExpanded::dump(const std::string fileName) const
725 {
726 #ifdef USE_NETCDF
727 const int ldims=2+DataTypes::maxRank;
728 const NcDim* ncdims[ldims];
729 NcVar *var, *ids;
730 int rank = getRank();
731 int type= getFunctionSpace().getTypeCode();
732 int ndims =0;
733 long dims[ldims];
734 const double* d_ptr=&(m_data[0]);
735 const DataTypes::ShapeType& shape = getShape();
736 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
737 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
738 #ifdef PASO_MPI
739 MPI_Status status;
740 #endif
741
742 #ifdef PASO_MPI
743 /* Serialize NetCDF I/O */
744 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
745 #endif
746
747 // netCDF error handler
748 NcError err(NcError::verbose_nonfatal);
749 // Create the file.
750 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
751 NcFile dataFile(newFileName, NcFile::Replace);
752 // check if writing was successful
753 if (!dataFile.is_valid())
754 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
755 if (!dataFile.add_att("type_id",2) )
756 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
757 if (!dataFile.add_att("rank",rank) )
758 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
759 if (!dataFile.add_att("function_space_type",type))
760 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
761 ndims=rank+2;
762 if ( rank >0 ) {
763 dims[0]=shape[0];
764 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
765 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
766 }
767 if ( rank >1 ) {
768 dims[1]=shape[1];
769 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
770 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
771 }
772 if ( rank >2 ) {
773 dims[2]=shape[2];
774 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
775 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
776 }
777 if ( rank >3 ) {
778 dims[3]=shape[3];
779 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
780 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
781 }
782 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
783 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
784 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
785 dims[rank+1]=getFunctionSpace().getNumSamples();
786 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
787 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
788
789 if (getFunctionSpace().getNumSamples()>0)
790 {
791
792 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
793 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
794 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
795 if (! (ids->put(ids_p,dims[rank+1])) )
796 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
797 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
798 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
799 if (! (var->put(d_ptr,dims)) )
800 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
801 }
802 #ifdef PASO_MPI
803 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
804 #endif
805 #else
806 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
807 #endif
808 }
809
810 void
811 DataExpanded::setTaggedValue(int tagKey,
812 const DataTypes::ShapeType& pointshape,
813 const DataTypes::ValueType& value,
814 int dataOffset)
815 {
816 int numSamples = getNumSamples();
817 int numDataPointsPerSample = getNumDPPSample();
818 int sampleNo,dataPointNo, i;
819 DataTypes::ValueType::size_type n = getNoValues();
820 double* p;
821 const double* in=&value[0+dataOffset];
822
823 if (value.size() != n) {
824 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
825 }
826
827 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
828 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
829 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
830 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
831 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
832 for (i=0; i<n ;++i) p[i]=in[i];
833 }
834 }
835 }
836 }
837
838
839 void
840 DataExpanded::reorderByReferenceIDs(int *reference_ids)
841 {
842 int numSamples = getNumSamples();
843 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
844 int sampleNo, sampleNo2,i;
845 double* p,*p2;
846 register double rtmp;
847 FunctionSpace fs=getFunctionSpace();
848
849 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
850 const int id_in=reference_ids[sampleNo];
851 const int id=fs.getReferenceIDOfSample(sampleNo);
852 if (id!=id_in) {
853 bool matched=false;
854 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
855 if (id == reference_ids[sampleNo2]) {
856 p=&(m_data[getPointOffset(sampleNo,0)]);
857 p2=&(m_data[getPointOffset(sampleNo2,0)]);
858 for (i=0; i<n ;i++) {
859 rtmp=p[i];
860 p[i]=p2[i];
861 p2[i]=rtmp;
862 }
863 reference_ids[sampleNo]=id;
864 reference_ids[sampleNo2]=id_in;
865 matched=true;
866 break;
867 }
868 }
869 if (! matched) {
870 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
871 }
872 }
873 }
874 }
875
876 DataTypes::ValueType&
877 DataExpanded::getVector()
878 {
879 return m_data.getData();
880 }
881
882 const DataTypes::ValueType&
883 DataExpanded::getVector() const
884 {
885 return m_data.getData();
886 }
887
888 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26