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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1388 - (show annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 11 months ago) by trankine
File size: 27725 byte(s)
And get the *(&(*&(* name right
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(getPointOffset(i,j),value);
221 }
222 }
223 }
224
225 void
226 DataExpanded::copy(const boost::python::numeric::array& value)
227 {
228
229 // extract the shape of the numarray
230 DataArrayView::ShapeType tempShape;
231 for (int i=0; i < value.getrank(); i++) {
232 tempShape.push_back(extract<int>(value.getshape()[i]));
233 }
234
235 // get the space for the data vector
236 int len = DataArrayView::noValues(tempShape);
237 DataVector temp_data(len, 0.0, len);
238 DataArrayView temp_dataView(temp_data, tempShape);
239 temp_dataView.copy(value);
240
241 //
242 // check the input shape matches this shape
243 if (!getPointDataView().checkShape(temp_dataView.getShape())) {
244 throw DataException(getPointDataView().createShapeErrorMessage(
245 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
246 temp_dataView.getShape()));
247 }
248 //
249 // now copy over the data
250 copy(temp_dataView);
251 }
252
253 void
254 DataExpanded::initialise(const DataArrayView::ShapeType& shape,
255 int noSamples,
256 int noDataPointsPerSample)
257 {
258 //
259 // resize data array to the required size
260 m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
261 //
262 // create the data view of the data array
263 DataArrayView temp(m_data.getData(),shape);
264 setPointDataView(temp);
265 }
266
267 string
268 DataExpanded::toString() const
269 {
270 stringstream temp;
271 FunctionSpace fs=getFunctionSpace();
272 //
273 // create a temporary view as the offset will be changed
274 DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
275 for (int i=0;i<m_data.getNumRows();i++) {
276 for (int j=0;j<m_data.getNumCols();j++) {
277 tempView.setOffset(getPointOffset(i,j));
278 stringstream suffix;
279 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
280 temp << tempView.toString(suffix.str());
281 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
282 temp << endl;
283 }
284 }
285 }
286 return temp.str();
287 }
288
289 DataArrayView::ValueType::size_type
290 DataExpanded::getPointOffset(int sampleNo,
291 int dataPointNo) const
292 {
293 return m_data.index(sampleNo,dataPointNo);
294 }
295
296 DataArrayView
297 DataExpanded::getDataPoint(int sampleNo,
298 int dataPointNo)
299 {
300 DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));
301 return temp;
302 }
303
304 DataArrayView::ValueType::size_type
305 DataExpanded::getLength() const
306 {
307 return m_data.size();
308 }
309
310 int
311 DataExpanded::archiveData(ofstream& archiveFile,
312 const DataArrayView::ValueType::size_type noValues) const
313 {
314 return(m_data.archiveData(archiveFile, noValues));
315 }
316
317 int
318 DataExpanded::extractData(ifstream& archiveFile,
319 const DataArrayView::ValueType::size_type noValues)
320 {
321 return(m_data.extractData(archiveFile, noValues));
322 }
323
324 void
325 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
326 //
327 // Get the number of samples and data-points per sample.
328 int numSamples = getNumSamples();
329 int numDataPointsPerSample = getNumDPPSample();
330 int dataPointRank = getPointDataView().getRank();
331 ShapeType dataPointShape = getPointDataView().getShape();
332 if (numSamples*numDataPointsPerSample > 0) {
333 //TODO: global error handling
334 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
335 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
336 }
337 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
338 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
339 }
340 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
341 if (dataPointRank==0) {
342 dataPointView()=value;
343 } else if (dataPointRank==1) {
344 for (int i=0; i<dataPointShape[0]; i++) {
345 dataPointView(i)=value;
346 }
347 } else if (dataPointRank==2) {
348 for (int i=0; i<dataPointShape[0]; i++) {
349 for (int j=0; j<dataPointShape[1]; j++) {
350 dataPointView(i,j)=value;
351 }
352 }
353 } else if (dataPointRank==3) {
354 for (int i=0; i<dataPointShape[0]; i++) {
355 for (int j=0; j<dataPointShape[1]; j++) {
356 for (int k=0; k<dataPointShape[2]; k++) {
357 dataPointView(i,j,k)=value;
358 }
359 }
360 }
361 } else if (dataPointRank==4) {
362 for (int i=0; i<dataPointShape[0]; i++) {
363 for (int j=0; j<dataPointShape[1]; j++) {
364 for (int k=0; k<dataPointShape[2]; k++) {
365 for (int l=0; l<dataPointShape[3]; l++) {
366 dataPointView(i,j,k,l)=value;
367 }
368 }
369 }
370 }
371 }
372 }
373 }
374 void
375 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
376 //
377 // Get the number of samples and data-points per sample.
378 int numSamples = getNumSamples();
379 int numDataPointsPerSample = getNumDPPSample();
380 int dataPointRank = getPointDataView().getRank();
381 ShapeType dataPointShape = getPointDataView().getShape();
382 //
383 // check rank:
384 if (value.getrank()!=dataPointRank)
385 throw DataException("Rank of numarray does not match Data object rank");
386 if (numSamples*numDataPointsPerSample > 0) {
387 //TODO: global error handling
388 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
389 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
390 }
391 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
392 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
393 }
394 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
395 if (dataPointRank==0) {
396 dataPointView()=extract<double>(value[0]);
397 } else if (dataPointRank==1) {
398 for (int i=0; i<dataPointShape[0]; i++) {
399 dataPointView(i)=extract<double>(value[i]);
400 }
401 } else if (dataPointRank==2) {
402 for (int i=0; i<dataPointShape[0]; i++) {
403 for (int j=0; j<dataPointShape[1]; j++) {
404 dataPointView(i,j)=extract<double>(value[i][j]);
405 }
406 }
407 } else if (dataPointRank==3) {
408 for (int i=0; i<dataPointShape[0]; i++) {
409 for (int j=0; j<dataPointShape[1]; j++) {
410 for (int k=0; k<dataPointShape[2]; k++) {
411 dataPointView(i,j,k)=extract<double>(value[i][j][k]);
412 }
413 }
414 }
415 } else if (dataPointRank==4) {
416 for (int i=0; i<dataPointShape[0]; i++) {
417 for (int j=0; j<dataPointShape[1]; j++) {
418 for (int k=0; k<dataPointShape[2]; k++) {
419 for (int l=0; l<dataPointShape[3]; l++) {
420 dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]);
421 }
422 }
423 }
424 }
425 }
426 }
427 }
428 void
429 DataExpanded::copyAll(const boost::python::numeric::array& value) {
430 //
431 // Get the number of samples and data-points per sample.
432 int numSamples = getNumSamples();
433 int numDataPointsPerSample = getNumDPPSample();
434 int dataPointRank = getPointDataView().getRank();
435 ShapeType dataPointShape = getPointDataView().getShape();
436 //
437 // check rank:
438 if (value.getrank()!=dataPointRank+1)
439 throw DataException("Rank of numarray does not match Data object rank");
440 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
441 throw DataException("leading dimension of numarray is too small");
442 //
443 int dataPoint = 0;
444 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
445 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
446 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
447 if (dataPointRank==0) {
448 dataPointView()=extract<double>(value[dataPoint]);
449 } else if (dataPointRank==1) {
450 for (int i=0; i<dataPointShape[0]; i++) {
451 dataPointView(i)=extract<double>(value[dataPoint][i]);
452 }
453 } else if (dataPointRank==2) {
454 for (int i=0; i<dataPointShape[0]; i++) {
455 for (int j=0; j<dataPointShape[1]; j++) {
456 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
457 }
458 }
459 } else if (dataPointRank==3) {
460 for (int i=0; i<dataPointShape[0]; i++) {
461 for (int j=0; j<dataPointShape[1]; j++) {
462 for (int k=0; k<dataPointShape[2]; k++) {
463 dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
464 }
465 }
466 }
467 } else if (dataPointRank==4) {
468 for (int i=0; i<dataPointShape[0]; i++) {
469 for (int j=0; j<dataPointShape[1]; j++) {
470 for (int k=0; k<dataPointShape[2]; k++) {
471 for (int l=0; l<dataPointShape[3]; l++) {
472 dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
473 }
474 }
475 }
476 }
477 }
478 dataPoint++;
479 }
480 }
481 }
482 void
483 DataExpanded::symmetric(DataAbstract* ev)
484 {
485 int sampleNo,dataPointNo;
486 int numSamples = getNumSamples();
487 int numDataPointsPerSample = getNumDPPSample();
488 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
489 if (temp_ev==0) {
490 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
491 }
492 DataArrayView& thisView=getPointDataView();
493 DataArrayView& evView=ev->getPointDataView();
494 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
495 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
496 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
497 DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
498 evView,ev->getPointOffset(sampleNo,dataPointNo));
499 }
500 }
501 }
502 void
503 DataExpanded::nonsymmetric(DataAbstract* ev)
504 {
505 int sampleNo,dataPointNo;
506 int numSamples = getNumSamples();
507 int numDataPointsPerSample = getNumDPPSample();
508 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
509 if (temp_ev==0) {
510 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
511 }
512 DataArrayView& thisView=getPointDataView();
513 DataArrayView& evView=ev->getPointDataView();
514 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
515 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
516 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
517 DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
518 evView,ev->getPointOffset(sampleNo,dataPointNo));
519 }
520 }
521 }
522 void
523 DataExpanded::trace(DataAbstract* ev, int axis_offset)
524 {
525 int sampleNo,dataPointNo;
526 int numSamples = getNumSamples();
527 int numDataPointsPerSample = getNumDPPSample();
528 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
529 if (temp_ev==0) {
530 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
531 }
532 DataArrayView& thisView=getPointDataView();
533 DataArrayView& evView=ev->getPointDataView();
534 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
535 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
536 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
537 DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
538 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
539 }
540 }
541 }
542
543 void
544 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
545 {
546 int sampleNo,dataPointNo;
547 int numSamples = getNumSamples();
548 int numDataPointsPerSample = getNumDPPSample();
549 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
550 if (temp_ev==0) {
551 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
552 }
553 DataArrayView& thisView=getPointDataView();
554 DataArrayView& evView=ev->getPointDataView();
555 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
556 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
557 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
558 DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
559 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
560 }
561 }
562 }
563
564 void
565 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
566 {
567 int sampleNo,dataPointNo;
568 int numSamples = getNumSamples();
569 int numDataPointsPerSample = getNumDPPSample();
570 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
571 if (temp_ev==0) {
572 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
573 }
574 DataArrayView& thisView=getPointDataView();
575 DataArrayView& evView=ev->getPointDataView();
576 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
577 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
578 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
579 DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
580 evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
581 }
582 }
583 }
584 void
585 DataExpanded::eigenvalues(DataAbstract* ev)
586 {
587 int sampleNo,dataPointNo;
588 int numSamples = getNumSamples();
589 int numDataPointsPerSample = getNumDPPSample();
590 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
591 if (temp_ev==0) {
592 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
593 }
594 DataArrayView& thisView=getPointDataView();
595 DataArrayView& evView=ev->getPointDataView();
596 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
597 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
598 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
599 DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
600 evView,ev->getPointOffset(sampleNo,dataPointNo));
601 }
602 }
603 }
604 void
605 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
606 {
607 int numSamples = getNumSamples();
608 int numDataPointsPerSample = getNumDPPSample();
609 int sampleNo,dataPointNo;
610 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
611 if (temp_ev==0) {
612 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
613 }
614 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
615 if (temp_V==0) {
616 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
617 }
618 DataArrayView& thisView=getPointDataView();
619 DataArrayView& evView=ev->getPointDataView();
620 DataArrayView& VView=V->getPointDataView();
621 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
622 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
623 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
624 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
625 evView,ev->getPointOffset(sampleNo,dataPointNo),
626 VView,V->getPointOffset(sampleNo,dataPointNo),
627 tol);
628 }
629 }
630 }
631
632 void
633 DataExpanded::setToZero(){
634 int numSamples = getNumSamples();
635 int numDataPointsPerSample = getNumDPPSample();
636 DataArrayView& thisView=getPointDataView();
637 DataArrayView::ValueType::size_type n = thisView.noValues();
638 double* p;
639 int sampleNo,dataPointNo, i;
640 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
641 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
642 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
643 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
644 for (int i=0; i<n ;++i) p[i]=0.;
645 }
646 }
647 }
648
649
650 void
651 DataExpanded::dump(const std::string fileName) const
652 {
653 #ifdef PASO_MPI
654 throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");
655 #endif
656 #ifdef USE_NETCDF
657 const int ldims=2+DataArrayView::maxRank;
658 const NcDim* ncdims[ldims];
659 NcVar *var, *ids;
660 int rank = getPointDataView().getRank();
661 int type= getFunctionSpace().getTypeCode();
662 int ndims =0;
663 long dims[ldims];
664 const double* d_ptr=&(m_data[0]);
665 DataArrayView::ShapeType shape = getPointDataView().getShape();
666
667 // netCDF error handler
668 NcError err(NcError::verbose_nonfatal);
669 // Create the file.
670 NcFile dataFile(fileName.c_str(), NcFile::Replace);
671 // check if writing was successful
672 if (!dataFile.is_valid())
673 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
674 if (!dataFile.add_att("type_id",2) )
675 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
676 if (!dataFile.add_att("rank",rank) )
677 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
678 if (!dataFile.add_att("function_space_type",type))
679 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
680 ndims=rank+2;
681 if ( rank >0 ) {
682 dims[0]=shape[0];
683 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
684 throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");
685 }
686 if ( rank >1 ) {
687 dims[1]=shape[1];
688 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
689 throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");
690 }
691 if ( rank >2 ) {
692 dims[2]=shape[2];
693 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
694 throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");
695 }
696 if ( rank >3 ) {
697 dims[3]=shape[3];
698 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
699 throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");
700 }
701 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
702 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
703 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
704 dims[rank+1]=getFunctionSpace().getNumSamples();
705 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
706 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
707
708 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
709 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
710 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
711 if (! (ids->put(ids_p,dims[rank+1])) )
712 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
713
714 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
715 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
716 if (! (var->put(d_ptr,dims)) )
717 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
718 #else
719 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
720 #endif
721 }
722
723 void
724 DataExpanded::setTaggedValue(int tagKey,
725 const DataArrayView& value)
726 {
727 int numSamples = getNumSamples();
728 int numDataPointsPerSample = getNumDPPSample();
729 int sampleNo,dataPointNo, i;
730 DataArrayView& thisView=getPointDataView();
731 DataArrayView::ValueType::size_type n = thisView.noValues();
732 double* p,*in=&(value.getData()[0]);
733
734 if (value.noValues() != n) {
735 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
736 }
737
738 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
739 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
740 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
741 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
742 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
743 for (int i=0; i<n ;++i) p[i]=in[i];
744 }
745 }
746 }
747 }
748
749
750 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26