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

Annotation of /trunk/escript/src/DataArrayView.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 100 - (hide annotations)
Wed Dec 15 03:48:48 2004 UTC (14 years, 10 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/DataArrayView.cpp
File size: 19310 byte(s)
*** empty log message ***

1 jgs 82 // $Id$
2    
3     /*
4     ******************************************************************************
5     * *
6     * COPYRIGHT ACcESS 2004 - All Rights Reserved *
7     * *
8     * This software is the property of ACcESS. No part of this code *
9     * may be copied in any form or by any means without the expressed written *
10     * consent of ACcESS. Copying, use or modification of this software *
11     * by any unauthorised person is illegal unless that person has a software *
12     * license agreement with ACcESS. *
13     * *
14     ******************************************************************************
15     */
16    
17     #include "escript/Data/DataArrayView.h"
18     #include "escript/Data/DataException.h"
19    
20     #include <sstream>
21     #include <vector>
22     #include <iostream>
23     #include <boost/python/extract.hpp>
24     #include <boost/python/object.hpp>
25    
26     using namespace std;
27     using namespace boost::python;
28    
29     namespace escript {
30    
31     DataArrayView::DataArrayView():
32     m_shape(),
33     m_offset(),
34     m_noValues()
35     {
36     m_data=0;
37     }
38    
39     DataArrayView::DataArrayView(ValueType& data,
40     const ShapeType& viewShape,
41     int offset):
42     m_data(&data),
43     m_shape(viewShape),
44     m_offset(offset),
45     m_noValues(noValues(viewShape))
46     {
47     //
48     // check the shape and throw an exception if an illegal shape is specified
49     if (viewShape.size() > m_maxRank) {
50     stringstream temp;
51     temp << "Error- Couldn't construct DataArrayView, invalid rank."
52     << "Input data rank: " << viewShape.size()
53     << " Maximum rank allowed for DataArrayView: "
54     << m_maxRank;
55     throw DataException(temp.str());
56     }
57 jgs 100 if (m_data->size()<(noValues()+m_offset)) {
58 jgs 82 stringstream temp;
59     temp << "Error- Couldn't construct DataArrayView, insufficient storage."
60     << " Shape requires: " << noValues()+m_offset << " values."
61     << " Data values available: " << m_data->size();
62     throw DataException(temp.str());
63     }
64     }
65    
66     DataArrayView::DataArrayView(const DataArrayView& other):
67     m_data(other.m_data),
68     m_offset(other.m_offset),
69     m_shape(other.m_shape),
70     m_noValues(other.m_noValues)
71     {
72     }
73    
74     bool
75     DataArrayView::isEmpty() const
76     {
77     return(m_data==0);
78     }
79    
80     void
81     DataArrayView::copy(const boost::python::numeric::array& value)
82     {
83     DataArrayView::ShapeType tempShape;
84     for (int i=0; i<value.getrank(); ++i) {
85     tempShape.push_back(extract<int>(value.getshape()[i]));
86     }
87     EsysAssert((!isEmpty()&&checkShape(tempShape)),
88     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
89    
90     if (value.getrank()==0) {
91     (*this)()=extract<double>(value[value.getshape()]);
92     } else if (value.getrank()==1) {
93     for (ValueType::size_type i=0;i<tempShape[0];++i) {
94     (*this)(i)=extract<double>(value[i]);
95     }
96     } else if (value.getrank()==2) {
97     for (ValueType::size_type i=0;i<tempShape[0];++i) {
98     for (ValueType::size_type j=0;j<tempShape[1];++j) {
99     (*this)(i,j)=extract<double>(value[i][j]);
100     }
101     }
102     } else if (value.getrank()==3) {
103     for (ValueType::size_type i=0;i<tempShape[0];++i) {
104     for (ValueType::size_type j=0;j<tempShape[1];++j) {
105     for (ValueType::size_type k=0;k<tempShape[2];++k) {
106     (*this)(i,j,k)=extract<double>(value[i][j][k]);
107     }
108     }
109     }
110     } else if (value.getrank()==4) {
111     for (ValueType::size_type i=0;i<tempShape[0];++i) {
112     for (ValueType::size_type j=0;j<tempShape[1];++j) {
113     for (ValueType::size_type k=0;k<tempShape[2];++k) {
114     for (ValueType::size_type m=0;m<tempShape[3];++m) {
115     (*this)(i,j,k,m)=extract<double>(value[i][j][k][m]);
116     }
117     }
118     }
119     }
120     }
121     }
122 jgs 100
123 jgs 82 void
124     DataArrayView::copy(const DataArrayView& other)
125     {
126     copy(m_offset,other);
127     }
128    
129     void
130     DataArrayView::copy(ValueType::size_type offset,
131     const DataArrayView& other)
132     {
133     EsysAssert((!isEmpty()&&checkShape(other.getShape())),
134     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
135     memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
136     }
137    
138     void
139     DataArrayView::copy(ValueType::size_type thisOffset,
140     const DataArrayView& other,
141     ValueType::size_type otherOffset)
142     {
143     EsysAssert((!isEmpty()&&checkShape(other.getShape())),
144     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
145     memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
146     }
147    
148     void
149     DataArrayView::copy(ValueType::size_type offset,
150     ValueType::value_type value)
151     {
152     //
153     // fill the entire view with the single value
154     EsysAssert(!isEmpty(),"Error - View is empty");
155     ValueType temp(noValues(),value);
156     memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
157     }
158    
159     int
160     DataArrayView::getRank() const
161     {
162     return m_shape.size();
163     }
164    
165     const DataArrayView::ShapeType&
166     DataArrayView::getShape() const
167     {
168     return m_shape;
169     }
170    
171     int
172     DataArrayView::noValues(const ShapeType& shape)
173     {
174     ShapeType::const_iterator i;
175     //
176     // An empty shape vector means rank 0 which contains 1 value
177     int noValues=1;
178     for (i=shape.begin();i!=shape.end();++i) {
179     noValues*=(*i);
180     }
181     return noValues;
182     }
183 jgs 100
184 jgs 82 int
185     DataArrayView::noValues() const
186     {
187     return m_noValues;
188     }
189    
190     bool
191     DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
192     {
193     return (m_shape==other);
194     }
195    
196     string
197     DataArrayView::createShapeErrorMessage(const string& messagePrefix,
198     const DataArrayView::ShapeType& other) const
199     {
200     stringstream temp;
201     temp << messagePrefix
202     << " This shape: " << shapeToString(m_shape)
203     << " Other shape: " << shapeToString(other);
204     return temp.str();
205     }
206    
207     void
208     DataArrayView::setOffset(ValueType::size_type offset)
209     {
210 jgs 100 EsysAssert((m_data->size()>=(noValues()+offset)), "Invalid offset");
211 jgs 82 m_offset=offset;
212     }
213    
214     DataArrayView::ShapeType
215     DataArrayView::getResultSliceShape(const RegionType& region)
216     {
217     RegionType::const_iterator i;
218     DataArrayView::ShapeType result;
219 jgs 100 for (i=region.begin();i!=region.end();++i) {
220     int dimSize=(i->second-i->first);
221 jgs 82 if (dimSize!=0) {
222     result.push_back(dimSize);
223     }
224     }
225     return result;
226     }
227    
228     void
229     DataArrayView::copySlice(const DataArrayView& other,
230     const RegionType& region)
231     {
232     //
233 jgs 100 // version of slice that uses the default offset
234 jgs 82 copySlice(m_offset,other,other.m_offset,region);
235     }
236    
237     void
238     DataArrayView::copySlice(ValueType::size_type thisOffset,
239     const DataArrayView& other,
240     ValueType::size_type otherOffset,
241     const RegionType& region)
242     {
243     //
244 jgs 100 // Assume region is okay, Other asserts will detect range errors
245     EsysAssert((other.getRank()==region.size()),"Error - Invalid slice region.");
246     EsysAssert((!isEmpty()&&checkShape(getResultSliceShape(region))),
247     createShapeErrorMessage("Error - Couldn't copy slice due to shape mismatch.",getResultSliceShape(region)));
248 jgs 82 //
249 jgs 100 // take a local copy of the region and modify for the case
250     // of 2:2 being equivalent to 2:3. The only difference in meaning being in
251     // determining the output shape.
252     RegionType localRegion=region;
253     for (int i=0;i<localRegion.size();++i) {
254     if (localRegion[i].first==localRegion[i].second) {
255     ++(localRegion[i].second);
256 jgs 82 }
257     }
258     int numCopy=0;
259 jgs 100 switch (localRegion.size()) {
260 jgs 82 case 0:
261     (*m_data)[thisOffset]=(*other.m_data)[otherOffset];
262     break;
263     case 1:
264 jgs 100 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
265 jgs 82 (*m_data)[numCopy+thisOffset]=(*other.m_data)[i+otherOffset];
266 jgs 100 ++numCopy;
267 jgs 82 }
268     break;
269     case 2:
270 jgs 100 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
271     for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
272 jgs 82 (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j)+otherOffset];
273 jgs 100 ++numCopy;
274 jgs 82 }
275     }
276     break;
277     case 3:
278 jgs 100 for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
279     for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
280     for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
281 jgs 82 (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k)+otherOffset];
282 jgs 100 ++numCopy;
283 jgs 82 }
284     }
285     }
286     break;
287     case 4:
288 jgs 100 for (int m=localRegion[2].first;m<localRegion[2].second;++m) {
289     for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
290     for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
291     for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
292     (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k,m)+otherOffset];
293     ++numCopy;
294 jgs 82 }
295     }
296     }
297     }
298     break;
299     default:
300     stringstream mess;
301 jgs 100 mess << "Error - (copySlice) Invalid rank.";
302 jgs 82 throw DataException(mess.str());
303     }
304     }
305 jgs 100
306 jgs 82 void
307     DataArrayView::copySliceFrom(const DataArrayView& other,
308     const RegionType& region)
309     {
310     //
311     // version of slice that uses the default offset
312     copySliceFrom(m_offset,other,other.m_offset,region);
313     }
314    
315     void
316     DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
317     const DataArrayView& other,
318     ValueType::size_type otherOffset,
319     const RegionType& region)
320     {
321     //
322 jgs 100 // Assume region is okay, Other asserts will detect range errors
323     EsysAssert((getRank()==region.size()),"Error - Invalid slice region.");
324     EsysAssert((!isEmpty()&&other.checkShape(getResultSliceShape(region))),
325     createShapeErrorMessage("Error - Couldn't copy slice due to shape mismatch.",
326     getResultSliceShape(region)));
327 jgs 82 //
328 jgs 100 // take a local copy of the region and modify for the case
329     // of 2:2 being equivalent to 2:3. The only difference in meaning being in
330     // determining the output shape.
331     RegionType localRegion=region;
332     for (int i=0;i<localRegion.size();++i) {
333     if (localRegion[i].first==localRegion[i].second) {
334     ++(localRegion[i].second);
335 jgs 82 }
336     }
337     int numCopy=0;
338 jgs 100 switch (localRegion.size()) {
339 jgs 82 case 0:
340     (*m_data)[thisOffset]=(*other.m_data)[otherOffset];
341     break;
342     case 1:
343 jgs 100 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
344 jgs 82 (*m_data)[i+thisOffset]=(*other.m_data)[numCopy+otherOffset];
345 jgs 100 ++numCopy;
346 jgs 82 }
347     break;
348     case 2:
349 jgs 100 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
350     for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
351 jgs 82 (*m_data)[relIndex(i,j)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
352 jgs 100 ++numCopy;
353 jgs 82 }
354     }
355     break;
356     case 3:
357 jgs 100 for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
358     for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
359     for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
360 jgs 82 (*m_data)[relIndex(i,j,k)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
361 jgs 100 ++numCopy;
362 jgs 82 }
363     }
364     }
365     break;
366     case 4:
367 jgs 100 for (int m=localRegion[2].first;m<localRegion[2].second;++m) {
368     for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
369     for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
370     for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
371     (*m_data)[relIndex(i,j,k,m)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
372     ++numCopy;
373 jgs 82 }
374     }
375     }
376     }
377     break;
378     default:
379     stringstream mess;
380 jgs 100 mess << "Error - (copySliceFrom) Invalid rank: " << localRegion.size() << " for region.";
381 jgs 82 throw DataException(mess.str());
382     }
383     }
384 jgs 100
385 jgs 82 DataArrayView::ShapeType
386     DataArrayView::determineResultShape(const DataArrayView& left,
387     const DataArrayView& right)
388     {
389     DataArrayView::ShapeType temp;
390 jgs 100 for (int i=0; i<(left.getRank()-1); ++i) {
391 jgs 82 temp.push_back(left.getShape()[i]);
392     }
393 jgs 100 for (int i=1; i<right.getRank(); ++i) {
394 jgs 82 temp.push_back(right.getShape()[i]);
395     }
396     return temp;
397     }
398    
399     string
400     DataArrayView::toString(const string& suffix) const
401     {
402     EsysAssert(!isEmpty(),"Error - View is empty");
403     stringstream temp;
404     string finalSuffix=suffix;
405     if (suffix.length() > 0) {
406     finalSuffix+=" ";
407     }
408     switch (getRank()) {
409     case 0:
410     temp << finalSuffix << (*this)();
411     break;
412     case 1:
413 jgs 100 for (int i=0;i<getShape()[0];++i) {
414 jgs 82 temp << "(" << i << ") " << finalSuffix << (*this)(i);
415     if (i!=(getShape()[0]-1)) {
416     temp << endl;
417     }
418     }
419     break;
420     case 2:
421 jgs 100 for (int i=0;i<getShape()[0];++i) {
422     for (int j=0;j<getShape()[1];++j) {
423 jgs 82 temp << "(" << i << "," << j << ") " << finalSuffix << (*this)(i,j);
424     if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
425     temp << endl;
426     }
427     }
428     }
429     break;
430     case 3:
431 jgs 100 for (int i=0;i<getShape()[0];++i) {
432     for (int j=0;j<getShape()[1];++j) {
433     for (int k=0;k<getShape()[2];++k) {
434     temp << "("
435     << i << ","
436     << j << ","
437     << k << ") "
438     << finalSuffix << (*this)(i,j,k);
439     //
440     // don't put an endl after the last value
441     if (!(i==(getShape()[0]-1) &&
442     j==(getShape()[1]-1) &&
443     k==(getShape()[2]-1))) {
444 jgs 82 temp << endl;
445     }
446     }
447     }
448     }
449     break;
450     case 4:
451 jgs 100 // break;
452 jgs 82 default:
453     stringstream mess;
454     mess << "Error - (toString) Invalid rank: " << getRank();
455     throw DataException(mess.str());
456     }
457     return temp.str();
458     }
459    
460     string
461     DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)
462     {
463     stringstream temp;
464     temp << "(";
465 jgs 100 for (int i=0;i<shape.size();++i) {
466 jgs 82 temp << shape[i];
467 jgs 100 if (shape.size()-i != 1) {
468 jgs 82 temp << ",";
469     }
470     }
471     temp << ")";
472     return temp.str();
473     }
474    
475     void
476     DataArrayView::matMult(const DataArrayView& left,
477     const DataArrayView& right,
478     DataArrayView& result)
479     {
480     if (left.getRank()==0 || right.getRank()==0) {
481     stringstream temp;
482     temp << "Error - (matMult) Invalid for rank 0 objects.";
483     throw DataException(temp.str());
484     }
485    
486     if (left.getShape()[left.getRank()-1]!=right.getShape()[0]) {
487     stringstream temp;
488     temp << "Error - (matMult) Dimension: " << left.getRank()
489     << ", size: " << left.getShape()[left.getRank()-1]
490     << " of LHS and dimension: 1, size: " << right.getShape()[0]
491     << " of RHS don't match.";
492     throw DataException(temp.str());
493     }
494     int outputRank=left.getRank()+right.getRank()-2;
495    
496     if (outputRank < 0) {
497     stringstream temp;
498     temp << "Error - (matMult) Left and right operands cannot be multiplied "
499     << "as they have incompatable rank.";
500     throw DataException(temp.str());
501     }
502    
503     if (outputRank != result.getRank()) {
504     stringstream temp;
505     temp << "Error - (matMult) Rank of result array is: "
506     << result.getRank()
507     << " it must be: " << outputRank;
508     throw DataException(temp.str());
509     }
510    
511     for (int i=0; i<(left.getRank()-1); ++i) {
512     if (left.getShape()[i]!=result.getShape()[i]) {
513     stringstream temp;
514     temp << "Error - (matMult) Dimension: " << i
515     << " of LHS and output array don't match.";
516     throw DataException(temp.str());
517     }
518     }
519    
520     for (int i=1; i<right.getRank(); ++i) {
521     if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
522     stringstream temp;
523     temp << "Error - (matMult) Dimension: " << i
524     << ", size: " << right.getShape()[i]
525     << " of RHS and dimension: " << i+left.getRank()-1
526     << ", size: " << result.getShape()[i+left.getRank()-1]
527     << " of output array don't match.";
528     throw DataException(temp.str());
529     }
530     }
531    
532     switch (left.getRank()) {
533     case 1:
534     switch (right.getRank()) {
535     case 1:
536     result()=0;
537     for (int i=0;i<left.getShape()[0];++i) {
538     result()+=left(i)*right(i);
539     }
540     break;
541     case 2:
542     for (int i=0;i<result.getShape()[0];++i) {
543     result(i)=0;
544     for (int j=0;j<right.getShape()[0];++j) {
545     result(i)+=left(j)*right(j,i);
546     }
547     }
548     break;
549     default:
550     stringstream temp;
551     temp << "Error - (matMult) Invalid rank. Programming error.";
552     throw DataException(temp.str());
553     break;
554     }
555     break;
556     case 2:
557     switch (right.getRank()) {
558     case 1:
559     result()=0;
560     for (int i=0;i<left.getShape()[0];++i) {
561     result(i)=0;
562     for (int j=0;j<left.getShape()[1];++j) {
563     result(i)+=left(i,j)*right(i);
564     }
565     }
566     break;
567     case 2:
568     for (int i=0;i<result.getShape()[0];++i) {
569     for (int j=0;j<result.getShape()[1];++j) {
570     result(i,j)=0;
571     for (int jR=0;jR<right.getShape()[0];++jR) {
572     result(i,j)+=left(i,jR)*right(jR,j);
573     }
574     }
575     }
576     break;
577     default:
578     stringstream temp;
579     temp << "Error - (matMult) Invalid rank. Programming error.";
580     throw DataException(temp.str());
581     break;
582     }
583     break;
584     default:
585     stringstream temp;
586     temp << "Error - (matMult) Not supported for rank: " << left.getRank();
587     throw DataException(temp.str());
588     break;
589     }
590     }
591    
592 jgs 100 DataArrayView::RegionType
593     DataArrayView::getSliceRegion(const boost::python::object& key) const
594     {
595     int slice_len=0,i=0,rank=getRank(),out_len=0;
596     extract<tuple> key_tuple(key);
597     if (key_tuple.check()) {
598     slice_len=extract<int>(key.attr("__len__")());
599     } else {
600     slice_len=1;
601     }
602     if (slice_len>=rank) {
603     out_len=slice_len;
604     } else {
605     out_len=rank;
606     }
607     DataArrayView::RegionType out(out_len);
608     if (key_tuple.check()) {
609     for (i=0;i<slice_len;i++) out[i]=getSliceRange(getShape()[i],key[i]);
610     } else {
611     out[0]=getSliceRange(getShape()[0],key);
612     }
613     for (i=slice_len;i<rank;i++) out[i]=std::pair<int,int>(0,getShape()[i]);
614     // throw DataException("Error - number of slices is bigger then the rank.");
615     return out;
616     }
617    
618     std::pair<int,int>
619     getSliceRange(const int s,
620     const boost::python::object& key)
621     {
622     int s1=s;
623     int s0=0;
624     /* if the key is an int we set start=end=key */
625     extract<int> slice_int(key);
626     if (slice_int.check()) {
627     s0=slice_int();
628     s1=s0;
629     } else {
630     extract<int> step(key.attr("step"));
631     if (step.check() && step()!=1) {
632     throw DataException("Error - Data does not support increments in slicing ");
633     } else {
634     extract<int> start(key.attr("start"));
635     if (start.check()) s0=start();
636     extract<int> stop(key.attr("stop"));
637     if (stop.check()) s1=stop();
638     }
639     }
640     return std::pair<int,int>(s0,s1);
641     }
642    
643 jgs 82 bool operator==(const DataArrayView& left, const DataArrayView& right)
644     {
645     //
646     // The views are considered equal if they have the same shape and each value
647     // of the view is the same.
648     bool result=(left.m_shape==right.m_shape);
649     if (result) {
650     for (int i=0;i<left.noValues();++i) {
651     result=result && (left.getData(i)==right.getData(i));
652     }
653     }
654     return result;
655     }
656    
657     bool operator!=(const DataArrayView& left, const DataArrayView& right)
658     {
659     return (!operator==(left,right));
660     }
661    
662     } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26