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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26