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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1796 - (show annotations)
Wed Sep 17 01:45:46 2008 UTC (11 years, 1 month ago) by jfenwick
File size: 33970 byte(s)
Merged noarrayview branch onto trunk.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26