Warning, /include/Geant4/tools/lina/MATCOM is written in an unsupported language. File is not indexed.
0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003
0004 #ifndef tools_MATCOM
0005 #define tools_MATCOM
0006
0007 #include "../eqT"
0008
0009 #include <cstddef> //size_t
0010
0011 // common code to class mat and nmat.
0012
0013 #define TOOLS_MATCOM \
0014 protected:\
0015 static T zero() {return T();}\
0016 static T one() {return T(1);}\
0017 static T minus_one() {return T(-1);}\
0018 static T two() {return T(2);}\
0019 public:\
0020 typedef T elem_t;\
0021 typedef unsigned int size_type;\
0022 public:\
0023 unsigned int rows() const {return dimension();}\
0024 unsigned int cols() const {return dimension();}\
0025 \
0026 void set_value(unsigned int aR,unsigned int aC,const T& a_value) { \
0027 m_vec[aR + aC * dimension()] = a_value;\
0028 }\
0029 \
0030 const T& value(unsigned int aR,unsigned int aC) const { \
0031 return m_vec[aR + aC * dimension()];\
0032 }\
0033 \
0034 T value(unsigned int aR,unsigned int aC) { \
0035 return m_vec[aR + aC * dimension()];\
0036 }\
0037 \
0038 void set_matrix(const TOOLS_MAT_CLASS& a_m){ /*optimization.*/\
0039 _copy(a_m.m_vec);\
0040 }\
0041 \
0042 void set_constant(const T& a_v){\
0043 for(unsigned int i=0;i<dim2();i++) m_vec[i] = a_v;\
0044 }\
0045 void set_zero(){\
0046 set_constant(zero());\
0047 }\
0048 void set_identity() {\
0049 set_zero();\
0050 for(unsigned int i=0;i<dimension();i++) m_vec[i+i*dimension()] = one();\
0051 }\
0052 void set_diagonal(const T& a_s) {\
0053 set_zero();\
0054 for(unsigned int i=0;i<dimension();i++) m_vec[i+i*dimension()] = a_s;\
0055 }\
0056 typedef T (*func)(const T&);\
0057 void apply_func(func a_func) {\
0058 T* pos = m_vec;\
0059 for(unsigned int i=0;i<dim2();i++,pos++) *pos = a_func(*pos);\
0060 }\
0061 template <class RANDOM>\
0062 void set_random(RANDOM& a_random) {\
0063 for(unsigned int i=0;i<dim2();i++) m_vec[i] = a_random.shoot();\
0064 }\
0065 template <class RANDOM>\
0066 void set_symmetric_random(RANDOM& a_random) {\
0067 unsigned int _D = dimension();\
0068 {for(unsigned int r=0;r<_D;r++) set_value(r,r,a_random.shoot());}\
0069 T rd;\
0070 {for(unsigned int r=0;r<_D;r++) {\
0071 for(unsigned int c=(r+1);c<_D;c++) {\
0072 rd = a_random.shoot();\
0073 set_value(r,c,rd);\
0074 set_value(c,r,rd);\
0075 }\
0076 }}\
0077 }\
0078 template <class RANDOM>\
0079 void set_antisymmetric_random(RANDOM& a_random) {\
0080 unsigned int _D = dimension();\
0081 {for(unsigned int r=0;r<_D;r++) set_value(r,r,zero());}\
0082 T rd;\
0083 {for(unsigned int r=0;r<_D;r++) {\
0084 for(unsigned int c=(r+1);c<_D;c++) {\
0085 rd = a_random.shoot();\
0086 set_value(r,c,rd);\
0087 set_value(c,r,minus_one()*rd);\
0088 }\
0089 }}\
0090 }\
0091 public:\
0092 template <class ARRAY>\
0093 bool mul_array(ARRAY& a_array,T a_tmp[]) const {\
0094 /* a_array = this *= a_array */\
0095 unsigned int _dim = dimension();\
0096 T* pos = a_tmp;\
0097 for(unsigned int r=0;r<_dim;r++,pos++) {\
0098 *pos = T();\
0099 for(unsigned int c=0;c<_dim;c++) *pos += m_vec[r+c*_dim]*a_array[c];\
0100 }\
0101 {for(unsigned int i=0;i<_dim;i++) a_array[i] = a_tmp[i];}\
0102 return true;\
0103 }\
0104 template <class VEC>\
0105 bool mul_vec(VEC& a_vec,T a_tmp[]) const {\
0106 /* a_vec = this *= a_vec */\
0107 unsigned int _dim = dimension();\
0108 if(a_vec.dimension()!=_dim) return false;\
0109 T* pos = a_tmp;\
0110 for(unsigned int r=0;r<_dim;r++,pos++) {\
0111 *pos = T();\
0112 for(unsigned int c=0;c<_dim;c++) *pos += m_vec[r+c*_dim]*a_vec[c];\
0113 }\
0114 {for(unsigned int i=0;i<_dim;i++) a_vec[i] = a_tmp[i];}\
0115 return true;\
0116 }\
0117 template <class VEC>\
0118 bool mul_vec(VEC& a_vec) const {\
0119 T* _tmp = new T[dimension()];\
0120 bool status = mul_vec(a_vec,_tmp);\
0121 delete [] _tmp;\
0122 return status;\
0123 }\
0124 \
0125 bool mul_array(T a_vec[],T a_tmp[]) const {\
0126 /* a_vec = this *= a_vec */\
0127 unsigned int _dim = dimension();\
0128 T* pos = a_tmp;\
0129 for(unsigned int r=0;r<_dim;r++,pos++) {\
0130 *pos = T();\
0131 for(unsigned int c=0;c<_dim;c++) *pos += m_vec[r+c*_dim]*a_vec[c];\
0132 }\
0133 {for(unsigned int i=0;i<_dim;i++) a_vec[i] = a_tmp[i];}\
0134 return true;\
0135 }\
0136 bool mul_array(T a_vec[]) const {\
0137 T* _tmp = new T[dimension()];\
0138 bool status = mul_array(a_vec,_tmp);\
0139 delete [] _tmp;\
0140 return status;\
0141 }\
0142 \
0143 void mul_mtx(const TOOLS_MAT_CLASS& a_m) {\
0144 _mul_mtx(a_m.m_vec);\
0145 }\
0146 void mul_mtx(const TOOLS_MAT_CLASS& a_m,T a_tmp[]) {\
0147 _mul_mtx(a_m.m_vec,a_tmp);\
0148 }\
0149 void left_mul_mtx(const TOOLS_MAT_CLASS& a_m) { \
0150 /* this = a_m * this :*/\
0151 _left_mul_mtx(a_m.m_vec);\
0152 }\
0153 bool equal(const TOOLS_MAT_CLASS& a_m) const {\
0154 if(&a_m==this) return true;\
0155 for(unsigned int i=0;i<dim2();i++) {\
0156 if(m_vec[i]!=a_m.m_vec[i]) return false;\
0157 }\
0158 return true;\
0159 }\
0160 \
0161 bool equal(const TOOLS_MAT_CLASS& a_m,const T& a_prec) const {\
0162 if(&a_m==this) return true;\
0163 T* tp = (T*)m_vec;\
0164 T* mp = (T*)a_m.m_vec;\
0165 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
0166 T diff = (*tp) - (*mp);\
0167 if(diff<zero()) diff *= minus_one();\
0168 if(diff>=a_prec) return false;\
0169 }\
0170 return true;\
0171 }\
0172 \
0173 template <class PREC>\
0174 bool equal_prec(const TOOLS_MAT_CLASS& a_m,const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0175 if(&a_m==this) return true;\
0176 T* tp = (T*)m_vec;\
0177 T* mp = (T*)a_m.m_vec;\
0178 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
0179 T diff = (*tp) - (*mp);\
0180 if(a_fabs(diff)>=a_prec) return false;\
0181 }\
0182 return true;\
0183 }\
0184 \
0185 void mx_diff(const TOOLS_MAT_CLASS& a_m,T& a_mx_diff) const {\
0186 T* tp = (T*)m_vec;\
0187 T* mp = (T*)a_m.m_vec;\
0188 a_mx_diff = (*tp) - (*mp);\
0189 if(a_mx_diff<zero()) a_mx_diff *= minus_one();\
0190 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
0191 T diff = (*tp) - (*mp);\
0192 if(diff<zero()) diff *= minus_one();\
0193 a_mx_diff = (diff>a_mx_diff?diff:a_mx_diff);\
0194 }\
0195 }\
0196 \
0197 bool to_rm_is_proportional(const TOOLS_MAT_CLASS& a_m,const T& a_prec,T& a_factor) const {\
0198 if(&a_m==this) {a_factor=one();return true;}\
0199 /* If true, then : a_m = a_factor * this.*/\
0200 a_factor = zero();\
0201 T* tp = (T*)m_vec;\
0202 T* mp = (T*)a_m.m_vec;\
0203 bool first = true;\
0204 for(unsigned int i=0;i<dim2();i++,tp++,mp++) {\
0205 if( ((*tp)==zero()) && ((*mp)==zero())) {\
0206 continue;\
0207 } else if( ((*tp)==zero()) && ((*mp)!=zero())) {\
0208 return false;\
0209 } else {\
0210 if(first) {\
0211 a_factor = (*mp)/(*tp);\
0212 first = false;\
0213 } else {\
0214 T diff = (*tp)*a_factor - (*mp);\
0215 if(diff<zero()) diff *= minus_one();\
0216 if(diff>=a_prec) return false;\
0217 }\
0218 }\
0219 }\
0220 return true;\
0221 }\
0222 \
0223 public:\
0224 bool is_proportional(const TOOLS_MAT_CLASS& a_right,T& a_factor) const {\
0225 /* If true, then : a_right = a_factor * this. a_factor could be zero.*/\
0226 if(this==&a_right) {a_factor=T(1);return true;}\
0227 a_factor = zero();\
0228 if(dimension()!=a_right.dimension()) return false;\
0229 T* lp = (T*)m_vec;\
0230 T* rp = (T*)a_right.m_vec;\
0231 bool first = true;\
0232 size_t _data_size = data_size();\
0233 for(size_t i=0;i<_data_size;i++,lp++,rp++) {\
0234 if(*lp==zero()) {\
0235 if(*rp==zero()) continue;\
0236 return false;\
0237 }\
0238 if(first) {\
0239 a_factor = (*rp)/(*lp);\
0240 first = false;\
0241 continue;\
0242 }\
0243 if((*lp)*a_factor!=(*rp)) return false;\
0244 }\
0245 return true;\
0246 }\
0247 \
0248 template <class PREC>\
0249 bool is_proportional_prec(const TOOLS_MAT_CLASS& a_right,const PREC& a_prec,PREC(*a_fabs)(const T&),T& a_factor) const {\
0250 /* If true, then : a_right = a_factor * this. a_factor could be zero.*/\
0251 if(this==&a_right) {a_factor=T(1);return true;}\
0252 a_factor = zero();\
0253 if(dimension()!=a_right.dimension()) return false;\
0254 T* lp = (T*)m_vec;\
0255 T* rp = (T*)a_right.m_vec;\
0256 bool first = true;\
0257 size_t _data_size = data_size();\
0258 for(size_t i=0;i<_data_size;i++,lp++,rp++) {\
0259 if(is_zero(*lp,a_prec,a_fabs)) {\
0260 if(is_zero(*rp,a_prec,a_fabs)) continue;\
0261 return false;\
0262 }\
0263 if(first) {\
0264 a_factor = (*rp)/(*lp);\
0265 first = false;\
0266 continue;\
0267 }\
0268 if(!numbers_are_equal((*lp)*a_factor,*rp,a_prec,a_fabs)) return false;\
0269 }\
0270 return true;\
0271 }\
0272 \
0273 bool is_diagonal() const {\
0274 unsigned int _D = dimension();\
0275 for(unsigned int r=0;r<_D;r++) {\
0276 for(unsigned int c=0;c<_D;c++) {\
0277 if(c!=r) {if(value(r,c)!=zero()) return false;}\
0278 }\
0279 }\
0280 return true;\
0281 }\
0282 \
0283 template <class PREC>\
0284 bool is_diagonal_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0285 unsigned int _D = dimension();\
0286 for(unsigned int r=0;r<_D;r++) {\
0287 for(unsigned int c=0;c<_D;c++) {\
0288 if(c!=r) {if(!is_zero(value(r,c),a_prec,a_fabs)) return false;}\
0289 }\
0290 }\
0291 return true;\
0292 }\
0293 \
0294 bool is_identity() const {\
0295 unsigned int _D = dimension();\
0296 {for(unsigned int r=0;r<_D;r++) {\
0297 if(value(r,r)!=one()) return false;\
0298 }}\
0299 {for(unsigned int r=0;r<_D;r++) {\
0300 for(unsigned int c=0;c<_D;c++) {\
0301 if(c!=r) {if(value(r,c)!=zero()) return false;}\
0302 }\
0303 }}\
0304 return true;\
0305 }\
0306 \
0307 template <class PREC>\
0308 bool is_identity_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0309 unsigned int _D = dimension();\
0310 {for(unsigned int r=0;r<_D;r++) {\
0311 if(!numbers_are_equal(value(r,r),one(),a_prec,a_fabs)) return false;\
0312 }}\
0313 {for(unsigned int r=0;r<_D;r++) {\
0314 for(unsigned int c=0;c<_D;c++) {\
0315 if(c!=r) {if(!is_zero(value(r,c),a_prec,a_fabs)) return false;}\
0316 }\
0317 }}\
0318 return true;\
0319 }\
0320 \
0321 const T* data() const {return m_vec;}\
0322 unsigned int size() const {return dim2();}\
0323 unsigned int data_size() const {return dim2();} /*for mathz*/\
0324 \
0325 T trace() const {\
0326 T _value = zero();\
0327 unsigned int _D = dimension();\
0328 for(unsigned int c=0;c<_D;c++) _value += m_vec[c+c*_D];\
0329 return _value;\
0330 }\
0331 \
0332 void transpose() {\
0333 unsigned int _D = dimension();\
0334 for(unsigned int r=0;r<_D;r++) {\
0335 for(unsigned int c=(r+1);c<_D;c++) {\
0336 T vrc = value(r,c);\
0337 T vcr = value(c,r);\
0338 set_value(r,c,vcr);\
0339 set_value(c,r,vrc);\
0340 }\
0341 }\
0342 }\
0343 \
0344 void multiply(const T& a_T) {\
0345 for(unsigned int i=0;i<dim2();i++) m_vec[i] *= a_T;\
0346 }\
0347 \
0348 bool is_symmetric() const {\
0349 unsigned int _D = dimension();\
0350 for(unsigned int r=0;r<_D;r++) {\
0351 for(unsigned int c=(r+1);c<_D;c++) {\
0352 if(value(r,c)!=value(c,r)) return false;\
0353 }\
0354 }\
0355 return true;\
0356 }\
0357 \
0358 template <class PREC>\
0359 bool is_symmetric_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0360 unsigned int _D = dimension();\
0361 for(unsigned int r=0;r<_D;r++) {\
0362 for(unsigned int c=(r+1);c<_D;c++) {\
0363 T diff = value(r,c)-value(c,r);\
0364 if(a_fabs(diff)>=a_prec) return false;\
0365 }\
0366 }\
0367 return true;\
0368 }\
0369 \
0370 bool is_antisymmetric() const {\
0371 unsigned int _D = dimension();\
0372 {for(unsigned int r=0;r<_D;r++) {\
0373 if(value(r,r)!=zero()) return false;\
0374 }}\
0375 for(unsigned int r=0;r<_D;r++) {\
0376 for(unsigned int c=(r+1);c<_D;c++) {\
0377 if(value(r,c)!=minus_one()*value(c,r)) return false;\
0378 }\
0379 }\
0380 return true;\
0381 }\
0382 \
0383 template <class PREC>\
0384 bool is_antisymmetric_prec(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0385 unsigned int _D = dimension();\
0386 {for(unsigned int r=0;r<_D;r++) {\
0387 if(a_fabs(value(r,r))>=a_prec) return false;\
0388 }}\
0389 for(unsigned int r=0;r<_D;r++) {\
0390 for(unsigned int c=(r+1);c<_D;c++) {\
0391 T diff = value(r,c)-minus_one()*value(c,r);\
0392 if(a_fabs(diff)>=a_prec) return false;\
0393 }\
0394 }\
0395 return true;\
0396 }\
0397 \
0398 void symmetric_part(TOOLS_MAT_CLASS& a_res) const {\
0399 a_res = *this;\
0400 a_res.transpose();\
0401 a_res += *this;\
0402 a_res.multiply(one()/two());\
0403 }\
0404 \
0405 void antisymmetric_part(TOOLS_MAT_CLASS& a_res) const {\
0406 a_res = *this;\
0407 a_res.transpose();\
0408 a_res.multiply(minus_one());\
0409 a_res += *this;\
0410 a_res.multiply(one()/two());\
0411 }\
0412 \
0413 template <class PREC>\
0414 bool is_block_UL_DR(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0415 /* Look if even dim matrix is of the form : |X 0|\
0416 |0 Y|*/\
0417 unsigned int _D = dimension();\
0418 unsigned int _D_2 = _D/2;\
0419 if(2*_D_2!=_D) return false;\
0420 for(unsigned int r=0;r<_D_2;r++) {\
0421 for(unsigned int c=_D_2;c<_D;c++) {\
0422 if(a_fabs(value(r,c))>=a_prec) return false;\
0423 }\
0424 }\
0425 for(unsigned int r=_D_2;r<_D;r++) {\
0426 for(unsigned int c=0;c<_D_2;c++) {\
0427 if(a_fabs(value(r,c))>=a_prec) return false;\
0428 }\
0429 }\
0430 return true;\
0431 }\
0432 \
0433 template <class PREC>\
0434 bool is_block_UR_DL(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0435 /* Look if even dim matrix is of the form : |0 X|\
0436 |Y 0|*/\
0437 unsigned int _D = dimension();\
0438 unsigned int _D_2 = _D/2;\
0439 if(2*_D_2!=_D) return false;\
0440 for(unsigned int r=0;r<_D_2;r++) {\
0441 for(unsigned int c=0;c<_D_2;c++) {\
0442 if(a_fabs(value(r,c))>=a_prec) return false;\
0443 }\
0444 }\
0445 for(unsigned int r=_D_2;r<_D;r++) {\
0446 for(unsigned int c=_D_2;c<_D;c++) {\
0447 if(a_fabs(value(r,c))>=a_prec) return false;\
0448 }\
0449 }\
0450 return true;\
0451 }\
0452 \
0453 template <class PREC>\
0454 bool is_decomplex(const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0455 /* Look if even dim matrix is of the form : | X Y|\
0456 |-Y X|*/\
0457 unsigned int _D = dimension();\
0458 unsigned int _D_2 = _D/2;\
0459 if(2*_D_2!=_D) return false;\
0460 for(unsigned int r=0;r<_D_2;r++) {\
0461 for(unsigned int c=0;c<_D_2;c++) {\
0462 if(a_fabs(value(r,c)-value(r+_D_2,c+_D_2))>=a_prec) return false;\
0463 }\
0464 }\
0465 for(unsigned int r=0;r<_D_2;r++) {\
0466 for(unsigned int c=_D_2;c<_D;c++) {\
0467 if(a_fabs(value(r,c)+value(r+_D_2,c-_D_2))>=a_prec) return false;\
0468 }\
0469 }\
0470 return true;\
0471 }\
0472 \
0473 template <class PREC>\
0474 T determinant_prec(unsigned int a_tmp_rs[],unsigned int a_tmp_cs[], /*[rord=dim-1]*/\
0475 const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0476 unsigned int ord = dimension();\
0477 if(ord==0) {\
0478 return zero();\
0479 } else if(ord==1) {\
0480 return *m_vec;\
0481 } else if(ord==2) {\
0482 T v00 = *m_vec;\
0483 T v10 = *(m_vec+1);\
0484 T v01 = *(m_vec+2);\
0485 T v11 = *(m_vec+3);\
0486 return (v00 * v11 - v10 * v01);\
0487 } else if(ord==3) {\
0488 /* 00 01 02 \
0489 10 11 12 \
0490 20 21 22 \
0491 */\
0492 T v00 = *m_vec;\
0493 T v10 = *(m_vec+1);\
0494 T v20 = *(m_vec+2);\
0495 T v01 = *(m_vec+3);\
0496 T v11 = *(m_vec+4);\
0497 T v21 = *(m_vec+5);\
0498 T v02 = *(m_vec+6);\
0499 T v12 = *(m_vec+7);\
0500 T v22 = *(m_vec+8);\
0501 T cof_00 = v11 * v22 - v21 * v12;\
0502 T cof_10 = v01 * v22 - v21 * v02;\
0503 T cof_20 = v01 * v12 - v11 * v02;\
0504 return (v00*cof_00-v10*cof_10+v20*cof_20);\
0505 }\
0506 \
0507 unsigned int rord = ord-1;\
0508 \
0509 T v_rc;\
0510 \
0511 T det = zero();\
0512 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs[i] = i+1;}}\
0513 unsigned int c = 0;\
0514 \
0515 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs[i] = i+1;}}\
0516 bool sg = true; /*c=0+r=0*/\
0517 for(unsigned int r=0;r<ord;r++) {\
0518 if(r>=1) a_tmp_rs[r-1] = r-1;\
0519 v_rc = value(r,c);\
0520 if(!is_zero(v_rc,a_prec,a_fabs)) {\
0521 T subdet = sub_determinant_prec(rord,a_tmp_rs,a_tmp_cs,a_prec,a_fabs);\
0522 if(sg) \
0523 det += v_rc * subdet;\
0524 else\
0525 det -= v_rc * subdet;\
0526 }\
0527 sg = sg?false:true;\
0528 }\
0529 \
0530 return det;\
0531 }\
0532 \
0533 T determinant(unsigned int a_tmp_rs[],unsigned int a_tmp_cs[]) const { /*[rord=dim-1]*/ \
0534 return determinant_prec<double>(a_tmp_rs,a_tmp_cs,0,zero_fabs);\
0535 }\
0536 \
0537 T determinant() const {\
0538 unsigned int ord = dimension();\
0539 if(ord==0) {\
0540 return zero();\
0541 } else if(ord==1) {\
0542 return *m_vec;\
0543 } else if(ord==2) {\
0544 T v00 = *m_vec;\
0545 T v10 = *(m_vec+1);\
0546 T v01 = *(m_vec+2);\
0547 T v11 = *(m_vec+3);\
0548 return (v00 * v11 - v10 * v01);\
0549 } else if(ord==3) {\
0550 /* 00 01 02 \
0551 10 11 12 \
0552 20 21 22 \
0553 */\
0554 T v00 = *m_vec;\
0555 T v10 = *(m_vec+1);\
0556 T v20 = *(m_vec+2);\
0557 T v01 = *(m_vec+3);\
0558 T v11 = *(m_vec+4);\
0559 T v21 = *(m_vec+5);\
0560 T v02 = *(m_vec+6);\
0561 T v12 = *(m_vec+7);\
0562 T v22 = *(m_vec+8);\
0563 T cof_00 = v11 * v22 - v21 * v12;\
0564 T cof_10 = v01 * v22 - v21 * v02;\
0565 T cof_20 = v01 * v12 - v11 * v02;\
0566 return (v00*cof_00-v10*cof_10+v20*cof_20);\
0567 }\
0568 unsigned int rord = ord-1;\
0569 unsigned int* rs = new unsigned int[rord];\
0570 unsigned int* cs = new unsigned int[rord];\
0571 T det = determinant(rs,cs);\
0572 delete [] rs;\
0573 delete [] cs;\
0574 return det;\
0575 }\
0576 \
0577 template <class PREC>\
0578 bool invert_prec(TOOLS_MAT_CLASS& a_res,\
0579 unsigned int a_tmp_rs[],unsigned int a_tmp_cs[], /*[rord=dim-1]*/\
0580 const PREC& a_prec,PREC(*a_fabs)(const T&) /*for det=?zero*/ \
0581 ) const { \
0582 unsigned int ord = dimension();\
0583 if(ord==0) return true;\
0584 \
0585 if(ord==1) {\
0586 T det = value(0,0);\
0587 if(is_zero(det,a_prec,a_fabs)) return false;\
0588 a_res.set_value(0,0,one()/det);\
0589 return true;\
0590 } else if(ord==2) {\
0591 T v00 = *m_vec;\
0592 T v10 = *(m_vec+1);\
0593 T v01 = *(m_vec+2);\
0594 T v11 = *(m_vec+3);\
0595 T det = (v00 * v11 - v10 * v01);\
0596 if(is_zero(det,a_prec,a_fabs)) return false;\
0597 a_res.set_value(0,0,v11/det);\
0598 a_res.set_value(1,1,v00/det);\
0599 a_res.set_value(0,1,minus_one()*v01/det);\
0600 a_res.set_value(1,0,minus_one()*v10/det);\
0601 return true;\
0602 } else if(ord==3) {\
0603 /* 00 01 02 \
0604 10 11 12 \
0605 20 21 22 \
0606 */\
0607 T v00 = *m_vec;\
0608 T v10 = *(m_vec+1);\
0609 T v20 = *(m_vec+2);\
0610 T v01 = *(m_vec+3);\
0611 T v11 = *(m_vec+4);\
0612 T v21 = *(m_vec+5);\
0613 T v02 = *(m_vec+6);\
0614 T v12 = *(m_vec+7);\
0615 T v22 = *(m_vec+8);\
0616 T cof_00 = v11 * v22 - v21 * v12;\
0617 T cof_10 = v01 * v22 - v21 * v02;\
0618 T cof_20 = v01 * v12 - v11 * v02;\
0619 T det = (v00*cof_00-v10*cof_10+v20*cof_20);\
0620 if(is_zero(det,a_prec,a_fabs)) return false;\
0621 T cof_01 = v10 * v22 - v20 * v12;\
0622 T cof_11 = v00 * v22 - v20 * v02;\
0623 T cof_21 = v00 * v12 - v10 * v02;\
0624 T cof_02 = v10 * v21 - v20 * v11;\
0625 T cof_12 = v00 * v21 - v20 * v01;\
0626 T cof_22 = v00 * v11 - v10 * v01;\
0627 a_res.set_value(0,0,cof_00/det);\
0628 a_res.set_value(1,0,minus_one()*cof_01/det);\
0629 a_res.set_value(2,0,cof_02/det);\
0630 a_res.set_value(0,1,minus_one()*cof_10/det);\
0631 a_res.set_value(1,1,cof_11/det);\
0632 a_res.set_value(2,1,minus_one()*cof_12/det);\
0633 a_res.set_value(0,2,cof_20/det);\
0634 a_res.set_value(1,2,minus_one()*cof_21/det);\
0635 a_res.set_value(2,2,cof_22/det);\
0636 return true;\
0637 }\
0638 \
0639 /*Generic invertion method.*/\
0640 \
0641 unsigned int rord = ord-1;\
0642 \
0643 /* Get det with r = 0;*/\
0644 T det = zero();\
0645 T subdet;\
0646 {\
0647 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs[i] = i+1;}}\
0648 unsigned int r = 0;\
0649 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs[i] = i+1;}}\
0650 bool sg = true; /*r=0+c=0*/\
0651 for(unsigned int c=0;c<ord;c++) {\
0652 if(c>=1) a_tmp_cs[c-1] = c-1;\
0653 subdet = sub_determinant_prec(rord,a_tmp_rs,a_tmp_cs,a_prec,a_fabs);\
0654 if(sg) {\
0655 det += value(r,c) * subdet;\
0656 a_res.set_value(c,r,subdet);\
0657 sg = false;\
0658 } else {\
0659 det += value(r,c) * subdet * minus_one();\
0660 a_res.set_value(c,r,subdet * minus_one());\
0661 sg = true;\
0662 }\
0663 }}\
0664 \
0665 if(is_zero(det,a_prec,a_fabs)) return false;\
0666 \
0667 {for(unsigned int c=0;c<ord;c++) {\
0668 a_res.set_value(c,0,a_res.value(c,0)/det);\
0669 }}\
0670 \
0671 {for(unsigned int i=0;i<rord;i++) {a_tmp_rs[i] = i+1;}}\
0672 bool sgr = false; /*r=1+c=0*/\
0673 for(unsigned int r=1;r<ord;r++) {\
0674 if(r>=1) a_tmp_rs[r-1] = r-1;\
0675 {for(unsigned int i=0;i<rord;i++) {a_tmp_cs[i] = i+1;}}\
0676 bool sg = sgr;\
0677 for(unsigned int c=0;c<ord;c++) {\
0678 if(c>=1) a_tmp_cs[c-1] = c-1;\
0679 subdet = sub_determinant_prec(rord,a_tmp_rs,a_tmp_cs,a_prec,a_fabs);\
0680 if(sg) {\
0681 a_res.set_value(c,r,subdet/det);\
0682 sg = false;\
0683 } else {\
0684 a_res.set_value(c,r,(subdet * minus_one())/det);\
0685 sg = true;\
0686 }\
0687 }\
0688 sgr = sgr?false:true;\
0689 }\
0690 \
0691 return true;\
0692 }\
0693 \
0694 template <class PREC>\
0695 bool invert_prec(TOOLS_MAT_CLASS& a_res,const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0696 unsigned int ord = dimension();\
0697 if(ord==0) return true;\
0698 \
0699 if(ord==1) {\
0700 T det = value(0,0);\
0701 if(is_zero(det,a_prec,a_fabs)) return false;\
0702 a_res.set_value(0,0,one()/det);\
0703 return true;\
0704 } else if(ord==2) {\
0705 T v00 = *m_vec;\
0706 T v10 = *(m_vec+1);\
0707 T v01 = *(m_vec+2);\
0708 T v11 = *(m_vec+3);\
0709 T det = (v00 * v11 - v10 * v01);\
0710 if(is_zero(det,a_prec,a_fabs)) return false;\
0711 a_res.set_value(0,0,v11/det);\
0712 a_res.set_value(1,1,v00/det);\
0713 a_res.set_value(0,1,minus_one()*v01/det);\
0714 a_res.set_value(1,0,minus_one()*v10/det);\
0715 return true;\
0716 } else if(ord==3) {\
0717 /* 00 01 02 \
0718 10 11 12 \
0719 20 21 22 \
0720 */\
0721 T v00 = *m_vec;\
0722 T v10 = *(m_vec+1);\
0723 T v20 = *(m_vec+2);\
0724 T v01 = *(m_vec+3);\
0725 T v11 = *(m_vec+4);\
0726 T v21 = *(m_vec+5);\
0727 T v02 = *(m_vec+6);\
0728 T v12 = *(m_vec+7);\
0729 T v22 = *(m_vec+8);\
0730 T cof_00 = v11 * v22 - v21 * v12;\
0731 T cof_10 = v01 * v22 - v21 * v02;\
0732 T cof_20 = v01 * v12 - v11 * v02;\
0733 T det = (v00*cof_00-v10*cof_10+v20*cof_20);\
0734 if(is_zero(det,a_prec,a_fabs)) return false;\
0735 T cof_01 = v10 * v22 - v20 * v12;\
0736 T cof_11 = v00 * v22 - v20 * v02;\
0737 T cof_21 = v00 * v12 - v10 * v02;\
0738 T cof_02 = v10 * v21 - v20 * v11;\
0739 T cof_12 = v00 * v21 - v20 * v01;\
0740 T cof_22 = v00 * v11 - v10 * v01;\
0741 a_res.set_value(0,0,cof_00/det);\
0742 a_res.set_value(1,0,minus_one()*cof_01/det);\
0743 a_res.set_value(2,0,cof_02/det);\
0744 a_res.set_value(0,1,minus_one()*cof_10/det);\
0745 a_res.set_value(1,1,cof_11/det);\
0746 a_res.set_value(2,1,minus_one()*cof_12/det);\
0747 a_res.set_value(0,2,cof_20/det);\
0748 a_res.set_value(1,2,minus_one()*cof_21/det);\
0749 a_res.set_value(2,2,cof_22/det);\
0750 return true;\
0751 }\
0752 \
0753 unsigned int rord = ord-1;\
0754 unsigned int* cs = new unsigned int[rord];\
0755 unsigned int* rs = new unsigned int[rord];\
0756 bool status = invert_prec(a_res,rs,cs,a_prec,a_fabs);\
0757 delete [] cs;\
0758 delete [] rs;\
0759 return status;\
0760 }\
0761 \
0762 bool invert(TOOLS_MAT_CLASS& a_res,unsigned int a_tmp_rs[],unsigned int a_tmp_cs[]) const { /*[rord=dim-1]*/ \
0763 return invert_prec<double>(a_res,a_tmp_rs,a_tmp_cs,0,zero_fabs);\
0764 }\
0765 \
0766 bool invert(TOOLS_MAT_CLASS& a_res) const {\
0767 return invert_prec<double>(a_res,0,zero_fabs);\
0768 }\
0769 \
0770 void power(unsigned int a_n,TOOLS_MAT_CLASS& a_res) const {\
0771 a_res.set_identity();\
0772 T* _tmp = new T[dim2()];\
0773 for(unsigned int i=0;i<a_n;i++) {\
0774 a_res._mul_mtx(m_vec,_tmp);\
0775 }\
0776 delete [] _tmp;\
0777 }\
0778 \
0779 void exp(unsigned int a_le,TOOLS_MAT_CLASS& a_res,TOOLS_MAT_CLASS& a_tmp,T a_tmp_2[]) const { /*OPTIMIZATION*/\
0780 /* result = I + M + M**2/2! + M**3/3! + .... */\
0781 a_res.set_identity();\
0782 a_tmp.set_identity();\
0783 for(unsigned int i=1;i<=a_le;i++) {\
0784 a_tmp._mul_mtx(m_vec,a_tmp_2);\
0785 a_tmp.multiply(one()/T(i));\
0786 a_res += a_tmp;\
0787 }\
0788 }\
0789 \
0790 void exp(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
0791 /* result = I + M + M**2/2! + M**3/3! + .... */\
0792 TOOLS_MAT_CLASS tmp(*this);\
0793 tmp.set_identity();\
0794 T* tmp_2 = new T[dim2()];\
0795 exp(a_le,a_res,tmp,tmp_2);\
0796 delete [] tmp_2;\
0797 }\
0798 \
0799 void cosh(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
0800 /* result = I + M**2/2! + M**4/4! + .... */\
0801 a_res.set_identity();\
0802 TOOLS_MAT_CLASS M_2(*this);\
0803 M_2._mul_mtx(m_vec);\
0804 TOOLS_MAT_CLASS tmp(*this);\
0805 tmp.set_identity();\
0806 T* _tmp = new T[dim2()];\
0807 for(unsigned int i=1;i<=a_le;i+=2) {\
0808 tmp._mul_mtx(M_2.m_vec,_tmp);\
0809 tmp.multiply(one()/T(i+0)); \
0810 tmp.multiply(one()/T(i+1)); \
0811 a_res += tmp;\
0812 }\
0813 delete [] _tmp;\
0814 }\
0815 \
0816 void sinh(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
0817 /* result = M + M**3/3! + .... */\
0818 a_res._copy(m_vec);\
0819 TOOLS_MAT_CLASS M_2(*this);\
0820 M_2._mul_mtx(m_vec);\
0821 TOOLS_MAT_CLASS tmp(*this);\
0822 tmp._copy(m_vec);\
0823 T* _tmp = new T[dim2()];\
0824 for(unsigned int i=1;i<=a_le;i+=2) {\
0825 tmp._mul_mtx(M_2.m_vec,_tmp);\
0826 tmp.multiply(one()/T(i+1)); \
0827 tmp.multiply(one()/T(i+2)); \
0828 a_res += tmp;\
0829 }\
0830 delete [] _tmp;\
0831 }\
0832 \
0833 void log(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
0834 /* result = (M-I) - (M-I)**2/2 + (M-I)**3/3 +... */\
0835 /* WARNING : touchy, it may not converge ! */\
0836 a_res.set_zero();\
0837 \
0838 TOOLS_MAT_CLASS M_I;\
0839 M_I.set_identity();\
0840 M_I.multiply(minus_one());\
0841 M_I._add_mtx(m_vec);\
0842 \
0843 TOOLS_MAT_CLASS M_Ip(M_I);\
0844 T fact = minus_one();\
0845 \
0846 TOOLS_MAT_CLASS tmp;\
0847 \
0848 T* _tmp = new T[dim2()];\
0849 for(unsigned int i=0;i<=a_le;i++) {\
0850 fact *= minus_one(); \
0851 tmp = M_Ip;\
0852 tmp.multiply(fact/T(i+1)); \
0853 a_res += tmp;\
0854 M_Ip._mul_mtx(M_I.m_vec,_tmp);\
0855 }\
0856 delete [] _tmp;\
0857 }\
0858 \
0859 void omega(unsigned int a_le,TOOLS_MAT_CLASS& a_res,TOOLS_MAT_CLASS& a_tmp,T a_tmp_2[]) const { /*OPTIMIZATION*/\
0860 /* result = I + M/2! + M**2/3! + M**3/4! + .... */\
0861 a_res.set_identity();\
0862 a_tmp.set_identity();\
0863 for(unsigned int i=1;i<=a_le;i++) {\
0864 a_tmp._mul_mtx(m_vec,a_tmp_2);\
0865 a_tmp.multiply(one()/T(i+1));\
0866 a_res += a_tmp;\
0867 }\
0868 }\
0869 \
0870 void omega(unsigned int a_le,TOOLS_MAT_CLASS& a_res) const {\
0871 /* result = I + M/2! + M**2/3! + M**3/4! + .... */\
0872 TOOLS_MAT_CLASS tmp(*this);\
0873 tmp.set_identity();\
0874 T* tmp_2 = new T[dim2()];\
0875 omega(a_le,a_res,tmp,tmp_2);\
0876 delete [] tmp_2;\
0877 }\
0878 \
0879 template <class MAT>\
0880 bool copy(const MAT& a_from) {\
0881 /*for exa from a double matrix to a symbol matrix*/\
0882 unsigned int _D = dimension();\
0883 if(a_from.dimension()!=_D) return false;\
0884 for(unsigned int r=0;r<_D;r++) {\
0885 for(unsigned int c=0;c<_D;c++) {\
0886 set_value(r,c,a_from.value(r,c));\
0887 }\
0888 }\
0889 return true;\
0890 }\
0891 public: /*operators*/\
0892 T operator()(unsigned int a_r,unsigned int a_c) const {\
0893 /*WARNING : no check on a_r,a_c.*/\
0894 return m_vec[a_r + a_c * dimension()];\
0895 }\
0896 \
0897 T& operator[](size_t a_index) { /*for tools/sg/sf_vec*/\
0898 /*WARNING : no check on a_index.*/\
0899 return m_vec[a_index];\
0900 }\
0901 const T& operator[](size_t a_index) const {\
0902 /*WARNING : no check on a_index.*/\
0903 return m_vec[a_index];\
0904 }\
0905 bool operator==(const TOOLS_MAT_CLASS& a_array) const {\
0906 return equal(a_array);\
0907 }\
0908 bool operator!=(const TOOLS_MAT_CLASS& a_array) const {\
0909 return !operator==(a_array);\
0910 }\
0911 TOOLS_MAT_CLASS& operator*=(const TOOLS_MAT_CLASS& a_m) {\
0912 _mul_mtx(a_m.m_vec);\
0913 return *this;\
0914 }\
0915 TOOLS_MAT_CLASS& operator+=(const TOOLS_MAT_CLASS& a_m) {\
0916 _add_mtx(a_m.m_vec);\
0917 return *this;\
0918 }\
0919 TOOLS_MAT_CLASS& operator-=(const TOOLS_MAT_CLASS& a_m) {\
0920 _sub_mtx(a_m.m_vec);\
0921 return *this;\
0922 }\
0923 TOOLS_MAT_CLASS& operator*=(const T& a_fac) {\
0924 for(unsigned int i=0;i<dim2();i++) m_vec[i] *= a_fac;\
0925 return *this;\
0926 }\
0927 \
0928 TOOLS_MAT_CLASS operator*(const T& a_fac) {\
0929 TOOLS_MAT_CLASS res;\
0930 res.operator*=(a_fac);\
0931 return res;\
0932 }\
0933 protected:\
0934 void _copy(const T a_m[]) {\
0935 T* tp = (T*)m_vec;T* ap = (T*)a_m;\
0936 for(unsigned int i=0;i<dim2();i++,tp++,ap++) *tp = *ap;\
0937 }\
0938 \
0939 void _add_mtx(const T a_m[]) { /* this = this + a_m, */\
0940 T* tp = (T*)m_vec;T* ap = (T*)a_m;\
0941 for(unsigned int i=0;i<dim2();i++,tp++,ap++) *tp += *ap;\
0942 }\
0943 void _sub_mtx(const T a_m[]) { /* this = this - a_m, */\
0944 T* tp = (T*)m_vec;T* ap = (T*)a_m;\
0945 for(unsigned int i=0;i<dim2();i++,tp++,ap++) *tp -= *ap;\
0946 }\
0947 \
0948 void _mul_mtx(const T a_m[],T a_tmp[]) { /*OPTIMIZATION*/\
0949 /* this = this * a_m */\
0950 typedef T* Tp;\
0951 Tp tpos,ttpos,rpos,apos,mpos,aapos;\
0952 T _value;\
0953 unsigned int r,c,i;\
0954 \
0955 unsigned int _D = dimension();\
0956 \
0957 tpos = a_tmp;\
0958 for(r=0;r<_D;r++,tpos++) {\
0959 ttpos = tpos;\
0960 rpos = m_vec+r;\
0961 apos = (T*)a_m;\
0962 for(c=0;c<_D;c++,ttpos+=_D,apos+=_D) {\
0963 _value = zero();\
0964 mpos = rpos;\
0965 aapos = apos;\
0966 for(i=0;i<_D;i++,mpos+=_D,aapos++) _value += (*mpos) * (*aapos);\
0967 *ttpos = _value;\
0968 }\
0969 }\
0970 _copy(a_tmp);\
0971 }\
0972 \
0973 void _mul_mtx(const T a_m[]) {\
0974 T* _tmp = new T[dim2()];\
0975 _mul_mtx(a_m,_tmp);\
0976 delete [] _tmp;\
0977 }\
0978 \
0979 void _left_mul_mtx(const T a_m[]) {\
0980 /* this = a_m * this */\
0981 unsigned int _D = dimension();\
0982 T* _tmp = new T[dim2()];\
0983 for(unsigned int r=0;r<_D;r++) {\
0984 for(unsigned int c=0;c<_D;c++) {\
0985 T _value = zero();\
0986 for(unsigned int i=0;i<_D;i++) {\
0987 _value += (*(a_m+r+i*_D)) * (*(m_vec+i+c*_D)); /*optimize.*/\
0988 }\
0989 *(_tmp+r+c*_D) = _value;\
0990 }\
0991 }\
0992 _copy(_tmp);\
0993 delete [] _tmp;\
0994 }\
0995 \
0996 template <class PREC>\
0997 T sub_determinant_prec(unsigned int a_ord,unsigned int aRs[],unsigned int aCs[],\
0998 const PREC& a_prec,PREC(*a_fabs)(const T&)) const {\
0999 /*WARNING : to optimize, we do not check the content of aRs, aCs.*/\
1000 unsigned int ord = a_ord;\
1001 if(ord==0) return zero();\
1002 else if(ord==1) return value(aRs[0],aCs[0]);\
1003 else if(ord==2) {\
1004 /*return (value(aRs[0],aCs[0]) * value(aRs[1],aCs[1]) -\
1005 value(aRs[1],aCs[0]) * value(aRs[0],aCs[1])); \
1006 Optimize the upper :*/\
1007 \
1008 unsigned int _ord = dimension();\
1009 \
1010 return ( (*(m_vec+aRs[0]+aCs[0]*_ord)) * (*(m_vec+aRs[1]+aCs[1]*_ord)) -\
1011 (*(m_vec+aRs[1]+aCs[0]*_ord)) * (*(m_vec+aRs[0]+aCs[1]*_ord)) );\
1012 \
1013 } else if(ord==3) {\
1014 /* 00 01 02 \
1015 10 11 12 \
1016 20 21 22 \
1017 */\
1018 unsigned int _ord = dimension();\
1019 \
1020 T v00 = *(m_vec+aRs[0]+aCs[0]*_ord);\
1021 T v10 = *(m_vec+aRs[1]+aCs[0]*_ord);\
1022 T v20 = *(m_vec+aRs[2]+aCs[0]*_ord);\
1023 T v01 = *(m_vec+aRs[0]+aCs[1]*_ord);\
1024 T v11 = *(m_vec+aRs[1]+aCs[1]*_ord);\
1025 T v21 = *(m_vec+aRs[2]+aCs[1]*_ord);\
1026 T v02 = *(m_vec+aRs[0]+aCs[2]*_ord);\
1027 T v12 = *(m_vec+aRs[1]+aCs[2]*_ord);\
1028 T v22 = *(m_vec+aRs[2]+aCs[2]*_ord);\
1029 T cof_00 = v11 * v22 - v21 * v12;\
1030 T cof_10 = v01 * v22 - v21 * v02;\
1031 T cof_20 = v01 * v12 - v11 * v02;\
1032 return (v00*cof_00-v10*cof_10+v20*cof_20);\
1033 }\
1034 \
1035 unsigned int rord = ord-1;\
1036 unsigned int* cs = new unsigned int[rord];\
1037 unsigned int* rs = new unsigned int[rord];\
1038 \
1039 T v_rc;\
1040 \
1041 T det = zero();\
1042 {for(unsigned int i=0;i<rord;i++) {cs[i] = aCs[i+1];}}\
1043 unsigned int c = 0;\
1044 /*if(c>=1) cs[c-1] = c-1;*/\
1045 \
1046 {for(unsigned int i=0;i<rord;i++) {rs[i] = aRs[i+1];}}\
1047 bool sg = true; /*c=0+r=0*/\
1048 for(unsigned int r=0;r<ord;r++) {\
1049 if(r>=1) rs[r-1] = aRs[r-1];\
1050 v_rc = value(aRs[r],aCs[c]);\
1051 if(!is_zero(v_rc,a_prec,a_fabs)) {\
1052 T subdet = sub_determinant_prec(rord,rs,cs,a_prec,a_fabs);\
1053 if(sg)\
1054 det += v_rc * subdet;\
1055 else\
1056 det -= v_rc * subdet;\
1057 }\
1058 sg = sg?false:true;\
1059 }\
1060 \
1061 delete [] cs;\
1062 delete [] rs;\
1063 \
1064 return det;\
1065 }\
1066 \
1067 static double zero_fabs(const T& a_number) {return a_number==zero()?0:1000000;}
1068
1069 #endif
