Warning, /include/Geant4/tools/lina/mat 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_mat
0005 #define tools_mat
0006
0007 #ifdef TOOLS_MEM
0008 #include "../mem"
0009 #endif
0010
0011 #include "MATCOM"
0012
0013 //#include <cstring> //memcpy
0014
0015 //#define TOOLS_MAT_NEW
0016 namespace tools {
0017
0018 template <class T,unsigned int D>
0019 class mat {
0020 static const unsigned int _D2 = D*D;
0021 unsigned int dim2() const {return _D2;}
0022 #define TOOLS_MAT_CLASS mat
0023 TOOLS_MATCOM
0024 #undef TOOLS_MAT_CLASS
0025 private:
0026 #ifdef TOOLS_MEM
0027 static const std::string& s_class() {
0028 static const std::string s_v("tools::mat");
0029 return s_v;
0030 }
0031 #endif
0032 public:
0033 mat(
0034 #ifdef TOOLS_MEM
0035 bool a_inc = true
0036 #endif
0037 ) {
0038 #ifdef TOOLS_MEM
0039 if(a_inc) mem::increment(s_class().c_str());
0040 #endif
0041 #ifdef TOOLS_MAT_NEW
0042 m_vec = new T[D*D];
0043 #endif
0044 for(unsigned int i=0;i<_D2;i++) m_vec[i] = zero();
0045 }
0046 virtual ~mat() {
0047 #ifdef TOOLS_MAT_NEW
0048 delete [] m_vec;
0049 #endif
0050 #ifdef TOOLS_MEM
0051 mem::decrement(s_class().c_str());
0052 #endif
0053 }
0054 public:
0055 mat(const mat& a_from) {
0056 #ifdef TOOLS_MEM
0057 mem::increment(s_class().c_str());
0058 #endif
0059 #ifdef TOOLS_MAT_NEW
0060 m_vec = new T[D*D];
0061 #endif
0062 _copy(a_from.m_vec);
0063 }
0064 mat& operator=(const mat& a_from){
0065 if(&a_from==this) return *this;
0066 _copy(a_from.m_vec);
0067 return *this;
0068 }
0069 public:
0070 mat(const T a_v[]){
0071 #ifdef TOOLS_MEM
0072 mem::increment(s_class().c_str());
0073 #endif
0074 #ifdef TOOLS_MAT_NEW
0075 m_vec = new T[D*D];
0076 #endif
0077 _copy(a_v);
0078 }
0079 public:
0080 unsigned int dimension() const {return D;}
0081 protected:
0082 #ifdef TOOLS_MAT_NEW
0083 T* m_vec;
0084 #else
0085 T m_vec[D*D];
0086 #endif
0087
0088 private:static void check_instantiation() {mat<float,2> dummy;}
0089 };
0090
0091 template <class T>
0092 class nmat {
0093 unsigned int dim2() const {return m_D2;}
0094 #define TOOLS_MAT_CLASS nmat
0095 TOOLS_MATCOM
0096 #undef TOOLS_MAT_CLASS
0097 private:
0098 #ifdef TOOLS_MEM
0099 static const std::string& s_class() {
0100 static const std::string s_v("tools::nmat");
0101 return s_v;
0102 }
0103 #endif
0104 public:
0105 nmat(unsigned int a_D):m_D(a_D),m_D2(a_D*a_D),m_vec(0) {
0106 #ifdef TOOLS_MEM
0107 mem::increment(s_class().c_str());
0108 #endif
0109 m_vec = new T[m_D2];
0110 for(unsigned int i=0;i<m_D2;i++) m_vec[i] = zero();
0111 }
0112 virtual ~nmat() {
0113 delete [] m_vec;
0114 #ifdef TOOLS_MEM
0115 mem::decrement(s_class().c_str());
0116 #endif
0117 }
0118 public:
0119 nmat(const nmat& a_from)
0120 :m_D(a_from.m_D),m_D2(a_from.m_D2),m_vec(0)
0121 {
0122 #ifdef TOOLS_MEM
0123 mem::increment(s_class().c_str());
0124 #endif
0125 m_vec = new T[m_D2];
0126 _copy(a_from.m_vec);
0127 }
0128 nmat& operator=(const nmat& a_from){
0129 if(&a_from==this) return *this;
0130 if(a_from.m_D!=m_D) {
0131 m_D = a_from.m_D;
0132 m_D2 = a_from.m_D2;
0133 delete [] m_vec;
0134 m_vec = new T[m_D2];
0135 }
0136 _copy(a_from.m_vec);
0137 return *this;
0138 }
0139 public:
0140 nmat(unsigned int a_D,const T a_v[])
0141 :m_D(a_D),m_D2(a_D*a_D),m_vec(0)
0142 {
0143 #ifdef TOOLS_MEM
0144 mem::increment(s_class().c_str());
0145 #endif
0146 m_vec = new T[m_D2];
0147 _copy(a_v);
0148 }
0149 public:
0150 unsigned int dimension() const {return m_D;}
0151 protected:
0152 unsigned int m_D;
0153 unsigned int m_D2;
0154 T* m_vec;
0155 private:static void check_instantiation() {nmat<float> dummy(2);}
0156 };
0157
0158 template <class T,unsigned int D>
0159 inline nmat<T> copy(const mat<T,D>& a_from) {
0160 unsigned int D2 = D*D;
0161 nmat<T> v(D);
0162 for(unsigned int i=0;i<D2;i++) v[i] = a_from[i];
0163 return v;
0164 }
0165
0166 template <class VECTOR>
0167 inline void multiply(VECTOR& a_vec,const typename VECTOR::value_type& a_mat) {
0168 //typedef typename VECTOR::size_type sz_t;
0169 //sz_t number = a_vec.size();
0170 //for(sz_t index=0;index<number;index++) a_vec[index] *= a_mat;
0171 typedef typename VECTOR::iterator it_t;
0172 for(it_t it=a_vec.begin();it!=a_vec.end();++it) *it *= a_mat;
0173 }
0174
0175 template <class VECTOR>
0176 inline void multiply(VECTOR& a_vec,const typename VECTOR::value_type::elem_t& a_value) {
0177 typedef typename VECTOR::iterator it_t;
0178 for(it_t it=a_vec.begin();it!=a_vec.end();++it) (*it).multiply(a_value);
0179 }
0180
0181 ///////////////////////////////////////////////////////////////////////////////////////
0182 /// related to complex numbers : //////////////////////////////////////////////////////
0183 ///////////////////////////////////////////////////////////////////////////////////////
0184
0185 template <class MAT>
0186 inline void conjugate(MAT& a_m,typename MAT::elem_t (*a_conj)(const typename MAT::elem_t&)) {
0187 typedef typename MAT::elem_t T;
0188 T* pos = const_cast<T*>(a_m.data());
0189 unsigned int D2 = a_m.dimension()*a_m.dimension();
0190 for(unsigned int i=0;i<D2;i++,pos++) {
0191 *pos = a_conj(*pos); //T = std::complex<>
0192 }
0193 }
0194
0195 template <class MAT>
0196 inline bool is_real(MAT& a_m,typename MAT::elem_t::value_type (*a_imag)(const typename MAT::elem_t&)) {
0197 typedef typename MAT::elem_t T;
0198 T* pos = const_cast<T*>(a_m.data());
0199 unsigned int D2 = a_m.dimension()*a_m.dimension();
0200 for(unsigned int i=0;i<D2;i++,pos++) {if(a_imag(*pos)) return false;}
0201 return true;
0202 }
0203
0204 template <class MAT,class PREC>
0205 inline bool is_real_prec(MAT& a_m,typename MAT::elem_t::value_type (*a_imag)(const typename MAT::elem_t&),
0206 const PREC& a_prec,PREC(*a_fabs)(const typename MAT::elem_t::value_type&)) {\
0207 typedef typename MAT::elem_t T;
0208 T* pos = const_cast<T*>(a_m.data());
0209 unsigned int D2 = a_m.dimension()*a_m.dimension();
0210 for(unsigned int i=0;i<D2;i++,pos++) {if(a_fabs(a_imag(*pos))>=a_prec) return false;}
0211 return true;
0212 }
0213
0214 template <class MAT>
0215 inline bool is_imag(MAT& a_m,typename MAT::elem_t::value_type (*a_real)(const typename MAT::elem_t&)) {
0216 typedef typename MAT::elem_t T;
0217 T* pos = const_cast<T*>(a_m.data());
0218 unsigned int D2 = a_m.dimension()*a_m.dimension();
0219 for(unsigned int i=0;i<D2;i++,pos++) {if(a_real(*pos)) return false;}
0220 return true;
0221 }
0222
0223 template <class CMAT,class RMAT>
0224 inline bool to_real(const CMAT& a_c,RMAT& a_r,typename CMAT::elem_t::value_type (*a_real)(const typename CMAT::elem_t&)) {
0225 if(a_r.dimension()!=a_c.dimension()) return false;
0226 typedef typename CMAT::elem_t CT;
0227 const CT* cpos = a_c.data();
0228 typedef typename RMAT::elem_t RT;
0229 RT* rpos = const_cast<RT*>(a_r.data());
0230 unsigned int D2 = a_c.dimension()*a_c.dimension();
0231 for(unsigned int i=0;i<D2;i++,cpos++,rpos++) *rpos = a_real(*cpos);
0232 return true;
0233 }
0234
0235 template <class RMAT,class CMAT>
0236 inline bool to_complex(const RMAT& a_r,CMAT& a_c) {
0237 if(a_c.dimension()!=a_r.dimension()) return false;
0238 typedef typename RMAT::elem_t RT;
0239 const RT* rpos = a_r.data();
0240 typedef typename CMAT::elem_t CT;
0241 CT* cpos = const_cast<CT*>(a_c.data());
0242 unsigned int D2 = a_r.dimension()*a_r.dimension();
0243 for(unsigned int i=0;i<D2;i++,rpos++,cpos++) *cpos = CT(*rpos);
0244 return true;
0245 }
0246
0247 template <class MAT>
0248 inline void dagger(MAT& a_m,typename MAT::elem_t (*a_conj)(const typename MAT::elem_t&)) {
0249 conjugate<MAT>(a_m,a_conj);
0250 a_m.transpose();
0251 }
0252
0253 template <class CMAT,class RMAT>
0254 inline bool decomplex(const CMAT& a_c,RMAT& a_r,
0255 typename CMAT::elem_t::value_type (*a_real)(const typename CMAT::elem_t&),
0256 typename CMAT::elem_t::value_type (*a_imag)(const typename CMAT::elem_t&)) {
0257 // CMAT = X+iY
0258 // RMAT = | X Y |
0259 // | -Y X |
0260 typedef typename CMAT::elem_t CT; //std::complex<double>
0261 typedef typename RMAT::elem_t RT; //double
0262
0263 unsigned int cdim = a_c.dimension();
0264 if(a_r.dimension()!=2*cdim) {a_r.set_zero();return false;}
0265 RT value;unsigned int r,c;
0266 for(r=0;r<cdim;r++) {
0267 for(c=0;c<cdim;c++) {
0268 const CT& cvalue = a_c.value(r,c);
0269 value = a_real(cvalue);
0270 a_r.set_value(r,c,value);
0271 a_r.set_value(r+cdim,c+cdim,value);
0272 value = a_imag(cvalue);
0273 a_r.set_value(r,c+cdim,value);
0274 a_r.set_value(r+cdim,c,-value);
0275 }
0276 }
0277 return true;
0278 }
0279
0280 template <class VEC_CMAT,class VEC_RMAT>
0281 inline bool decomplex(
0282 const VEC_CMAT& a_vc,VEC_RMAT& a_vr
0283 ,typename VEC_CMAT::value_type::elem_t::value_type (*a_real)(const typename VEC_CMAT::value_type::elem_t&)
0284 ,typename VEC_CMAT::value_type::elem_t::value_type (*a_imag)(const typename VEC_CMAT::value_type::elem_t&)
0285 ) {
0286 // CMAT = X+iY
0287 // RMAT = | X Y |
0288 // | -Y X |
0289 typedef typename VEC_CMAT::size_type sz_t;
0290 sz_t number = a_vc.size();
0291 a_vr.resize(number);
0292 for(sz_t index=0;index<number;index++) {
0293 if(!decomplex(a_vc[index],a_vr[index],a_real,a_imag)) {
0294 a_vr.clear();
0295 return false;
0296 }
0297 }
0298 return true;
0299 }
0300
0301 ///////////////////////////////////////////////////////////////////////////////////////
0302 ///////////////////////////////////////////////////////////////////////////////////////
0303 ///////////////////////////////////////////////////////////////////////////////////////
0304
0305 //for sf, mf :
0306 //template <class T,unsigned int D>
0307 //inline const T* get_data(const mat<T,D>& a_v) {return a_v.data();}
0308
0309 template <class MAT>
0310 inline MAT commutator(const MAT& a1,const MAT& a2) {
0311 return a1*a2-a2*a1;
0312 }
0313
0314 template <class MAT>
0315 inline MAT anticommutator(const MAT& a1,const MAT& a2) {
0316 return a1*a2+a2*a1;
0317 }
0318
0319 template <class MAT>
0320 inline void commutator(const MAT& a1,const MAT& a2,MAT& a_tmp,MAT& a_res) {
0321 a_res = a1;
0322 a_res *= a2;
0323 a_tmp = a2;
0324 a_tmp *= a1;
0325 a_res -= a_tmp;
0326 }
0327
0328 template <class MAT,class T>
0329 inline void commutator(const MAT& a1,const MAT& a2,MAT& a_tmp,T a_vtmp[],MAT& a_res) {
0330 a_res = a1;
0331 a_res.mul_mtx(a2,a_vtmp);
0332 a_tmp = a2;
0333 a_tmp.mul_mtx(a1,a_vtmp);
0334 a_res -= a_tmp;
0335 }
0336
0337 template <class MAT>
0338 inline void anticommutator(const MAT& a1,const MAT& a2,MAT& a_tmp,MAT& a_res) {
0339 a_res = a1;
0340 a_res *= a2;
0341 a_tmp = a2;
0342 a_tmp *= a1;
0343 a_res += a_tmp;
0344 }
0345
0346 template <class MAT,class T>
0347 inline void anticommutator(const MAT& a1,const MAT& a2,MAT& a_tmp,T a_vtmp[],MAT& a_res) {
0348 a_res = a1;
0349 a_res.mul_mtx(a2,a_vtmp);
0350 a_tmp = a2;
0351 a_tmp.mul_mtx(a1,a_vtmp);
0352 a_res += a_tmp;
0353 }
0354
0355 template <class T,unsigned int D>
0356 inline bool commutator_equal(const mat<T,D>& a_1,const mat<T,D>& a_2,const mat<T,D>& a_c,const T& a_prec) {
0357 unsigned int order = D;
0358 const T* p1 = a_1.data();
0359 const T* p2 = a_2.data();
0360 const T* pc = a_c.data();
0361 for(unsigned int r=0;r<order;r++) {
0362 for(unsigned int c=0;c<order;c++) {
0363 T _12 = T();
0364 {for(unsigned int i=0;i<order;i++) {
0365 _12 += (*(p1+r+i*order)) * (*(p2+i+c*order));
0366 }}
0367 T _21 = T();
0368 {for(unsigned int i=0;i<order;i++) {
0369 _21 += (*(p2+r+i*order)) * (*(p1+i+c*order));
0370 }}
0371 T diff = (_12-_21) - *(pc+r+c*order);
0372 if(diff<T()) diff *= T(-1);
0373 if(diff>=a_prec) return false;
0374 }
0375 }
0376 return true;
0377 }
0378
0379 template <class T,unsigned int D>
0380 inline bool anticommutator_equal(const mat<T,D>& a_1,const mat<T,D>& a_2,const mat<T,D>& a_c,const T& a_prec) {
0381 unsigned int order = D;
0382 const T* p1 = a_1.data();
0383 const T* p2 = a_2.data();
0384 const T* pc = a_c.data();
0385 for(unsigned int r=0;r<order;r++) {
0386 for(unsigned int c=0;c<order;c++) {
0387 T _12 = T();
0388 {for(unsigned int i=0;i<order;i++) {
0389 _12 += (*(p1+r+i*order)) * (*(p2+i+c*order));
0390 }}
0391 T _21 = T();
0392 {for(unsigned int i=0;i<order;i++) {
0393 _21 += (*(p2+r+i*order)) * (*(p1+i+c*order));
0394 }}
0395 T diff = (_12+_21) - *(pc+r+c*order);
0396 if(diff<T()) diff *= T(-1);
0397 if(diff>=a_prec) return false;
0398 }
0399 }
0400 return true;
0401 }
0402
0403 template <class T,unsigned int D>
0404 inline mat<T,D> operator-(const mat<T,D>& a1,const mat<T,D>& a2) {
0405 mat<T,D> res(a1);
0406 res -= a2;
0407 return res;
0408 }
0409 template <class T,unsigned int D>
0410 inline mat<T,D> operator+(const mat<T,D>& a1,const mat<T,D>& a2) {
0411 mat<T,D> res(a1);
0412 res += a2;
0413 return res;
0414 }
0415 template <class T,unsigned int D>
0416 inline mat<T,D> operator*(const mat<T,D>& a1,const mat<T,D>& a2) {
0417 mat<T,D> res(a1);
0418 res *= a2;
0419 return res;
0420 }
0421 template <class T,unsigned int D>
0422 inline mat<T,D> operator*(const T& a_fac,const mat<T,D>& a_m) {
0423 mat<T,D> res(a_m);
0424 res *= a_fac;
0425 return res;
0426 }
0427 /*
0428 template <class T,unsigned int D>
0429 inline mat<T,D> operator*(const mat<T,D>& a_m,const T& a_fac) {
0430 mat<T,D> res(a_m);
0431 res *= a_fac;
0432 return res;
0433 }
0434 */
0435
0436 template <class T>
0437 inline nmat<T> operator-(const nmat<T>& a1,const nmat<T>& a2) {
0438 nmat<T> res(a1);
0439 res -= a2;
0440 return res;
0441 }
0442 template <class T>
0443 inline nmat<T> operator+(const nmat<T>& a1,const nmat<T>& a2) {
0444 nmat<T> res(a1);
0445 res += a2;
0446 return res;
0447 }
0448 template <class T>
0449 inline nmat<T> operator*(const nmat<T>& a1,const nmat<T>& a2) {
0450 nmat<T> res(a1);
0451 res *= a2;
0452 return res;
0453 }
0454 template <class T>
0455 inline nmat<T> operator*(const T& a_fac,const nmat<T>& a_m) {
0456 nmat<T> res(a_m);
0457 res *= a_fac;
0458 return res;
0459 }
0460 /*
0461 template <class T>
0462 inline nmat<T> operator*(const nmat<T>& a_m,const T& a_fac) {
0463 nmat<T> res(a_m);
0464 res *= a_fac;
0465 return res;
0466 }
0467 */
0468
0469 template <class MAT,class REAL>
0470 inline bool mat_fabs(const MAT& a_in,MAT& a_ou,REAL(*a_fabs)(const typename MAT::elem_t&)) {
0471 if(a_in.dimension()!=a_ou.dimension()) {a_ou.set_zero();return false;}
0472 typedef typename MAT::elem_t T;
0473 T* in_pos = const_cast<T*>(a_in.data());
0474 T* ou_pos = const_cast<T*>(a_ou.data());
0475 unsigned int D2 = a_in.dimension()*a_in.dimension();
0476 for(unsigned int i=0;i<D2;i++,in_pos++,ou_pos++) {*ou_pos = a_fabs(*in_pos);}
0477 return true;
0478 }
0479
0480 }
0481
0482 /*
0483 #include <cstdarg>
0484
0485 namespace tools {
0486
0487 template <class MAT>
0488 inline void matrix_va_set(MAT& a_m,size_t a_dummy,...) {
0489 typedef typename MAT::elem_t T;
0490 va_list args;
0491 va_start(args,a_dummy);
0492 unsigned int D = a_m.dimension();
0493 for(unsigned int r=0;r<D;r++) {
0494 for(unsigned int c=0;c<D;c++) {
0495 a_m.set_value(r,c,va_arg(args,T));
0496 }
0497 }
0498 va_end(args);
0499 }
0500
0501 }
0502 */
0503
0504 namespace tools {
0505
0506 #define TOOLS_MELEM const typename MAT::elem_t&
0507
0508 ////////////////////////////////////////////////
0509 /// specific D=2 ///////////////////////////////
0510 ////////////////////////////////////////////////
0511 template <class MAT>
0512 inline void matrix_set(MAT& a_m
0513 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01
0514 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11
0515 ){
0516 //a_<R><C>
0517 //vec[R + C * 2];
0518 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0519 vec[0] = a_00;vec[2] = a_01;
0520 vec[1] = a_10;vec[3] = a_11;
0521 }
0522
0523 ////////////////////////////////////////////////
0524 /// specific D=3 ///////////////////////////////
0525 ////////////////////////////////////////////////
0526 template <class MAT>
0527 inline void matrix_set(MAT& a_m
0528 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02 //1 row
0529 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12 //2 row
0530 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22 //3 row
0531 ){
0532 //a_<R><C>
0533 //vec[R + C * 3];
0534 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0535 vec[0] = a_00;vec[3] = a_01;vec[6] = a_02;
0536 vec[1] = a_10;vec[4] = a_11;vec[7] = a_12;
0537 vec[2] = a_20;vec[5] = a_21;vec[8] = a_22;
0538 }
0539
0540 ////////////////////////////////////////////////
0541 /// specific D=4 ///////////////////////////////
0542 ////////////////////////////////////////////////
0543 template <class MAT>
0544 inline void matrix_set(MAT& a_m
0545 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02,TOOLS_MELEM a_03 //1 row
0546 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12,TOOLS_MELEM a_13 //2 row
0547 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22,TOOLS_MELEM a_23 //3 row
0548 ,TOOLS_MELEM a_30,TOOLS_MELEM a_31,TOOLS_MELEM a_32,TOOLS_MELEM a_33 //4 row
0549 ){
0550 //a_<R><C>
0551 //vec[R + C * 4];
0552 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0553 vec[0] = a_00;vec[4] = a_01;vec[ 8] = a_02;vec[12] = a_03;
0554 vec[1] = a_10;vec[5] = a_11;vec[ 9] = a_12;vec[13] = a_13;
0555 vec[2] = a_20;vec[6] = a_21;vec[10] = a_22;vec[14] = a_23;
0556 vec[3] = a_30;vec[7] = a_31;vec[11] = a_32;vec[15] = a_33;
0557 }
0558
0559 ////////////////////////////////////////////////
0560 /// specific D=5 ///////////////////////////////
0561 ////////////////////////////////////////////////
0562 template <class MAT>
0563 inline void matrix_set(MAT& a_m
0564 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02,TOOLS_MELEM a_03,TOOLS_MELEM a_04 //1 row
0565 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12,TOOLS_MELEM a_13,TOOLS_MELEM a_14 //2 row
0566 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22,TOOLS_MELEM a_23,TOOLS_MELEM a_24 //3 row
0567 ,TOOLS_MELEM a_30,TOOLS_MELEM a_31,TOOLS_MELEM a_32,TOOLS_MELEM a_33,TOOLS_MELEM a_34 //4 row
0568 ,TOOLS_MELEM a_40,TOOLS_MELEM a_41,TOOLS_MELEM a_42,TOOLS_MELEM a_43,TOOLS_MELEM a_44 //5 row
0569 ){
0570 //a_<R><C>
0571 //vec[R + C * 5];
0572 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0573 vec[0] = a_00;vec[5] = a_01;vec[10] = a_02;vec[15] = a_03;vec[20] = a_04;
0574 vec[1] = a_10;vec[6] = a_11;vec[11] = a_12;vec[16] = a_13;vec[21] = a_14;
0575 vec[2] = a_20;vec[7] = a_21;vec[12] = a_22;vec[17] = a_23;vec[22] = a_24;
0576 vec[3] = a_30;vec[8] = a_31;vec[13] = a_32;vec[18] = a_33;vec[23] = a_34;
0577 vec[4] = a_40;vec[9] = a_41;vec[14] = a_42;vec[19] = a_43;vec[24] = a_44;
0578 }
0579
0580 ////////////////////////////////////////////////
0581 /// specific D=6 ///////////////////////////////
0582 ////////////////////////////////////////////////
0583 template <class MAT>
0584 inline void matrix_set(MAT& a_m
0585 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02,TOOLS_MELEM a_03,TOOLS_MELEM a_04,TOOLS_MELEM a_05 //1 row
0586 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12,TOOLS_MELEM a_13,TOOLS_MELEM a_14,TOOLS_MELEM a_15 //2 row
0587 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22,TOOLS_MELEM a_23,TOOLS_MELEM a_24,TOOLS_MELEM a_25 //3 row
0588 ,TOOLS_MELEM a_30,TOOLS_MELEM a_31,TOOLS_MELEM a_32,TOOLS_MELEM a_33,TOOLS_MELEM a_34,TOOLS_MELEM a_35 //4 row
0589 ,TOOLS_MELEM a_40,TOOLS_MELEM a_41,TOOLS_MELEM a_42,TOOLS_MELEM a_43,TOOLS_MELEM a_44,TOOLS_MELEM a_45 //5 row
0590 ,TOOLS_MELEM a_50,TOOLS_MELEM a_51,TOOLS_MELEM a_52,TOOLS_MELEM a_53,TOOLS_MELEM a_54,TOOLS_MELEM a_55 //6 row
0591 ){
0592 //a_<R><C>
0593 //vec[R + C * 6];
0594 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0595 vec[0] = a_00;vec[ 6] = a_01;vec[12] = a_02;vec[18] = a_03;vec[24] = a_04;vec[30] = a_05;
0596 vec[1] = a_10;vec[ 7] = a_11;vec[13] = a_12;vec[19] = a_13;vec[25] = a_14;vec[31] = a_15;
0597 vec[2] = a_20;vec[ 8] = a_21;vec[14] = a_22;vec[20] = a_23;vec[26] = a_24;vec[32] = a_25;
0598 vec[3] = a_30;vec[ 9] = a_31;vec[15] = a_32;vec[21] = a_33;vec[27] = a_34;vec[33] = a_35;
0599 vec[4] = a_40;vec[10] = a_41;vec[16] = a_42;vec[22] = a_43;vec[28] = a_44;vec[34] = a_45;
0600 vec[5] = a_50;vec[11] = a_51;vec[17] = a_52;vec[23] = a_53;vec[29] = a_54;vec[35] = a_55;
0601 }
0602
0603 ////////////////////////////////////////////////
0604 /// specific D=10 //////////////////////////////
0605 ////////////////////////////////////////////////
0606 template <class MAT>
0607 inline void matrix_set(MAT& a_m
0608 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02,TOOLS_MELEM a_03,TOOLS_MELEM a_04,TOOLS_MELEM a_05,TOOLS_MELEM a_06,TOOLS_MELEM a_07,TOOLS_MELEM a_08,TOOLS_MELEM a_09 //1 row
0609 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12,TOOLS_MELEM a_13,TOOLS_MELEM a_14,TOOLS_MELEM a_15,TOOLS_MELEM a_16,TOOLS_MELEM a_17,TOOLS_MELEM a_18,TOOLS_MELEM a_19 //2 row
0610 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22,TOOLS_MELEM a_23,TOOLS_MELEM a_24,TOOLS_MELEM a_25,TOOLS_MELEM a_26,TOOLS_MELEM a_27,TOOLS_MELEM a_28,TOOLS_MELEM a_29 //3 row
0611 ,TOOLS_MELEM a_30,TOOLS_MELEM a_31,TOOLS_MELEM a_32,TOOLS_MELEM a_33,TOOLS_MELEM a_34,TOOLS_MELEM a_35,TOOLS_MELEM a_36,TOOLS_MELEM a_37,TOOLS_MELEM a_38,TOOLS_MELEM a_39 //4 row
0612 ,TOOLS_MELEM a_40,TOOLS_MELEM a_41,TOOLS_MELEM a_42,TOOLS_MELEM a_43,TOOLS_MELEM a_44,TOOLS_MELEM a_45,TOOLS_MELEM a_46,TOOLS_MELEM a_47,TOOLS_MELEM a_48,TOOLS_MELEM a_49 //5 row
0613 ,TOOLS_MELEM a_50,TOOLS_MELEM a_51,TOOLS_MELEM a_52,TOOLS_MELEM a_53,TOOLS_MELEM a_54,TOOLS_MELEM a_55,TOOLS_MELEM a_56,TOOLS_MELEM a_57,TOOLS_MELEM a_58,TOOLS_MELEM a_59 //6 row
0614 ,TOOLS_MELEM a_60,TOOLS_MELEM a_61,TOOLS_MELEM a_62,TOOLS_MELEM a_63,TOOLS_MELEM a_64,TOOLS_MELEM a_65,TOOLS_MELEM a_66,TOOLS_MELEM a_67,TOOLS_MELEM a_68,TOOLS_MELEM a_69 //7 row
0615 ,TOOLS_MELEM a_70,TOOLS_MELEM a_71,TOOLS_MELEM a_72,TOOLS_MELEM a_73,TOOLS_MELEM a_74,TOOLS_MELEM a_75,TOOLS_MELEM a_76,TOOLS_MELEM a_77,TOOLS_MELEM a_78,TOOLS_MELEM a_79 //8 row
0616 ,TOOLS_MELEM a_80,TOOLS_MELEM a_81,TOOLS_MELEM a_82,TOOLS_MELEM a_83,TOOLS_MELEM a_84,TOOLS_MELEM a_85,TOOLS_MELEM a_86,TOOLS_MELEM a_87,TOOLS_MELEM a_88,TOOLS_MELEM a_89 //9 row
0617 ,TOOLS_MELEM a_90,TOOLS_MELEM a_91,TOOLS_MELEM a_92,TOOLS_MELEM a_93,TOOLS_MELEM a_94,TOOLS_MELEM a_95,TOOLS_MELEM a_96,TOOLS_MELEM a_97,TOOLS_MELEM a_98,TOOLS_MELEM a_99 //10 row
0618
0619 ){
0620 //a_<R><C>
0621 //vec[R + C * 10];
0622 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0623 vec[0] = a_00;vec[10] = a_01;vec[20] = a_02;vec[30] = a_03;vec[40] = a_04;vec[50] = a_05;vec[60] = a_06;vec[70] = a_07;vec[80] = a_08;vec[90] = a_09;
0624 vec[1] = a_10;vec[11] = a_11;vec[21] = a_12;vec[31] = a_13;vec[41] = a_14;vec[51] = a_15;vec[61] = a_16;vec[71] = a_17;vec[81] = a_18;vec[91] = a_19;
0625 vec[2] = a_20;vec[12] = a_21;vec[22] = a_22;vec[32] = a_23;vec[42] = a_24;vec[52] = a_25;vec[62] = a_26;vec[72] = a_27;vec[82] = a_28;vec[92] = a_29;
0626 vec[3] = a_30;vec[13] = a_31;vec[23] = a_32;vec[33] = a_33;vec[43] = a_34;vec[53] = a_35;vec[63] = a_36;vec[73] = a_37;vec[83] = a_38;vec[93] = a_39;
0627 vec[4] = a_40;vec[14] = a_41;vec[24] = a_42;vec[34] = a_43;vec[44] = a_44;vec[54] = a_45;vec[64] = a_46;vec[74] = a_47;vec[84] = a_48;vec[94] = a_49;
0628 vec[5] = a_50;vec[15] = a_51;vec[25] = a_52;vec[35] = a_53;vec[45] = a_54;vec[55] = a_55;vec[65] = a_56;vec[75] = a_57;vec[85] = a_58;vec[95] = a_59;
0629 vec[6] = a_60;vec[16] = a_61;vec[26] = a_62;vec[36] = a_63;vec[46] = a_64;vec[56] = a_65;vec[66] = a_66;vec[76] = a_67;vec[86] = a_68;vec[96] = a_69;
0630 vec[7] = a_70;vec[17] = a_71;vec[27] = a_72;vec[37] = a_73;vec[47] = a_74;vec[57] = a_75;vec[67] = a_76;vec[77] = a_77;vec[87] = a_78;vec[97] = a_79;
0631 vec[8] = a_80;vec[18] = a_81;vec[28] = a_82;vec[38] = a_83;vec[48] = a_84;vec[58] = a_85;vec[68] = a_86;vec[78] = a_87;vec[88] = a_88;vec[98] = a_89;
0632 vec[9] = a_90;vec[19] = a_91;vec[29] = a_92;vec[39] = a_93;vec[49] = a_94;vec[59] = a_95;vec[69] = a_96;vec[79] = a_97;vec[89] = a_98;vec[99] = a_99;
0633 }
0634
0635 #undef TOOLS_MELEM
0636
0637 }
0638
0639 ////////////////////////////////////////////////
0640 ////////////////////////////////////////////////
0641 ////////////////////////////////////////////////
0642 #include <ostream>
0643
0644 namespace tools {
0645
0646 //NOTE : print is a Python keyword.
0647 template <class MAT>
0648 inline void dump(std::ostream& a_out,const std::string& aCMT,const MAT& a_matrix) {
0649 if(aCMT.size()) a_out << aCMT << std::endl;
0650 unsigned int D = a_matrix.dimension();
0651 for(unsigned int r=0;r<D;r++) {
0652 for(unsigned int c=0;c<D;c++) {
0653 a_out << " " << a_matrix.value(r,c);
0654 }
0655 a_out << std::endl;
0656 }
0657 }
0658
0659 template <class MAT>
0660 inline bool check_invert(const MAT& a_matrix,std::ostream& a_out) {
0661 MAT I;I.set_identity();
0662 MAT tmp;
0663 if(!a_matrix.invert(tmp)) return false;
0664 tmp.mul_mtx(a_matrix);
0665 if(!tmp.equal(I)) {
0666 dump(a_out,"problem with inv of :",a_matrix);
0667 return false;
0668 }
0669 return true;
0670 }
0671
0672 }
0673
0674 #endif
