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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 576 - (show annotations)
Fri Mar 3 08:28:42 2006 UTC (13 years, 7 months ago) by gross
File size: 27767 byte(s)
some steps towards eigenvalue and eigenvector calculation
1 // $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 "DataArrayView.h"
18 #include "DataException.h"
19
20 #include <sstream>
21
22 #include <boost/python/extract.hpp>
23
24 using namespace std;
25 using namespace boost::python;
26
27 namespace escript {
28
29 DataArrayView::DataArrayView():
30 m_shape(),
31 m_offset(0),
32 m_noValues(0)
33 {
34 m_data=0;
35 }
36
37 DataArrayView::DataArrayView(ValueType& data,
38 const ShapeType& viewShape,
39 int offset):
40 m_data(&data),
41 m_shape(viewShape),
42 m_offset(offset),
43 m_noValues(noValues(viewShape))
44 {
45 //
46 // check the shape rank and size and throw an exception if a
47 // shape incompatible with the given data array is specified
48 if (viewShape.size() > m_maxRank) {
49 stringstream temp;
50 temp << "Error- Couldn't construct DataArrayView, invalid rank."
51 << "Input data rank: " << viewShape.size() << " "
52 << "Maximum rank allowed for DataArrayView: " << m_maxRank;
53 throw DataException(temp.str());
54 }
55 if (m_data->size() < (m_offset+noValues())) {
56 stringstream temp;
57 temp << "Error- Couldn't construct DataArrayView, insufficient data values. "
58 << "Shape requires: " << noValues()+m_offset << " values."
59 << "Data values available: " << m_data->size();
60 throw DataException(temp.str());
61 }
62 }
63
64 DataArrayView::DataArrayView(const DataArrayView& other):
65 m_data(other.m_data),
66 m_offset(other.m_offset),
67 m_shape(other.m_shape),
68 m_noValues(other.m_noValues)
69 {
70 }
71
72 bool
73 DataArrayView::isEmpty() const
74 {
75 return (m_data==0);
76 }
77
78 void
79 DataArrayView::copy(const boost::python::numeric::array& value)
80 {
81 ShapeType tempShape;
82 for (int i=0; i<value.getrank(); i++) {
83 tempShape.push_back(extract<int>(value.getshape()[i]));
84 }
85
86 EsysAssert((!isEmpty()&&checkShape(tempShape)),
87 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
88
89 if (value.getrank()==0) {
90 (*this)()=extract<double>(value[value.getshape()]);
91 } else if (value.getrank()==1) {
92 for (ValueType::size_type i=0;i<tempShape[0];i++) {
93 (*this)(i)=extract<double>(value[i]);
94 }
95 } else if (value.getrank()==2) {
96 for (ValueType::size_type i=0;i<tempShape[0];i++) {
97 for (ValueType::size_type j=0;j<tempShape[1];j++) {
98 (*this)(i,j)=extract<double>(value[i][j]);
99 }
100 }
101 } else if (value.getrank()==3) {
102 for (ValueType::size_type i=0;i<tempShape[0];i++) {
103 for (ValueType::size_type j=0;j<tempShape[1];j++) {
104 for (ValueType::size_type k=0;k<tempShape[2];k++) {
105 (*this)(i,j,k)=extract<double>(value[i][j][k]);
106 }
107 }
108 }
109 } else if (value.getrank()==4) {
110 for (ValueType::size_type i=0;i<tempShape[0];i++) {
111 for (ValueType::size_type j=0;j<tempShape[1];j++) {
112 for (ValueType::size_type k=0;k<tempShape[2];k++) {
113 for (ValueType::size_type l=0;l<tempShape[3];l++) {
114 (*this)(i,j,k,l)=extract<double>(value[i][j][k][l]);
115 }
116 }
117 }
118 }
119 }
120 }
121
122 void
123 DataArrayView::copy(const DataArrayView& other)
124 {
125 copy(m_offset,other);
126 }
127
128 void
129 DataArrayView::copy(ValueType::size_type offset,
130 const DataArrayView& other)
131 {
132 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),
133 "Error - Couldn't copy due to insufficient storage.");
134 EsysAssert((checkShape(other.getShape())),
135 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
136 if (checkOffset(offset)) {
137 memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
138 } else {
139 throw DataException("Error - invalid offset specified.");
140 }
141 }
142
143 void
144 DataArrayView::copy(ValueType::size_type thisOffset,
145 const DataArrayView& other,
146 ValueType::size_type otherOffset)
147 {
148 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),
149 "Error - Couldn't copy due to insufficient storage.");
150 EsysAssert((checkShape(other.getShape())),
151 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
152 if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {
153 memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
154 } else {
155 throw DataException("Error - invalid offset specified.");
156 }
157 }
158
159 void
160 DataArrayView::copy(ValueType::size_type offset,
161 ValueType::value_type value)
162 {
163 EsysAssert((!isEmpty()&&checkOffset(offset)),
164 "Error - Couldn't copy due to insufficient storage.");
165 if (checkOffset(offset)) {
166 vector<double> temp(noValues(),value);
167 memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
168 } else {
169 throw DataException("Error - invalid offset specified.");
170 }
171 }
172
173 int
174 DataArrayView::getRank() const
175 {
176 return m_shape.size();
177 }
178
179 const
180 DataArrayView::ShapeType&
181 DataArrayView::getShape() const
182 {
183 return m_shape;
184 }
185
186 int
187 DataArrayView::noValues(const ShapeType& shape)
188 {
189 ShapeType::const_iterator i;
190 //
191 // An empty shape vector means rank 0 which contains 1 value
192 int noValues=1;
193 for (i=shape.begin();i!=shape.end();i++) {
194 noValues*=(*i);
195 }
196 return noValues;
197 }
198
199 int
200 DataArrayView::noValues(const RegionLoopRangeType& region)
201 {
202 //
203 // An empty region vector means rank 0 which contains 1 value
204 int noValues=1;
205 for (int i=0;i<region.size();i++) {
206 noValues*=region[i].second-region[i].first;
207 }
208 return noValues;
209 }
210
211 int
212 DataArrayView::noValues() const
213 {
214 return m_noValues;
215 }
216
217 bool
218 DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
219 {
220 return (m_shape==other);
221 }
222
223 string
224 DataArrayView::createShapeErrorMessage(const string& messagePrefix,
225 const DataArrayView::ShapeType& other) const
226 {
227 stringstream temp;
228 temp << messagePrefix
229 << " This shape: " << shapeToString(m_shape)
230 << " Other shape: " << shapeToString(other);
231 return temp.str();
232 }
233
234 DataArrayView::ValueType::size_type
235 DataArrayView::getOffset() const
236 {
237 return m_offset;
238 }
239
240 void
241 DataArrayView::setOffset(ValueType::size_type offset)
242 {
243 EsysAssert((checkOffset(offset)), "Error - Invalid offset.");
244 if (checkOffset(offset)) {
245 m_offset=offset;
246 } else {
247 throw DataException("Error - invalid offset specified.");
248 }
249 }
250
251 void
252 DataArrayView::incrOffset()
253 {
254 EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");
255 if (checkOffset(m_offset+noValues())) {
256 m_offset=m_offset+noValues();
257 } else {
258 throw DataException("Error - Cannot increment offset.");
259 }
260 }
261
262 bool
263 DataArrayView::checkOffset() const
264 {
265 return checkOffset(m_offset);
266 }
267
268 bool
269 DataArrayView::checkOffset(ValueType::size_type offset) const
270 {
271 return (m_data->size() >= (offset+noValues()));
272 }
273
274 DataArrayView::ValueType&
275 DataArrayView::getData() const
276 {
277 EsysAssert(!isEmpty(),"Error - View is empty");
278 return *m_data;
279 }
280
281 DataArrayView::ValueType::reference
282 DataArrayView::getData(ValueType::size_type i) const
283 {
284 //
285 // don't do any index checking so allow one past the end of the
286 // vector providing the equivalent of end()
287 return (*m_data)[m_offset+i];
288 }
289
290 DataArrayView::ShapeType
291 DataArrayView::getResultSliceShape(const RegionType& region)
292 {
293 int dimSize;
294 RegionType::const_iterator i;
295 ShapeType result;
296 for (i=region.begin();i!=region.end();i++) {
297 dimSize=((i->second)-(i->first));
298 if (dimSize!=0) {
299 result.push_back(dimSize);
300 }
301 }
302 return result;
303 }
304
305 DataArrayView::RegionType
306 DataArrayView::getSliceRegion(const boost::python::object& key) const
307 {
308 int slice_rank, i;
309 int this_rank=getRank();
310 RegionType out(this_rank);
311 /* allow for case where key is singular eg: [1], this implies we
312 want to generate a rank-1 dimension object, as opposed to eg: [1,2]
313 which implies we want to take a rank dimensional object with one
314 dimension of size 1 */
315 extract<tuple> key_tuple(key);
316 if (key_tuple.check()) {
317 slice_rank=extract<int> (key.attr("__len__")());
318 /* ensure slice is correctly dimensioned */
319 if (slice_rank>this_rank) {
320 throw DataException("Error - rank of slices does not match rank of slicee");
321 } else {
322 /* calculate values for slice range */
323 for (i=0;i<slice_rank;i++) {
324 out[i]=getSliceRange(key[i],getShape()[i]);
325 }
326 }
327 } else {
328 slice_rank=1;
329 if (slice_rank>this_rank) {
330 throw DataException("Error - rank of slices does not match rank of slicee");
331 } else {
332 out[0]=getSliceRange(key,getShape()[0]);
333 }
334 }
335 for (i=slice_rank;i<this_rank;i++) {
336 out[i]=std::pair<int,int>(0,getShape()[i]);
337 }
338 return out;
339 }
340
341 std::pair<int,int>
342 getSliceRange(const boost::python::object& key,
343 const int shape)
344 {
345 /* default slice range is range of entire shape dimension */
346 int s0=0, s1=shape;;
347 extract<int> slice_int(key);
348 if (slice_int.check()) {
349 /* if the key is a single int set start=key and end=key */
350 /* in this case, we want to return a rank-1 dimension object from
351 this object, taken from a single index value for one of this
352 object's dimensions */
353 s0=slice_int();
354 s1=s0;
355 } else {
356 /* if key is a pair extract begin and end values */
357 extract<int> step(key.attr("step"));
358 if (step.check() && step()!=1) {
359 throw DataException("Error - Data does not support increments in slicing ");
360 } else {
361 extract<int> start(key.attr("start"));
362 if (start.check()) {
363 s0=start();
364 }
365 extract<int> stop(key.attr("stop"));
366 if (stop.check()) {
367 s1=stop();
368 }
369 }
370 }
371 if (s0 < 0)
372 throw DataException("Error - slice index out of range.");
373 if (s0 == s1 and s1 >= shape)
374 throw DataException("Error - slice index out of range.");
375 if (s0 != s1 and s1>shape)
376 throw DataException("Error - slice index out of range.");
377 if (s0 > s1)
378 throw DataException("Error - lower index must less or equal upper index.");
379 return std::pair<int,int>(s0,s1);
380 }
381
382 DataArrayView::RegionLoopRangeType
383 getSliceRegionLoopRange(const DataArrayView::RegionType& region)
384 {
385 DataArrayView::RegionLoopRangeType region_loop_range(region.size());
386 for (int i=0;i<region.size();i++) {
387 if (region[i].first==region[i].second) {
388 region_loop_range[i].first=region[i].first;
389 region_loop_range[i].second=region[i].second+1;
390 } else {
391 region_loop_range[i].first=region[i].first;
392 region_loop_range[i].second=region[i].second;
393 }
394 }
395 return region_loop_range;
396 }
397
398 void
399 DataArrayView::copySlice(const DataArrayView& other,
400 const RegionLoopRangeType& region)
401 {
402 //
403 // version of copySlice that uses the default offsets
404 copySlice(m_offset,other,other.m_offset,region);
405 }
406
407 void
408 DataArrayView::copySlice(ValueType::size_type thisOffset,
409 const DataArrayView& other,
410 ValueType::size_type otherOffset,
411 const RegionLoopRangeType& region)
412 {
413
414 //
415 // Make sure views are not empty
416
417 EsysAssert(!isEmpty(),
418 "Error - this view is empty.");
419 EsysAssert(!other.isEmpty(),
420 "Error - other view is empty.");
421
422 //
423 // Check the view to be sliced from is compatible with the region to be sliced,
424 // and that the region to be sliced is compatible with this view:
425
426 EsysAssert(checkOffset(thisOffset),
427 "Error - offset incompatible with this view.");
428 EsysAssert(otherOffset+noValues()<=other.getData().size(),
429 "Error - offset incompatible with other view.");
430
431 EsysAssert(other.getRank()==region.size(),
432 "Error - slice not same rank as view to be sliced from.");
433
434 EsysAssert(noValues()==noValues(getResultSliceShape(region)),
435 "Error - slice shape not compatible shape for this view.");
436
437 //
438 // copy the values in the specified region of the other view into this view
439
440 // the following loops cannot be parallelised due to the numCopy counter
441 int numCopy=0;
442
443 switch (region.size()) {
444 case 0:
445 /* this case should never be encountered,
446 as python will never pass us an empty region.
447 here for completeness only, allows slicing of a scalar */
448 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
449 numCopy++;
450 break;
451 case 1:
452 for (int i=region[0].first;i<region[0].second;i++) {
453 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];
454 numCopy++;
455 }
456 break;
457 case 2:
458 for (int j=region[1].first;j<region[1].second;j++) {
459 for (int i=region[0].first;i<region[0].second;i++) {
460 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];
461 numCopy++;
462 }
463 }
464 break;
465 case 3:
466 for (int k=region[2].first;k<region[2].second;k++) {
467 for (int j=region[1].first;j<region[1].second;j++) {
468 for (int i=region[0].first;i<region[0].second;i++) {
469 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
470 numCopy++;
471 }
472 }
473 }
474 break;
475 case 4:
476 for (int l=region[3].first;l<region[3].second;l++) {
477 for (int k=region[2].first;k<region[2].second;k++) {
478 for (int j=region[1].first;j<region[1].second;j++) {
479 for (int i=region[0].first;i<region[0].second;i++) {
480 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];
481 numCopy++;
482 }
483 }
484 }
485 }
486 break;
487 default:
488 stringstream mess;
489 mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
490 throw DataException(mess.str());
491 }
492 }
493
494 void
495 DataArrayView::copySliceFrom(const DataArrayView& other,
496 const RegionLoopRangeType& region_loop_range)
497 {
498 //
499 // version of copySliceFrom that uses the default offsets
500 copySliceFrom(m_offset,other,other.m_offset,region_loop_range);
501 }
502
503 void
504 DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
505 const DataArrayView& other,
506 ValueType::size_type otherOffset,
507 const RegionLoopRangeType& region)
508 {
509
510 //
511 // Make sure views are not empty
512
513 EsysAssert(!isEmpty(),
514 "Error - this view is empty.");
515 EsysAssert(!other.isEmpty(),
516 "Error - other view is empty.");
517
518 //
519 // Check this view is compatible with the region to be sliced,
520 // and that the region to be sliced is compatible with the other view:
521
522 EsysAssert(other.checkOffset(otherOffset),
523 "Error - offset incompatible with other view.");
524 EsysAssert(thisOffset+other.noValues()<=getData().size(),
525 "Error - offset incompatible with this view.");
526
527 EsysAssert(getRank()==region.size(),
528 "Error - slice not same rank as this view.");
529
530 EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),
531 "Error - slice shape not compatible shape for other view.");
532
533 //
534 // copy the values in the other view into the specified region of this view
535
536 // allow for case where other view is a scalar
537 if (other.getRank()==0) {
538
539 // the following loops cannot be parallelised due to the numCopy counter
540 int numCopy=0;
541
542 switch (region.size()) {
543 case 0:
544 /* this case should never be encountered,
545 as python will never pass us an empty region.
546 here for completeness only, allows slicing of a scalar */
547 (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
548 numCopy++;
549 break;
550 case 1:
551 for (int i=region[0].first;i<region[0].second;i++) {
552 (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];
553 numCopy++;
554 }
555 break;
556 case 2:
557 for (int j=region[1].first;j<region[1].second;j++) {
558 for (int i=region[0].first;i<region[0].second;i++) {
559 (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset];
560 numCopy++;
561 }
562 }
563 break;
564 case 3:
565 for (int k=region[2].first;k<region[2].second;k++) {
566 for (int j=region[1].first;j<region[1].second;j++) {
567 for (int i=region[0].first;i<region[0].second;i++) {
568 (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset];
569 numCopy++;
570 }
571 }
572 }
573 break;
574 case 4:
575 for (int l=region[3].first;l<region[3].second;l++) {
576 for (int k=region[2].first;k<region[2].second;k++) {
577 for (int j=region[1].first;j<region[1].second;j++) {
578 for (int i=region[0].first;i<region[0].second;i++) {
579 (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset];
580 numCopy++;
581 }
582 }
583 }
584 }
585 break;
586 default:
587 stringstream mess;
588 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
589 throw DataException(mess.str());
590 }
591
592 } else {
593
594 // the following loops cannot be parallelised due to the numCopy counter
595 int numCopy=0;
596
597 switch (region.size()) {
598 case 0:
599 /* this case should never be encountered,
600 as python will never pass us an empty region.
601 here for completeness only, allows slicing of a scalar */
602 (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
603 numCopy++;
604 break;
605 case 1:
606 for (int i=region[0].first;i<region[0].second;i++) {
607 (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset+numCopy];
608 numCopy++;
609 }
610 break;
611 case 2:
612 for (int j=region[1].first;j<region[1].second;j++) {
613 for (int i=region[0].first;i<region[0].second;i++) {
614 (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset+numCopy];
615 numCopy++;
616 }
617 }
618 break;
619 case 3:
620 for (int k=region[2].first;k<region[2].second;k++) {
621 for (int j=region[1].first;j<region[1].second;j++) {
622 for (int i=region[0].first;i<region[0].second;i++) {
623 (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset+numCopy];
624 numCopy++;
625 }
626 }
627 }
628 break;
629 case 4:
630 for (int l=region[3].first;l<region[3].second;l++) {
631 for (int k=region[2].first;k<region[2].second;k++) {
632 for (int j=region[1].first;j<region[1].second;j++) {
633 for (int i=region[0].first;i<region[0].second;i++) {
634 (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset+numCopy];
635 numCopy++;
636 }
637 }
638 }
639 }
640 break;
641 default:
642 stringstream mess;
643 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
644 throw DataException(mess.str());
645 }
646
647 }
648
649 }
650
651 string
652 DataArrayView::toString(const string& suffix) const
653 {
654 EsysAssert(!isEmpty(),"Error - View is empty.");
655 stringstream temp;
656 string finalSuffix=suffix;
657 if (suffix.length() > 0) {
658 finalSuffix+=" ";
659 }
660 switch (getRank()) {
661 case 0:
662 temp << finalSuffix << (*this)();
663 break;
664 case 1:
665 for (int i=0;i<getShape()[0];i++) {
666 temp << "(" << i << ") " << finalSuffix << (*this)(i);
667 if (i!=(getShape()[0]-1)) {
668 temp << endl;
669 }
670 }
671 break;
672 case 2:
673 for (int i=0;i<getShape()[0];i++) {
674 for (int j=0;j<getShape()[1];j++) {
675 temp << "(" << i << "," << j << ") " << finalSuffix << (*this)(i,j);
676 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
677 temp << endl;
678 }
679 }
680 }
681 break;
682 case 3:
683 for (int i=0;i<getShape()[0];i++) {
684 for (int j=0;j<getShape()[1];j++) {
685 for (int k=0;k<getShape()[2];k++) {
686 temp << "(" << i << "," << j << "," << k << ") " << finalSuffix << (*this)(i,j,k);
687 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1))) {
688 temp << endl;
689 }
690 }
691 }
692 }
693 break;
694 case 4:
695 for (int i=0;i<getShape()[0];i++) {
696 for (int j=0;j<getShape()[1];j++) {
697 for (int k=0;k<getShape()[2];k++) {
698 for (int l=0;l<getShape()[3];l++) {
699 temp << "(" << i << "," << j << "," << k << "," << l << ") " << finalSuffix << (*this)(i,j,k,l);
700 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1) && l==(getShape()[3]-1))) {
701 temp << endl;
702 }
703 }
704 }
705 }
706 }
707 break;
708 default:
709 stringstream mess;
710 mess << "Error - (toString) Invalid rank: " << getRank();
711 throw DataException(mess.str());
712 }
713 return temp.str();
714 }
715
716 string
717 DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)
718 {
719 stringstream temp;
720 temp << "(";
721 for (int i=0;i<shape.size();i++) {
722 temp << shape[i];
723 if (i < shape.size()-1) {
724 temp << ",";
725 }
726 }
727 temp << ")";
728 return temp.str();
729 }
730
731 void
732 DataArrayView::eigenvalues(const DataArrayView& in,
733 DataArrayView& ev)
734 {
735 EsysAssert(!in.isEmpty(),
736 "Error - in view is empty.");
737 EsysAssert(!ev.isEmpty(),
738 "Error - ev view is empty.");
739 EsysAssert(in.getRank()==2,
740 "Error - in has not rank 2.");
741 EsysAssert(in.getShape()[0]==in.getShape()[1],
742 "Error - in is not square");
743 EsysAssert(ev.getRank()==1,
744 "Error - ev has not rank 1.");
745 EsysAssert(ev.getShape()[0]==in.getShape()[0],
746 "Error - ev is too short");
747 EsysAssert(0<ev.getShape()[0] && ev.getShape()[0]<4,
748 "Error - dimension is restricted to 3.");
749 double ev0,ev1,ev2;
750 int s=in.getShape()[0];
751 if (s==1) {
752 eigenvalues1(in(0,0),&ev0);
753 ev(0)=ev0;
754
755 } else if (s==2) {
756 eigenvalues2(in(0,0),(in(0,1)+in(1,0))/2.,in(1,1),
757 &ev0,&ev1);
758 ev(0)=ev0;
759 ev(1)=ev1;
760
761 } else if (s==3) {
762 eigenvalues3(in(0,0),(in(0,1)+in(1,0))/2.,(in(0,2)+in(2,0))/2.,in(1,1),(in(2,1)+in(1,2))/2.,in(2,2),
763 &ev0,&ev1,&ev2);
764 ev(0)=ev0;
765 ev(1)=ev1;
766 ev(2)=ev2;
767
768 }
769 }
770
771 void
772 DataArrayView::matMult(const DataArrayView& left,
773 const DataArrayView& right,
774 DataArrayView& result)
775 {
776
777 if (left.getRank()==0 || right.getRank()==0) {
778 stringstream temp;
779 temp << "Error - (matMult) Invalid for rank 0 objects.";
780 throw DataException(temp.str());
781 }
782
783 if (left.getShape()[left.getRank()-1] != right.getShape()[0]) {
784 stringstream temp;
785 temp << "Error - (matMult) Dimension: " << left.getRank()
786 << ", size: " << left.getShape()[left.getRank()-1]
787 << " of LHS and dimension: 1, size: " << right.getShape()[0]
788 << " of RHS don't match.";
789 throw DataException(temp.str());
790 }
791
792 int outputRank = left.getRank()+right.getRank()-2;
793
794 if (outputRank < 0) {
795 stringstream temp;
796 temp << "Error - (matMult) LHS and RHS cannot be multiplied "
797 << "as they have incompatable rank.";
798 throw DataException(temp.str());
799 }
800
801 if (outputRank != result.getRank()) {
802 stringstream temp;
803 temp << "Error - (matMult) Rank of result array is: "
804 << result.getRank()
805 << " it must be: " << outputRank;
806 throw DataException(temp.str());
807 }
808
809 for (int i=0; i<(left.getRank()-1); i++) {
810 if (left.getShape()[i] != result.getShape()[i]) {
811 stringstream temp;
812 temp << "Error - (matMult) Dimension: " << i
813 << " of LHS and result array don't match.";
814 throw DataException(temp.str());
815 }
816 }
817
818 for (int i=1; i<right.getRank(); i++) {
819 if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
820 stringstream temp;
821 temp << "Error - (matMult) Dimension: " << i
822 << ", size: " << right.getShape()[i]
823 << " of RHS and dimension: " << i+left.getRank()-1
824 << ", size: " << result.getShape()[i+left.getRank()-1]
825 << " of result array don't match.";
826 throw DataException(temp.str());
827 }
828 }
829
830 switch (left.getRank()) {
831
832 case 1:
833 switch (right.getRank()) {
834 case 1:
835 result()=0;
836 for (int i=0;i<left.getShape()[0];i++) {
837 result()+=left(i)*right(i);
838 }
839 break;
840 case 2:
841 for (int i=0;i<result.getShape()[0];i++) {
842 result(i)=0;
843 for (int j=0;j<right.getShape()[0];j++) {
844 result(i)+=left(j)*right(j,i);
845 }
846 }
847 break;
848 default:
849 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
850 throw DataException(temp.str());
851 break;
852 }
853 break;
854
855 case 2:
856 switch (right.getRank()) {
857 case 1:
858 result()=0;
859 for (int i=0;i<left.getShape()[0];i++) {
860 result(i)=0;
861 for (int j=0;j<left.getShape()[1];j++) {
862 result(i)+=left(i,j)*right(i);
863 }
864 }
865 break;
866 case 2:
867 for (int i=0;i<result.getShape()[0];i++) {
868 for (int j=0;j<result.getShape()[1];j++) {
869 result(i,j)=0;
870 for (int jR=0;jR<right.getShape()[0];jR++) {
871 result(i,j)+=left(i,jR)*right(jR,j);
872 }
873 }
874 }
875 break;
876 default:
877 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
878 throw DataException(temp.str());
879 break;
880 }
881 break;
882
883 default:
884 stringstream temp; temp << "Error - (matMult) Not supported for rank: " << left.getRank();
885 throw DataException(temp.str());
886 break;
887 }
888
889 }
890
891 DataArrayView::ShapeType
892 DataArrayView::determineResultShape(const DataArrayView& left,
893 const DataArrayView& right)
894 {
895 ShapeType temp;
896 for (int i=0; i<(left.getRank()-1); i++) {
897 temp.push_back(left.getShape()[i]);
898 }
899 for (int i=1; i<right.getRank(); i++) {
900 temp.push_back(right.getShape()[i]);
901 }
902 return temp;
903 }
904
905 bool operator==(const DataArrayView& left, const DataArrayView& right)
906 {
907 //
908 // The views are considered equal if they have the same shape and each value
909 // of the view is the same.
910 bool result=(left.m_shape==right.m_shape);
911 if (result) {
912 for (int i=0;i<left.noValues();i++) {
913 result=result && (left.getData(i)==right.getData(i));
914 }
915 }
916 return result;
917 }
918
919 bool operator!=(const DataArrayView& left, const DataArrayView& right)
920 {
921 return (!operator==(left,right));
922 }
923
924 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26