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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1748 - (show annotations)
Wed Sep 3 06:10:39 2008 UTC (10 years, 11 months ago) by ksteube
File size: 29630 byte(s)
MPI parallelism for Data().dump and load.  Use multiple NetCDF
files, one file per MPI process

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26