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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (hide annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 10 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/DataArrayView.cpp
File size: 23197 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 102 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 102
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 102
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 102 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 jgs 102 int dimSize;
218 jgs 82 RegionType::const_iterator i;
219     DataArrayView::ShapeType result;
220 jgs 102 for (i=region.begin();i!=region.end();i++) {
221     dimSize=((i->second) - (i->first));
222 jgs 82 if (dimSize!=0) {
223     result.push_back(dimSize);
224     }
225     }
226     return result;
227     }
228    
229 jgs 102 DataArrayView::RegionType
230     DataArrayView::getSliceRegion(const boost::python::object& key) const
231     {
232     int slice_rank, out_rank, i;
233     int this_rank=getRank();
234     /* allow for case where key is singular eg: [1], this implies we
235     want to generate a rank-1 dimension object, as opposed to eg: [1,2]
236     which implies we want to take a rank dimensional object with one
237     dimension of size 1 */
238     extract<tuple> key_tuple(key);
239     if (key_tuple.check()) {
240     slice_rank=extract<int> (key.attr("__len__")());
241     } else {
242     slice_rank=1;
243     }
244     /* ensure slice is correctly dimensioned */
245     if (slice_rank>=this_rank) {
246     out_rank=slice_rank;
247     } else {
248     out_rank=this_rank;
249     }
250     DataArrayView::RegionType out(out_rank);
251     /* calculate values for slice range */
252     if (key_tuple.check()) {
253     for (i=0;i<slice_rank;i++) {
254     out[i]=getSliceRange(key[i],getShape()[i]);
255     }
256     } else {
257     out[0]=getSliceRange(key,getShape()[0]);
258     }
259     for (i=slice_rank;i<this_rank;i++) {
260     out[i]=std::pair<int,int>(0,getShape()[i]);
261     }
262     return out;
263     }
264    
265     std::pair<int,int>
266     getSliceRange(const boost::python::object& key,
267     const int shape)
268     {
269     /* default slice range is range of entire shape dimension */
270     int s0=0, s1=shape;;
271     extract<int> slice_int(key);
272     if (slice_int.check()) {
273     /* if the key is a single int set start=key and end=key */
274     /* in this case, we want to return a rank-1 dimension object from
275     this object, taken from a single index value for one of this
276     object's dimensions */
277     s0=slice_int();
278     s1=s0;
279     } else {
280     /* if key is a pair extract begin and end values */
281     extract<int> step(key.attr("step"));
282     if (step.check() && step()!=1) {
283     throw DataException("Error - Data does not support increments in slicing ");
284     } else {
285     extract<int> start(key.attr("start"));
286     if (start.check()) {
287     s0=start();
288     }
289     extract<int> stop(key.attr("stop"));
290     if (stop.check()) {
291     s1=stop();
292     }
293     }
294     }
295     return std::pair<int,int>(s0,s1);
296     }
297    
298 jgs 82 void
299     DataArrayView::copySlice(const DataArrayView& other,
300     const RegionType& region)
301     {
302     //
303 jgs 102 // version of copySlice that uses the default offsets
304 jgs 82 copySlice(m_offset,other,other.m_offset,region);
305     }
306    
307     void
308     DataArrayView::copySlice(ValueType::size_type thisOffset,
309     const DataArrayView& other,
310     ValueType::size_type otherOffset,
311     const RegionType& region)
312     {
313 jgs 102
314 jgs 82 //
315 jgs 102 // Make sure this object is not empty
316     EsysAssert(!isEmpty(),"Error - this data object is empty.");
317    
318 jgs 82 //
319 jgs 102 // Modify region to copy from in order to
320     // deal with the case where one range in the region contains identical indexes,
321     // eg: <<1,1><0,3><0,3>>
322     // This situation implies we want to copy from an object with rank greater than that of this
323     // object. eg: we want to copy the values from a two dimensional slice out of a three
324     // dimensional object into a two dimensional object.
325     // We do this by taking a slice from the other object where one dimension of
326     // the slice region is of size 1. So in the above example, we modify the above
327     // region like so: <<1,2><0,3><0,3>> and take this slice.
328    
329     RegionType slice_region=region;
330     for (int i=0;i<slice_region.size();i++) {
331     if (slice_region[i].first==slice_region[i].second) {
332     (slice_region[i].second)++;
333 jgs 82 }
334     }
335 jgs 102
336     //
337     // Check the object to be sliced from is compatible with the region to be sliced,
338     // and that the region to be sliced is compatible with this object:
339     // 1: the region to be sliced should be the same rank as the other object.
340     // 2: the other object should have at least as many datapoints at the region requires.
341     // 3: the region should have the same number of datapoints as this object
342    
343     DataArrayView::ShapeType slice_shape=getResultSliceShape(slice_region);
344     EsysAssert(other.getRank()==slice_region.size(),
345     "Error - slice rank incompatible with object to be sliced from.");
346     EsysAssert((noValues(slice_shape)<=other.noValues()),
347     "Error - slice shape incompatible with object to be sliced from.");
348     EsysAssert((noValues(slice_shape)==noValues()),
349     "Error - slice shape incompatible with this object.");
350    
351     //
352     // copy the values in the specified slice_region of the other object into this object
353    
354 jgs 82 int numCopy=0;
355 jgs 102 switch (slice_region.size()) {
356 jgs 82 case 0:
357 jgs 102 /* this case should never be encountered, as python will never pass us an empty slice */
358 jgs 82 (*m_data)[thisOffset]=(*other.m_data)[otherOffset];
359     break;
360     case 1:
361 jgs 102 for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
362 jgs 82 (*m_data)[numCopy+thisOffset]=(*other.m_data)[i+otherOffset];
363 jgs 102 numCopy++;
364 jgs 82 }
365     break;
366     case 2:
367 jgs 102 for (int j=slice_region[1].first;j<slice_region[1].second;j++) {
368     for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
369 jgs 82 (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j)+otherOffset];
370 jgs 102 numCopy++;
371 jgs 82 }
372     }
373     break;
374     case 3:
375 jgs 102 for (int k=slice_region[2].first;k<slice_region[2].second;k++) {
376     for (int j=slice_region[1].first;j<slice_region[1].second;j++) {
377     for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
378 jgs 82 (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k)+otherOffset];
379 jgs 102 numCopy++;
380 jgs 82 }
381     }
382     }
383     break;
384     case 4:
385 jgs 102 for (int l=slice_region[3].first;l<slice_region[3].second;l++) {
386     for (int k=slice_region[2].first;k<slice_region[2].second;k++) {
387     for (int j=slice_region[1].first;j<slice_region[1].second;j++) {
388     for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
389     (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k,l)+otherOffset];
390     numCopy++;
391 jgs 82 }
392     }
393     }
394     }
395     break;
396     default:
397     stringstream mess;
398 jgs 102 mess << "Error - (copySlice) Invalid rank: " << slice_region.size() << " for region.";
399 jgs 82 throw DataException(mess.str());
400     }
401     }
402 jgs 102
403 jgs 82 void
404     DataArrayView::copySliceFrom(const DataArrayView& other,
405     const RegionType& region)
406     {
407     //
408     // version of slice that uses the default offset
409     copySliceFrom(m_offset,other,other.m_offset,region);
410     }
411    
412     void
413     DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
414     const DataArrayView& other,
415     ValueType::size_type otherOffset,
416     const RegionType& region)
417     {
418 jgs 102
419 jgs 82 //
420 jgs 102 // Make sure this object is not empty
421     EsysAssert(!isEmpty(),"Error - this data object is empty.");
422    
423 jgs 82 //
424 jgs 102 // Modify region to copy in order to
425     // deal with the case where one range in the region contains identical idexes,
426     // eg: <<1,1><0,3><0,3>>
427     // This situation implies we want to copy from an object with rank less than that of this
428     // object. eg: we want to copy the values from a two dimensional object into a two dimensional
429     // slice of a three dimensional object.
430     // We do this by copying to a slice from this object where one dimension of
431     // the slice region is of size 1. So in the above example, we modify the above
432     // region like so: <<1,2><0,3><0,3>> and copy into this slice.
433    
434     RegionType slice_region=region;
435     for (int i=0;i<slice_region.size();i++) {
436     if (slice_region[i].first==slice_region[i].second) {
437     (slice_region[i].second)++;
438 jgs 82 }
439     }
440 jgs 102
441     //
442     // Check the object to be coped from is compatible with the slice region to be copied
443     // into, and that the region to be sliced into is compatible with this object:
444     // 1: the region to be sliced into should be the same rank as this object.
445     // 2: the other object should have exactly the same number of data points as the
446     // region to be sliced into requires.
447     // 3: the region can be a higher rank than the other object, provided 2 holds.
448     // 4: this object must have at least as many data points as the region requires.
449    
450     //
451     // Check region size and shape
452     DataArrayView::ShapeType slice_shape=getResultSliceShape(slice_region);
453     EsysAssert((getRank()==slice_region.size()),"Error - Invalid slice region.");
454     EsysAssert((noValues(slice_shape)==other.noValues()),
455     "Error - slice shape incompatible with object to be sliced from.");
456     EsysAssert(other.getRank()<=slice_region.size(),
457     "Error - slice rank incompatible with object to be sliced from.");
458     EsysAssert((noValues(slice_shape)<=noValues()),
459     "Error - slice shape incompatible with this object.");
460    
461     //
462     // copy the values in the other object into the specified slice_region of this object
463    
464 jgs 82 int numCopy=0;
465 jgs 102 switch (slice_region.size()) {
466 jgs 82 case 0:
467 jgs 102 /* this case should never be encountered, as python will never pass us an empty slice */
468 jgs 82 (*m_data)[thisOffset]=(*other.m_data)[otherOffset];
469     break;
470     case 1:
471 jgs 102 for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
472 jgs 82 (*m_data)[i+thisOffset]=(*other.m_data)[numCopy+otherOffset];
473 jgs 102 numCopy++;
474 jgs 82 }
475     break;
476     case 2:
477 jgs 102 for (int j=slice_region[1].first;j<slice_region[1].second;j++) {
478     for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
479 jgs 82 (*m_data)[relIndex(i,j)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
480 jgs 102 numCopy++;
481 jgs 82 }
482     }
483     break;
484     case 3:
485 jgs 102 for (int k=slice_region[2].first;k<slice_region[2].second;k++) {
486     for (int j=slice_region[1].first;j<slice_region[1].second;j++) {
487     for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
488 jgs 82 (*m_data)[relIndex(i,j,k)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
489 jgs 102 numCopy++;
490 jgs 82 }
491     }
492     }
493     break;
494     case 4:
495 jgs 102 for (int l=slice_region[3].first;l<slice_region[3].second;l++) {
496     for (int k=slice_region[2].first;k<slice_region[2].second;k++) {
497     for (int j=slice_region[1].first;j<slice_region[1].second;j++) {
498     for (int i=slice_region[0].first;i<slice_region[0].second;i++) {
499     (*m_data)[relIndex(i,j,k,l)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
500     numCopy++;
501 jgs 82 }
502     }
503     }
504     }
505     break;
506     default:
507     stringstream mess;
508 jgs 102 mess << "Error - (copySliceFrom) Invalid rank: " << slice_region.size() << " for region.";
509 jgs 82 throw DataException(mess.str());
510     }
511     }
512 jgs 102
513 jgs 82 DataArrayView::ShapeType
514     DataArrayView::determineResultShape(const DataArrayView& left,
515     const DataArrayView& right)
516     {
517     DataArrayView::ShapeType temp;
518 jgs 102 for (int i=0; i<(left.getRank()-1); i++) {
519 jgs 82 temp.push_back(left.getShape()[i]);
520     }
521 jgs 102 for (int i=1; i<right.getRank(); i++) {
522 jgs 82 temp.push_back(right.getShape()[i]);
523     }
524     return temp;
525     }
526    
527     string
528     DataArrayView::toString(const string& suffix) const
529     {
530     EsysAssert(!isEmpty(),"Error - View is empty");
531     stringstream temp;
532     string finalSuffix=suffix;
533     if (suffix.length() > 0) {
534     finalSuffix+=" ";
535     }
536     switch (getRank()) {
537     case 0:
538     temp << finalSuffix << (*this)();
539     break;
540     case 1:
541 jgs 102 for (int i=0;i<getShape()[0];i++) {
542 jgs 82 temp << "(" << i << ") " << finalSuffix << (*this)(i);
543     if (i!=(getShape()[0]-1)) {
544     temp << endl;
545     }
546     }
547     break;
548     case 2:
549 jgs 102 for (int i=0;i<getShape()[0];i++) {
550     for (int j=0;j<getShape()[1];j++) {
551 jgs 82 temp << "(" << i << "," << j << ") " << finalSuffix << (*this)(i,j);
552     if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
553     temp << endl;
554     }
555     }
556     }
557     break;
558     case 3:
559 jgs 102 for (int i=0;i<getShape()[0];i++) {
560     for (int j=0;j<getShape()[1];j++) {
561     for (int k=0;k<getShape()[2];k++) {
562     temp << "(" << i << "," << j << "," << k << ") " << finalSuffix << (*this)(i,j,k);
563     if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1))) {
564 jgs 82 temp << endl;
565     }
566     }
567     }
568     }
569     break;
570     case 4:
571 jgs 102 for (int i=0;i<getShape()[0];i++) {
572     for (int j=0;j<getShape()[1];j++) {
573     for (int k=0;k<getShape()[2];k++) {
574     for (int l=0;l<getShape()[3];l++) {
575     temp << "(" << i << "," << j << "," << k << "," << l << ") " << finalSuffix << (*this)(i,j,k,l);
576     if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1) && l==(getShape()[3]-1))) {
577     temp << endl;
578     }
579     }
580     }
581     }
582     }
583     break;
584 jgs 82 default:
585     stringstream mess;
586     mess << "Error - (toString) Invalid rank: " << getRank();
587     throw DataException(mess.str());
588     }
589     return temp.str();
590     }
591    
592     string
593     DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)
594     {
595     stringstream temp;
596     temp << "(";
597 jgs 102 for (int i=0;i<shape.size();i++) {
598 jgs 82 temp << shape[i];
599 jgs 102 if (i < shape.size()-1) {
600 jgs 82 temp << ",";
601     }
602     }
603     temp << ")";
604     return temp.str();
605     }
606    
607     void
608     DataArrayView::matMult(const DataArrayView& left,
609     const DataArrayView& right,
610     DataArrayView& result)
611     {
612     if (left.getRank()==0 || right.getRank()==0) {
613     stringstream temp;
614     temp << "Error - (matMult) Invalid for rank 0 objects.";
615     throw DataException(temp.str());
616     }
617    
618     if (left.getShape()[left.getRank()-1]!=right.getShape()[0]) {
619     stringstream temp;
620     temp << "Error - (matMult) Dimension: " << left.getRank()
621     << ", size: " << left.getShape()[left.getRank()-1]
622     << " of LHS and dimension: 1, size: " << right.getShape()[0]
623     << " of RHS don't match.";
624     throw DataException(temp.str());
625     }
626     int outputRank=left.getRank()+right.getRank()-2;
627    
628     if (outputRank < 0) {
629     stringstream temp;
630     temp << "Error - (matMult) Left and right operands cannot be multiplied "
631     << "as they have incompatable rank.";
632     throw DataException(temp.str());
633     }
634    
635     if (outputRank != result.getRank()) {
636     stringstream temp;
637     temp << "Error - (matMult) Rank of result array is: "
638     << result.getRank()
639     << " it must be: " << outputRank;
640     throw DataException(temp.str());
641     }
642    
643     for (int i=0; i<(left.getRank()-1); ++i) {
644     if (left.getShape()[i]!=result.getShape()[i]) {
645     stringstream temp;
646     temp << "Error - (matMult) Dimension: " << i
647     << " of LHS and output array don't match.";
648     throw DataException(temp.str());
649     }
650     }
651    
652     for (int i=1; i<right.getRank(); ++i) {
653     if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
654     stringstream temp;
655     temp << "Error - (matMult) Dimension: " << i
656     << ", size: " << right.getShape()[i]
657     << " of RHS and dimension: " << i+left.getRank()-1
658     << ", size: " << result.getShape()[i+left.getRank()-1]
659     << " of output array don't match.";
660     throw DataException(temp.str());
661     }
662     }
663    
664     switch (left.getRank()) {
665     case 1:
666     switch (right.getRank()) {
667     case 1:
668     result()=0;
669     for (int i=0;i<left.getShape()[0];++i) {
670     result()+=left(i)*right(i);
671     }
672     break;
673     case 2:
674     for (int i=0;i<result.getShape()[0];++i) {
675     result(i)=0;
676     for (int j=0;j<right.getShape()[0];++j) {
677     result(i)+=left(j)*right(j,i);
678     }
679     }
680     break;
681     default:
682     stringstream temp;
683     temp << "Error - (matMult) Invalid rank. Programming error.";
684     throw DataException(temp.str());
685     break;
686     }
687     break;
688     case 2:
689     switch (right.getRank()) {
690     case 1:
691     result()=0;
692     for (int i=0;i<left.getShape()[0];++i) {
693     result(i)=0;
694     for (int j=0;j<left.getShape()[1];++j) {
695     result(i)+=left(i,j)*right(i);
696     }
697     }
698     break;
699     case 2:
700     for (int i=0;i<result.getShape()[0];++i) {
701     for (int j=0;j<result.getShape()[1];++j) {
702     result(i,j)=0;
703     for (int jR=0;jR<right.getShape()[0];++jR) {
704     result(i,j)+=left(i,jR)*right(jR,j);
705     }
706     }
707     }
708     break;
709     default:
710     stringstream temp;
711     temp << "Error - (matMult) Invalid rank. Programming error.";
712     throw DataException(temp.str());
713     break;
714     }
715     break;
716     default:
717     stringstream temp;
718     temp << "Error - (matMult) Not supported for rank: " << left.getRank();
719     throw DataException(temp.str());
720     break;
721     }
722     }
723    
724     bool operator==(const DataArrayView& left, const DataArrayView& right)
725     {
726     //
727     // The views are considered equal if they have the same shape and each value
728     // of the view is the same.
729     bool result=(left.m_shape==right.m_shape);
730     if (result) {
731     for (int i=0;i<left.noValues();++i) {
732     result=result && (left.getData(i)==right.getData(i));
733     }
734     }
735     return result;
736     }
737    
738     bool operator!=(const DataArrayView& left, const DataArrayView& right)
739     {
740     return (!operator==(left,right));
741     }
742    
743     } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26