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