Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
SmallMatrix.hpp
1#ifndef SMALLMATRIX_hpp
2#define SMALLMATRIX_hpp
3
12
13namespace FEDD {
19template <class T>
20class SmallMatrix{
21
22public:
23 typedef unsigned UN;
24 typedef std::vector<T> Vec_T_Type;
25 typedef const std::vector<T> Vec_T_Const_Type;
26 typedef std::vector<Vec_T_Type> Vec2D_T_Type;
27
28 SmallMatrix();
29
30 SmallMatrix(int size);
31
32 SmallMatrix(int size, T v);
33
34 Vec_T_Type &operator[](int i);
35
36 Vec_T_Const_Type &operator[](int i) const;
37
38 SmallMatrix<T> &operator=(const SmallMatrix<T> &sm);
39
40 SmallMatrix<T> operator*(const SmallMatrix<T> &other);
41
42 SmallMatrix<T> operator+(const SmallMatrix<T> &other);
43
44 SmallMatrix<T> &operator*=(const SmallMatrix<T> &other);
45
46 SmallMatrix<T> &operator+=(const SmallMatrix<T> &other);
47
48 int multiply(const SmallMatrix<T> &bMat, SmallMatrix<T> &cMat) const; //this*B=C
49
50 int add(const SmallMatrix<T> &bMat, SmallMatrix<T> &cMat) const; //this+B=C
51
52 int innerProduct(const SmallMatrix<T> &bMat, T &res) const;
53
54 double innerProduct(const SmallMatrix<T> &bMat) const;
55
56 void trace(T &res) const;
57
58 void scale(T scalar);
59
60 int size() const;
61
62 void resize(UN size);
63
64 void print() const;
65
66 double computeInverse(SmallMatrix<T> &inverse) const;
67
68 double computeDet() const;
69
70 double computeScaling() const;
71
72private:
73 Vec2D_T_Type values_;
74 int size_;
75};
76
77
78
79
80
81template<class T>
82SmallMatrix<T>::SmallMatrix():
83values_(0){
84
85}
86template<class T>
87SmallMatrix<T>::SmallMatrix(int size):
88values_(size,Vec_T_Type(size)),
89size_(size){
90
91}
92template<class T>
93SmallMatrix<T>::SmallMatrix(int size, T v):
94values_( size, Vec_T_Type( size, v ) ),
95size_(size){
96
97}
98template<class T>
99typename SmallMatrix<T>::Vec_T_Type &SmallMatrix<T>::operator[](int i) {
100
101 return values_.at(i);
102}
103template<class T>
104typename SmallMatrix<T>::Vec_T_Const_Type &SmallMatrix<T>::operator[](int i) const{
105
106 return values_.at(i);
107}
108template<class T>
109int SmallMatrix<T>::multiply(const SmallMatrix<T> &bMat, SmallMatrix &cMat) const {
110
111 if (size_==2) {
112 cMat[0][0] = values_[0][0]*bMat[0][0] + values_[0][1]*bMat[1][0];
113 cMat[0][1] = values_[0][0]*bMat[0][1] + values_[0][1]*bMat[1][1];
114 cMat[1][0] = values_[1][0]*bMat[0][0] + values_[1][1]*bMat[1][0];
115 cMat[1][1] = values_[1][0]*bMat[0][1] + values_[1][1]*bMat[1][1];
116 }
117 else if(size_==3){
118 cMat[0][0] = values_[0][0]*bMat[0][0] + values_[0][1]*bMat[1][0] + values_[0][2]*bMat[2][0];
119 cMat[0][1] = values_[0][0]*bMat[0][1] + values_[0][1]*bMat[1][1] + values_[0][2]*bMat[2][1];
120 cMat[0][2] = values_[0][0]*bMat[0][2] + values_[0][1]*bMat[1][2] + values_[0][2]*bMat[2][2];
121
122 cMat[1][0] = values_[1][0]*bMat[0][0] + values_[1][1]*bMat[1][0] + values_[1][2]*bMat[2][0];
123 cMat[1][1] = values_[1][0]*bMat[0][1] + values_[1][1]*bMat[1][1] + values_[1][2]*bMat[2][1];
124 cMat[1][2] = values_[1][0]*bMat[0][2] + values_[1][1]*bMat[1][2] + values_[1][2]*bMat[2][2];
125
126 cMat[2][0] = values_[2][0]*bMat[0][0] + values_[2][1]*bMat[1][0] + values_[2][2]*bMat[2][0];
127 cMat[2][1] = values_[2][0]*bMat[0][1] + values_[2][1]*bMat[1][1] + values_[2][2]*bMat[2][1];
128 cMat[2][2] = values_[2][0]*bMat[0][2] + values_[2][1]*bMat[1][2] + values_[2][2]*bMat[2][2];
129 }
130 else
131 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "SmallMatrix wrong size!")
132
133 return 0;
134}
135
136template<class T>
137int SmallMatrix<T>::add(const SmallMatrix<T> &bMat, SmallMatrix &cMat) const {
138
139 for(int i=0; i< size_; i++){
140 for(int j=0; j< size_;j++){
141 cMat[i][j] = values_[i][j]+bMat[i][j];
142 }
143 }
144
145 /*if (size_==2) {
146 cMat[0][0] = values_[0][0]+bMat[0][0];
147 cMat[0][1] = values_[0][1]+bMat[0][1];
148 cMat[1][0] = values_[1][0]+bMat[1][0];
149 cMat[1][1] = values_[1][1]+bMat[1][1];
150 }
151 else if(size_==3){
152 cMat[0][0] = values_[0][0]+bMat[0][0];
153 cMat[0][1] = values_[0][1]+bMat[0][1];
154 cMat[0][2] = values_[0][2]+bMat[0][2];
155
156 cMat[1][0] = values_[1][0]+bMat[1][0];
157 cMat[1][1] = values_[1][1]+bMat[1][1];
158 cMat[1][2] = values_[1][2]+bMat[1][2];
159
160 cMat[2][0] = values_[2][0]+bMat[2][0];
161 cMat[2][1] = values_[2][1]+bMat[2][1];
162 cMat[2][2] = values_[2][2]+bMat[2][2];
163 }
164 else
165 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "SmallMatrix wrong size!")*/
166
167 return 0;
168}
169
170
171template<class T>
172int SmallMatrix<T>::innerProduct(const SmallMatrix<T> &bMat, T &res) const {
173 res = 0.;
174 if (size_==2) {
175 res = values_[0][0]*bMat[0][0] + values_[0][1]*bMat[0][1] +
176 values_[1][0]*bMat[1][0] + values_[1][1]*bMat[1][1];
177 }
178 else if(size_==3){
179 res = values_[0][0]*bMat[0][0] + values_[0][1]*bMat[0][1] + values_[0][2]*bMat[0][2] +
180 values_[1][0]*bMat[1][0] + values_[1][1]*bMat[1][1] + values_[1][2]*bMat[1][2] +
181 values_[2][0]*bMat[2][0] + values_[2][1]*bMat[2][1] + values_[2][2]*bMat[2][2];
182 }
183 else
184 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "SmallMatrix wrong size!")
185
186 return 0;
187}
188
189template<class T>
190double SmallMatrix<T>::innerProduct(const SmallMatrix<T> &bMat) const {
191
192 if (size_==2) {
193 return values_[0][0]*bMat[0][0] + values_[0][1]*bMat[0][1] +
194 values_[1][0]*bMat[1][0] + values_[1][1]*bMat[1][1];
195 }
196 else if(size_==3){
197 return values_[0][0]*bMat[0][0] + values_[0][1]*bMat[0][1] + values_[0][2]*bMat[0][2] +
198 values_[1][0]*bMat[1][0] + values_[1][1]*bMat[1][1] + values_[1][2]*bMat[1][2] +
199 values_[2][0]*bMat[2][0] + values_[2][1]*bMat[2][1] + values_[2][2]*bMat[2][2];
200 }
201 else
202 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "SmallMatrix wrong size!")
203
204 return 0.;
205}
206
207template<class T>
208void SmallMatrix<T>::trace(T &res) const
209{
210 // In res steht das Resultat
211 if(size_ == 2)
212 {
213 res = values_[0][0] + values_[1][1];
214 }
215 else if(size_ == 3)
216 {
217 res = values_[0][0] + values_[1][1] + values_[2][2];
218 }
219 else
220 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "SmallMatrix wrong size!")
221
222
223}
224
225template<class T>
226void SmallMatrix<T>::scale(T scalar){
227
228 for (int i=0; i<size_; i++) {
229 for (int j=0; j<size_; j++) {
230 values_[i][j] *= scalar;
231 }
232 }
233
234}
235
236template<class T>
237int SmallMatrix<T>::size() const{
238 return size_;
239}
240
241template<class T>
242void SmallMatrix<T>::resize(UN size){
243 for (int i=0; i<size_; i++)
244 values_[i].resize(size);
245
246 values_.resize(size,Vec_T_Type(size));
247 size_ = size;
248}
249
250template<class T>
251typename SmallMatrix<T>::SmallMatrix &SmallMatrix<T>::operator=(const SmallMatrix<T> &sm){
252
253 size_ = sm.size();
254 values_.resize(size_);
255
256 for (int i=0; i<sm.size(); i++) {
257 values_.at(i) = sm[i];
258 }
259
260 return *this;
261}
262
263template<class T>
264typename SmallMatrix<T>::SmallMatrix SmallMatrix<T>::operator*(const SmallMatrix<T> &other){
265
266 size_ = other.size();
267 SmallMatrix<T> result(size_);
268
269 this->multiply(other, result);
270
271 return result;
272}
273
274template<class T>
275typename SmallMatrix<T>::SmallMatrix SmallMatrix<T>::operator+(const SmallMatrix<T> &other){
276
277 size_ = other.size();
278 SmallMatrix<T> result(size_);
279
280 this->add(other, result);
281
282 return result;
283}
284
285template<class T>
286typename SmallMatrix<T>::SmallMatrix &SmallMatrix<T>::operator*=(const SmallMatrix<T> &other){
287
288 this->multiply(other, *this);
289
290 return *this;
291}
292
293template<class T>
294typename SmallMatrix<T>::SmallMatrix &SmallMatrix<T>::operator+=(const SmallMatrix<T> &other){
295
296 this->add(other, *this);
297
298 return *this;
299}
300
301template<class T>
302void SmallMatrix<T>::print() const {
303
304 for (int i=0; i<size_; i++) {
305 std::cout << "row "<< i << " values:";
306 for (int j=0; j<size_; j++) {
307 std::cout << " "<< values_[i][j];
308 }
309 std::cout << std::endl;
310 }
311}
312
313template<class T>
314double SmallMatrix<T>::computeInverse(SmallMatrix<T> &inverse) const {
315
316 T det = this->computeDet();
317
318 inverse.resize(size_);
319
320 if (size_==2) {
321 inverse[0][0] = values_[1][1] / det;
322 inverse[0][1] = (-values_[0][1]) / det;
323 inverse[1][0] = (-values_[1][0]) / det;
324 inverse[1][1] = values_[0][0] / det;
325 }
326 else if (size_==3) {
327 inverse[0][0] = (values_[1][1] * values_[2][2] - values_[1][2] * values_[2][1]) / det;
328 inverse[0][1] = (values_[0][2] * values_[2][1] - values_[0][1] * values_[2][2]) / det;
329 inverse[0][2] = (values_[0][1] * values_[1][2] - values_[0][2] * values_[1][1]) / det;
330
331 inverse[1][0] = (values_[1][2] * values_[2][0] - values_[1][0] * values_[2][2]) / det;
332 inverse[1][1] = (values_[0][0] * values_[2][2] - values_[0][2] * values_[2][0]) / det;
333 inverse[1][2] = (values_[0][2] * values_[1][0] - values_[0][0] * values_[1][2]) / det;
334
335 inverse[2][0] = (values_[1][0] * values_[2][1] - values_[1][1] * values_[2][0]) / det;
336 inverse[2][1] = (values_[0][1] * values_[2][0] - values_[0][0] * values_[2][1]) / det;
337 inverse[2][2] = (values_[0][0] * values_[1][1] - values_[0][1] * values_[1][0]) / det;
338
339 }
340 else
341 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Matrix size is not 2 or 3.")
342 return det;
343}
344
345template<class T>
346double SmallMatrix<T>::computeDet() const {
347
348 double det;
349 if (size_==2) {
350 det = values_[0][0] * values_[1][1] - values_[1][0] * values_[0][1];
351 }
352 else if(size_==3){
353 det = values_[0][0] * values_[1][1] * values_[2][2] +
354 values_[0][1] * values_[1][2] * values_[2][0] +
355 values_[0][2] * values_[1][0] * values_[2][1] -
356 values_[2][0] * values_[1][1] * values_[0][2] -
357 values_[2][1] * values_[1][2] * values_[0][0] -
358 values_[2][2] * values_[1][0] * values_[0][1] ;
359 }
360 else
361 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Matrix size is not 2 or 3.")
362
363
364 return det;
365}
366
367template<class T>
368double SmallMatrix<T>::computeScaling() const {
369
370 double scaling;
371 if (size_==2)
372 scaling = std::sqrt( values_[0][0] * values_[0][0] + values_[1][0] * values_[1][0] );
373 else if(size_==3){
374
375 double c1 = ( values_[1][0] * values_[2][1] - values_[2][0] * values_[1][1] );
376 double c2 = ( values_[2][0] * values_[0][1] - values_[0][0] * values_[2][1] );
377 double c3 = ( values_[0][0] * values_[1][1] - values_[1][0] * values_[0][1] );
378
379 scaling = std::sqrt( c1*c1 + c2*c2 + c3*c3 );
380 }
381 else
382 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Matrix size is not 2 or 3.")
383
384
385 return scaling;
386}
387
388 //Pseudo-Inverse: (A^T A)^-1 * A^T
389
390}
391#endif
FEDD::SmallMatrix
This class represents a templated small Matrix of type T. Primarily created for 2x2 and 3x3 matrices....
Definition SmallMatrix.hpp:20
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36