| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms. 2 // See the file tools.license for terms.
3 3
4 #ifndef tools_MATCOM 4 #ifndef tools_MATCOM
5 #define tools_MATCOM 5 #define tools_MATCOM
6 6
7 /* NOTE : bof, no big improvement. 7 /* NOTE : bof, no big improvement.
8 #include <cstring> //memcpy 8 #include <cstring> //memcpy
9 inline void vec_copy(double* a_to,const double 9 inline void vec_copy(double* a_to,const double* a_from,unsigned int a_number) {
10 ::memcpy(a_to,a_from,a_number*sizeof(double) 10 ::memcpy(a_to,a_from,a_number*sizeof(double));
11 } 11 }
12 12
13 template <class T> 13 template <class T>
14 inline void vec_copy(T* a_to,const T* a_from,u 14 inline void vec_copy(T* a_to,const T* a_from,unsigned int a_number) {
15 T* pto = (T*)a_to; 15 T* pto = (T*)a_to;
16 T* pfm = (T*)a_from; 16 T* pfm = (T*)a_from;
17 for(unsigned int i=0;i<a_number;i++,pto++,pf 17 for(unsigned int i=0;i<a_number;i++,pto++,pfm++) *pto = *pfm;
18 } 18 }
19 */ 19 */
20 20
21 /*NOTE : bof, no big improvement. 21 /*NOTE : bof, no big improvement.
22 #include <Accelerate/Accelerate.h> 22 #include <Accelerate/Accelerate.h>
23 23
24 inline void vec_add(const float* a_1,const flo 24 inline void vec_add(const float* a_1,const float* a_2,float* a_res,unsigned int a_number) {
25 ::vDSP_vadd(a_1,1,a_2,1,a_res,1,a_number); 25 ::vDSP_vadd(a_1,1,a_2,1,a_res,1,a_number);
26 } 26 }
27 inline void vec_sub(const float* a_1,const flo 27 inline void vec_sub(const float* a_1,const float* a_2,float* a_res,unsigned int a_number) {
28 ::vDSP_vsub(a_1,1,a_2,1,a_res,1,a_number); 28 ::vDSP_vsub(a_1,1,a_2,1,a_res,1,a_number);
29 } 29 }
30 */ 30 */
31 31
32 /* 32 /*
33 template <class T> 33 template <class T>
34 inline void vec_add(const T* a_1,const T* a_2, 34 inline void vec_add(const T* a_1,const T* a_2,T* a_res,unsigned int a_number) {
35 T* p1 = (T*)a_1; 35 T* p1 = (T*)a_1;
36 T* p2 = (T*)a_2; 36 T* p2 = (T*)a_2;
37 T* pr = (T*)a_res; 37 T* pr = (T*)a_res;
38 for(unsigned int i=0;i<a_number;i++,p1++,p2+ 38 for(unsigned int i=0;i<a_number;i++,p1++,p2++,pr++) *pr = *p1+*p2;
39 } 39 }
40 40
41 template <class T> 41 template <class T>
42 inline void vec_sub(const T* a_1,const T* a_2, 42 inline void vec_sub(const T* a_1,const T* a_2,T* a_res,unsigned int a_number) {
43 T* p1 = (T*)a_1; 43 T* p1 = (T*)a_1;
44 T* p2 = (T*)a_2; 44 T* p2 = (T*)a_2;
45 T* pr = (T*)a_res; 45 T* pr = (T*)a_res;
46 for(unsigned int i=0;i<a_number;i++,p1++,p2+ 46 for(unsigned int i=0;i<a_number;i++,p1++,p2++,pr++) *pr = *p1-*p2;
47 } 47 }
48 */ 48 */
49 49
50 #include "../eqT" 50 #include "../eqT"
51 51
52 #include <cstddef> //size_t 52 #include <cstddef> //size_t
53 53
54 // common code to class mat and nmat. 54 // common code to class mat and nmat.
55 55
56 #define TOOLS_MATCOM \ 56 #define TOOLS_MATCOM \
57 protected:\ 57 protected:\
58 static T zero() {return T();}\ 58 static T zero() {return T();}\
59 static T one() {return T(1);}\ 59 static T one() {return T(1);}\
60 static T minus_one() {return T(-1);}\ 60 static T minus_one() {return T(-1);}\
61 static T two() {return T(2);}\ 61 static T two() {return T(2);}\
62 public:\ 62 public:\
63 typedef T elem_t;\ 63 typedef T elem_t;\
64 typedef unsigned int size_type;\ 64 typedef unsigned int size_type;\
65 public:\ 65 public:\
66 unsigned int rows() const {return dimension( 66 unsigned int rows() const {return dimension();}\
67 unsigned int cols() const {return dimension( 67 unsigned int cols() const {return dimension();}\
68 \ 68 \
69 void set_value(unsigned int aR,unsigned int 69 void set_value(unsigned int aR,unsigned int aC,const T& a_value) { \
70 m_vec[aR + aC * dimension()] = a_value;\ 70 m_vec[aR + aC * dimension()] = a_value;\
71 }\ 71 }\
72 \ 72 \
73 const T& value(unsigned int aR,unsigned int 73 const T& value(unsigned int aR,unsigned int aC) const { \
74 return m_vec[aR + aC * dimension()];\ 74 return m_vec[aR + aC * dimension()];\
75 }\ 75 }\
76 \ 76 \
77 T value(unsigned int aR,unsigned int aC) { \ 77 T value(unsigned int aR,unsigned int aC) { \
78 return m_vec[aR + aC * dimension()];\ 78 return m_vec[aR + aC * dimension()];\
79 }\ 79 }\
80 \ 80 \
81 void set_matrix(const TOOLS_MAT_CLASS& a_m){ 81 void set_matrix(const TOOLS_MAT_CLASS& a_m){ /*optimization.*/\
82 _copy(a_m.m_vec);\ 82 _copy(a_m.m_vec);\
83 }\ 83 }\
84 \ 84 \
85 void set_constant(const T& a_v){\ 85 void set_constant(const T& a_v){\
86 for(unsigned int i=0;i<dim2();i++) m_vec[i 86 for(unsigned int i=0;i<dim2();i++) m_vec[i] = a_v;\
87 }\ 87 }\
88 void set_zero(){\ 88 void set_zero(){\
89 set_constant(zero());\ 89 set_constant(zero());\
90 }\ 90 }\
91 void set_identity() {\ 91 void set_identity() {\
92 set_zero();\ 92 set_zero();\
93 for(unsigned int i=0;i<dimension();i++) m_ 93 for(unsigned int i=0;i<dimension();i++) m_vec[i+i*dimension()] = one();\
94 }\ 94 }\
95 void set_diagonal(const T& a_s) {\ 95 void set_diagonal(const T& a_s) {\
96 set_zero();\ 96 set_zero();\
97 for(unsigned int i=0;i<dimension();i++) m_ 97 for(unsigned int i=0;i<dimension();i++) m_vec[i+i*dimension()] = a_s;\
98 }\ 98 }\
99 typedef T (*func)(const T&);\ 99 typedef T (*func)(const T&);\
100 void apply_func(func a_func) {\ 100 void apply_func(func a_func) {\
101 T* pos = m_vec;\ 101 T* pos = m_vec;\
102 for(unsigned int i=0;i<dim2();i++,pos++) * 102 for(unsigned int i=0;i<dim2();i++,pos++) *pos = a_func(*pos);\
103 }\ 103 }\
104 template <class RANDOM>\ 104 template <class RANDOM>\
105 void set_random(RANDOM& a_random) {\ 105 void set_random(RANDOM& a_random) {\
106 for(unsigned int i=0;i<dim2();i++) m_vec[i 106 for(unsigned int i=0;i<dim2();i++) m_vec[i] = a_random.shoot();\
107 }\ 107 }\
108 template <class RANDOM>\ 108 template <class RANDOM>\
109 void set_symmetric_random(RANDOM& a_random) 109 void set_symmetric_random(RANDOM& a_random) {\
110 unsigned int _D = dimension();\ 110 unsigned int _D = dimension();\
111 {for(unsigned int r=0;r<_D;r++) set_value(r 111 {for(unsigned int r=0;r<_D;r++) set_value(r,r,a_random.shoot());}\
112 T rd;\ 112 T rd;\
113 {for(unsigned int r=0;r<_D;r++) {\ 113 {for(unsigned int r=0;r<_D;r++) {\
114 for(unsigned int c=(r+1);c<_D;c++) {\ 114 for(unsigned int c=(r+1);c<_D;c++) {\
115 rd = a_random.shoot();\ 115 rd = a_random.shoot();\
116 set_value(r,c,rd);\ 116 set_value(r,c,rd);\
117 set_value(c,r,rd);\ 117 set_value(c,r,rd);\
118 }\ 118 }\
119 }}\ 119 }}\
120 }\ 120 }\
121 template <class RANDOM>\ 121 template <class RANDOM>\
122 void set_antisymmetric_random(RANDOM& a_rand 122 void set_antisymmetric_random(RANDOM& a_random) {\
123 unsigned int _D = dimension();\ 123 unsigned int _D = dimension();\
124 {for(unsigned int r=0;r<_D;r++) set_value(r 124 {for(unsigned int r=0;r<_D;r++) set_value(r,r,zero());}\
125 T rd;\ 125 T rd;\
126 {for(unsigned int r=0;r<_D;r++) {\ 126 {for(unsigned int r=0;r<_D;r++) {\
127 for(unsigned int c=(r+1);c<_D;c++) {\ 127 for(unsigned int c=(r+1);c<_D;c++) {\
128 rd = a_random.shoot();\ 128 rd = a_random.shoot();\
129 set_value(r,c,rd);\ 129 set_value(r,c,rd);\
130 set_value(c,r,minus_one()*rd);\ 130 set_value(c,r,minus_one()*rd);\
131 }\ 131 }\
132 }}\ 132 }}\
133 }\ 133 }\
134 public:\ 134 public:\
135 template <class ARRAY>\ 135 template <class ARRAY>\
136 bool mul_array(ARRAY& a_array,T a_tmp[]) con 136 bool mul_array(ARRAY& a_array,T a_tmp[]) const {\
137 /* a_array = this *= a_array */\ 137 /* a_array = this *= a_array */\
138 unsigned int _dim = dimension();\ 138 unsigned int _dim = dimension();\
139 T* pos = a_tmp;\ 139 T* pos = a_tmp;\
140 for(unsigned int r=0;r<_dim;r++,pos++) {\ 140 for(unsigned int r=0;r<_dim;r++,pos++) {\
141 *pos = T();\ 141 *pos = T();\
142 for(unsigned int c=0;c<_dim;c++) *pos += 142 for(unsigned int c=0;c<_dim;c++) *pos += m_vec[r+c*_dim]*a_array[c];\
143 }\ 143 }\
144 {for(unsigned int i=0;i<_dim;i++) a_array[i 144 {for(unsigned int i=0;i<_dim;i++) a_array[i] = a_tmp[i];}\
145 return true;\ 145 return true;\
146 }\ 146 }\
147 template <class VEC>\ 147 template <class VEC>\
148 bool mul_vec(VEC& a_vec,T a_tmp[]) const {\ 148 bool mul_vec(VEC& a_vec,T a_tmp[]) const {\
149 /* a_vec = this *= a_vec */\ 149 /* a_vec = this *= a_vec */\
150 unsigned int _dim = dimension();\ 150 unsigned int _dim = dimension();\
151 if(a_vec.dimension()!=_dim) return false;\ 151 if(a_vec.dimension()!=_dim) return false;\
152 T* pos = a_tmp;\ 152 T* pos = a_tmp;\
153 for(unsigned int r=0;r<_dim;r++,pos++) {\ 153 for(unsigned int r=0;r<_dim;r++,pos++) {\
154 *pos = T();\ 154 *pos = T();\
155 for(unsigned int c=0;c<_dim;c++) *pos += 155 for(unsigned int c=0;c<_dim;c++) *pos += m_vec[r+c*_dim]*a_vec[c];\
156 }\ 156 }\
157 {for(unsigned int i=0;i<_dim;i++) a_vec[i] 157 {for(unsigned int i=0;i<_dim;i++) a_vec[i] = a_tmp[i];}\
158 return true;\ 158 return true;\
159 }\ 159 }\
160 template <class VEC>\ 160 template <class VEC>\
161 bool mul_vec(VEC& a_vec) const {\ 161 bool mul_vec(VEC& a_vec) const {\
162 T* _tmp = new T[dimension()];\ 162 T* _tmp = new T[dimension()];\
163 bool status = mul_vec(a_vec,_tmp);\ 163 bool status = mul_vec(a_vec,_tmp);\
164 delete [] _tmp;\ 164 delete [] _tmp;\
165 return status;\ 165 return status;\
166 }\ 166 }\
167 \ 167 \
168 bool mul_array(T a_vec[],T a_tmp[]) const {\ 168 bool mul_array(T a_vec[],T a_tmp[]) const {\
169 /* a_vec = this *= a_vec */\ 169 /* a_vec = this *= a_vec */\
170 unsigned int _dim = dimension();\ 170 unsigned int _dim = dimension();\
171 T* pos = a_tmp;\ 171 T* pos = a_tmp;\
172 for(unsigned int r=0;r<_dim;r++,pos++) {\ 172 for(unsigned int r=0;r<_dim;r++,pos++) {\
173 *pos = T();\ 173 *pos = T();\
174 for(unsigned int c=0;c<_dim;c++) *pos += 174 for(unsigned int c=0;c<_dim;c++) *pos += m_vec[r+c*_dim]*a_vec[c];\
175 }\ 175 }\
176 {for(unsigned int i=0;i<_dim;i++) a_vec[i] 176 {for(unsigned int i=0;i<_dim;i++) a_vec[i] = a_tmp[i];}\
177 return true;\ 177 return true;\
178 }\ 178 }\
179 bool mul_array(T a_vec[]) const {\ 179 bool mul_array(T a_vec[]) const {\
180 T* _tmp = new T[dimension()];\ 180 T* _tmp = new T[dimension()];\
181 bool status = mul_array(a_vec,_tmp);\ 181 bool status = mul_array(a_vec,_tmp);\
182 delete [] _tmp;\ 182 delete [] _tmp;\
183 return status;\ 183 return status;\
184 }\ 184 }\
185 \ 185 \
186 void mul_mtx(const TOOLS_MAT_CLASS& a_m) {\ 186 void mul_mtx(const TOOLS_MAT_CLASS& a_m) {\
187 _mul_mtx(a_m.m_vec);\ 187 _mul_mtx(a_m.m_vec);\
188 }\ 188 }\
189 void mul_mtx(const TOOLS_MAT_CLASS& a_m,T a_ 189 void mul_mtx(const TOOLS_MAT_CLASS& a_m,T a_tmp[]) {\
190 _mul_mtx(a_m.m_vec,a_tmp);\ 190 _mul_mtx(a_m.m_vec,a_tmp);\
191 }\ 191 }\
192 void left_mul_mtx(const TOOLS_MAT_CLASS& a_m 192 void left_mul_mtx(const TOOLS_MAT_CLASS& a_m) { \
193 /* this = a_m * this :*/\ 193 /* this = a_m * this :*/\
194 _left_mul_mtx(a_m.m_vec);\ 194 _left_mul_mtx(a_m.m_vec);\
195 }\ 195 }\
196 bool equal(const TOOLS_MAT_CLASS& a_m) const 196 bool equal(const TOOLS_MAT_CLASS& a_m) const {\
197 if(&a_m==this) return true;\ 197 if(&a_m==this) return true;\
198 for(unsigned int i=0;i<dim2();i++) {\ 198 for(unsigned int i=0;i<dim2();i++) {\
199 if(m_vec[i]!=a_m.m_vec[i]) return false; 199 if(m_vec[i]!=a_m.m_vec[i]) return false;\
200 }\ 200 }\
201 return true;\ 201 return true;\
202 }\ 202 }\
203 \ 203 \
204 bool equal(const TOOLS_MAT_CLASS& a_m,const 204 bool equal(const TOOLS_MAT_CLASS& a_m,const T& a_prec) const {\
205 if(&a_m==this) return true;\ 205 if(&a_m==this) return true;\
206 T* tp = (T*)m_vec;\ 206 T* tp = (T*)m_vec;\
207 T* mp = (T*)a_m.m_vec;\ 207 T* mp = (T*)a_m.m_vec;\
208 for(unsigned int i=0;i<dim2();i++,tp++,mp+ 208 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
209 T diff = (*tp) - (*mp);\ 209 T diff = (*tp) - (*mp);\
210 if(diff<zero()) diff *= minus_one();\ 210 if(diff<zero()) diff *= minus_one();\
211 if(diff>=a_prec) return false;\ 211 if(diff>=a_prec) return false;\
212 }\ 212 }\
213 return true;\ 213 return true;\
214 }\ 214 }\
215 \ 215 \
216 template <class PREC>\ 216 template <class PREC>\
217 bool equal_prec(const TOOLS_MAT_CLASS& a_m,c 217 bool equal_prec(const TOOLS_MAT_CLASS& a_m,const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
218 if(&a_m==this) return true;\ 218 if(&a_m==this) return true;\
219 T* tp = (T*)m_vec;\ 219 T* tp = (T*)m_vec;\
220 T* mp = (T*)a_m.m_vec;\ 220 T* mp = (T*)a_m.m_vec;\
221 for(unsigned int i=0;i<dim2();i++,tp++,mp+ 221 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
222 T diff = (*tp) - (*mp);\ 222 T diff = (*tp) - (*mp);\
223 if(a_fabs(diff)>=a_prec) return false;\ 223 if(a_fabs(diff)>=a_prec) return false;\
224 }\ 224 }\
225 return true;\ 225 return true;\
226 }\ 226 }\
227 \ 227 \
228 void mx_diff(const TOOLS_MAT_CLASS& a_m,T& a 228 void mx_diff(const TOOLS_MAT_CLASS& a_m,T& a_mx_diff) const {\
229 T* tp = (T*)m_vec;\ 229 T* tp = (T*)m_vec;\
230 T* mp = (T*)a_m.m_vec;\ 230 T* mp = (T*)a_m.m_vec;\
231 a_mx_diff = (*tp) - (*mp);\ 231 a_mx_diff = (*tp) - (*mp);\
232 if(a_mx_diff<zero()) a_mx_diff *= minus_on 232 if(a_mx_diff<zero()) a_mx_diff *= minus_one();\
233 for(unsigned int i=0;i<dim2();i++,tp++,mp+ 233 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
234 T diff = (*tp) - (*mp);\ 234 T diff = (*tp) - (*mp);\
235 if(diff<zero()) diff *= minus_one();\ 235 if(diff<zero()) diff *= minus_one();\
236 a_mx_diff = (diff>a_mx_diff?diff:a_mx_di 236 a_mx_diff = (diff>a_mx_diff?diff:a_mx_diff);\
237 }\ 237 }\
238 }\ 238 }\
239 \ 239 \
240 bool to_rm_is_proportional(const TOOLS_MAT_C 240 bool to_rm_is_proportional(const TOOLS_MAT_CLASS& a_m,const T& a_prec,T& a_factor) const {\
241 if(&a_m==this) {a_factor=one();return true 241 if(&a_m==this) {a_factor=one();return true;}\
242 /* If true, then : a_m = a_factor * this.* 242 /* If true, then : a_m = a_factor * this.*/\
243 a_factor = zero();\ 243 a_factor = zero();\
244 T* tp = (T*)m_vec;\ 244 T* tp = (T*)m_vec;\
245 T* mp = (T*)a_m.m_vec;\ 245 T* mp = (T*)a_m.m_vec;\
246 bool first = true;\ 246 bool first = true;\
247 for(unsigned int i=0;i<dim2();i++,tp++,mp+ 247 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
248 if( ((*tp)==zero()) && ((*mp)==ze 248 if( ((*tp)==zero()) && ((*mp)==zero())) {\
249 continue;\ 249 continue;\
250 /*} else if( ((*tp)!=zero()) && ((*mp)== 250 /*} else if( ((*tp)!=zero()) && ((*mp)==zero())) {\
251 return false;*/\ 251 return false;*/\
252 } else if( ((*tp)==zero()) && ((*mp)!=ze 252 } else if( ((*tp)==zero()) && ((*mp)!=zero())) {\
253 return false;\ 253 return false;\
254 } else {\ 254 } else {\
255 if(first) {\ 255 if(first) {\
256 a_factor = (*mp)/(*tp);\ 256 a_factor = (*mp)/(*tp);\
257 first = false;\ 257 first = false;\
258 } else {\ 258 } else {\
259 T diff = (*tp)*a_factor - (*mp);\ 259 T diff = (*tp)*a_factor - (*mp);\
260 if(diff<zero()) diff *= minus_one(); 260 if(diff<zero()) diff *= minus_one();\
261 if(diff>=a_prec) return false;\ 261 if(diff>=a_prec) return false;\
262 }\ 262 }\
263 }\ 263 }\
264 }\ 264 }\
265 return true;\ 265 return true;\
266 }\ 266 }\
267 \ 267 \
268 public:\ 268 public:\
269 bool is_proportional(const TOOLS_MAT_CLASS& 269 bool is_proportional(const TOOLS_MAT_CLASS& a_right,T& a_factor) const {\
270 /* If true, then : a_right = a_factor * th 270 /* If true, then : a_right = a_factor * this. a_factor could be zero.*/\
271 if(this==&a_right) {a_factor=T(1);return t 271 if(this==&a_right) {a_factor=T(1);return true;}\
272 a_factor = zero();\ 272 a_factor = zero();\
273 if(dimension()!=a_right.dimension()) retur 273 if(dimension()!=a_right.dimension()) return false;\
274 T* lp = (T*)m_vec;\ 274 T* lp = (T*)m_vec;\
275 T* rp = (T*)a_right.m_vec;\ 275 T* rp = (T*)a_right.m_vec;\
276 bool first = true;\ 276 bool first = true;\
277 size_t _data_size = data_size();\ 277 size_t _data_size = data_size();\
278 for(size_t i=0;i<_data_size;i++,lp++,rp++) 278 for(size_t i=0;i<_data_size;i++,lp++,rp++) {\
279 if(*lp==zero()) {\ 279 if(*lp==zero()) {\
280 if(*rp==zero()) continue;\ 280 if(*rp==zero()) continue;\
281 return false;\ 281 return false;\
282 }\ 282 }\
283 if(first) {\ 283 if(first) {\
284 a_factor = (*rp)/(*lp);\ 284 a_factor = (*rp)/(*lp);\
285 first = false;\ 285 first = false;\
286 continue;\ 286 continue;\
287 }\ 287 }\
288 if((*lp)*a_factor!=(*rp)) return false;\ 288 if((*lp)*a_factor!=(*rp)) return false;\
289 }\ 289 }\
290 return true;\ 290 return true;\
291 }\ 291 }\
292 \ 292 \
293 template <class PREC>\ 293 template <class PREC>\
294 bool is_proportional_prec(const TOOLS_MAT_CL 294 bool is_proportional_prec(const TOOLS_MAT_CLASS& a_right,const PREC& a_prec,PREC(*a_fabs)(const T&),T& a_factor) const {\
295 /* If true, then : a_right = a_factor * th 295 /* If true, then : a_right = a_factor * this. a_factor could be zero.*/\
296 if(this==&a_right) {a_factor=T(1);return t 296 if(this==&a_right) {a_factor=T(1);return true;}\
297 a_factor = zero();\ 297 a_factor = zero();\
298 if(dimension()!=a_right.dimension()) retur 298 if(dimension()!=a_right.dimension()) return false;\
299 T* lp = (T*)m_vec;\ 299 T* lp = (T*)m_vec;\
300 T* rp = (T*)a_right.m_vec;\ 300 T* rp = (T*)a_right.m_vec;\
301 bool first = true;\ 301 bool first = true;\
302 size_t _data_size = data_size();\ 302 size_t _data_size = data_size();\
303 for(size_t i=0;i<_data_size;i++,lp++,rp++) 303 for(size_t i=0;i<_data_size;i++,lp++,rp++) {\
304 if(is_zero(*lp,a_prec,a_fabs)) {\ 304 if(is_zero(*lp,a_prec,a_fabs)) {\
305 if(is_zero(*rp,a_prec,a_fabs)) continu 305 if(is_zero(*rp,a_prec,a_fabs)) continue;\
306 return false;\ 306 return false;\
307 }\ 307 }\
308 if(first) {\ 308 if(first) {\
309 a_factor = (*rp)/(*lp);\ 309 a_factor = (*rp)/(*lp);\
310 first = false;\ 310 first = false;\
311 continue;\ 311 continue;\
312 }\ 312 }\
313 if(!numbers_are_equal((*lp)*a_factor,*rp 313 if(!numbers_are_equal((*lp)*a_factor,*rp,a_prec,a_fabs)) return false;\
314 }\ 314 }\
315 return true;\ 315 return true;\
316 }\ 316 }\
317 \ 317 \
318 bool is_diagonal() const {\ 318 bool is_diagonal() const {\
319 unsigned int _D = dimension();\ 319 unsigned int _D = dimension();\
320 for(unsigned int r=0;r<_D;r++) {\ 320 for(unsigned int r=0;r<_D;r++) {\
321 for(unsigned int c=0;c<_D;c++) {\ 321 for(unsigned int c=0;c<_D;c++) {\
322 if(c!=r) {if(value(r,c)!=zero()) retur 322 if(c!=r) {if(value(r,c)!=zero()) return false;}\
323 }\ 323 }\
324 }\ 324 }\
325 return true;\ 325 return true;\
326 }\ 326 }\
327 \ 327 \
328 template <class PREC>\ 328 template <class PREC>\
329 bool is_diagonal_prec(const PREC& a_prec,PRE 329 bool is_diagonal_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
330 unsigned int _D = dimension();\ 330 unsigned int _D = dimension();\
331 for(unsigned int r=0;r<_D;r++) {\ 331 for(unsigned int r=0;r<_D;r++) {\
332 for(unsigned int c=0;c<_D;c++) {\ 332 for(unsigned int c=0;c<_D;c++) {\
333 if(c!=r) {if(!is_zero(value(r,c),a_pre 333 if(c!=r) {if(!is_zero(value(r,c),a_prec,a_fabs)) return false;}\
334 }\ 334 }\
335 }\ 335 }\
336 return true;\ 336 return true;\
337 }\ 337 }\
338 \ 338 \
339 bool is_identity() const {\ 339 bool is_identity() const {\
340 unsigned int _D = dimension();\ 340 unsigned int _D = dimension();\
341 {for(unsigned int r=0;r<_D;r++) {\ 341 {for(unsigned int r=0;r<_D;r++) {\
342 if(value(r,r)!=one()) return false;\ 342 if(value(r,r)!=one()) return false;\
343 }}\ 343 }}\
344 {for(unsigned int r=0;r<_D;r++) {\ 344 {for(unsigned int r=0;r<_D;r++) {\
345 for(unsigned int c=0;c<_D;c++) {\ 345 for(unsigned int c=0;c<_D;c++) {\
346 if(c!=r) {if(value(r,c)!=zero()) retur 346 if(c!=r) {if(value(r,c)!=zero()) return false;}\
347 }\ 347 }\
348 }}\ 348 }}\
349 return true;\ 349 return true;\
350 }\ 350 }\
351 \ 351 \
352 template <class PREC>\ 352 template <class PREC>\
353 bool is_identity_prec(const PREC& a_prec,PRE 353 bool is_identity_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
354 unsigned int _D = dimension();\ 354 unsigned int _D = dimension();\
355 {for(unsigned int r=0;r<_D;r++) {\ 355 {for(unsigned int r=0;r<_D;r++) {\
356 if(!numbers_are_equal(value(r,r),one(),a 356 if(!numbers_are_equal(value(r,r),one(),a_prec,a_fabs)) return false;\
357 }}\ 357 }}\
358 {for(unsigned int r=0;r<_D;r++) {\ 358 {for(unsigned int r=0;r<_D;r++) {\
359 for(unsigned int c=0;c<_D;c++) {\ 359 for(unsigned int c=0;c<_D;c++) {\
360 if(c!=r) {if(!is_zero(value(r,c),a_pre 360 if(c!=r) {if(!is_zero(value(r,c),a_prec,a_fabs)) return false;}\
361 }\ 361 }\
362 }}\ 362 }}\
363 return true;\ 363 return true;\
364 }\ 364 }\
365 \ 365 \
366 const T* data() const {return m_vec;}\ 366 const T* data() const {return m_vec;}\
367 unsigned int size() const {return dim2();}\ 367 unsigned int size() const {return dim2();}\
368 unsigned int data_size() const {return dim2( 368 unsigned int data_size() const {return dim2();} /*for mathz*/\
369 \ 369 \
370 T trace() const {\ 370 T trace() const {\
371 T _value = zero();\ 371 T _value = zero();\
372 unsigned int _D = dimension();\ 372 unsigned int _D = dimension();\
373 for(unsigned int c=0;c<_D;c++) _value += m 373 for(unsigned int c=0;c<_D;c++) _value += m_vec[c+c*_D];\
374 return _value;\ 374 return _value;\
375 }\ 375 }\
376 \ 376 \
377 void transpose() {\ 377 void transpose() {\
378 unsigned int _D = dimension();\ 378 unsigned int _D = dimension();\
379 for(unsigned int r=0;r<_D;r++) {\ 379 for(unsigned int r=0;r<_D;r++) {\
380 for(unsigned int c=(r+1);c<_D;c++) {\ 380 for(unsigned int c=(r+1);c<_D;c++) {\
381 T vrc = value(r,c);\ 381 T vrc = value(r,c);\
382 T vcr = value(c,r);\ 382 T vcr = value(c,r);\
383 set_value(r,c,vcr);\ 383 set_value(r,c,vcr);\
384 set_value(c,r,vrc);\ 384 set_value(c,r,vrc);\
385 }\ 385 }\
386 }\ 386 }\
387 }\ 387 }\
388 \ 388 \
389 void multiply(const T& a_T) {\ 389 void multiply(const T& a_T) {\
390 for(unsigned int i=0;i<dim2();i++) m_vec[i 390 for(unsigned int i=0;i<dim2();i++) m_vec[i] *= a_T;\
391 }\ 391 }\
392 \ 392 \
393 bool is_symmetric() const {\ 393 bool is_symmetric() const {\
394 unsigned int _D = dimension();\ 394 unsigned int _D = dimension();\
395 for(unsigned int r=0;r<_D;r++) {\ 395 for(unsigned int r=0;r<_D;r++) {\
396 for(unsigned int c=(r+1);c<_D;c++) {\ 396 for(unsigned int c=(r+1);c<_D;c++) {\
397 if(value(r,c)!=value(c,r)) return fals 397 if(value(r,c)!=value(c,r)) return false;\
398 }\ 398 }\
399 }\ 399 }\
400 return true;\ 400 return true;\
401 }\ 401 }\
402 \ 402 \
403 template <class PREC>\ 403 template <class PREC>\
404 bool is_symmetric_prec(const PREC& a_prec,PR 404 bool is_symmetric_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
405 unsigned int _D = dimension();\ 405 unsigned int _D = dimension();\
406 for(unsigned int r=0;r<_D;r++) {\ 406 for(unsigned int r=0;r<_D;r++) {\
407 for(unsigned int c=(r+1);c<_D;c++) {\ 407 for(unsigned int c=(r+1);c<_D;c++) {\
408 T diff = value(r,c)-value(c,r);\ 408 T diff = value(r,c)-value(c,r);\
409 if(a_fabs(diff)>=a_prec) return false; 409 if(a_fabs(diff)>=a_prec) return false;\
410 }\ 410 }\
411 }\ 411 }\
412 return true;\ 412 return true;\
413 }\ 413 }\
414 \ 414 \
415 bool is_antisymmetric() const {\ 415 bool is_antisymmetric() const {\
416 unsigned int _D = dimension();\ 416 unsigned int _D = dimension();\
417 {for(unsigned int r=0;r<_D;r++) {\ 417 {for(unsigned int r=0;r<_D;r++) {\
418 if(value(r,r)!=zero()) return false;\ 418 if(value(r,r)!=zero()) return false;\
419 }}\ 419 }}\
420 for(unsigned int r=0;r<_D;r++) {\ 420 for(unsigned int r=0;r<_D;r++) {\
421 for(unsigned int c=(r+1);c<_D;c++) {\ 421 for(unsigned int c=(r+1);c<_D;c++) {\
422 if(value(r,c)!=minus_one()*value(c,r)) 422 if(value(r,c)!=minus_one()*value(c,r)) return false;\
423 }\ 423 }\
424 }\ 424 }\
425 return true;\ 425 return true;\
426 }\ 426 }\
427 \ 427 \
428 template <class PREC>\ 428 template <class PREC>\
429 bool is_antisymmetric_prec(const PREC& a_pre 429 bool is_antisymmetric_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
430 unsigned int _D = dimension();\ 430 unsigned int _D = dimension();\
431 {for(unsigned int r=0;r<_D;r++) {\ 431 {for(unsigned int r=0;r<_D;r++) {\
432 if(a_fabs(value(r,r))>=a_prec) return fa 432 if(a_fabs(value(r,r))>=a_prec) return false;\
433 }}\ 433 }}\
434 for(unsigned int r=0;r<_D;r++) {\ 434 for(unsigned int r=0;r<_D;r++) {\
435 for(unsigned int c=(r+1);c<_D;c++) {\ 435 for(unsigned int c=(r+1);c<_D;c++) {\
436 T diff = value(r,c)-minus_one()*value( 436 T diff = value(r,c)-minus_one()*value(c,r);\
437 if(a_fabs(diff)>=a_prec) return false; 437 if(a_fabs(diff)>=a_prec) return false;\
438 }\ 438 }\
439 }\ 439 }\
440 return true;\ 440 return true;\
441 }\ 441 }\
442 \ 442 \
443 void symmetric_part(TOOLS_MAT_CLASS& a_res) 443 void symmetric_part(TOOLS_MAT_CLASS& a_res) const {\
444 a_res = *this;\ 444 a_res = *this;\
445 a_res.transpose();\ 445 a_res.transpose();\
446 a_res += *this;\ 446 a_res += *this;\
447 a_res.multiply(one()/two());\ 447 a_res.multiply(one()/two());\
448 }\ 448 }\
449 \ 449 \
450 void antisymmetric_part(TOOLS_MAT_CLASS& a_r 450 void antisymmetric_part(TOOLS_MAT_CLASS& a_res) const {\
451 a_res = *this;\ 451 a_res = *this;\
452 a_res.transpose();\ 452 a_res.transpose();\
453 a_res.multiply(minus_one());\ 453 a_res.multiply(minus_one());\
454 a_res += *this;\ 454 a_res += *this;\
455 a_res.multiply(one()/two());\ 455 a_res.multiply(one()/two());\
456 }\ 456 }\
457 \ 457 \
458 template <class PREC>\ 458 template <class PREC>\
459 bool is_block_UL_DR(const PREC& a_prec,PREC( 459 bool is_block_UL_DR(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
460 /* Look if even dim matrix is of the form 460 /* Look if even dim matrix is of the form : |X 0|\
461 461 |0 Y|*/\
462 unsigned int _D = dimension();\ 462 unsigned int _D = dimension();\
463 unsigned int _D_2 = _D/2;\ 463 unsigned int _D_2 = _D/2;\
464 if(2*_D_2!=_D) return false;\ 464 if(2*_D_2!=_D) return false;\
465 for(unsigned int r=0;r<_D_2;r++) {\ 465 for(unsigned int r=0;r<_D_2;r++) {\
466 for(unsigned int c=_D_2;c<_D;c++) {\ 466 for(unsigned int c=_D_2;c<_D;c++) {\
467 if(a_fabs(value(r,c))>=a_prec) return 467 if(a_fabs(value(r,c))>=a_prec) return false;\
468 }\ 468 }\
469 }\ 469 }\
470 for(unsigned int r=_D_2;r<_D;r++) {\ 470 for(unsigned int r=_D_2;r<_D;r++) {\
471 for(unsigned int c=0;c<_D_2;c++) {\ 471 for(unsigned int c=0;c<_D_2;c++) {\
472 if(a_fabs(value(r,c))>=a_prec) return 472 if(a_fabs(value(r,c))>=a_prec) return false;\
473 }\ 473 }\
474 }\ 474 }\
475 return true;\ 475 return true;\
476 }\ 476 }\
477 \ 477 \
478 template <class PREC>\ 478 template <class PREC>\
479 bool is_block_UR_DL(const PREC& a_prec,PREC( 479 bool is_block_UR_DL(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
480 /* Look if even dim matrix is of the form 480 /* Look if even dim matrix is of the form : |0 X|\
481 481 |Y 0|*/\
482 unsigned int _D = dimension();\ 482 unsigned int _D = dimension();\
483 unsigned int _D_2 = _D/2;\ 483 unsigned int _D_2 = _D/2;\
484 if(2*_D_2!=_D) return false;\ 484 if(2*_D_2!=_D) return false;\
485 for(unsigned int r=0;r<_D_2;r++) {\ 485 for(unsigned int r=0;r<_D_2;r++) {\
486 for(unsigned int c=0;c<_D_2;c++) {\ 486 for(unsigned int c=0;c<_D_2;c++) {\
487 if(a_fabs(value(r,c))>=a_prec) return 487 if(a_fabs(value(r,c))>=a_prec) return false;\
488 }\ 488 }\
489 }\ 489 }\
490 for(unsigned int r=_D_2;r<_D;r++) {\ 490 for(unsigned int r=_D_2;r<_D;r++) {\
491 for(unsigned int c=_D_2;c<_D;c++) {\ 491 for(unsigned int c=_D_2;c<_D;c++) {\
492 if(a_fabs(value(r,c))>=a_prec) return 492 if(a_fabs(value(r,c))>=a_prec) return false;\
493 }\ 493 }\
494 }\ 494 }\
495 return true;\ 495 return true;\
496 }\ 496 }\
497 \ 497 \
498 template <class PREC>\ 498 template <class PREC>\
499 bool is_decomplex(const PREC& a_prec,PREC(*a 499 bool is_decomplex(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
500 /* Look if even dim matrix is of the form 500 /* Look if even dim matrix is of the form : | X Y|\
501 501 |-Y X|*/\
502 unsigned int _D = dimension();\ 502 unsigned int _D = dimension();\
503 unsigned int _D_2 = _D/2;\ 503 unsigned int _D_2 = _D/2;\
504 if(2*_D_2!=_D) return false;\ 504 if(2*_D_2!=_D) return false;\
505 for(unsigned int r=0;r<_D_2;r++) {\ 505 for(unsigned int r=0;r<_D_2;r++) {\
506 for(unsigned int c=0;c<_D_2;c++) {\ 506 for(unsigned int c=0;c<_D_2;c++) {\
507 if(a_fabs(value(r,c)-value(r+_D_2,c+_D 507 if(a_fabs(value(r,c)-value(r+_D_2,c+_D_2))>=a_prec) return false;\
508 }\ 508 }\
509 }\ 509 }\
510 for(unsigned int r=0;r<_D_2;r++) {\ 510 for(unsigned int r=0;r<_D_2;r++) {\
511 for(unsigned int c=_D_2;c<_D;c++) {\ 511 for(unsigned int c=_D_2;c<_D;c++) {\
512 if(a_fabs(value(r,c)+value(r+_D_2,c-_D 512 if(a_fabs(value(r,c)+value(r+_D_2,c-_D_2))>=a_prec) return false;\
513 }\ 513 }\
514 }\ 514 }\
515 return true;\ 515 return true;\
516 }\ 516 }\
517 \ 517 \
518 template <class PREC>\ 518 template <class PREC>\
519 T determinant_prec(unsigned int a_tmp_rs[],u 519 T determinant_prec(unsigned int a_tmp_rs[],unsigned int a_tmp_cs[], /*[rord=dim-1]*/\
520 const PREC& a_prec,PREC(* 520 const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
521 unsigned int ord = dimension();\ 521 unsigned int ord = dimension();\
522 if(ord==0) {\ 522 if(ord==0) {\
523 return zero();\ 523 return zero();\
524 } else if(ord==1) {\ 524 } else if(ord==1) {\
525 return *m_vec;\ 525 return *m_vec;\
526 } else if(ord==2) {\ 526 } else if(ord==2) {\
527 T v00 = *m_vec;\ 527 T v00 = *m_vec;\
528 T v10 = *(m_vec+1);\ 528 T v10 = *(m_vec+1);\
529 T v01 = *(m_vec+2);\ 529 T v01 = *(m_vec+2);\
530 T v11 = *(m_vec+3);\ 530 T v11 = *(m_vec+3);\
531 return (v00 * v11 - v10 * v01);\ 531 return (v00 * v11 - v10 * v01);\
532 } else if(ord==3) {\ 532 } else if(ord==3) {\
533 /* 00 01 02 \ 533 /* 00 01 02 \
534 10 11 12 \ 534 10 11 12 \
535 20 21 22 \ 535 20 21 22 \
536 */\ 536 */\
537 T v00 = *m_vec;\ 537 T v00 = *m_vec;\
538 T v10 = *(m_vec+1);\ 538 T v10 = *(m_vec+1);\
539 T v20 = *(m_vec+2);\ 539 T v20 = *(m_vec+2);\
540 T v01 = *(m_vec+3);\ 540 T v01 = *(m_vec+3);\
541 T v11 = *(m_vec+4);\ 541 T v11 = *(m_vec+4);\
542 T v21 = *(m_vec+5);\ 542 T v21 = *(m_vec+5);\
543 T v02 = *(m_vec+6);\ 543 T v02 = *(m_vec+6);\
544 T v12 = *(m_vec+7);\ 544 T v12 = *(m_vec+7);\
545 T v22 = *(m_vec+8);\ 545 T v22 = *(m_vec+8);\
546 T cof_00 = v11 * v22 - v21 * v12;\ 546 T cof_00 = v11 * v22 - v21 * v12;\
547 T cof_10 = v01 * v22 - v21 * v02;\ 547 T cof_10 = v01 * v22 - v21 * v02;\
548 T cof_20 = v01 * v12 - v11 * v02;\ 548 T cof_20 = v01 * v12 - v11 * v02;\
549 return (v00*cof_00-v10*cof_10+v20*cof_20 549 return (v00*cof_00-v10*cof_10+v20*cof_20);\
550 }\ 550 }\
551 \ 551 \
552 unsigned int rord = ord-1;\ 552 unsigned int rord = ord-1;\
553 \ 553 \
554 T v_rc;\ 554 T v_rc;\
555 \ 555 \
556 T det = zero();\ 556 T det = zero();\
557 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs 557 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs[i] = i+1;}}\
558 unsigned int c = 0;\ 558 unsigned int c = 0;\
559 \ 559 \
560 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs 560 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs[i] = i+1;}}\
561 bool sg = true; /*c=0+r=0*/\ 561 bool sg = true; /*c=0+r=0*/\
562 for(unsigned int r=0;r<ord;r++) {\ 562 for(unsigned int r=0;r<ord;r++) {\
563 if(r>=1) a_tmp_rs[r-1] = r-1;\ 563 if(r>=1) a_tmp_rs[r-1] = r-1;\
564 v_rc = value(r,c);\ 564 v_rc = value(r,c);\
565 if(!is_zero(v_rc,a_prec,a_fabs)) {\ 565 if(!is_zero(v_rc,a_prec,a_fabs)) {\
566 T subdet = sub_determinant_prec(rord,a 566 T subdet = sub_determinant_prec(rord,a_tmp_rs,a_tmp_cs,a_prec,a_fabs);\
567 if(sg) \ 567 if(sg) \
568 det += v_rc * subdet;\ 568 det += v_rc * subdet;\
569 else\ 569 else\
570 det -= v_rc * subdet;\ 570 det -= v_rc * subdet;\
571 }\ 571 }\
572 sg = sg?false:true;\ 572 sg = sg?false:true;\
573 }\ 573 }\
574 \ 574 \
575 return det;\ 575 return det;\
576 }\ 576 }\
577 \ 577 \
578 T determinant(unsigned int a_tmp_rs[],unsign 578 T determinant(unsigned int a_tmp_rs[],unsigned int a_tmp_cs[]) const { /*[rord=dim-1]*/ \
579 return determinant_prec<double>(a_tmp_rs,a 579 return determinant_prec<double>(a_tmp_rs,a_tmp_cs,0,zero_fabs);\
580 }\ 580 }\
581 \ 581 \
582 T determinant() const {\ 582 T determinant() const {\
583 unsigned int ord = dimension();\ 583 unsigned int ord = dimension();\
584 if(ord==0) {\ 584 if(ord==0) {\
585 return zero();\ 585 return zero();\
586 } else if(ord==1) {\ 586 } else if(ord==1) {\
587 return *m_vec;\ 587 return *m_vec;\
588 } else if(ord==2) {\ 588 } else if(ord==2) {\
589 T v00 = *m_vec;\ 589 T v00 = *m_vec;\
590 T v10 = *(m_vec+1);\ 590 T v10 = *(m_vec+1);\
591 T v01 = *(m_vec+2);\ 591 T v01 = *(m_vec+2);\
592 T v11 = *(m_vec+3);\ 592 T v11 = *(m_vec+3);\
593 return (v00 * v11 - v10 * v01);\ 593 return (v00 * v11 - v10 * v01);\
594 } else if(ord==3) {\ 594 } else if(ord==3) {\
595 /* 00 01 02 \ 595 /* 00 01 02 \
596 10 11 12 \ 596 10 11 12 \
597 20 21 22 \ 597 20 21 22 \
598 */\ 598 */\
599 T v00 = *m_vec;\ 599 T v00 = *m_vec;\
600 T v10 = *(m_vec+1);\ 600 T v10 = *(m_vec+1);\
601 T v20 = *(m_vec+2);\ 601 T v20 = *(m_vec+2);\
602 T v01 = *(m_vec+3);\ 602 T v01 = *(m_vec+3);\
603 T v11 = *(m_vec+4);\ 603 T v11 = *(m_vec+4);\
604 T v21 = *(m_vec+5);\ 604 T v21 = *(m_vec+5);\
605 T v02 = *(m_vec+6);\ 605 T v02 = *(m_vec+6);\
606 T v12 = *(m_vec+7);\ 606 T v12 = *(m_vec+7);\
607 T v22 = *(m_vec+8);\ 607 T v22 = *(m_vec+8);\
608 T cof_00 = v11 * v22 - v21 * v12;\ 608 T cof_00 = v11 * v22 - v21 * v12;\
609 T cof_10 = v01 * v22 - v21 * v02;\ 609 T cof_10 = v01 * v22 - v21 * v02;\
610 T cof_20 = v01 * v12 - v11 * v02;\ 610 T cof_20 = v01 * v12 - v11 * v02;\
611 return (v00*cof_00-v10*cof_10+v20*cof_20 611 return (v00*cof_00-v10*cof_10+v20*cof_20);\
612 }\ 612 }\
613 unsigned int rord = ord-1;\ 613 unsigned int rord = ord-1;\
614 unsigned int* rs = new unsigned int[rord]; 614 unsigned int* rs = new unsigned int[rord];\
615 unsigned int* cs = new unsigned int[rord]; 615 unsigned int* cs = new unsigned int[rord];\
616 T det = determinant(rs,cs);\ 616 T det = determinant(rs,cs);\
617 delete [] rs;\ 617 delete [] rs;\
618 delete [] cs;\ 618 delete [] cs;\
619 return det;\ 619 return det;\
620 }\ 620 }\
621 \ 621 \
622 template <class PREC>\ 622 template <class PREC>\
623 bool invert_prec(TOOLS_MAT_CLASS& a_res,\ 623 bool invert_prec(TOOLS_MAT_CLASS& a_res,\
624 unsigned int a_tmp_rs[],uns 624 unsigned int a_tmp_rs[],unsigned int a_tmp_cs[], /*[rord=dim-1]*/\
625 const PREC& a_prec,PREC(*a_ 625 const PREC& a_prec,PREC(*a_fabs)(const T&) /*for det=?zero*/ \
626 ) const { \ 626 ) const { \
627 unsigned int ord = dimension();\ 627 unsigned int ord = dimension();\
628 if(ord==0) return true;\ 628 if(ord==0) return true;\
629 \ 629 \
630 if(ord==1) {\ 630 if(ord==1) {\
631 T det = value(0,0);\ 631 T det = value(0,0);\
632 if(is_zero(det,a_prec,a_fabs)) return fa 632 if(is_zero(det,a_prec,a_fabs)) return false;\
633 a_res.set_value(0,0,one()/det);\ 633 a_res.set_value(0,0,one()/det);\
634 return true;\ 634 return true;\
635 } else if(ord==2) {\ 635 } else if(ord==2) {\
636 T v00 = *m_vec;\ 636 T v00 = *m_vec;\
637 T v10 = *(m_vec+1);\ 637 T v10 = *(m_vec+1);\
638 T v01 = *(m_vec+2);\ 638 T v01 = *(m_vec+2);\
639 T v11 = *(m_vec+3);\ 639 T v11 = *(m_vec+3);\
640 T det = (v00 * v11 - v10 * v01);\ 640 T det = (v00 * v11 - v10 * v01);\
641 if(is_zero(det,a_prec,a_fabs)) return fa 641 if(is_zero(det,a_prec,a_fabs)) return false;\
642 a_res.set_value(0,0,v11/det);\ 642 a_res.set_value(0,0,v11/det);\
643 a_res.set_value(1,1,v00/det);\ 643 a_res.set_value(1,1,v00/det);\
644 a_res.set_value(0,1,minus_one()*v01/det) 644 a_res.set_value(0,1,minus_one()*v01/det);\
645 a_res.set_value(1,0,minus_one()*v10/det) 645 a_res.set_value(1,0,minus_one()*v10/det);\
646 return true;\ 646 return true;\
647 } else if(ord==3) {\ 647 } else if(ord==3) {\
648 /* 00 01 02 \ 648 /* 00 01 02 \
649 10 11 12 \ 649 10 11 12 \
650 20 21 22 \ 650 20 21 22 \
651 */\ 651 */\
652 T v00 = *m_vec;\ 652 T v00 = *m_vec;\
653 T v10 = *(m_vec+1);\ 653 T v10 = *(m_vec+1);\
654 T v20 = *(m_vec+2);\ 654 T v20 = *(m_vec+2);\
655 T v01 = *(m_vec+3);\ 655 T v01 = *(m_vec+3);\
656 T v11 = *(m_vec+4);\ 656 T v11 = *(m_vec+4);\
657 T v21 = *(m_vec+5);\ 657 T v21 = *(m_vec+5);\
658 T v02 = *(m_vec+6);\ 658 T v02 = *(m_vec+6);\
659 T v12 = *(m_vec+7);\ 659 T v12 = *(m_vec+7);\
660 T v22 = *(m_vec+8);\ 660 T v22 = *(m_vec+8);\
661 T cof_00 = v11 * v22 - v21 * v12;\ 661 T cof_00 = v11 * v22 - v21 * v12;\
662 T cof_10 = v01 * v22 - v21 * v02;\ 662 T cof_10 = v01 * v22 - v21 * v02;\
663 T cof_20 = v01 * v12 - v11 * v02;\ 663 T cof_20 = v01 * v12 - v11 * v02;\
664 T det = (v00*cof_00-v10*cof_10+v20*cof_2 664 T det = (v00*cof_00-v10*cof_10+v20*cof_20);\
665 if(is_zero(det,a_prec,a_fabs)) return fa 665 if(is_zero(det,a_prec,a_fabs)) return false;\
666 T cof_01 = v10 * v22 - v20 * v12;\ 666 T cof_01 = v10 * v22 - v20 * v12;\
667 T cof_11 = v00 * v22 - v20 * v02;\ 667 T cof_11 = v00 * v22 - v20 * v02;\
668 T cof_21 = v00 * v12 - v10 * v02;\ 668 T cof_21 = v00 * v12 - v10 * v02;\
669 T cof_02 = v10 * v21 - v20 * v11;\ 669 T cof_02 = v10 * v21 - v20 * v11;\
670 T cof_12 = v00 * v21 - v20 * v01;\ 670 T cof_12 = v00 * v21 - v20 * v01;\
671 T cof_22 = v00 * v11 - v10 * v01;\ 671 T cof_22 = v00 * v11 - v10 * v01;\
672 a_res.set_value(0,0,cof_00/det);\ 672 a_res.set_value(0,0,cof_00/det);\
673 a_res.set_value(1,0,minus_one()*cof_01/d 673 a_res.set_value(1,0,minus_one()*cof_01/det);\
674 a_res.set_value(2,0,cof_02/det);\ 674 a_res.set_value(2,0,cof_02/det);\
675 a_res.set_value(0,1,minus_one()*cof_10/d 675 a_res.set_value(0,1,minus_one()*cof_10/det);\
676 a_res.set_value(1,1,cof_11/det);\ 676 a_res.set_value(1,1,cof_11/det);\
677 a_res.set_value(2,1,minus_one()*cof_12/d 677 a_res.set_value(2,1,minus_one()*cof_12/det);\
678 a_res.set_value(0,2,cof_20/det);\ 678 a_res.set_value(0,2,cof_20/det);\
679 a_res.set_value(1,2,minus_one()*cof_21/d 679 a_res.set_value(1,2,minus_one()*cof_21/det);\
680 a_res.set_value(2,2,cof_22/det);\ 680 a_res.set_value(2,2,cof_22/det);\
681 return true;\ 681 return true;\
682 }\ 682 }\
683 \ 683 \
684 /*Generic invertion method.*/\ 684 /*Generic invertion method.*/\
685 \ 685 \
686 unsigned int rord = ord-1;\ 686 unsigned int rord = ord-1;\
687 \ 687 \
688 /* Get det with r = 0;*/\ 688 /* Get det with r = 0;*/\
689 T det = zero();\ 689 T det = zero();\
690 T subdet;\ 690 T subdet;\
691 {\ 691 {\
692 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs 692 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs[i] = i+1;}}\
693 unsigned int r = 0;\ 693 unsigned int r = 0;\
694 /*if(r>=1) a_tmp_rs[r-1] = r-1;*/\ 694 /*if(r>=1) a_tmp_rs[r-1] = r-1;*/\
695 \ 695 \
696 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs 696 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs[i] = i+1;}}\
697 bool sg = true; /*r=0+c=0*/\ 697 bool sg = true; /*r=0+c=0*/\
698 for(unsigned int c=0;c<ord;c++) {\ 698 for(unsigned int c=0;c<ord;c++) {\
699 if(c>=1) a_tmp_cs[c-1] = c-1;\ 699 if(c>=1) a_tmp_cs[c-1] = c-1;\
700 subdet = sub_determinant_prec(rord,a_tmp 700 subdet = sub_determinant_prec(rord,a_tmp_rs,a_tmp_cs,a_prec,a_fabs);\
701 if(sg) {\ 701 if(sg) {\
702 det += value(r,c) * subdet;\ 702 det += value(r,c) * subdet;\
703 a_res.set_value(c,r,subdet);\ 703 a_res.set_value(c,r,subdet);\
704 sg = false;\ 704 sg = false;\
705 } else {\ 705 } else {\
706 det += value(r,c) * subdet * minus_one 706 det += value(r,c) * subdet * minus_one();\
707 a_res.set_value(c,r,subdet * minus_one 707 a_res.set_value(c,r,subdet * minus_one());\
708 sg = true;\ 708 sg = true;\
709 }\ 709 }\
710 }}\ 710 }}\
711 \ 711 \
712 if(is_zero(det,a_prec,a_fabs)) return false 712 if(is_zero(det,a_prec,a_fabs)) return false;\
713 \ 713 \
714 {for(unsigned int c=0;c<ord;c++) {\ 714 {for(unsigned int c=0;c<ord;c++) {\
715 a_res.set_value(c,0,a_res.value(c,0)/det 715 a_res.set_value(c,0,a_res.value(c,0)/det);\
716 }}\ 716 }}\
717 \ 717 \
718 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs 718 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs[i] = i+1;}}\
719 bool sgr = false; /*r=1+c=0*/\ 719 bool sgr = false; /*r=1+c=0*/\
720 for(unsigned int r=1;r<ord;r++) {\ 720 for(unsigned int r=1;r<ord;r++) {\
721 if(r>=1) a_tmp_rs[r-1] = r-1;\ 721 if(r>=1) a_tmp_rs[r-1] = r-1;\
722 {for(unsigned int i=0;i<rord;i++) {a_tmp_ 722 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs[i] = i+1;}}\
723 bool sg = sgr;\ 723 bool sg = sgr;\
724 for(unsigned int c=0;c<ord;c++) {\ 724 for(unsigned int c=0;c<ord;c++) {\
725 if(c>=1) a_tmp_cs[c-1] = c-1;\ 725 if(c>=1) a_tmp_cs[c-1] = c-1;\
726 subdet = sub_determinant_prec(rord,a_t 726 subdet = sub_determinant_prec(rord,a_tmp_rs,a_tmp_cs,a_prec,a_fabs);\
727 if(sg) {\ 727 if(sg) {\
728 a_res.set_value(c,r,subdet/det);\ 728 a_res.set_value(c,r,subdet/det);\
729 sg = false;\ 729 sg = false;\
730 } else {\ 730 } else {\
731 a_res.set_value(c,r,(subdet * minus_ 731 a_res.set_value(c,r,(subdet * minus_one())/det);\
732 sg = true;\ 732 sg = true;\
733 }\ 733 }\
734 }\ 734 }\
735 sgr = sgr?false:true;\ 735 sgr = sgr?false:true;\
736 }\ 736 }\
737 \ 737 \
738 return true;\ 738 return true;\
739 }\ 739 }\
740 \ 740 \
741 template <class PREC>\ 741 template <class PREC>\
742 bool invert_prec(TOOLS_MAT_CLASS& a_res,cons 742 bool invert_prec(TOOLS_MAT_CLASS& a_res,const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
743 unsigned int ord = dimension();\ 743 unsigned int ord = dimension();\
744 if(ord==0) return true;\ 744 if(ord==0) return true;\
745 \ 745 \
746 if(ord==1) {\ 746 if(ord==1) {\
747 T det = value(0,0);\ 747 T det = value(0,0);\
748 if(is_zero(det,a_prec,a_fabs)) return fa 748 if(is_zero(det,a_prec,a_fabs)) return false;\
749 a_res.set_value(0,0,one()/det);\ 749 a_res.set_value(0,0,one()/det);\
750 return true;\ 750 return true;\
751 } else if(ord==2) {\ 751 } else if(ord==2) {\
752 T v00 = *m_vec;\ 752 T v00 = *m_vec;\
753 T v10 = *(m_vec+1);\ 753 T v10 = *(m_vec+1);\
754 T v01 = *(m_vec+2);\ 754 T v01 = *(m_vec+2);\
755 T v11 = *(m_vec+3);\ 755 T v11 = *(m_vec+3);\
756 T det = (v00 * v11 - v10 * v01);\ 756 T det = (v00 * v11 - v10 * v01);\
757 if(is_zero(det,a_prec,a_fabs)) return fa 757 if(is_zero(det,a_prec,a_fabs)) return false;\
758 a_res.set_value(0,0,v11/det);\ 758 a_res.set_value(0,0,v11/det);\
759 a_res.set_value(1,1,v00/det);\ 759 a_res.set_value(1,1,v00/det);\
760 a_res.set_value(0,1,minus_one()*v01/det) 760 a_res.set_value(0,1,minus_one()*v01/det);\
761 a_res.set_value(1,0,minus_one()*v10/det) 761 a_res.set_value(1,0,minus_one()*v10/det);\
762 return true;\ 762 return true;\
763 } else if(ord==3) {\ 763 } else if(ord==3) {\
764 /* 00 01 02 \ 764 /* 00 01 02 \
765 10 11 12 \ 765 10 11 12 \
766 20 21 22 \ 766 20 21 22 \
767 */\ 767 */\
768 T v00 = *m_vec;\ 768 T v00 = *m_vec;\
769 T v10 = *(m_vec+1);\ 769 T v10 = *(m_vec+1);\
770 T v20 = *(m_vec+2);\ 770 T v20 = *(m_vec+2);\
771 T v01 = *(m_vec+3);\ 771 T v01 = *(m_vec+3);\
772 T v11 = *(m_vec+4);\ 772 T v11 = *(m_vec+4);\
773 T v21 = *(m_vec+5);\ 773 T v21 = *(m_vec+5);\
774 T v02 = *(m_vec+6);\ 774 T v02 = *(m_vec+6);\
775 T v12 = *(m_vec+7);\ 775 T v12 = *(m_vec+7);\
776 T v22 = *(m_vec+8);\ 776 T v22 = *(m_vec+8);\
777 T cof_00 = v11 * v22 - v21 * v12;\ 777 T cof_00 = v11 * v22 - v21 * v12;\
778 T cof_10 = v01 * v22 - v21 * v02;\ 778 T cof_10 = v01 * v22 - v21 * v02;\
779 T cof_20 = v01 * v12 - v11 * v02;\ 779 T cof_20 = v01 * v12 - v11 * v02;\
780 T det = (v00*cof_00-v10*cof_10+v20*cof_2 780 T det = (v00*cof_00-v10*cof_10+v20*cof_20);\
781 if(is_zero(det,a_prec,a_fabs)) return fa 781 if(is_zero(det,a_prec,a_fabs)) return false;\
782 T cof_01 = v10 * v22 - v20 * v12;\ 782 T cof_01 = v10 * v22 - v20 * v12;\
783 T cof_11 = v00 * v22 - v20 * v02;\ 783 T cof_11 = v00 * v22 - v20 * v02;\
784 T cof_21 = v00 * v12 - v10 * v02;\ 784 T cof_21 = v00 * v12 - v10 * v02;\
785 T cof_02 = v10 * v21 - v20 * v11;\ 785 T cof_02 = v10 * v21 - v20 * v11;\
786 T cof_12 = v00 * v21 - v20 * v01;\ 786 T cof_12 = v00 * v21 - v20 * v01;\
787 T cof_22 = v00 * v11 - v10 * v01;\ 787 T cof_22 = v00 * v11 - v10 * v01;\
788 a_res.set_value(0,0,cof_00/det);\ 788 a_res.set_value(0,0,cof_00/det);\
789 a_res.set_value(1,0,minus_one()*cof_01/d 789 a_res.set_value(1,0,minus_one()*cof_01/det);\
790 a_res.set_value(2,0,cof_02/det);\ 790 a_res.set_value(2,0,cof_02/det);\
791 a_res.set_value(0,1,minus_one()*cof_10/d 791 a_res.set_value(0,1,minus_one()*cof_10/det);\
792 a_res.set_value(1,1,cof_11/det);\ 792 a_res.set_value(1,1,cof_11/det);\
793 a_res.set_value(2,1,minus_one()*cof_12/d 793 a_res.set_value(2,1,minus_one()*cof_12/det);\
794 a_res.set_value(0,2,cof_20/det);\ 794 a_res.set_value(0,2,cof_20/det);\
795 a_res.set_value(1,2,minus_one()*cof_21/d 795 a_res.set_value(1,2,minus_one()*cof_21/det);\
796 a_res.set_value(2,2,cof_22/det);\ 796 a_res.set_value(2,2,cof_22/det);\
797 return true;\ 797 return true;\
798 }\ 798 }\
799 \ 799 \
800 unsigned int rord = ord-1;\ 800 unsigned int rord = ord-1;\
801 unsigned int* cs = new unsigned int[rord]; 801 unsigned int* cs = new unsigned int[rord];\
802 unsigned int* rs = new unsigned int[rord]; 802 unsigned int* rs = new unsigned int[rord];\
803 bool status = invert_prec(a_res,rs,cs,a_pr 803 bool status = invert_prec(a_res,rs,cs,a_prec,a_fabs);\
804 delete [] cs;\ 804 delete [] cs;\
805 delete [] rs;\ 805 delete [] rs;\
806 return status;\ 806 return status;\
807 }\ 807 }\
808 \ 808 \
809 bool invert(TOOLS_MAT_CLASS& a_res,unsigned 809 bool invert(TOOLS_MAT_CLASS& a_res,unsigned int a_tmp_rs[],unsigned int a_tmp_cs[]) const { /*[rord=dim-1]*/ \
810 return invert_prec<double>(a_res,a_tmp_rs, 810 return invert_prec<double>(a_res,a_tmp_rs,a_tmp_cs,0,zero_fabs);\
811 }\ 811 }\
812 \ 812 \
813 bool invert(TOOLS_MAT_CLASS& a_res) const {\ 813 bool invert(TOOLS_MAT_CLASS& a_res) const {\
814 return invert_prec<double>(a_res,0,zero_fa 814 return invert_prec<double>(a_res,0,zero_fabs);\
815 }\ 815 }\
816 \ 816 \
817 void power(unsigned int a_n,TOOLS_MAT_CLASS& 817 void power(unsigned int a_n,TOOLS_MAT_CLASS& a_res) const {\
818 a_res.set_identity();\ 818 a_res.set_identity();\
819 T* _tmp = new T[dim2()];\ 819 T* _tmp = new T[dim2()];\
820 for(unsigned int i=0;i<a_n;i++) {\ 820 for(unsigned int i=0;i<a_n;i++) {\
821 a_res._mul_mtx(m_vec,_tmp);\ 821 a_res._mul_mtx(m_vec,_tmp);\
822 }\ 822 }\
823 delete [] _tmp;\ 823 delete [] _tmp;\
824 }\ 824 }\
825 \ 825 \
826 void exp(unsigned int a_le,TOOLS_MAT_CLASS& 826 void exp(unsigned int a_le,TOOLS_MAT_CLASS& a_res,TOOLS_MAT_CLASS& a_tmp,T a_tmp_2[]) const { /*OPTIMIZATION*/\
827 /* result = I + M + M**2/2! + M**3/3! + .. 827 /* result = I + M + M**2/2! + M**3/3! + .... */\
828 a_res.set_identity();\ 828 a_res.set_identity();\
829 a_tmp.set_identity();\ 829 a_tmp.set_identity();\
830 for(unsigned int i=1;i<=a_le;i++) {\ 830 for(unsigned int i=1;i<=a_le;i++) {\
831 a_tmp._mul_mtx(m_vec,a_tmp_2);\ 831 a_tmp._mul_mtx(m_vec,a_tmp_2);\
832 a_tmp.multiply(one()/T(i));\ 832 a_tmp.multiply(one()/T(i));\
833 a_res += a_tmp;\ 833 a_res += a_tmp;\
834 }\ 834 }\
835 }\ 835 }\
836 \ 836 \
837 void exp(unsigned int a_le,TOOLS_MAT_CLASS& 837 void exp(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
838 /* result = I + M + M**2/2! + M**3/3! + .. 838 /* result = I + M + M**2/2! + M**3/3! + .... */\
839 TOOLS_MAT_CLASS tmp(*this);\ 839 TOOLS_MAT_CLASS tmp(*this);\
840 tmp.set_identity();\ 840 tmp.set_identity();\
841 T* tmp_2 = new T[dim2()];\ 841 T* tmp_2 = new T[dim2()];\
842 exp(a_le,a_res,tmp,tmp_2);\ 842 exp(a_le,a_res,tmp,tmp_2);\
843 delete [] tmp_2;\ 843 delete [] tmp_2;\
844 }\ 844 }\
845 \ 845 \
846 void cosh(unsigned int a_le,TOOLS_MAT_CLASS& 846 void cosh(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
847 /* result = I + M**2/2! + M**4/4! + .... * 847 /* result = I + M**2/2! + M**4/4! + .... */\
848 a_res.set_identity();\ 848 a_res.set_identity();\
849 TOOLS_MAT_CLASS M_2(*this);\ 849 TOOLS_MAT_CLASS M_2(*this);\
850 M_2._mul_mtx(m_vec);\ 850 M_2._mul_mtx(m_vec);\
851 TOOLS_MAT_CLASS tmp(*this);\ 851 TOOLS_MAT_CLASS tmp(*this);\
852 tmp.set_identity();\ 852 tmp.set_identity();\
853 T* _tmp = new T[dim2()];\ 853 T* _tmp = new T[dim2()];\
854 for(unsigned int i=1;i<=a_le;i+=2) {\ 854 for(unsigned int i=1;i<=a_le;i+=2) {\
855 tmp._mul_mtx(M_2.m_vec,_tmp);\ 855 tmp._mul_mtx(M_2.m_vec,_tmp);\
856 tmp.multiply(one()/T(i+0)); \ 856 tmp.multiply(one()/T(i+0)); \
857 tmp.multiply(one()/T(i+1)); \ 857 tmp.multiply(one()/T(i+1)); \
858 a_res += tmp;\ 858 a_res += tmp;\
859 }\ 859 }\
860 delete [] _tmp;\ 860 delete [] _tmp;\
861 }\ 861 }\
862 \ 862 \
863 void sinh(unsigned int a_le,TOOLS_MAT_CLASS& 863 void sinh(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
864 /* result = M + M**3/3! + .... */\ 864 /* result = M + M**3/3! + .... */\
865 a_res._copy(m_vec);\ 865 a_res._copy(m_vec);\
866 TOOLS_MAT_CLASS M_2(*this);\ 866 TOOLS_MAT_CLASS M_2(*this);\
867 M_2._mul_mtx(m_vec);\ 867 M_2._mul_mtx(m_vec);\
868 TOOLS_MAT_CLASS tmp(*this);\ 868 TOOLS_MAT_CLASS tmp(*this);\
869 tmp._copy(m_vec);\ 869 tmp._copy(m_vec);\
870 T* _tmp = new T[dim2()];\ 870 T* _tmp = new T[dim2()];\
871 for(unsigned int i=1;i<=a_le;i+=2) {\ 871 for(unsigned int i=1;i<=a_le;i+=2) {\
872 tmp._mul_mtx(M_2.m_vec,_tmp);\ 872 tmp._mul_mtx(M_2.m_vec,_tmp);\
873 tmp.multiply(one()/T(i+1)); \ 873 tmp.multiply(one()/T(i+1)); \
874 tmp.multiply(one()/T(i+2)); \ 874 tmp.multiply(one()/T(i+2)); \
875 a_res += tmp;\ 875 a_res += tmp;\
876 }\ 876 }\
877 delete [] _tmp;\ 877 delete [] _tmp;\
878 }\ 878 }\
879 \ 879 \
880 void log(unsigned int a_le,TOOLS_MAT_CLASS& 880 void log(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
881 /* result = (M-I) - (M-I)**2/2 + (M-I)**3/ 881 /* result = (M-I) - (M-I)**2/2 + (M-I)**3/3 +... */\
882 /* WARNING : touchy, it may not converge ! 882 /* WARNING : touchy, it may not converge ! */\
883 a_res.set_zero();\ 883 a_res.set_zero();\
884 \ 884 \
885 TOOLS_MAT_CLASS M_I;\ 885 TOOLS_MAT_CLASS M_I;\
886 M_I.set_identity();\ 886 M_I.set_identity();\
887 M_I.multiply(minus_one());\ 887 M_I.multiply(minus_one());\
888 M_I._add_mtx(m_vec);\ 888 M_I._add_mtx(m_vec);\
889 \ 889 \
890 TOOLS_MAT_CLASS M_Ip(M_I);\ 890 TOOLS_MAT_CLASS M_Ip(M_I);\
891 T fact = minus_one();\ 891 T fact = minus_one();\
892 \ 892 \
893 TOOLS_MAT_CLASS tmp;\ 893 TOOLS_MAT_CLASS tmp;\
894 \ 894 \
895 T* _tmp = new T[dim2()];\ 895 T* _tmp = new T[dim2()];\
896 for(unsigned int i=0;i<=a_le;i++) {\ 896 for(unsigned int i=0;i<=a_le;i++) {\
897 fact *= minus_one(); \ 897 fact *= minus_one(); \
898 tmp = M_Ip;\ 898 tmp = M_Ip;\
899 tmp.multiply(fact/T(i+1)); \ 899 tmp.multiply(fact/T(i+1)); \
900 a_res += tmp;\ 900 a_res += tmp;\
901 M_Ip._mul_mtx(M_I.m_vec,_tmp);\ 901 M_Ip._mul_mtx(M_I.m_vec,_tmp);\
902 }\ 902 }\
903 delete [] _tmp;\ 903 delete [] _tmp;\
904 }\ 904 }\
905 \ 905 \
906 void omega(unsigned int a_le,TOOLS_MAT_CLASS 906 void omega(unsigned int a_le,TOOLS_MAT_CLASS& a_res,TOOLS_MAT_CLASS& a_tmp,T a_tmp_2[]) const { /*OPTIMIZATION*/\
907 /* result = I + M/2! + M**2/3! + M**3/4! + 907 /* result = I + M/2! + M**2/3! + M**3/4! + .... */\
908 a_res.set_identity();\ 908 a_res.set_identity();\
909 a_tmp.set_identity();\ 909 a_tmp.set_identity();\
910 for(unsigned int i=1;i<=a_le;i++) {\ 910 for(unsigned int i=1;i<=a_le;i++) {\
911 a_tmp._mul_mtx(m_vec,a_tmp_2);\ 911 a_tmp._mul_mtx(m_vec,a_tmp_2);\
912 a_tmp.multiply(one()/T(i+1));\ 912 a_tmp.multiply(one()/T(i+1));\
913 a_res += a_tmp;\ 913 a_res += a_tmp;\
914 }\ 914 }\
915 }\ 915 }\
916 \ 916 \
917 void omega(unsigned int a_le,TOOLS_MAT_CLASS 917 void omega(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
918 /* result = I + M/2! + M**2/3! + M**3/4! + 918 /* result = I + M/2! + M**2/3! + M**3/4! + .... */\
919 TOOLS_MAT_CLASS tmp(*this);\ 919 TOOLS_MAT_CLASS tmp(*this);\
920 tmp.set_identity();\ 920 tmp.set_identity();\
921 T* tmp_2 = new T[dim2()];\ 921 T* tmp_2 = new T[dim2()];\
922 omega(a_le,a_res,tmp,tmp_2);\ 922 omega(a_le,a_res,tmp,tmp_2);\
923 delete [] tmp_2;\ 923 delete [] tmp_2;\
924 }\ 924 }\
925 \ 925 \
926 template <class MAT>\ 926 template <class MAT>\
927 bool copy(const MAT& a_from) {\ 927 bool copy(const MAT& a_from) {\
928 /*for exa from a double matrix to a symbol 928 /*for exa from a double matrix to a symbol matrix*/\
929 unsigned int _D = dimension();\ 929 unsigned int _D = dimension();\
930 if(a_from.dimension()!=_D) return false;\ 930 if(a_from.dimension()!=_D) return false;\
931 for(unsigned int r=0;r<_D;r++) {\ 931 for(unsigned int r=0;r<_D;r++) {\
932 for(unsigned int c=0;c<_D;c++) {\ 932 for(unsigned int c=0;c<_D;c++) {\
933 set_value(r,c,a_from.value(r,c));\ 933 set_value(r,c,a_from.value(r,c));\
934 }\ 934 }\
935 }\ 935 }\
936 return true;\ 936 return true;\
937 }\ 937 }\
938 public: /*operators*/\ 938 public: /*operators*/\
939 T operator()(unsigned int a_r,unsigned int a 939 T operator()(unsigned int a_r,unsigned int a_c) const {\
940 /*WARNING : no check on a_r,a_c.*/\ 940 /*WARNING : no check on a_r,a_c.*/\
941 return m_vec[a_r + a_c * dimension()];\ 941 return m_vec[a_r + a_c * dimension()];\
942 }\ 942 }\
943 \ 943 \
944 T& operator[](size_t a_index) { /*for tools/ 944 T& operator[](size_t a_index) { /*for tools/sg/sf_vec*/\
945 /*WARNING : no check on a_index.*/\ 945 /*WARNING : no check on a_index.*/\
946 return m_vec[a_index];\ 946 return m_vec[a_index];\
947 }\ 947 }\
948 const T& operator[](size_t a_index) const {\ 948 const T& operator[](size_t a_index) const {\
949 /*WARNING : no check on a_index.*/\ 949 /*WARNING : no check on a_index.*/\
950 return m_vec[a_index];\ 950 return m_vec[a_index];\
951 }\ 951 }\
952 bool operator==(const TOOLS_MAT_CLASS& a_arr 952 bool operator==(const TOOLS_MAT_CLASS& a_array) const {\
953 return equal(a_array);\ 953 return equal(a_array);\
954 }\ 954 }\
955 bool operator!=(const TOOLS_MAT_CLASS& a_arr 955 bool operator!=(const TOOLS_MAT_CLASS& a_array) const {\
956 return !operator==(a_array);\ 956 return !operator==(a_array);\
957 }\ 957 }\
958 TOOLS_MAT_CLASS& operator*=(const TOOLS_MAT_ 958 TOOLS_MAT_CLASS& operator*=(const TOOLS_MAT_CLASS& a_m) {\
959 _mul_mtx(a_m.m_vec);\ 959 _mul_mtx(a_m.m_vec);\
960 return *this;\ 960 return *this;\
961 }\ 961 }\
962 TOOLS_MAT_CLASS& operator+=(const TOOLS_MAT_ 962 TOOLS_MAT_CLASS& operator+=(const TOOLS_MAT_CLASS& a_m) {\
963 _add_mtx(a_m.m_vec);\ 963 _add_mtx(a_m.m_vec);\
964 return *this;\ 964 return *this;\
965 }\ 965 }\
966 TOOLS_MAT_CLASS& operator-=(const TOOLS_MAT_ 966 TOOLS_MAT_CLASS& operator-=(const TOOLS_MAT_CLASS& a_m) {\
967 _sub_mtx(a_m.m_vec);\ 967 _sub_mtx(a_m.m_vec);\
968 return *this;\ 968 return *this;\
969 }\ 969 }\
970 TOOLS_MAT_CLASS& operator*=(const T& a_fac) 970 TOOLS_MAT_CLASS& operator*=(const T& a_fac) {\
971 for(unsigned int i=0;i<dim2();i++) m_vec[i 971 for(unsigned int i=0;i<dim2();i++) m_vec[i] *= a_fac;\
972 return *this;\ 972 return *this;\
973 }\ 973 }\
974 \ 974 \
975 TOOLS_MAT_CLASS operator*(const T& a_fac) {\ 975 TOOLS_MAT_CLASS operator*(const T& a_fac) {\
976 TOOLS_MAT_CLASS res;\ 976 TOOLS_MAT_CLASS res;\
977 res.operator*=(a_fac);\ 977 res.operator*=(a_fac);\
978 return res;\ 978 return res;\
979 }\ 979 }\
980 protected:\ 980 protected:\
981 void _copy(const T a_m[]) {\ 981 void _copy(const T a_m[]) {\
982 T* tp = (T*)m_vec;T* ap = (T*)a_m;\ 982 T* tp = (T*)m_vec;T* ap = (T*)a_m;\
983 for(unsigned int i=0;i<dim2();i++,tp++,ap+ 983 for(unsigned int i=0;i<dim2();i++,tp++,ap++) *tp = *ap;\
984 /*{for(unsigned int i=0;i<dim2();i++) m_vec 984 /*{for(unsigned int i=0;i<dim2();i++) m_vec[i] = a_m[i];}*/\
985 /* memcpy does not work with std::complex<> 985 /* memcpy does not work with std::complex<> and mat<symbol,4> see tools/tests/symbolic.cpp */\
986 /*vec_copy(m_vec,a_m,dim2());*/\ 986 /*vec_copy(m_vec,a_m,dim2());*/\
987 }\ 987 }\
988 \ 988 \
989 void _add_mtx(const T a_m[]) { /* this = th 989 void _add_mtx(const T a_m[]) { /* this = this + a_m, */\
990 T* tp = (T*)m_vec;T* ap = (T*)a_m;\ 990 T* tp = (T*)m_vec;T* ap = (T*)a_m;\
991 for(unsigned int i=0;i<dim2();i++,tp++,ap+ 991 for(unsigned int i=0;i<dim2();i++,tp++,ap++) *tp += *ap;\
992 /*for(unsigned int i=0;i<dim2();i++) m_vec 992 /*for(unsigned int i=0;i<dim2();i++) m_vec[i] += a_m[i];*/\
993 /*vec_add(m_vec,a_m,m_vec,dim2());*/\ 993 /*vec_add(m_vec,a_m,m_vec,dim2());*/\
994 }\ 994 }\
995 void _sub_mtx(const T a_m[]) { /* this = th 995 void _sub_mtx(const T a_m[]) { /* this = this - a_m, */\
996 T* tp = (T*)m_vec;T* ap = (T*)a_m;\ 996 T* tp = (T*)m_vec;T* ap = (T*)a_m;\
997 for(unsigned int i=0;i<dim2();i++,tp++,ap+ 997 for(unsigned int i=0;i<dim2();i++,tp++,ap++) *tp -= *ap;\
998 /*for(unsigned int i=0;i<dim2();i++) m_vec 998 /*for(unsigned int i=0;i<dim2();i++) m_vec[i] -= a_m[i];*/\
999 /*vec_sub(m_vec,a_m,m_vec,dim2());*/\ 999 /*vec_sub(m_vec,a_m,m_vec,dim2());*/\
1000 }\ 1000 }\
1001 \ 1001 \
1002 /*\ 1002 /*\
1003 void _mul_mtx(const T a_m[],T a_tmp[]) {\ 1003 void _mul_mtx(const T a_m[],T a_tmp[]) {\
1004 unsigned int ord = dimension();\ 1004 unsigned int ord = dimension();\
1005 for(unsigned int r=0;r<ord;r++) {\ 1005 for(unsigned int r=0;r<ord;r++) {\
1006 for(unsigned int c=0;c<ord;c++) {\ 1006 for(unsigned int c=0;c<ord;c++) {\
1007 T _value = zero();\ 1007 T _value = zero();\
1008 for(unsigned int i=0;i<ord;i++) {\ 1008 for(unsigned int i=0;i<ord;i++) {\
1009 _value += (*(m_vec+r+i*ord)) * (*(a 1009 _value += (*(m_vec+r+i*ord)) * (*(a_m+i+c*ord)); //optimize.\
1010 }\ 1010 }\
1011 *(a_tmp+r+c*ord) = _value;\ 1011 *(a_tmp+r+c*ord) = _value;\
1012 }\ 1012 }\
1013 }\ 1013 }\
1014 _copy(a_tmp);\ 1014 _copy(a_tmp);\
1015 }\ 1015 }\
1016 */\ 1016 */\
1017 void _mul_mtx(const T a_m[],T a_tmp[]) { /* 1017 void _mul_mtx(const T a_m[],T a_tmp[]) { /*OPTIMIZATION*/\
1018 /* this = this * a_m */\ 1018 /* this = this * a_m */\
1019 typedef T* Tp;\ 1019 typedef T* Tp;\
1020 Tp tpos,ttpos,rpos,apos,mpos,aapos;\ 1020 Tp tpos,ttpos,rpos,apos,mpos,aapos;\
1021 T _value;\ 1021 T _value;\
1022 unsigned int r,c,i;\ 1022 unsigned int r,c,i;\
1023 \ 1023 \
1024 unsigned int _D = dimension();\ 1024 unsigned int _D = dimension();\
1025 \ 1025 \
1026 tpos = a_tmp;\ 1026 tpos = a_tmp;\
1027 for(r=0;r<_D;r++,tpos++) {\ 1027 for(r=0;r<_D;r++,tpos++) {\
1028 ttpos = tpos;\ 1028 ttpos = tpos;\
1029 rpos = m_vec+r;\ 1029 rpos = m_vec+r;\
1030 apos = (T*)a_m;\ 1030 apos = (T*)a_m;\
1031 for(c=0;c<_D;c++,ttpos+=_D,apos+=_D) {\ 1031 for(c=0;c<_D;c++,ttpos+=_D,apos+=_D) {\
1032 _value = zero();\ 1032 _value = zero();\
1033 mpos = rpos;\ 1033 mpos = rpos;\
1034 aapos = apos;\ 1034 aapos = apos;\
1035 for(i=0;i<_D;i++,mpos+=_D,aapos++) _v 1035 for(i=0;i<_D;i++,mpos+=_D,aapos++) _value += (*mpos) * (*aapos);\
1036 *ttpos = _value;\ 1036 *ttpos = _value;\
1037 }\ 1037 }\
1038 }\ 1038 }\
1039 _copy(a_tmp);\ 1039 _copy(a_tmp);\
1040 }\ 1040 }\
1041 \ 1041 \
1042 void _mul_mtx(const T a_m[]) {\ 1042 void _mul_mtx(const T a_m[]) {\
1043 T* _tmp = new T[dim2()];\ 1043 T* _tmp = new T[dim2()];\
1044 _mul_mtx(a_m,_tmp);\ 1044 _mul_mtx(a_m,_tmp);\
1045 delete [] _tmp;\ 1045 delete [] _tmp;\
1046 }\ 1046 }\
1047 \ 1047 \
1048 void _left_mul_mtx(const T a_m[]) {\ 1048 void _left_mul_mtx(const T a_m[]) {\
1049 /* this = a_m * this */\ 1049 /* this = a_m * this */\
1050 unsigned int _D = dimension();\ 1050 unsigned int _D = dimension();\
1051 T* _tmp = new T[dim2()];\ 1051 T* _tmp = new T[dim2()];\
1052 for(unsigned int r=0;r<_D;r++) {\ 1052 for(unsigned int r=0;r<_D;r++) {\
1053 for(unsigned int c=0;c<_D;c++) {\ 1053 for(unsigned int c=0;c<_D;c++) {\
1054 T _value = zero();\ 1054 T _value = zero();\
1055 for(unsigned int i=0;i<_D;i++) {\ 1055 for(unsigned int i=0;i<_D;i++) {\
1056 _value += (*(a_m+r+i*_D)) * (*(m_ve 1056 _value += (*(a_m+r+i*_D)) * (*(m_vec+i+c*_D)); /*optimize.*/\
1057 }\ 1057 }\
1058 *(_tmp+r+c*_D) = _value;\ 1058 *(_tmp+r+c*_D) = _value;\
1059 }\ 1059 }\
1060 }\ 1060 }\
1061 _copy(_tmp);\ 1061 _copy(_tmp);\
1062 delete [] _tmp;\ 1062 delete [] _tmp;\
1063 }\ 1063 }\
1064 \ 1064 \
1065 template <class PREC>\ 1065 template <class PREC>\
1066 T sub_determinant_prec(unsigned int a_ord,u 1066 T sub_determinant_prec(unsigned int a_ord,unsigned int aRs[],unsigned int aCs[],\
1067 const PREC& a_prec,P 1067 const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
1068 /*WARNING : to optimize, we do not check 1068 /*WARNING : to optimize, we do not check the content of aRs, aCs.*/\
1069 unsigned int ord = a_ord;\ 1069 unsigned int ord = a_ord;\
1070 if(ord==0) return zero();\ 1070 if(ord==0) return zero();\
1071 else if(ord==1) return value(aRs[0],aCs[0 1071 else if(ord==1) return value(aRs[0],aCs[0]);\
1072 else if(ord==2) {\ 1072 else if(ord==2) {\
1073 /*return (value(aRs[0],aCs[0]) * value( 1073 /*return (value(aRs[0],aCs[0]) * value(aRs[1],aCs[1]) -\
1074 value(aRs[1],aCs[0]) * value(a 1074 value(aRs[1],aCs[0]) * value(aRs[0],aCs[1])); \
1075 Optimize the upper :*/\ 1075 Optimize the upper :*/\
1076 \ 1076 \
1077 unsigned int _ord = dimension();\ 1077 unsigned int _ord = dimension();\
1078 \ 1078 \
1079 return ( (*(m_vec+aRs[0]+aCs[0]*_ord)) 1079 return ( (*(m_vec+aRs[0]+aCs[0]*_ord)) * (*(m_vec+aRs[1]+aCs[1]*_ord)) -\
1080 (*(m_vec+aRs[1]+aCs[0]*_ord)) 1080 (*(m_vec+aRs[1]+aCs[0]*_ord)) * (*(m_vec+aRs[0]+aCs[1]*_ord)) );\
1081 \ 1081 \
1082 } else if(ord==3) {\ 1082 } else if(ord==3) {\
1083 /* 00 01 02 \ 1083 /* 00 01 02 \
1084 10 11 12 \ 1084 10 11 12 \
1085 20 21 22 \ 1085 20 21 22 \
1086 */\ 1086 */\
1087 unsigned int _ord = dimension();\ 1087 unsigned int _ord = dimension();\
1088 \ 1088 \
1089 T v00 = *(m_vec+aRs[0]+aCs[0]*_ord);\ 1089 T v00 = *(m_vec+aRs[0]+aCs[0]*_ord);\
1090 T v10 = *(m_vec+aRs[1]+aCs[0]*_ord);\ 1090 T v10 = *(m_vec+aRs[1]+aCs[0]*_ord);\
1091 T v20 = *(m_vec+aRs[2]+aCs[0]*_ord);\ 1091 T v20 = *(m_vec+aRs[2]+aCs[0]*_ord);\
1092 T v01 = *(m_vec+aRs[0]+aCs[1]*_ord);\ 1092 T v01 = *(m_vec+aRs[0]+aCs[1]*_ord);\
1093 T v11 = *(m_vec+aRs[1]+aCs[1]*_ord);\ 1093 T v11 = *(m_vec+aRs[1]+aCs[1]*_ord);\
1094 T v21 = *(m_vec+aRs[2]+aCs[1]*_ord);\ 1094 T v21 = *(m_vec+aRs[2]+aCs[1]*_ord);\
1095 T v02 = *(m_vec+aRs[0]+aCs[2]*_ord);\ 1095 T v02 = *(m_vec+aRs[0]+aCs[2]*_ord);\
1096 T v12 = *(m_vec+aRs[1]+aCs[2]*_ord);\ 1096 T v12 = *(m_vec+aRs[1]+aCs[2]*_ord);\
1097 T v22 = *(m_vec+aRs[2]+aCs[2]*_ord);\ 1097 T v22 = *(m_vec+aRs[2]+aCs[2]*_ord);\
1098 T cof_00 = v11 * v22 - v21 * v12;\ 1098 T cof_00 = v11 * v22 - v21 * v12;\
1099 T cof_10 = v01 * v22 - v21 * v02;\ 1099 T cof_10 = v01 * v22 - v21 * v02;\
1100 T cof_20 = v01 * v12 - v11 * v02;\ 1100 T cof_20 = v01 * v12 - v11 * v02;\
1101 return (v00*cof_00-v10*cof_10+v20*cof_2 1101 return (v00*cof_00-v10*cof_10+v20*cof_20);\
1102 }\ 1102 }\
1103 \ 1103 \
1104 unsigned int rord = ord-1;\ 1104 unsigned int rord = ord-1;\
1105 unsigned int* cs = new unsigned int[rord] 1105 unsigned int* cs = new unsigned int[rord];\
1106 unsigned int* rs = new unsigned int[rord] 1106 unsigned int* rs = new unsigned int[rord];\
1107 \ 1107 \
1108 T v_rc;\ 1108 T v_rc;\
1109 \ 1109 \
1110 T det = zero();\ 1110 T det = zero();\
1111 {for(unsigned int i=0;i<rord;i++) {cs[i] = 1111 {for(unsigned int i=0;i<rord;i++) {cs[i] = aCs[i+1];}}\
1112 unsigned int c = 0;\ 1112 unsigned int c = 0;\
1113 /*if(c>=1) cs[c-1] = c-1;*/\ 1113 /*if(c>=1) cs[c-1] = c-1;*/\
1114 \ 1114 \
1115 {for(unsigned int i=0;i<rord;i++) {rs[i] = 1115 {for(unsigned int i=0;i<rord;i++) {rs[i] = aRs[i+1];}}\
1116 bool sg = true; /*c=0+r=0*/\ 1116 bool sg = true; /*c=0+r=0*/\
1117 for(unsigned int r=0;r<ord;r++) {\ 1117 for(unsigned int r=0;r<ord;r++) {\
1118 if(r>=1) rs[r-1] = aRs[r-1];\ 1118 if(r>=1) rs[r-1] = aRs[r-1];\
1119 v_rc = value(aRs[r],aCs[c]);\ 1119 v_rc = value(aRs[r],aCs[c]);\
1120 if(!is_zero(v_rc,a_prec,a_fabs)) {\ 1120 if(!is_zero(v_rc,a_prec,a_fabs)) {\
1121 T subdet = sub_determinant_prec(rord, 1121 T subdet = sub_determinant_prec(rord,rs,cs,a_prec,a_fabs);\
1122 if(sg)\ 1122 if(sg)\
1123 det += v_rc * subdet;\ 1123 det += v_rc * subdet;\
1124 else\ 1124 else\
1125 det -= v_rc * subdet;\ 1125 det -= v_rc * subdet;\
1126 }\ 1126 }\
1127 sg = sg?false:true;\ 1127 sg = sg?false:true;\
1128 }\ 1128 }\
1129 \ 1129 \
1130 delete [] cs;\ 1130 delete [] cs;\
1131 delete [] rs;\ 1131 delete [] rs;\
1132 \ 1132 \
1133 return det;\ 1133 return det;\
1134 }\ 1134 }\
1135 \ 1135 \
1136 static double zero_fabs(const T& a_number) 1136 static double zero_fabs(const T& a_number) {return a_number==zero()?0:1000000;}
1137 1137
1138 #endif 1138 #endif