1 #ifndef TNTSUPP_H
2 #define TNTSUPP_H
3
4 #include <iostream>
5 #include "tnt/subscript.h"
6 #include "tnt/tnt.h"
7 #include "tnt/vec.h"
8 #include "tnt/fmat.h"
9 #include "tnt/region1d.h"
10 #include "tnt/region2d.h"
11 #include "tnt/transv.h"
12 #include "tnt/lu.h"
13
14 //using namespace std;
15 //using namespace TNT;
16
17
18 namespace TNT
19 {
20
21 typedef Vector<double> DVector;
22 typedef Vector<int> IVector;
23 typedef Fortran_Matrix<double> DMatrix;
24
25 template <class T>
min(T a,T b)26 inline T min(T a, T b)
27 {
28 return (a < b ? a : b);
29 }
30
31 template <class T>
max(T a,T b)32 inline T max(T a, T b)
33 {
34 return (a > b ? a : b);
35 }
36
37 //typedef Region1D<DVector> subDVector;
38 //typedef Region2D<DMatrix> subDMatrix;
39
40 //matrix operation on Region2D;
41
42 template <class T>
43 Fortran_Matrix<T> operator+(const_Region2D<Fortran_Matrix<T> > &A,
44 const_Region2D<Fortran_Matrix<T> > &B) {
45 Subscript M = A.num_rows();
46 Subscript N = A.num_cols();
47
48 assert(M==B.num_rows());
49 assert(N==B.num_cols());
50
51 Fortran_Matrix<T> ans(M, N);
52 for (Subscript i = 1; i <= M; i++)
53 for (Subscript j = 1; j <= N; j++)
54 ans(i, j) = A(i, j) + B(i, j);
55 return ans;
56 }
57
58 template <class T>
59 Fortran_Matrix<T> operator+(const Region2D<Fortran_Matrix<T> > &A,
60 const Region2D<Fortran_Matrix<T> > &B) {
61 Subscript M = A.num_rows();
62 Subscript N = A.num_cols();
63
64 assert(M==B.num_rows());
65 assert(N==B.num_cols());
66
67 Fortran_Matrix<T> ans(M, N);
68 for (Subscript i = 1; i <= M; i++)
69 for (Subscript j = 1; j <= N; j++)
70 ans(i, j) = A(i, j) + B(i, j);
71 return ans;
72 }
73
74 template <class T>
75 Fortran_Matrix<T> operator-(const_Region2D<Fortran_Matrix<T> > &A,
76 const_Region2D<Fortran_Matrix<T> > &B) {
77 Subscript M = A.num_rows();
78 Subscript N = A.num_cols();
79
80 assert(M==B.num_rows());
81 assert(N==B.num_cols());
82
83 Fortran_Matrix<T> ans(M, N);
84 for (Subscript i = 1; i <= M; i++)
85 for (Subscript j = 1; j <= N; j++)
86 ans(i, j) = A(i, j) - B(i, j);
87 return ans;
88 }
89 template <class T>
90 Fortran_Matrix<T> operator-(const Region2D<Fortran_Matrix<T> > &A,
91 const Region2D<Fortran_Matrix<T> > &B) {
92 Subscript M = A.num_rows();
93 Subscript N = A.num_cols();
94
95 assert(M==B.num_rows());
96 assert(N==B.num_cols());
97
98 Fortran_Matrix<T> ans(M, N);
99 for (Subscript i = 1; i <= M; i++)
100 for (Subscript j = 1; j <= N; j++)
101 ans(i, j) = A(i, j) - B(i, j);
102 return ans;
103 }
104
105 template <class T>
mult_element(const Region2D<Fortran_Matrix<T>> & A,const Region2D<Fortran_Matrix<T>> & B)106 Fortran_Matrix<T> mult_element(const Region2D<Fortran_Matrix<T> > &A,
107 const Region2D<Fortran_Matrix<T> > &B)
108 {
109 Subscript M = A.num_rows();
110 Subscript N = A.num_cols();
111
112 assert(M==B.num_rows());
113 assert(N==B.num_cols());
114
115 Fortran_Matrix<T> tmp(M,N);
116 Subscript i,j;
117
118 for (i=1; i<=M; i++)
119 for (j=1; j<=N; j++)
120 tmp(i,j) = A(i,j) * B(i,j);
121
122 return tmp;
123 }
124
125
126 template <class T>
transpose(const Region2D<Fortran_Matrix<T>> & A)127 Fortran_Matrix<T> transpose(const Region2D<Fortran_Matrix<T> > &A)
128 {
129 Subscript M = A.num_rows();
130 Subscript N = A.num_cols();
131
132 Fortran_Matrix<T> S(N,M);
133 Subscript i, j;
134
135 for (i=1; i<=M; i++)
136 for (j=1; j<=N; j++)
137 S(j,i) = A(i,j);
138
139 return S;
140 }
141
142
143 template <class T>
matmult(const Region2D<Fortran_Matrix<T>> & A,const Region2D<Fortran_Matrix<T>> & B)144 inline Fortran_Matrix<T> matmult(const Region2D<Fortran_Matrix<T> > &A,
145 const Region2D<Fortran_Matrix<T> > &B)
146 {
147
148 #ifdef TNT_BOUNDS_CHECK
149 assert(A.num_cols() == B.num_rows());
150 #endif
151
152 Subscript M = A.num_rows();
153 Subscript N = A.num_cols();
154 Subscript K = B.num_cols();
155
156 Fortran_Matrix<T> tmp(M,K);
157 T sum;
158
159 for (Subscript i=1; i<=M; i++)
160 for (Subscript k=1; k<=K; k++)
161 {
162 sum = 0;
163 for (Subscript j=1; j<=N; j++)
164 sum = sum + A(i,j) * B(j,k);
165
166 tmp(i,k) = sum;
167 }
168
169 return tmp;
170 }
171
172 template <class T>
173 inline Fortran_Matrix<T> operator*(const Region2D<Fortran_Matrix<T> > &A,
174 const Region2D<Fortran_Matrix<T> > &B)
175 {
176 return matmult(A,B);
177 }
178
179 template <class T>
matmult(Fortran_Matrix<T> & C,const Region2D<Fortran_Matrix<T>> & A,const Region2D<Fortran_Matrix<T>> & B)180 inline int matmult(Fortran_Matrix<T>& C, const Region2D<Fortran_Matrix<T> >&A,
181 const Region2D<Fortran_Matrix<T> >&B)
182 {
183
184 assert(A.num_cols() == B.num_rows());
185
186 Subscript M = A.num_rows();
187 Subscript N = A.num_cols();
188 Subscript K = B.num_cols();
189
190 C.newsize(M,K); // adjust shape of C, if necessary
191
192
193 T sum;
194
195 const T* row_i;
196 const T* col_k;
197
198 for (Subscript i=1; i<=M; i++)
199 {
200 for (Subscript k=1; k<=K; k++)
201 {
202 row_i = &A(i,1);
203 col_k = &B(1,k);
204 sum = 0;
205 for (Subscript j=1; j<=N; j++)
206 {
207 sum += *row_i * *col_k;
208 row_i += M;
209 col_k ++;
210 }
211
212 C(i,k) = sum;
213 }
214
215 }
216
217 return 0;
218 }
219
220
221 template <class T>
matmult(const Region2D<Fortran_Matrix<T>> & A,const Vector<T> & x)222 Vector<T> matmult(const Region2D<Fortran_Matrix<T> > &A, const Vector<T> &x)
223 {
224
225 #ifdef TNT_BOUNDS_CHECK
226 assert(A.num_cols() == x.dim());
227 #endif
228
229 Subscript M = A.num_rows();
230 Subscript N = A.num_cols();
231
232 Vector<T> tmp(M);
233 T sum;
234
235 for (Subscript i=1; i<=M; i++)
236 {
237 sum = 0;
238 for (Subscript j=1; j<=N; j++)
239 sum = sum + A(i,j) * x(j);
240
241 tmp(i) = sum;
242 }
243
244 return tmp;
245 }
246
247 template <class T>
248 inline Vector<T> operator*(const Region2D<Fortran_Matrix<T> > &A, const Vector<T> &x)
249 {
250 return matmult(A,x);
251 }
252
253 template <class T>
254 inline Fortran_Matrix<T> operator*(const Region2D<Fortran_Matrix<T> > &A,
255 const T &x)
256 {
257 Subscript M = A.num_rows();
258 Subscript N = A.num_cols();
259
260 Subscript MN = M*N;
261
262 Fortran_Matrix<T> res(M,N);
263 const T* a = A.begin();
264 T* t = res.begin();
265 T* tend = res.end();
266
267 for (t=res.begin(); t < tend; t++, a++)
268 *t = *a * x;
269
270 return res;
271 }
272
273
274 //convert Region2D to matrix or vector
275 template <class T>
asMat(const Region2D<Fortran_Matrix<T>> & A)276 Fortran_Matrix<T> asMat(const Region2D<Fortran_Matrix<T> > &A) {
277 Subscript m = A.num_rows(), n = A.num_cols();
278 Fortran_Matrix<T> ans(m, n);
279 for (Subscript i = 1; i <= m; i++)
280 for (Subscript j = 1; j <= n; j++)
281 ans(i, j) = A(i, j);
282 return ans;
283 }
284
285 template <class T>
asMat(const_Region2D<Fortran_Matrix<T>> & A)286 Fortran_Matrix<T> asMat(const_Region2D<Fortran_Matrix<T> > &A) {
287 Subscript m = A.num_rows(), n = A.num_cols();
288 Fortran_Matrix<T> ans(m, n);
289 for (Subscript i = 1; i <= m; i++)
290 for (Subscript j = 1; j <= n; j++)
291 ans(i, j) = A(i, j);
292 return ans;
293 }
294
295 template <class T>
asVec(const Region2D<Fortran_Matrix<T>> & A)296 Vector<T> asVec(const Region2D<Fortran_Matrix<T> > &A) {
297 // A is 1 row or 1 col
298 Subscript m = A.num_rows(), n = A.num_cols();
299 if (m == 1) {
300 Vector<T> ans(n);
301 for (Subscript i = 1; i <= n; i++) ans(i) = A(1,i);
302 return ans;
303 }
304 else {
305 Vector<T> ans(m);
306 for (Subscript i = 1; i <= m; i++) ans(i) = A(i,1);
307 return ans;
308 }
309 }
310
311 template <class T>
asVec(const_Region2D<Fortran_Matrix<T>> A)312 Vector<T> asVec(const_Region2D<Fortran_Matrix<T> > A) {
313 // A is 1 row or 1 col
314 Subscript m = A.num_rows(), n = A.num_cols();
315 if (m == 1) {
316 Vector<T> ans(n);
317 for (Subscript i = 1; i <= n; i++) ans(i) = A(1,i);
318 return ans;
319 }
320 else {
321 Vector<T> ans(m);
322 for (Subscript i = 1; i <= m; i++) ans(i) = A(i,1);
323 return ans;
324 }
325 }
326
327 //convert vector to matrix
328 template <class T>
asRowMat(const Vector<T> & v)329 Fortran_Matrix<T> asRowMat(const Vector<T> &v) {
330 Subscript n = v.size();
331 Fortran_Matrix<T> ans(1,n);
332 for (Subscript i = 1; i <= n; i++) ans(1,i) = v(i);
333 return ans;
334 }
335
336 template <class T>
asColMat(const Vector<T> & v)337 Fortran_Matrix<T> asColMat(const Vector<T> &v) {
338 Subscript n = v.size();
339 Fortran_Matrix<T> ans(n,1);
340 for (Subscript i = 1; i <= n; i++) ans(i,1) = v(i);
341 return ans;
342 }
343
344
345 //scalar multiplication
346 template <class T>
347 inline Vector<T> operator*(const Vector<T> &v, const T &x) {
348 Subscript m = v.size();
349 Vector<T> ans(m);
350 for (Subscript i = 1; i <= m; i++)
351 ans(i) = v(i) * x;
352 return ans;
353 }
354
355 template <class T>
356 inline Vector<T> operator*(const T &x, const Vector<T> &v) {
357 return v * x;
358 }
359
360 template <class T>
361 inline Fortran_Matrix<T> operator*(const T &x, const Fortran_Matrix<T> &A) {
362 return A * x;
363 }
364
365 // utilities:
366
367 template <class T>
MatRow(Fortran_Matrix<T> & x,int i)368 Region2D<Fortran_Matrix<T> > MatRow(Fortran_Matrix<T> &x, int i) {
369 int n = x.num_cols();
370 Index1D I(i,i), J(1,n);
371 return x(I,J);
372 }
373
374 template <class T>
MatCol(Fortran_Matrix<T> & x,int i)375 Region2D<Fortran_Matrix<T> > MatCol(Fortran_Matrix<T> &x, int i) {
376 int m = x.num_rows();
377 Index1D I(1,m), J(i,i);
378 return x(I,J);
379 }
380
381 template <class T>
MatRows(Fortran_Matrix<T> & x,const Index1D & I)382 Region2D<Fortran_Matrix<T> > MatRows(Fortran_Matrix<T> &x, const Index1D &I) {
383 int n = x.num_cols(); Index1D J(1,n);
384 return x(I,J);
385 }
386
387 template <class T>
MatCols(Fortran_Matrix<T> & x,const Index1D & J)388 Region2D<Fortran_Matrix<T> > MatCols(Fortran_Matrix<T> &x, const Index1D &J) {
389 int m = x.num_rows(); Index1D I(1,m);
390 return x(I,J);
391 }
392
393 template <class T>
VecSubs(Vector<T> & x,const Index1D & I)394 Region1D<Vector<T> > VecSubs(Vector<T> &x, const Index1D &I) {
395 return Region1D<Vector<T> >(x, I);
396 }
397
398 template <class T>
asVec(const Region1D<Vector<T>> & x)399 Vector<T> asVec(const Region1D<Vector<T> > &x) {
400 Vector<T> ans(x.dim());
401 for (int i = 1; i <= ans.size(); i++) ans(i) = x(i);
402 return ans;
403 }
404
405 // transp(A) * inv(B) * C
406 template <class Matrix, class T>
matmult(const Transpose_View<Matrix> & A,const Fortran_Matrix<T> & B)407 Fortran_Matrix<T> matmult(
408 const Transpose_View<Matrix> & A,
409 const Fortran_Matrix<T> &B)
410 {
411 Subscript M = A.num_rows();
412 Subscript N = A.num_cols();
413
414 assert(B.num_rows() == N);
415 Subscript L = B.num_cols();
416
417 Fortran_Matrix<T> x(M,L);
418 Subscript i, j, k;
419 T tmp = 0;
420
421 for (i=1; i<=M; i++) {
422 for (j=1; j<=L; j++) {
423 tmp = 0;
424 for (k = 1; k <= N; k++) tmp += A(i,k) * B(k,j);
425 x(i,j) = tmp;
426 }
427 }
428
429 return x;
430 }
431
432 template <class Matrix, class T>
433 inline Fortran_Matrix<T> operator*(const Transpose_View<Matrix> & A, const Fortran_Matrix<T> &B)
434 {
435 return matmult(A,B);
436 }
437
438 //crossprod
439 template <class T>
outerprod(const Vector<T> & v)440 Fortran_Matrix<T> outerprod(const Vector<T> &v) {
441 int n = v.size();
442 Fortran_Matrix<T> ans(n,n);
443 for (int i = 1; i <= n; i++)
444 for (int j = 1; j <= n; j++)
445 ans(i,j) = v(i) * v(j);
446 return ans;
447 }
448
449 template <class T>
outerprod(const Vector<T> & v1,const Vector<T> & v2)450 Fortran_Matrix<T> outerprod(const Vector<T> &v1, const Vector<T> &v2) {
451 int m = v1.size(), n = v2.size();
452 Fortran_Matrix<T> ans(m,n);
453 for (int i = 1; i <= m; i++)
454 for (int j = 1; j <= n; j++)
455 ans(i,j) = v1(i) * v2(j);
456 return ans;
457 }
458
459 template <class T>
fmax(const Vector<T> & v)460 T fmax(const Vector<T> & v) {
461 T ans = v(1);
462 for (int i = 1; i <= v.dim(); i++)
463 if (ans < v(i)) ans = v(i);
464 return ans;
465 }
466
467 template <class T>
fmin(const Vector<T> & v)468 T fmin(const Vector<T> & v) {
469 T ans = v(1);
470 for (int i = 1; i <= v.dim(); i++)
471 if (ans > v(i)) ans = v(i);
472 return ans;
473 }
474
475 template <class T>
fmax(const Fortran_Matrix<T> & m)476 T fmax(const Fortran_Matrix<T> &m) {
477 T ans = m(1, 1);
478 for (int i = 1; i <= m.dim(1); i++)
479 for (int j = 1; j <= m.dim(2); j++)
480 if (ans < m(i, j)) ans = m(i, j);
481 return ans;
482 }
483
484 template <class T>
fmin(const Fortran_Matrix<T> & m)485 T fmin(const Fortran_Matrix<T> &m) {
486 T ans = m(1, 1);
487 for (int i = 1; i <= m.dim(1); i++)
488 for (int j = 1; j <= m.dim(2); j++)
489 if (ans > m(i, j)) ans = m(i, j);
490 return ans;
491 }
492 template <class T>
sum(const Vector<T> & v)493 T sum(const Vector<T> &v) {
494 T ans = 0;
495 for (int i = 1; i <= v.dim(); i++)
496 ans += v(i);
497 return ans;
498 }
499
500 template <class T>
sum(const Fortran_Matrix<T> & m)501 T sum(const Fortran_Matrix<T> &m) {
502 T ans = 0;
503 for (int i = 1; i <= m.dim(1); i++)
504 for (int j = 1; j <= m.dim(2); j++)
505 ans += m(i, j);
506 return ans;
507 }
508
509 } //namespace TNT
510
511 #endif //TNTSUPP_H
512