/[escript]/branches/subworld2/escriptcore/src/DataExpanded.cpp
ViewVC logotype

Contents of /branches/subworld2/escriptcore/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5464 - (show annotations)
Sat Feb 14 00:25:03 2015 UTC (4 years, 3 months ago) by jfenwick
Original Path: trunk/escriptcore/src/DataExpanded.cpp
File size: 29556 byte(s)
Fix for python compile warnings
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2015 by University of Queensland
5 * http://www.uq.edu.au
6 *
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 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #define ESNEEDPYTHON
18 #include "esysUtils/first.h"
19
20 #include "Data.h"
21 #include "DataExpanded.h"
22 #include "DataException.h"
23 #include "DataConstant.h"
24 #include "DataTagged.h"
25 #include <limits>
26
27 #include "esysUtils/Esys_MPI.h"
28
29 #ifdef USE_NETCDF
30 #include <netcdfcpp.h>
31 #endif
32
33 #include <boost/python/extract.hpp>
34 #include "DataMaths.h"
35
36 #include "esysUtils/EsysRandom.h"
37
38 using namespace std;
39 using namespace boost::python;
40 using namespace boost;
41 using namespace escript::DataTypes;
42
43
44 // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
45
46 #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; abort();throw DataException(ss.str());}
47
48 namespace escript {
49
50 DataExpanded::DataExpanded(const WrappedArray& value,
51 const FunctionSpace& what)
52 : parent(what,value.getShape())
53 {
54 //
55 // initialise the data array for this object
56 initialise(what.getNumSamples(),what.getNumDPPSample());
57 //
58 // copy the given value to every data point
59 copy(value);
60 }
61
62 DataExpanded::DataExpanded(const DataExpanded& other)
63 : parent(other.getFunctionSpace(), other.getShape()),
64 m_data(other.m_data)
65 {
66 }
67
68 DataExpanded::DataExpanded(const DataConstant& other)
69 : parent(other.getFunctionSpace(), other.getShape())
70 {
71 //
72 // initialise the data array for this object
73 initialise(other.getNumSamples(),other.getNumDPPSample());
74 //
75 // DataConstant only has one value, copy this to every data point
76 copy(other);
77 }
78
79 DataExpanded::DataExpanded(const DataTagged& other)
80 : parent(other.getFunctionSpace(), other.getShape())
81 {
82 //
83 // initialise the data array for this object
84 initialise(other.getNumSamples(),other.getNumDPPSample());
85 //
86 // for each data point in this object, extract and copy the corresponding data
87 // value from the given DataTagged object
88 int i,j;
89 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
90 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
91 #pragma omp parallel for private(i,j) schedule(static)
92 for (i=0;i<numRows;i++) {
93 for (j=0;j<numCols;j++) {
94 try {
95 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
96 other.getVectorRO(),
97 other.getPointOffset(i,j));
98 }
99 catch (std::exception& e) {
100 cout << e.what() << endl;
101 }
102 }
103 }
104 }
105
106 DataExpanded::DataExpanded(const DataExpanded& other,
107 const DataTypes::RegionType& region)
108 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
109 {
110 //
111 // get the shape of the slice
112 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
113 //
114 // initialise this Data object to the shape of the slice
115 initialise(other.getNumSamples(),other.getNumDPPSample());
116 //
117 // copy the data
118 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
119 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
120 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
121 int i,j;
122 #pragma omp parallel for private(i,j) schedule(static)
123 for (i=0;i<numRows;i++) {
124 for (j=0;j<numCols;j++) {
125 try {
126 DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
127 other.getVectorRO(),
128 other.getShape(),
129 other.getPointOffset(i,j),
130 region_loop_range);
131 }
132 catch (std::exception& e) {
133 cout << e.what() << endl;
134 }
135 }
136 }
137 }
138
139 DataExpanded::DataExpanded(const FunctionSpace& what,
140 const DataTypes::ShapeType &shape,
141 const DataTypes::ValueType &data)
142 : parent(what,shape)
143 {
144 EsysAssert(data.size()%getNoValues()==0,
145 "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
146
147 if (data.size()==getNoValues())
148 {
149 ValueType& vec=m_data.getData();
150 //
151 // create the view of the data
152 initialise(what.getNumSamples(),what.getNumDPPSample());
153 // now we copy this value to all elements
154 for (int i=0;i<getLength();)
155 {
156 for (unsigned int j=0;j<getNoValues();++j,++i)
157 {
158 vec[i]=data[j];
159 }
160 }
161 }
162 else
163 {
164 //
165 // copy the data in the correct format
166 m_data.getData()=data;
167 }
168
169
170 }
171
172 DataExpanded::DataExpanded(const FunctionSpace& what,
173 const DataTypes::ShapeType &shape,
174 const double v)
175 : parent(what,shape)
176 {
177 ValueType& vec=m_data.getData();
178 //
179 // create the view of the data
180 initialise(what.getNumSamples(),what.getNumDPPSample());
181 // now we copy this value to all elements
182 const int L=getLength();
183 int i;
184 #pragma omp parallel for schedule(static) private(i)
185 for (i=0;i<L;++i)
186 {
187 vec[i]=v;
188 }
189 }
190
191
192
193 DataExpanded::~DataExpanded()
194 {
195 }
196
197 DataAbstract*
198 DataExpanded::deepCopy()
199 {
200 return new DataExpanded(*this);
201 }
202
203
204 DataAbstract*
205 DataExpanded::getSlice(const DataTypes::RegionType& region) const
206 {
207 return new DataExpanded(*this,region);
208 }
209
210 void
211 DataExpanded::setSlice(const DataAbstract* value,
212 const DataTypes::RegionType& region)
213 {
214 const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
215 if (tempDataExp==0) {
216 throw DataException("Programming error - casting to DataExpanded.");
217 }
218 CHECK_FOR_EX_WRITE
219 //
220 // get shape of slice
221 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
222 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
223 //
224 // check shape
225 if (getRank()!=region.size()) {
226 throw DataException("Error - Invalid slice region.");
227 }
228 if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
229 throw DataException (DataTypes::createShapeErrorMessage(
230 "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
231 }
232 //
233 // copy the data from the slice into this object
234 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
235 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
236 int i, j;
237 ValueType& vec=getVectorRW();
238 const ShapeType& mshape=getShape();
239 const ValueType& tVec=tempDataExp->getVectorRO();
240 const ShapeType& tShape=tempDataExp->getShape();
241 #pragma omp parallel for private(i,j) schedule(static)
242 for (i=0;i<numRows;i++) {
243 for (j=0;j<numCols;j++) {
244 DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
245 tVec,
246 tShape,
247 tempDataExp->getPointOffset(i,j),
248 region_loop_range);
249
250 }
251 }
252 }
253
254 void
255 DataExpanded::copy(const DataConstant& value)
256 {
257 EsysAssert((checkShape(getShape(), value.getShape())),
258 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
259
260 //
261 // copy a single value to every data point in this object
262 int nRows=m_data.getNumRows();
263 int nCols=m_data.getNumCols();
264 int i,j;
265 #pragma omp parallel for private(i,j) schedule(static)
266 for (i=0;i<nRows;i++) {
267 for (j=0;j<nCols;j++) {
268 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(), value.getVectorRO(), 0);
269 }
270 }
271 }
272
273 void
274 DataExpanded::copy(const WrappedArray& value)
275 {
276 // check the input shape matches this shape
277 if (!DataTypes::checkShape(getShape(),value.getShape())) {
278 throw DataException(DataTypes::createShapeErrorMessage(
279 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
280 value.getShape(),getShape()));
281 }
282 getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
283 }
284
285
286 void
287 DataExpanded::initialise(int noSamples,
288 int noDataPointsPerSample)
289 {
290 if (noSamples==0) //retain the default empty object
291 {
292 return;
293 }
294 //
295 // resize data array to the required size
296 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
297 }
298 bool
299 DataExpanded::hasNaN() const
300 {
301 bool haveNaN=false;
302 const ValueType& v=m_data.getData();
303 #pragma omp parallel for
304 for (ValueType::size_type i=0;i<v.size();++i)
305 {
306 if (nancheck(v[i]))
307 {
308 #pragma omp critical
309 {
310 haveNaN=true;
311 }
312 }
313 }
314 return haveNaN;
315 }
316
317 void
318 DataExpanded::replaceNaN(double value)
319 {
320 #pragma omp parallel for
321 for (ValueType::size_type i=0;i<m_data.size();++i)
322 {
323 if (nancheck(m_data[i]))
324 {
325 m_data[i] = value;
326 }
327 }
328 }
329
330
331 string
332 DataExpanded::toString() const
333 {
334 stringstream temp;
335 FunctionSpace fs=getFunctionSpace();
336
337 int offset=0;
338 for (int i=0;i<m_data.getNumRows();i++) {
339 for (int j=0;j<m_data.getNumCols();j++) {
340 offset=getPointOffset(i,j);
341 stringstream suffix;
342 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
343 temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
344 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
345 temp << endl;
346 }
347 }
348 }
349 string result=temp.str();
350 if (result.empty())
351 {
352 return "(data contains no samples)\n";
353 }
354 return temp.str();
355 }
356
357 DataTypes::ValueType::size_type
358 DataExpanded::getPointOffset(int sampleNo,
359 int dataPointNo) const
360 {
361 return m_data.index(sampleNo,dataPointNo);
362 }
363
364 DataTypes::ValueType::size_type
365 DataExpanded::getPointOffset(int sampleNo,
366 int dataPointNo)
367 {
368 return m_data.index(sampleNo,dataPointNo);
369 }
370
371 DataTypes::ValueType::size_type
372 DataExpanded::getLength() const
373 {
374 return m_data.size();
375 }
376
377
378
379 void
380 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
381 CHECK_FOR_EX_WRITE
382 //
383 // Get the number of samples and data-points per sample.
384 int numSamples = getNumSamples();
385 int numDataPointsPerSample = getNumDPPSample();
386 int dataPointRank = getRank();
387 ShapeType dataPointShape = getShape();
388 if (numSamples*numDataPointsPerSample > 0) {
389 //TODO: global error handling
390 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
391 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
392 }
393 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
394 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
395 }
396 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
397 ValueType& vec=getVectorRW();
398 if (dataPointRank==0) {
399 vec[offset]=value;
400 } else if (dataPointRank==1) {
401 for (int i=0; i<dataPointShape[0]; i++) {
402 vec[offset+i]=value;
403 }
404 } else if (dataPointRank==2) {
405 for (int i=0; i<dataPointShape[0]; i++) {
406 for (int j=0; j<dataPointShape[1]; j++) {
407 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
408 }
409 }
410 } else if (dataPointRank==3) {
411 for (int i=0; i<dataPointShape[0]; i++) {
412 for (int j=0; j<dataPointShape[1]; j++) {
413 for (int k=0; k<dataPointShape[2]; k++) {
414 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
415 }
416 }
417 }
418 } else if (dataPointRank==4) {
419 for (int i=0; i<dataPointShape[0]; i++) {
420 for (int j=0; j<dataPointShape[1]; j++) {
421 for (int k=0; k<dataPointShape[2]; k++) {
422 for (int l=0; l<dataPointShape[3]; l++) {
423 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
424 }
425 }
426 }
427 }
428 }
429 }
430 }
431
432 void
433 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
434 CHECK_FOR_EX_WRITE
435 //
436 // Get the number of samples and data-points per sample.
437 int numSamples = getNumSamples();
438 int numDataPointsPerSample = getNumDPPSample();
439 //
440 // check rank:
441 if (value.getRank()!=getRank())
442 throw DataException("Rank of value does not match Data object rank");
443 if (numSamples*numDataPointsPerSample > 0) {
444 //TODO: global error handling
445 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
446 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
447 }
448 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
449 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
450 }
451 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
452 ValueType& vec=getVectorRW();
453 vec.copyFromArrayToOffset(value,offset,1);
454 }
455 }
456
457 void
458 DataExpanded::symmetric(DataAbstract* ev)
459 {
460 int sampleNo,dataPointNo;
461 int numSamples = getNumSamples();
462 int numDataPointsPerSample = getNumDPPSample();
463 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
464 if (temp_ev==0) {
465 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (probably a programming error).");
466 }
467 const ValueType& vec=getVectorRO();
468 const ShapeType& shape=getShape();
469 ValueType& evVec=temp_ev->getVectorRW();
470 const ShapeType& evShape=temp_ev->getShape();
471 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
472 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
473 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
474 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
475 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
476 }
477 }
478 }
479 void
480 DataExpanded::nonsymmetric(DataAbstract* ev)
481 {
482 int sampleNo,dataPointNo;
483 int numSamples = getNumSamples();
484 int numDataPointsPerSample = getNumDPPSample();
485 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
486 if (temp_ev==0) {
487 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (probably a programming error).");
488 }
489 const ValueType& vec=getVectorRO();
490 const ShapeType& shape=getShape();
491 ValueType& evVec=temp_ev->getVectorRW();
492 const ShapeType& evShape=temp_ev->getShape();
493 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
494 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
495 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
496 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
497 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
498 }
499 }
500 }
501 void
502 DataExpanded::trace(DataAbstract* ev, int axis_offset)
503 {
504 int sampleNo,dataPointNo;
505 int numSamples = getNumSamples();
506 int numDataPointsPerSample = getNumDPPSample();
507 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
508 if (temp_ev==0) {
509 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (probably a programming error).");
510 }
511 const ValueType& vec=getVectorRO();
512 const ShapeType& shape=getShape();
513 ValueType& evVec=temp_ev->getVectorRW();
514 const ShapeType& evShape=temp_ev->getShape();
515 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
516 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
517 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
518 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
519 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
520 }
521 }
522 }
523
524 void
525 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
526 {
527 int sampleNo,dataPointNo;
528 int numSamples = getNumSamples();
529 int numDataPointsPerSample = getNumDPPSample();
530 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
531 if (temp_ev==0) {
532 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (probably a programming error).");
533 }
534 const ValueType& vec=getVectorRO();
535 const ShapeType& shape=getShape();
536 ValueType& evVec=temp_ev->getVectorRW();
537 const ShapeType& evShape=temp_ev->getShape();
538 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
539 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
540 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
541 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
542 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
543 }
544 }
545 }
546
547 void
548 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
549 {
550 int sampleNo,dataPointNo;
551 int numSamples = getNumSamples();
552 int numDataPointsPerSample = getNumDPPSample();
553 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
554 if (temp_ev==0) {
555 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (probably a programming error).");
556 }
557 const ValueType& vec=getVectorRO();
558 const ShapeType& shape=getShape();
559 ValueType& evVec=temp_ev->getVectorRW();
560 const ShapeType& evShape=temp_ev->getShape();
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 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
565 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
566 }
567 }
568 }
569 void
570 DataExpanded::eigenvalues(DataAbstract* ev)
571 {
572 int sampleNo,dataPointNo;
573 int numSamples = getNumSamples();
574 int numDataPointsPerSample = getNumDPPSample();
575 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
576 if (temp_ev==0) {
577 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (probably a programming error).");
578 }
579 const ValueType& vec=getVectorRO();
580 const ShapeType& shape=getShape();
581 ValueType& evVec=temp_ev->getVectorRW();
582 const ShapeType& evShape=temp_ev->getShape();
583 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
585 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
586 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
587 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
588 }
589 }
590 }
591 void
592 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
593 {
594 int numSamples = getNumSamples();
595 int numDataPointsPerSample = getNumDPPSample();
596 int sampleNo,dataPointNo;
597 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
598 if (temp_ev==0) {
599 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (probably a programming error).");
600 }
601 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
602 if (temp_V==0) {
603 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (probably a programming error).");
604 }
605 const ValueType& vec=getVectorRO();
606 const ShapeType& shape=getShape();
607 ValueType& evVec=temp_ev->getVectorRW();
608 const ShapeType& evShape=temp_ev->getShape();
609 ValueType& VVec=temp_V->getVectorRW();
610 const ShapeType& VShape=temp_V->getShape();
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 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
615 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
616 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
617 }
618 }
619 }
620
621
622 int
623 DataExpanded::matrixInverse(DataAbstract* out) const
624 {
625 DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
626 if (temp==0)
627 {
628 throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (probably a programming error).");
629 }
630
631 if (getRank()!=2)
632 {
633 throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
634
635 }
636 int sampleNo;
637 const int numdpps=getNumDPPSample();
638 const int numSamples = getNumSamples();
639 const ValueType& vec=m_data.getData();
640 int errcode=0;
641 #pragma omp parallel private(sampleNo)
642 {
643 int errorcode=0;
644 LapackInverseHelper h(getShape()[0]);
645 #pragma omp for schedule(static)
646 for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
647 {
648 // not sure I like all those virtual calls to getPointOffset
649 DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
650 int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
651 if (res>errorcode)
652 {
653 errorcode=res;
654 #pragma omp critical
655 {
656 errcode=errorcode; // I'm not especially concerned which error gets reported as long as one is
657 }
658 }
659 }
660 }
661 return errcode;
662 }
663
664 void
665 DataExpanded::setToZero(){
666 // TODO: Surely there is a more efficient way to do this????
667 // Why is there no memset here? Parallel issues?
668 // A: This ensures that memory is touched by the correct thread.
669 CHECK_FOR_EX_WRITE
670 int numSamples = getNumSamples();
671 int numDataPointsPerSample = getNumDPPSample();
672 DataTypes::ValueType::size_type n = getNoValues();
673 double* p;
674 int sampleNo,dataPointNo, i;
675 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
676 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
677 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
678 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
679 for (i=0; i<n ;++i) p[i]=0.;
680 }
681 }
682 }
683
684 void
685 DataExpanded::dump(const std::string fileName) const
686 {
687 #ifdef USE_NETCDF
688 const int ldims=2+DataTypes::maxRank;
689 const NcDim* ncdims[ldims];
690 NcVar *var, *ids;
691 int rank = getRank();
692 int type= getFunctionSpace().getTypeCode();
693 int ndims =0;
694 long dims[ldims];
695 const double* d_ptr=&(m_data[0]);
696 const DataTypes::ShapeType& shape = getShape();
697 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
698 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
699 #ifdef ESYS_MPI
700 MPI_Status status;
701 #endif
702
703 #ifdef ESYS_MPI
704 /* Serialize NetCDF I/O */
705 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
706 #endif
707
708 // netCDF error handler
709 NcError err(NcError::verbose_nonfatal);
710 // Create the file.
711 std::string newFileName(esysUtils::appendRankToFileName(fileName,
712 mpi_num, mpi_iam));
713 NcFile dataFile(newFileName.c_str(), NcFile::Replace);
714 // check if writing was successful
715 if (!dataFile.is_valid())
716 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
717 if (!dataFile.add_att("type_id",2) )
718 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
719 if (!dataFile.add_att("rank",rank) )
720 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
721 if (!dataFile.add_att("function_space_type",type))
722 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
723 ndims=rank+2;
724 if ( rank >0 ) {
725 dims[0]=shape[0];
726 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
727 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
728 }
729 if ( rank >1 ) {
730 dims[1]=shape[1];
731 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
732 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
733 }
734 if ( rank >2 ) {
735 dims[2]=shape[2];
736 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
737 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
738 }
739 if ( rank >3 ) {
740 dims[3]=shape[3];
741 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
742 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
743 }
744 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
745 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
746 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
747 dims[rank+1]=getFunctionSpace().getNumSamples();
748 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
749 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
750
751 if (getFunctionSpace().getNumSamples()>0)
752 {
753
754 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
755 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
756 const dim_t* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
757 if (! (ids->put(ids_p,dims[rank+1])) )
758 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
759 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
760 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
761 if (! (var->put(d_ptr,dims)) )
762 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
763 }
764 #ifdef ESYS_MPI
765 if (mpi_iam<mpi_num-1)
766 MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
767 #endif
768 #else
769 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
770 #endif
771 }
772
773 void
774 DataExpanded::setTaggedValue(int tagKey,
775 const DataTypes::ShapeType& pointshape,
776 const DataTypes::ValueType& value,
777 int dataOffset)
778 {
779 CHECK_FOR_EX_WRITE
780 int numSamples = getNumSamples();
781 int numDataPointsPerSample = getNumDPPSample();
782 int sampleNo,dataPointNo, i;
783 DataTypes::ValueType::size_type n = getNoValues();
784 double* p;
785 const double* in=&value[0+dataOffset];
786
787 if (value.size() != n) {
788 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
789 }
790
791 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
792 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
793 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
794 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
795 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
796 for (i=0; i<n ;++i) p[i]=in[i];
797 }
798 }
799 }
800 }
801
802
803 void
804 DataExpanded::reorderByReferenceIDs(int *reference_ids)
805 {
806 CHECK_FOR_EX_WRITE
807 int numSamples = getNumSamples();
808 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
809 int sampleNo, sampleNo2,i;
810 double* p,*p2;
811 register double rtmp;
812 FunctionSpace fs=getFunctionSpace();
813
814 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
815 const int id_in=reference_ids[sampleNo];
816 const int id=fs.getReferenceIDOfSample(sampleNo);
817 if (id!=id_in) {
818 bool matched=false;
819 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
820 if (id == reference_ids[sampleNo2]) {
821 p=&(m_data[getPointOffset(sampleNo,0)]);
822 p2=&(m_data[getPointOffset(sampleNo2,0)]);
823 for (i=0; i<n ;i++) {
824 rtmp=p[i];
825 p[i]=p2[i];
826 p2[i]=rtmp;
827 }
828 reference_ids[sampleNo]=id;
829 reference_ids[sampleNo2]=id_in;
830 matched=true;
831 break;
832 }
833 }
834 if (! matched) {
835 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
836 }
837 }
838 }
839 }
840
841 DataTypes::ValueType&
842 DataExpanded::getVectorRW()
843 {
844 CHECK_FOR_EX_WRITE
845 return m_data.getData();
846 }
847
848 const DataTypes::ValueType&
849 DataExpanded::getVectorRO() const
850 {
851 return m_data.getData();
852 }
853
854 // void DataExpanded::randomFill(long seed) {
855 // CHECK_FOR_EX_WRITE
856 //
857 // DataVector& dv=getVectorRW();
858 // const size_t dvsize=dv.size();
859 //
860 // esysUtils::randomFillArray(seed, &(dv[0]), dvsize);
861 //
862 // }
863
864 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26