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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26