/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataExpanded.cpp
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1781 - (show annotations)
Thu Sep 11 05:03:14 2008 UTC (10 years, 7 months ago) by jfenwick
File size: 34344 byte(s)
Branch commit

Merged changes from trunk version 1695 up to and including version 1779.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26