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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 4 months ago) by ksteube
File size: 26448 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26