matrix.h
1.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#ifndef RTS_MATRIX_H
#define RTS_MATRIX_H
//#include "rts/vector.h"
#include <string.h>
#include <iostream>
#include <stim/math/vector.h>
#include <stim/math/vec3.h>
#include <stim/cuda/cudatools/callable.h>
namespace stim{
template <class T, int N>
struct matrix
{
//the matrix will be stored in column-major order (compatible with OpenGL)
T M[N*N];
CUDA_CALLABLE matrix()
{
for(int r=0; r<N; r++)
for(int c=0; c<N; c++)
if(r == c)
(*this)(r, c) = 1;
else
(*this)(r, c) = 0;
}
CUDA_CALLABLE matrix(T rhs[N*N])
{
memcpy(M,rhs, sizeof(T)*N*N);
}
CUDA_CALLABLE matrix<T,N> set(T rhs[N*N])
{
memcpy(M, rhs, sizeof(T)*N*N);
return *this;
}
CUDA_CALLABLE T& operator()(int row, int col)
{
return M[col * N + row];
}
CUDA_CALLABLE matrix<T, N> operator=(T rhs)
{
int Nsq = N*N;
for(int i=0; i<Nsq; i++)
M[i] = rhs;
return *this;
}
template<typename Y>
vec<Y> operator*(vec<Y> rhs){
unsigned int M = rhs.size();
vec<Y> result;
result.resize(M);
for(int r=0; r<M; r++)
for(int c=0; c<M; c++)
result[r] += (*this)(r, c) * rhs[c];
return result;
}
template<typename Y>
CUDA_CALLABLE vec3<Y> operator*(vec3<Y> rhs){
vec3<Y> result = 0;
for(int r=0; r<3; r++)
for(int c=0; c<3; c++)
result[r] += (*this)(r, c) * rhs[c];
return result;
}
std::string toStr()
{
std::stringstream ss;
for(int r = 0; r < N; r++)
{
ss << "| ";
for(int c=0; c<N; c++)
{
ss << (*this)(r, c) << " ";
}
ss << "|" << std::endl;
}
return ss.str();
}
};
} //end namespace rts
template <typename T, int N>
std::ostream& operator<<(std::ostream& os, stim::matrix<T, N> M)
{
os<<M.toStr();
return os;
}
//#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
//template<class T, int N> using rtsMatrix = rts::matrix<T, N>;
//#endif
#endif