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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 967 - (show annotations)
Tue Feb 13 09:40:12 2007 UTC (12 years, 7 months ago) by gross
File size: 25761 byte(s)
dump and load of expanded data via netCDF added. some test are still missing.
1 // $Id$
2 /*
3 ************************************************************
4 * Copyright 2006 by ACcESS MNRF *
5 * *
6 * http://www.access.edu.au *
7 * Primary Business: Queensland, Australia *
8 * Licensed under the Open Software License version 3.0 *
9 * http://www.opensource.org/licenses/osl-3.0.php *
10 * *
11 ************************************************************
12 */
13
14 #include "DataExpanded.h"
15 #include "DataException.h"
16 #include "DataConstant.h"
17 #include "DataTagged.h"
18 #include <netcdfcpp.h>
19
20 #include <boost/python/extract.hpp>
21
22 using namespace std;
23 using namespace boost::python;
24 using namespace boost;
25
26 namespace escript {
27
28 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
29 const FunctionSpace& what)
30 : DataAbstract(what)
31 {
32 DataArrayView::ShapeType tempShape;
33 //
34 // extract the shape of the python numarray
35 for (int i=0; i<value.getrank(); i++) {
36 tempShape.push_back(extract<int>(value.getshape()[i]));
37 }
38 //
39 // initialise the data array for this object
40 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
41 //
42 // copy the given value to every data point
43 copy(value);
44 }
45
46 DataExpanded::DataExpanded(const DataExpanded& other)
47 : DataAbstract(other.getFunctionSpace()),
48 m_data(other.m_data)
49 {
50 //
51 // create the view for the data
52 DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
53 setPointDataView(temp);
54 }
55
56 DataExpanded::DataExpanded(const DataConstant& other)
57 : DataAbstract(other.getFunctionSpace())
58 {
59 //
60 // initialise the data array for this object
61 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
62 //
63 // DataConstant only has one value, copy this to every data point
64 copy(other.getPointDataView());
65 }
66
67 DataExpanded::DataExpanded(const DataTagged& other)
68 : DataAbstract(other.getFunctionSpace())
69 {
70 //
71 // initialise the data array for this object
72 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
73 //
74 // for each data point in this object, extract and copy the corresponding data
75 // value from the given DataTagged object
76 int i,j;
77 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
78 DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
79 #pragma omp parallel for private(i,j) schedule(static)
80 for (i=0;i<numRows;i++) {
81 for (j=0;j<numCols;j++) {
82 try {
83 getPointDataView().copy(getPointOffset(i,j),
84 other.getPointDataView(),
85 other.getPointOffset(i,j));
86 }
87 catch (std::exception& e) {
88 cout << e.what() << endl;
89 }
90 }
91 }
92 }
93
94 DataExpanded::DataExpanded(const DataExpanded& other,
95 const DataArrayView::RegionType& region)
96 : DataAbstract(other.getFunctionSpace())
97 {
98 //
99 // get the shape of the slice
100 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
101 //
102 // initialise this Data object to the shape of the slice
103 initialise(shape,other.getNumSamples(),other.getNumDPPSample());
104 //
105 // copy the data
106 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
107 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
108 DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
109 int i,j;
110 #pragma omp parallel for private(i,j) schedule(static)
111 for (i=0;i<numRows;i++) {
112 for (j=0;j<numCols;j++) {
113 try {
114 getPointDataView().copySlice(getPointOffset(i,j),
115 other.getPointDataView(),
116 other.getPointOffset(i,j),
117 region_loop_range);
118 }
119 catch (std::exception& e) {
120 cout << e.what() << endl;
121 }
122 }
123 }
124 }
125
126 DataExpanded::DataExpanded(const DataArrayView& value,
127 const FunctionSpace& what)
128 : DataAbstract(what)
129 {
130 //
131 // get the shape of the given data value
132 DataArrayView::ShapeType tempShape=value.getShape();
133 //
134 // initialise this Data object to the shape of the given data value
135 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
136 //
137 // copy the given value to every data point
138 copy(value);
139 }
140
141 DataExpanded::DataExpanded(const FunctionSpace& what,
142 const DataArrayView::ShapeType &shape,
143 const DataArrayView::ValueType &data)
144 : DataAbstract(what)
145 {
146 //
147 // create the view of the data
148 initialise(shape,what.getNumSamples(),what.getNumDPPSample());
149 //
150 // copy the data in the correct format
151 m_data.getData()=data;
152 }
153
154 DataExpanded::~DataExpanded()
155 {
156 }
157
158 DataAbstract*
159 DataExpanded::getSlice(const DataArrayView::RegionType& region) const
160 {
161 return new DataExpanded(*this,region);
162 }
163
164 void
165 DataExpanded::setSlice(const DataAbstract* value,
166 const DataArrayView::RegionType& region)
167 {
168 const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
169 if (tempDataExp==0) {
170 throw DataException("Programming error - casting to DataExpanded.");
171 }
172 //
173 // get shape of slice
174 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
175 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
176 //
177 // check shape
178 if (getPointDataView().getRank()!=region.size()) {
179 throw DataException("Error - Invalid slice region.");
180 }
181 if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
182 throw DataException (value->getPointDataView().createShapeErrorMessage(
183 "Error - Couldn't copy slice due to shape mismatch.",shape));
184 }
185 //
186 // copy the data from the slice into this object
187 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
188 DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
189 int i, j;
190 #pragma omp parallel for private(i,j) schedule(static)
191 for (i=0;i<numRows;i++) {
192 for (j=0;j<numCols;j++) {
193 getPointDataView().copySliceFrom(getPointOffset(i,j),
194 tempDataExp->getPointDataView(),
195 tempDataExp->getPointOffset(i,j),
196 region_loop_range);
197 }
198 }
199 }
200
201 void
202 DataExpanded::copy(const DataArrayView& value)
203 {
204 //
205 // copy a single value to every data point in this object
206 int nRows=m_data.getNumRows();
207 int nCols=m_data.getNumCols();
208 int i,j;
209 #pragma omp parallel for private(i,j) schedule(static)
210 for (i=0;i<nRows;i++) {
211 for (j=0;j<nCols;j++) {
212 // NOTE: An exception may be thown from this call if
213 // DOASSERT is on which of course will play
214 // havoc with the omp threads. Run single threaded
215 // if using DOASSERT.
216 getPointDataView().copy(m_data.index(i,j),value);
217 }
218 }
219 }
220
221 void
222 DataExpanded::copy(const boost::python::numeric::array& value)
223 {
224 //
225 // first convert the numarray into a DataArray object
226 DataArray temp(value);
227 //
228 // check the input shape matches this shape
229 if (!getPointDataView().checkShape(temp.getView().getShape())) {
230 throw DataException(getPointDataView().createShapeErrorMessage(
231 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
232 temp.getView().getShape()));
233 }
234 //
235 // now copy over the data
236 copy(temp.getView());
237 }
238
239 void
240 DataExpanded::initialise(const DataArrayView::ShapeType& shape,
241 int noSamples,
242 int noDataPointsPerSample)
243 {
244 //
245 // resize data array to the required size
246 m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
247 //
248 // create the data view of the data array
249 DataArrayView temp(m_data.getData(),shape);
250 setPointDataView(temp);
251 }
252
253 string
254 DataExpanded::toString() const
255 {
256 stringstream temp;
257 FunctionSpace fs=getFunctionSpace();
258 //
259 // create a temporary view as the offset will be changed
260 DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
261 for (int i=0;i<m_data.getNumRows();i++) {
262 for (int j=0;j<m_data.getNumCols();j++) {
263 tempView.setOffset(m_data.index(i,j));
264 stringstream suffix;
265 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
266 temp << tempView.toString(suffix.str());
267 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
268 temp << endl;
269 }
270 }
271 }
272 return temp.str();
273 }
274
275 DataArrayView::ValueType::size_type
276 DataExpanded::getPointOffset(int sampleNo,
277 int dataPointNo) const
278 {
279 return m_data.index(sampleNo,dataPointNo);
280 }
281
282 DataArrayView
283 DataExpanded::getDataPoint(int sampleNo,
284 int dataPointNo)
285 {
286 DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));
287 return temp;
288 }
289
290 DataArrayView::ValueType::size_type
291 DataExpanded::getLength() const
292 {
293 return m_data.size();
294 }
295
296 int
297 DataExpanded::archiveData(ofstream& archiveFile,
298 const DataArrayView::ValueType::size_type noValues) const
299 {
300 return(m_data.archiveData(archiveFile, noValues));
301 }
302
303 int
304 DataExpanded::extractData(ifstream& archiveFile,
305 const DataArrayView::ValueType::size_type noValues)
306 {
307 return(m_data.extractData(archiveFile, noValues));
308 }
309
310 void
311 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
312 //
313 // Get the number of samples and data-points per sample.
314 int numSamples = getNumSamples();
315 int numDataPointsPerSample = getNumDPPSample();
316 int dataPointRank = getPointDataView().getRank();
317 ShapeType dataPointShape = getPointDataView().getShape();
318 if (numSamples*numDataPointsPerSample > 0) {
319 //TODO: global error handling
320 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
321 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
322 }
323 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
324 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
325 }
326 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
327 if (dataPointRank==0) {
328 dataPointView()=value;
329 } else if (dataPointRank==1) {
330 for (int i=0; i<dataPointShape[0]; i++) {
331 dataPointView(i)=value;
332 }
333 } else if (dataPointRank==2) {
334 for (int i=0; i<dataPointShape[0]; i++) {
335 for (int j=0; j<dataPointShape[1]; j++) {
336 dataPointView(i,j)=value;
337 }
338 }
339 } else if (dataPointRank==3) {
340 for (int i=0; i<dataPointShape[0]; i++) {
341 for (int j=0; j<dataPointShape[1]; j++) {
342 for (int k=0; k<dataPointShape[2]; k++) {
343 dataPointView(i,j,k)=value;
344 }
345 }
346 }
347 } else if (dataPointRank==4) {
348 for (int i=0; i<dataPointShape[0]; i++) {
349 for (int j=0; j<dataPointShape[1]; j++) {
350 for (int k=0; k<dataPointShape[2]; k++) {
351 for (int l=0; l<dataPointShape[3]; l++) {
352 dataPointView(i,j,k,l)=value;
353 }
354 }
355 }
356 }
357 }
358 }
359 }
360 void
361 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
362 //
363 // Get the number of samples and data-points per sample.
364 int numSamples = getNumSamples();
365 int numDataPointsPerSample = getNumDPPSample();
366 int dataPointRank = getPointDataView().getRank();
367 ShapeType dataPointShape = getPointDataView().getShape();
368 //
369 // check rank:
370 if (value.getrank()!=dataPointRank)
371 throw DataException("Rank of numarray does not match Data object rank");
372 if (numSamples*numDataPointsPerSample > 0) {
373 //TODO: global error handling
374 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
375 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
376 }
377 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
378 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
379 }
380 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
381 if (dataPointRank==0) {
382 dataPointView()=extract<double>(value[0]);
383 } else if (dataPointRank==1) {
384 for (int i=0; i<dataPointShape[0]; i++) {
385 dataPointView(i)=extract<double>(value[i]);
386 }
387 } else if (dataPointRank==2) {
388 for (int i=0; i<dataPointShape[0]; i++) {
389 for (int j=0; j<dataPointShape[1]; j++) {
390 dataPointView(i,j)=extract<double>(value[i][j]);
391 }
392 }
393 } else if (dataPointRank==3) {
394 for (int i=0; i<dataPointShape[0]; i++) {
395 for (int j=0; j<dataPointShape[1]; j++) {
396 for (int k=0; k<dataPointShape[2]; k++) {
397 dataPointView(i,j,k)=extract<double>(value[i][j][k]);
398 }
399 }
400 }
401 } else if (dataPointRank==4) {
402 for (int i=0; i<dataPointShape[0]; i++) {
403 for (int j=0; j<dataPointShape[1]; j++) {
404 for (int k=0; k<dataPointShape[2]; k++) {
405 for (int l=0; l<dataPointShape[3]; l++) {
406 dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]);
407 }
408 }
409 }
410 }
411 }
412 }
413 }
414 void
415 DataExpanded::copyAll(const boost::python::numeric::array& value) {
416 //
417 // Get the number of samples and data-points per sample.
418 int numSamples = getNumSamples();
419 int numDataPointsPerSample = getNumDPPSample();
420 int dataPointRank = getPointDataView().getRank();
421 ShapeType dataPointShape = getPointDataView().getShape();
422 //
423 // check rank:
424 if (value.getrank()!=dataPointRank+1)
425 throw DataException("Rank of numarray does not match Data object rank");
426 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
427 throw DataException("leading dimension of numarray is too small");
428 //
429 int dataPoint = 0;
430 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
431 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
432 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
433 if (dataPointRank==0) {
434 dataPointView()=extract<double>(value[dataPoint]);
435 } else if (dataPointRank==1) {
436 for (int i=0; i<dataPointShape[0]; i++) {
437 dataPointView(i)=extract<double>(value[dataPoint][i]);
438 }
439 } else if (dataPointRank==2) {
440 for (int i=0; i<dataPointShape[0]; i++) {
441 for (int j=0; j<dataPointShape[1]; j++) {
442 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
443 }
444 }
445 } else if (dataPointRank==3) {
446 for (int i=0; i<dataPointShape[0]; i++) {
447 for (int j=0; j<dataPointShape[1]; j++) {
448 for (int k=0; k<dataPointShape[2]; k++) {
449 dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
450 }
451 }
452 }
453 } else if (dataPointRank==4) {
454 for (int i=0; i<dataPointShape[0]; i++) {
455 for (int j=0; j<dataPointShape[1]; j++) {
456 for (int k=0; k<dataPointShape[2]; k++) {
457 for (int l=0; l<dataPointShape[3]; l++) {
458 dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
459 }
460 }
461 }
462 }
463 }
464 dataPoint++;
465 }
466 }
467 }
468 void
469 DataExpanded::symmetric(DataAbstract* ev)
470 {
471 int sampleNo,dataPointNo;
472 int numSamples = getNumSamples();
473 int numDataPointsPerSample = getNumDPPSample();
474 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
475 if (temp_ev==0) {
476 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
477 }
478 DataArrayView& thisView=getPointDataView();
479 DataArrayView& evView=ev->getPointDataView();
480 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
481 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
482 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
483 DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
484 evView,ev->getPointOffset(sampleNo,dataPointNo));
485 }
486 }
487 }
488 void
489 DataExpanded::nonsymmetric(DataAbstract* ev)
490 {
491 int sampleNo,dataPointNo;
492 int numSamples = getNumSamples();
493 int numDataPointsPerSample = getNumDPPSample();
494 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
495 if (temp_ev==0) {
496 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
497 }
498 DataArrayView& thisView=getPointDataView();
499 DataArrayView& evView=ev->getPointDataView();
500 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
501 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503 DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
504 evView,ev->getPointOffset(sampleNo,dataPointNo));
505 }
506 }
507 }
508 void
509 DataExpanded::trace(DataAbstract* ev, int axis_offset)
510 {
511 int sampleNo,dataPointNo;
512 int numSamples = getNumSamples();
513 int numDataPointsPerSample = getNumDPPSample();
514 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
515 if (temp_ev==0) {
516 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
517 }
518 DataArrayView& thisView=getPointDataView();
519 DataArrayView& evView=ev->getPointDataView();
520 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
521 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
522 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
523 DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
524 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
525 }
526 }
527 }
528
529 void
530 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
531 {
532 int sampleNo,dataPointNo;
533 int numSamples = getNumSamples();
534 int numDataPointsPerSample = getNumDPPSample();
535 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
536 if (temp_ev==0) {
537 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
538 }
539 DataArrayView& thisView=getPointDataView();
540 DataArrayView& evView=ev->getPointDataView();
541 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
542 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
543 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
544 DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
545 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
546 }
547 }
548 }
549
550 void
551 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
552 {
553 int sampleNo,dataPointNo;
554 int numSamples = getNumSamples();
555 int numDataPointsPerSample = getNumDPPSample();
556 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
557 if (temp_ev==0) {
558 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
559 }
560 DataArrayView& thisView=getPointDataView();
561 DataArrayView& evView=ev->getPointDataView();
562 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
563 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
564 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
565 DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
566 evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
567 }
568 }
569 }
570 void
571 DataExpanded::eigenvalues(DataAbstract* ev)
572 {
573 int sampleNo,dataPointNo;
574 int numSamples = getNumSamples();
575 int numDataPointsPerSample = getNumDPPSample();
576 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
577 if (temp_ev==0) {
578 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
579 }
580 DataArrayView& thisView=getPointDataView();
581 DataArrayView& evView=ev->getPointDataView();
582 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
583 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
584 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
585 DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
586 evView,ev->getPointOffset(sampleNo,dataPointNo));
587 }
588 }
589 }
590 void
591 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
592 {
593 int numSamples = getNumSamples();
594 int numDataPointsPerSample = getNumDPPSample();
595 int sampleNo,dataPointNo;
596 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
597 if (temp_ev==0) {
598 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
599 }
600 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
601 if (temp_V==0) {
602 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
603 }
604 DataArrayView& thisView=getPointDataView();
605 DataArrayView& evView=ev->getPointDataView();
606 DataArrayView& VView=V->getPointDataView();
607 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
608 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
609 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
610 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
611 evView,ev->getPointOffset(sampleNo,dataPointNo),
612 VView,V->getPointOffset(sampleNo,dataPointNo),
613 tol);
614 }
615 }
616 }
617
618 void
619 DataExpanded::dump(const std::string fileName) const
620 {
621 #ifdef PASO_MPI
622 throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");
623 #endif
624 const int ldims=2+DataArrayView::maxRank;
625 const NcDim* ncdims[ldims];
626 NcVar *var, *ids;
627 int rank = getPointDataView().getRank();
628 int type= getFunctionSpace().getTypeCode();
629 int ndims =0;
630 long dims[ldims];
631 DataArrayView::ShapeType shape = getPointDataView().getShape();
632
633 // netCDF error handler
634 NcError err(NcError::verbose_nonfatal);
635 // Create the file.
636 NcFile dataFile(fileName.c_str(), NcFile::Replace);
637 // check if writing was successful
638 if (!dataFile.is_valid())
639 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
640 if (!dataFile.add_att("type","expanded") )
641 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
642 if (!dataFile.add_att("rank",rank) )
643 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
644 if (!dataFile.add_att("function_space_type",type))
645 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
646 ndims=rank+2;
647 if ( rank >0 ) {
648 dims[0]=shape[0];
649 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
650 throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");
651 }
652 if ( rank >1 ) {
653 dims[1]=shape[1];
654 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
655 throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");
656 }
657 if ( rank >2 ) {
658 dims[2]=shape[2];
659 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
660 throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");
661 }
662 if ( rank >3 ) {
663 dims[3]=shape[3];
664 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
665 throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");
666 }
667 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
668 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
669 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
670 dims[rank+1]=getFunctionSpace().getNumSamples();
671 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
672 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
673
674 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
675 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
676 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
677 if (! (ids->put(ids_p,dims[rank+1])) )
678 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
679
680 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
681 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
682 if (! (var->put(&m_data[0],dims)) )
683 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
684 }
685
686
687 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26