/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataExpanded.cpp
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1724 - (show annotations)
Mon Aug 25 05:38:57 2008 UTC (11 years ago) by jfenwick
File size: 30222 byte(s)
Branch commit

Moved createShapeErrorMessage() into DataTypes.h
Modified functions in DataAlgorithm.h to use non-DataArrayView accesses.

Added getVector() to each of DataTagged, DataConstant, DataExpanded - This method returns 
the underlying DataVector by reference/constant reference. Note that this method does not 
exist in DataAbstract so (at the momement) in order to pull the data from something you 
need to know what you are looking at. (Lower level access is still possible via double* 
though).

DataTagged now has a getOffsetForTag method and a getDefaultOffset method.

DataMaths.h - A new file containing the reductionOps etc from DataArrayView (but without 
requiring DAV).
This file requires significant commenting improvements.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26