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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1781 - (show annotations)
Thu Sep 11 05:03:14 2008 UTC (10 years, 11 months ago) by jfenwick
File size: 5423 byte(s)
Branch commit

Merged changes from trunk version 1695 up to and including version 1779.


1
2 /* $Id: DataMaths.cpp 1697 2008-08-11 06:29:54Z jfenwick $ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #include "DataTypes.h"
17 #include "DataMaths.h"
18 #include <sstream>
19
20 namespace escript
21 {
22 namespace DataMaths
23 {
24
25 void
26 matMult(const DataTypes::ValueType& left,
27 const DataTypes::ShapeType& leftShape,
28 DataTypes::ValueType::size_type leftOffset,
29 const DataTypes::ValueType& right,
30 const DataTypes::ShapeType& rightShape,
31 DataTypes::ValueType::size_type rightOffset,
32 DataTypes::ValueType& result,
33 const DataTypes::ShapeType& resultShape)
34 {
35 using namespace escript::DataTypes;
36 using namespace std;
37
38 int leftRank=getRank(leftShape);
39 int rightRank=getRank(rightShape);
40 int resultRank=getRank(resultShape);
41 if (leftRank==0 || rightRank==0) {
42 stringstream temp;
43 temp << "Error - (matMult) Invalid for rank 0 objects.";
44 throw DataException(temp.str());
45 }
46
47 if (leftShape[leftRank-1] != rightShape[0]) {
48 stringstream temp;
49 temp << "Error - (matMult) Dimension: " << leftRank
50 << ", size: " << leftShape[leftRank-1]
51 << " of LHS and dimension: 1, size: " << rightShape[0]
52 << " of RHS don't match.";
53 throw DataException(temp.str());
54 }
55
56 int outputRank = leftRank+rightRank-2;
57
58 if (outputRank < 0) {
59 stringstream temp;
60 temp << "Error - (matMult) LHS and RHS cannot be multiplied "
61 << "as they have incompatable rank.";
62 throw DataException(temp.str());
63 }
64
65 if (outputRank != resultRank) {
66 stringstream temp;
67 temp << "Error - (matMult) Rank of result array is: "
68 << resultRank
69 << " it must be: " << outputRank;
70 throw DataException(temp.str());
71 }
72
73 for (int i=0; i<(leftRank-1); i++) {
74 if (leftShape[i] != resultShape[i]) {
75 stringstream temp;
76 temp << "Error - (matMult) Dimension: " << i
77 << " of LHS and result array don't match.";
78 throw DataException(temp.str());
79 }
80 }
81
82 for (int i=1; i<rightRank; i++) {
83 if (rightShape[i] != resultShape[i+leftRank-2]) {
84 stringstream temp;
85 temp << "Error - (matMult) Dimension: " << i
86 << ", size: " << rightShape[i]
87 << " of RHS and dimension: " << i+leftRank-1
88 << ", size: " << resultShape[i+leftRank-1]
89 << " of result array don't match.";
90 throw DataException(temp.str());
91 }
92 }
93
94 switch (leftRank) {
95
96 case 1:
97 switch (rightRank) {
98 case 1:
99 result[0]=0;
100 for (int i=0;i<leftShape[0];i++) {
101 result[0]+=left[i+leftOffset]*right[i+rightOffset];
102 }
103 break;
104 case 2:
105 for (int i=0;i<resultShape[0];i++) {
106 result[i]=0;
107 for (int j=0;j<rightShape[0];j++) {
108 result[i]+=left[j+leftOffset]*right[getRelIndex(rightShape,j,i)+rightOffset];
109 }
110 }
111 break;
112 default:
113 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
114 throw DataException(temp.str());
115 break;
116 }
117 break;
118
119 case 2:
120 switch (rightRank) {
121 case 1:
122 result[0]=0;
123 for (int i=0;i<leftShape[0];i++) {
124 result[i]=0;
125 for (int j=0;j<leftShape[1];j++) {
126 result[i]+=left[leftOffset+getRelIndex(leftShape,i,j)]*right[i+rightOffset];
127 }
128 }
129 break;
130 case 2:
131 for (int i=0;i<resultShape[0];i++) {
132 for (int j=0;j<resultShape[1];j++) {
133 result[getRelIndex(resultShape,i,j)]=0;
134 for (int jR=0;jR<rightShape[0];jR++) {
135 result[getRelIndex(resultShape,i,j)]+=left[leftOffset+getRelIndex(leftShape,i,jR)]*right[rightOffset+getRelIndex(rightShape,jR,j)];
136 }
137 }
138 }
139 break;
140 default:
141 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
142 throw DataException(temp.str());
143 break;
144 }
145 break;
146
147 default:
148 stringstream temp; temp << "Error - (matMult) Not supported for rank: " << leftRank;
149 throw DataException(temp.str());
150 break;
151 }
152
153 }
154
155
156 DataTypes::ShapeType
157 determineResultShape(const DataTypes::ShapeType& left,
158 const DataTypes::ShapeType& right)
159 {
160 DataTypes::ShapeType result;
161 for (int i=0; i<(DataTypes::getRank(left)-1); i++) {
162 result.push_back(left[i]);
163 }
164 for (int i=1; i<DataTypes::getRank(right); i++) {
165 result.push_back(right[i]);
166 }
167 return result;
168 }
169
170
171 } // end namespace
172 } // end namespace
173

  ViewVC Help
Powered by ViewVC 1.1.26