/[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 1773 - (show annotations)
Tue Sep 9 02:52:26 2008 UTC (13 years, 1 month ago) by jfenwick
File size: 33678 byte(s)
Branch commit

Fixed some bugs.
Code now passes all_tests.
Next step is to test under OMP and MPI


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26