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