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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2008 - (show annotations)
Mon Nov 10 08:59:14 2008 UTC (10 years, 9 months ago) by phornby
File size: 31849 byte(s)
Signed/unsigned comparisons eliminated (just loop variables, easy fix).
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 : parent(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 : parent(other.getFunctionSpace(), other.getShape()),
51 m_data(other.m_data)
52 {
53 }
54
55 DataExpanded::DataExpanded(const DataConstant& other)
56 : parent(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 : parent(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 : parent(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 : parent(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 (unsigned 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::getPointOffset(int sampleNo,
365 int dataPointNo)
366 {
367 return m_data.index(sampleNo,dataPointNo);
368 }
369
370 DataTypes::ValueType::size_type
371 DataExpanded::getLength() const
372 {
373 return m_data.size();
374 }
375
376
377
378 void
379 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
380 //
381 // Get the number of samples and data-points per sample.
382 int numSamples = getNumSamples();
383 int numDataPointsPerSample = getNumDPPSample();
384 int dataPointRank = getRank();
385 ShapeType dataPointShape = getShape();
386 if (numSamples*numDataPointsPerSample > 0) {
387 //TODO: global error handling
388 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
389 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
390 }
391 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
392 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
393 }
394 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
395 ValueType& vec=getVector();
396 if (dataPointRank==0) {
397 vec[offset]=value;
398 } else if (dataPointRank==1) {
399 for (int i=0; i<dataPointShape[0]; i++) {
400 vec[offset+i]=value;
401 }
402 } else if (dataPointRank==2) {
403 for (int i=0; i<dataPointShape[0]; i++) {
404 for (int j=0; j<dataPointShape[1]; j++) {
405 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
406 }
407 }
408 } else if (dataPointRank==3) {
409 for (int i=0; i<dataPointShape[0]; i++) {
410 for (int j=0; j<dataPointShape[1]; j++) {
411 for (int k=0; k<dataPointShape[2]; k++) {
412 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
413 }
414 }
415 }
416 } else if (dataPointRank==4) {
417 for (int i=0; i<dataPointShape[0]; i++) {
418 for (int j=0; j<dataPointShape[1]; j++) {
419 for (int k=0; k<dataPointShape[2]; k++) {
420 for (int l=0; l<dataPointShape[3]; l++) {
421 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
422 }
423 }
424 }
425 }
426 }
427 }
428 }
429 void
430 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
431 //
432 // Get the number of samples and data-points per sample.
433 int numSamples = getNumSamples();
434 int numDataPointsPerSample = getNumDPPSample();
435 int dataPointRank = getRank();
436 const ShapeType& shape = getShape();
437 //
438 // check rank:
439 if (value.getrank()!=dataPointRank)
440 throw DataException("Rank of numarray does not match Data object rank");
441 if (numSamples*numDataPointsPerSample > 0) {
442 //TODO: global error handling
443 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
444 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
445 }
446 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
447 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
448 }
449 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
450 ValueType& vec=getVector();
451 if (dataPointRank==0) {
452 vec[offset]=extract<double>(value[0]);
453 } else if (dataPointRank==1) {
454 for (int i=0; i<shape[0]; i++) {
455 vec[offset+i]=extract<double>(value[i]);
456 }
457 } else if (dataPointRank==2) {
458 for (int i=0; i<shape[0]; i++) {
459 for (int j=0; j<shape[1]; j++) {
460 vec[offset+getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
461 }
462 }
463 } else if (dataPointRank==3) {
464 for (int i=0; i<shape[0]; i++) {
465 for (int j=0; j<shape[1]; j++) {
466 for (int k=0; k<shape[2]; k++) {
467 vec[offset+getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
468 }
469 }
470 }
471 } else if (dataPointRank==4) {
472 for (int i=0; i<shape[0]; i++) {
473 for (int j=0; j<shape[1]; j++) {
474 for (int k=0; k<shape[2]; k++) {
475 for (int l=0; l<shape[3]; l++) {
476 vec[offset+getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
477 }
478 }
479 }
480 }
481 }
482 }
483 }
484 void
485 DataExpanded::copyAll(const boost::python::numeric::array& value) {
486 //
487 // Get the number of samples and data-points per sample.
488 int numSamples = getNumSamples();
489 int numDataPointsPerSample = getNumDPPSample();
490 int dataPointRank = getRank();
491 const ShapeType& dataPointShape = getShape();
492 //
493 // check rank:
494 if (value.getrank()!=dataPointRank+1)
495 throw DataException("Rank of numarray does not match Data object rank");
496 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
497 throw DataException("leading dimension of numarray is too small");
498 //
499 ValueType& vec=getVector();
500 int dataPoint = 0;
501 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503 ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
504 if (dataPointRank==0) {
505 vec[offset]=extract<double>(value[dataPoint]);
506 } else if (dataPointRank==1) {
507 for (int i=0; i<dataPointShape[0]; i++) {
508 vec[offset+i]=extract<double>(value[dataPoint][i]);
509 }
510 } else if (dataPointRank==2) {
511 for (int i=0; i<dataPointShape[0]; i++) {
512 for (int j=0; j<dataPointShape[1]; j++) {
513 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
514 }
515 }
516 } else if (dataPointRank==3) {
517 for (int i=0; i<dataPointShape[0]; i++) {
518 for (int j=0; j<dataPointShape[1]; j++) {
519 for (int k=0; k<dataPointShape[2]; k++) {
520 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
521 }
522 }
523 }
524 } else if (dataPointRank==4) {
525 for (int i=0; i<dataPointShape[0]; i++) {
526 for (int j=0; j<dataPointShape[1]; j++) {
527 for (int k=0; k<dataPointShape[2]; k++) {
528 for (int l=0; l<dataPointShape[3]; l++) {
529 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
530 }
531 }
532 }
533 }
534 }
535 dataPoint++;
536 }
537 }
538 }
539 void
540 DataExpanded::symmetric(DataAbstract* ev)
541 {
542 int sampleNo,dataPointNo;
543 int numSamples = getNumSamples();
544 int numDataPointsPerSample = getNumDPPSample();
545 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
546 if (temp_ev==0) {
547 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
548 }
549 ValueType& vec=getVector();
550 const ShapeType& shape=getShape();
551 ValueType& evVec=temp_ev->getVector();
552 const ShapeType& evShape=temp_ev->getShape();
553 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
554 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
555 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
556 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
557 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
558 }
559 }
560 }
561 void
562 DataExpanded::nonsymmetric(DataAbstract* ev)
563 {
564 int sampleNo,dataPointNo;
565 int numSamples = getNumSamples();
566 int numDataPointsPerSample = getNumDPPSample();
567 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
568 if (temp_ev==0) {
569 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
570 }
571 ValueType& vec=getVector();
572 const ShapeType& shape=getShape();
573 ValueType& evVec=temp_ev->getVector();
574 const ShapeType& evShape=temp_ev->getShape();
575 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
576 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
577 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
578 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
579 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
580 }
581 }
582 }
583 void
584 DataExpanded::trace(DataAbstract* ev, int axis_offset)
585 {
586 int sampleNo,dataPointNo;
587 int numSamples = getNumSamples();
588 int numDataPointsPerSample = getNumDPPSample();
589 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
590 if (temp_ev==0) {
591 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
592 }
593 ValueType& vec=getVector();
594 const ShapeType& shape=getShape();
595 ValueType& evVec=temp_ev->getVector();
596 const ShapeType& evShape=temp_ev->getShape();
597 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
598 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
599 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
600 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
601 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
602 }
603 }
604 }
605
606 void
607 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
608 {
609 int sampleNo,dataPointNo;
610 int numSamples = getNumSamples();
611 int numDataPointsPerSample = getNumDPPSample();
612 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
613 if (temp_ev==0) {
614 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
615 }
616 ValueType& vec=getVector();
617 const ShapeType& shape=getShape();
618 ValueType& evVec=temp_ev->getVector();
619 const ShapeType& evShape=temp_ev->getShape();
620 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
621 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
622 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
623 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
624 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
625 }
626 }
627 }
628
629 void
630 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
631 {
632 int sampleNo,dataPointNo;
633 int numSamples = getNumSamples();
634 int numDataPointsPerSample = getNumDPPSample();
635 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
636 if (temp_ev==0) {
637 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
638 }
639 ValueType& vec=getVector();
640 const ShapeType& shape=getShape();
641 ValueType& evVec=temp_ev->getVector();
642 const ShapeType& evShape=temp_ev->getShape();
643 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
644 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
645 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
646 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
647 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
648 }
649 }
650 }
651 void
652 DataExpanded::eigenvalues(DataAbstract* ev)
653 {
654 int sampleNo,dataPointNo;
655 int numSamples = getNumSamples();
656 int numDataPointsPerSample = getNumDPPSample();
657 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
658 if (temp_ev==0) {
659 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
660 }
661 ValueType& vec=getVector();
662 const ShapeType& shape=getShape();
663 ValueType& evVec=temp_ev->getVector();
664 const ShapeType& evShape=temp_ev->getShape();
665 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
666 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
667 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
668 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
669 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
670 }
671 }
672 }
673 void
674 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
675 {
676 int numSamples = getNumSamples();
677 int numDataPointsPerSample = getNumDPPSample();
678 int sampleNo,dataPointNo;
679 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
680 if (temp_ev==0) {
681 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
682 }
683 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
684 if (temp_V==0) {
685 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
686 }
687 ValueType& vec=getVector();
688 const ShapeType& shape=getShape();
689 ValueType& evVec=temp_ev->getVector();
690 const ShapeType& evShape=temp_ev->getShape();
691 ValueType& VVec=temp_V->getVector();
692 const ShapeType& VShape=temp_V->getShape();
693 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
694 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
695 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
696 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
697 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
698 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
699 }
700 }
701 }
702
703 void
704 DataExpanded::setToZero(){
705 // TODO: Surely there is a more efficient way to do this????
706 // Why is there no memset here? Parallel issues?
707 int numSamples = getNumSamples();
708 int numDataPointsPerSample = getNumDPPSample();
709 DataTypes::ValueType::size_type n = getNoValues();
710 double* p;
711 int sampleNo,dataPointNo, i;
712 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
713 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
714 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
715 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
716 for (i=0; i<n ;++i) p[i]=0.;
717 }
718 }
719 }
720
721 /* Append MPI rank to file name if multiple MPI processes */
722 char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
723 /* Make plenty of room for the mpi_rank number and terminating '\0' */
724 char *newFileName = (char *)malloc(strlen(fileName)+20);
725 strncpy(newFileName, fileName, strlen(fileName)+1);
726 if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
727 return(newFileName);
728 }
729
730 void
731 DataExpanded::dump(const std::string fileName) const
732 {
733 #ifdef USE_NETCDF
734 const int ldims=2+DataTypes::maxRank;
735 const NcDim* ncdims[ldims];
736 NcVar *var, *ids;
737 int rank = getRank();
738 int type= getFunctionSpace().getTypeCode();
739 int ndims =0;
740 long dims[ldims];
741 const double* d_ptr=&(m_data[0]);
742 const DataTypes::ShapeType& shape = getShape();
743 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
744 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
745 #ifdef PASO_MPI
746 MPI_Status status;
747 #endif
748
749 #ifdef PASO_MPI
750 /* Serialize NetCDF I/O */
751 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
752 #endif
753
754 // netCDF error handler
755 NcError err(NcError::verbose_nonfatal);
756 // Create the file.
757 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
758 NcFile dataFile(newFileName, NcFile::Replace);
759 // check if writing was successful
760 if (!dataFile.is_valid())
761 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
762 if (!dataFile.add_att("type_id",2) )
763 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
764 if (!dataFile.add_att("rank",rank) )
765 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
766 if (!dataFile.add_att("function_space_type",type))
767 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
768 ndims=rank+2;
769 if ( rank >0 ) {
770 dims[0]=shape[0];
771 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
772 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
773 }
774 if ( rank >1 ) {
775 dims[1]=shape[1];
776 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
777 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
778 }
779 if ( rank >2 ) {
780 dims[2]=shape[2];
781 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
782 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
783 }
784 if ( rank >3 ) {
785 dims[3]=shape[3];
786 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
787 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
788 }
789 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
790 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
791 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
792 dims[rank+1]=getFunctionSpace().getNumSamples();
793 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
794 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
795
796 if (getFunctionSpace().getNumSamples()>0)
797 {
798
799 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
800 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
801 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
802 if (! (ids->put(ids_p,dims[rank+1])) )
803 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
804 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
805 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
806 if (! (var->put(d_ptr,dims)) )
807 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
808 }
809 #ifdef PASO_MPI
810 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
811 #endif
812 #else
813 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
814 #endif
815 }
816
817 void
818 DataExpanded::setTaggedValue(int tagKey,
819 const DataTypes::ShapeType& pointshape,
820 const DataTypes::ValueType& value,
821 int dataOffset)
822 {
823 int numSamples = getNumSamples();
824 int numDataPointsPerSample = getNumDPPSample();
825 int sampleNo,dataPointNo, i;
826 DataTypes::ValueType::size_type n = getNoValues();
827 double* p;
828 const double* in=&value[0+dataOffset];
829
830 if (value.size() != n) {
831 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
832 }
833
834 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
835 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
836 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
837 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
838 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
839 for (i=0; i<n ;++i) p[i]=in[i];
840 }
841 }
842 }
843 }
844
845
846 void
847 DataExpanded::reorderByReferenceIDs(int *reference_ids)
848 {
849 int numSamples = getNumSamples();
850 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
851 int sampleNo, sampleNo2,i;
852 double* p,*p2;
853 register double rtmp;
854 FunctionSpace fs=getFunctionSpace();
855
856 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
857 const int id_in=reference_ids[sampleNo];
858 const int id=fs.getReferenceIDOfSample(sampleNo);
859 if (id!=id_in) {
860 bool matched=false;
861 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
862 if (id == reference_ids[sampleNo2]) {
863 p=&(m_data[getPointOffset(sampleNo,0)]);
864 p2=&(m_data[getPointOffset(sampleNo2,0)]);
865 for (i=0; i<n ;i++) {
866 rtmp=p[i];
867 p[i]=p2[i];
868 p2[i]=rtmp;
869 }
870 reference_ids[sampleNo]=id;
871 reference_ids[sampleNo2]=id_in;
872 matched=true;
873 break;
874 }
875 }
876 if (! matched) {
877 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
878 }
879 }
880 }
881 }
882
883 DataTypes::ValueType&
884 DataExpanded::getVector()
885 {
886 return m_data.getData();
887 }
888
889 const DataTypes::ValueType&
890 DataExpanded::getVector() const
891 {
892 return m_data.getData();
893 }
894
895 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26