7041ed378573409e95cb08e8ba517e25a4006720
[blender.git] / intern / smoke / intern / tnt / tnt_array2d_utils.h
1 /*
2 *
3 * Template Numerical Toolkit (TNT)
4 *
5 * Mathematical and Computational Sciences Division
6 * National Institute of Technology,
7 * Gaithersburg, MD USA
8 *
9 *
10 * This software was developed at the National Institute of Standards and
11 * Technology (NIST) by employees of the Federal Government in the course
12 * of their official duties. Pursuant to title 17 Section 105 of the
13 * United States Code, this software is not subject to copyright protection
14 * and is in the public domain. NIST assumes no responsibility whatsoever for
15 * its use by other parties, and makes no guarantees, expressed or implied,
16 * about its quality, reliability, or any other characteristic.
17 *
18 */
19
20
21 #ifndef TNT_ARRAY2D_UTILS_H
22 #define TNT_ARRAY2D_UTILS_H
23
24 #include <cstdlib>
25 #include <cassert>
26
27 namespace TNT
28 {
29
30
31 template <class T>
32 std::ostream& operator<<(std::ostream &s, const Array2D<T> &A)
33 {
34     int M=A.dim1();
35     int N=A.dim2();
36
37     s << M << " " << N << "\n";
38
39     for (int i=0; i<M; i++)
40     {
41         for (int j=0; j<N; j++)
42         {
43             s << A[i][j] << " ";
44         }
45         s << "\n";
46     }
47
48
49     return s;
50 }
51
52 template <class T>
53 std::istream& operator>>(std::istream &s, Array2D<T> &A)
54 {
55
56     int M, N;
57
58     s >> M >> N;
59
60         Array2D<T> B(M,N);
61
62     for (int i=0; i<M; i++)
63         for (int j=0; j<N; j++)
64         {
65             s >>  B[i][j];
66         }
67
68         A = B;
69     return s;
70 }
71
72
73 template <class T>
74 Array2D<T> operator+(const Array2D<T> &A, const Array2D<T> &B)
75 {
76         int m = A.dim1();
77         int n = A.dim2();
78
79         if (B.dim1() != m ||  B.dim2() != n )
80                 return Array2D<T>();
81
82         else
83         {
84                 Array2D<T> C(m,n);
85
86                 for (int i=0; i<m; i++)
87                 {
88                         for (int j=0; j<n; j++)
89                                 C[i][j] = A[i][j] + B[i][j];
90                 }
91                 return C;
92         }
93 }
94
95 template <class T>
96 Array2D<T> operator-(const Array2D<T> &A, const Array2D<T> &B)
97 {
98         int m = A.dim1();
99         int n = A.dim2();
100
101         if (B.dim1() != m ||  B.dim2() != n )
102                 return Array2D<T>();
103
104         else
105         {
106                 Array2D<T> C(m,n);
107
108                 for (int i=0; i<m; i++)
109                 {
110                         for (int j=0; j<n; j++)
111                                 C[i][j] = A[i][j] - B[i][j];
112                 }
113                 return C;
114         }
115 }
116
117
118 template <class T>
119 Array2D<T> operator*(const Array2D<T> &A, const Array2D<T> &B)
120 {
121         int m = A.dim1();
122         int n = A.dim2();
123
124         if (B.dim1() != m ||  B.dim2() != n )
125                 return Array2D<T>();
126
127         else
128         {
129                 Array2D<T> C(m,n);
130
131                 for (int i=0; i<m; i++)
132                 {
133                         for (int j=0; j<n; j++)
134                                 C[i][j] = A[i][j] * B[i][j];
135                 }
136                 return C;
137         }
138 }
139
140
141
142
143 template <class T>
144 Array2D<T> operator/(const Array2D<T> &A, const Array2D<T> &B)
145 {
146         int m = A.dim1();
147         int n = A.dim2();
148
149         if (B.dim1() != m ||  B.dim2() != n )
150                 return Array2D<T>();
151
152         else
153         {
154                 Array2D<T> C(m,n);
155
156                 for (int i=0; i<m; i++)
157                 {
158                         for (int j=0; j<n; j++)
159                                 C[i][j] = A[i][j] / B[i][j];
160                 }
161                 return C;
162         }
163 }
164
165
166
167
168
169 template <class T>
170 Array2D<T>&  operator+=(Array2D<T> &A, const Array2D<T> &B)
171 {
172         int m = A.dim1();
173         int n = A.dim2();
174
175         if (B.dim1() == m ||  B.dim2() == n )
176         {
177                 for (int i=0; i<m; i++)
178                 {
179                         for (int j=0; j<n; j++)
180                                 A[i][j] += B[i][j];
181                 }
182         }
183         return A;
184 }
185
186
187
188 template <class T>
189 Array2D<T>&  operator-=(Array2D<T> &A, const Array2D<T> &B)
190 {
191         int m = A.dim1();
192         int n = A.dim2();
193
194         if (B.dim1() == m ||  B.dim2() == n )
195         {
196                 for (int i=0; i<m; i++)
197                 {
198                         for (int j=0; j<n; j++)
199                                 A[i][j] -= B[i][j];
200                 }
201         }
202         return A;
203 }
204
205
206
207 template <class T>
208 Array2D<T>&  operator*=(Array2D<T> &A, const Array2D<T> &B)
209 {
210         int m = A.dim1();
211         int n = A.dim2();
212
213         if (B.dim1() == m ||  B.dim2() == n )
214         {
215                 for (int i=0; i<m; i++)
216                 {
217                         for (int j=0; j<n; j++)
218                                 A[i][j] *= B[i][j];
219                 }
220         }
221         return A;
222 }
223
224
225
226
227
228 template <class T>
229 Array2D<T>&  operator/=(Array2D<T> &A, const Array2D<T> &B)
230 {
231         int m = A.dim1();
232         int n = A.dim2();
233
234         if (B.dim1() == m ||  B.dim2() == n )
235         {
236                 for (int i=0; i<m; i++)
237                 {
238                         for (int j=0; j<n; j++)
239                                 A[i][j] /= B[i][j];
240                 }
241         }
242         return A;
243 }
244
245 /**
246     Matrix Multiply:  compute C = A*B, where C[i][j]
247     is the dot-product of row i of A and column j of B.
248
249
250     @param A an (m x n) array
251     @param B an (n x k) array
252     @return the (m x k) array A*B, or a null array (0x0)
253         if the matrices are non-conformant (i.e. the number
254         of columns of A are different than the number of rows of B.)
255
256
257 */
258 template <class T>
259 Array2D<T> matmult(const Array2D<T> &A, const Array2D<T> &B)
260 {
261     if (A.dim2() != B.dim1())
262         return Array2D<T>();
263
264     int M = A.dim1();
265     int N = A.dim2();
266     int K = B.dim2();
267
268     Array2D<T> C(M,K);
269
270     for (int i=0; i<M; i++)
271         for (int j=0; j<K; j++)
272         {
273             T sum = 0;
274
275             for (int k=0; k<N; k++)
276                 sum += A[i][k] * B [k][j];
277
278             C[i][j] = sum;
279         }
280
281     return C;
282
283 }
284
285 } // namespace TNT
286
287 #endif