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 #include "MATCOM"
0008
0009 namespace tools {
0010
0011 template <class T,unsigned int D>
0012 class mat {
0013 static const unsigned int _D2 = D*D;
0014 unsigned int dim2() const {return _D2;}
0015 #define TOOLS_MAT_CLASS mat
0016 TOOLS_MATCOM
0017 #undef TOOLS_MAT_CLASS
0018 private:
0019 public:
0020 mat() {
0021 #ifdef TOOLS_MAT_NEW
0022 m_vec = new T[D*D];
0023 #endif
0024 for(unsigned int i=0;i<_D2;i++) m_vec[i] = zero();
0025 }
0026 virtual ~mat() {
0027 #ifdef TOOLS_MAT_NEW
0028 delete [] m_vec;
0029 #endif
0030 }
0031 public:
0032 mat(const mat& a_from) {
0033 #ifdef TOOLS_MAT_NEW
0034 m_vec = new T[D*D];
0035 #endif
0036 _copy(a_from.m_vec);
0037 }
0038 mat& operator=(const mat& a_from){
0039 if(&a_from==this) return *this;
0040 _copy(a_from.m_vec);
0041 return *this;
0042 }
0043 public:
0044 mat(const T a_v[]){
0045 #ifdef TOOLS_MAT_NEW
0046 m_vec = new T[D*D];
0047 #endif
0048 _copy(a_v);
0049 }
0050 public:
0051 unsigned int dimension() const {return D;}
0052 protected:
0053 #ifdef TOOLS_MAT_NEW
0054 T* m_vec;
0055 #else
0056 T m_vec[D*D];
0057 #endif
0058
0059 private:static void check_instantiation() {mat<float,2> dummy;}
0060 };
0061
0062 template <class T>
0063 class nmat {
0064 unsigned int dim2() const {return m_D2;}
0065 #define TOOLS_MAT_CLASS nmat
0066 TOOLS_MATCOM
0067 #undef TOOLS_MAT_CLASS
0068 private:
0069 public:
0070 nmat(unsigned int a_D):m_D(a_D),m_D2(a_D*a_D),m_vec(0) {
0071 m_vec = new T[m_D2];
0072 for(unsigned int i=0;i<m_D2;i++) m_vec[i] = zero();
0073 }
0074 virtual ~nmat() {
0075 delete [] m_vec;
0076 }
0077 public:
0078 nmat(const nmat& a_from)
0079 :m_D(a_from.m_D),m_D2(a_from.m_D2),m_vec(0)
0080 {
0081 m_vec = new T[m_D2];
0082 _copy(a_from.m_vec);
0083 }
0084 nmat& operator=(const nmat& a_from){
0085 if(&a_from==this) return *this;
0086 if(a_from.m_D!=m_D) {
0087 m_D = a_from.m_D;
0088 m_D2 = a_from.m_D2;
0089 delete [] m_vec;
0090 m_vec = new T[m_D2];
0091 }
0092 _copy(a_from.m_vec);
0093 return *this;
0094 }
0095 public:
0096 nmat(unsigned int a_D,const T a_v[])
0097 :m_D(a_D),m_D2(a_D*a_D),m_vec(0)
0098 {
0099 m_vec = new T[m_D2];
0100 _copy(a_v);
0101 }
0102 public:
0103 unsigned int dimension() const {return m_D;}
0104 protected:
0105 unsigned int m_D;
0106 unsigned int m_D2;
0107 T* m_vec;
0108 private:static void check_instantiation() {nmat<float> dummy(2);}
0109 };
0110
0111 template <class T,unsigned int D>
0112 inline nmat<T> copy(const mat<T,D>& a_from) {
0113 unsigned int D2 = D*D;
0114 nmat<T> v(D);
0115 for(unsigned int i=0;i<D2;i++) v[i] = a_from[i];
0116 return v;
0117 }
0118
0119 template <class VECTOR>
0120 inline void multiply(VECTOR& a_vec,const typename VECTOR::value_type& a_mat) {
0121 typedef typename VECTOR::iterator it_t;
0122 for(it_t it=a_vec.begin();it!=a_vec.end();++it) *it *= a_mat;
0123 }
0124
0125 template <class VECTOR>
0126 inline void multiply(VECTOR& a_vec,const typename VECTOR::value_type::elem_t& a_value) {
0127 typedef typename VECTOR::iterator it_t;
0128 for(it_t it=a_vec.begin();it!=a_vec.end();++it) (*it).multiply(a_value);
0129 }
0130
0131 ///////////////////////////////////////////////////////////////////////////////////////
0132 /// related to complex numbers : //////////////////////////////////////////////////////
0133 ///////////////////////////////////////////////////////////////////////////////////////
0134
0135 template <class MAT>
0136 inline void conjugate(MAT& a_m,typename MAT::elem_t (*a_conj)(const typename MAT::elem_t&)) {
0137 typedef typename MAT::elem_t T;
0138 T* pos = const_cast<T*>(a_m.data());
0139 unsigned int D2 = a_m.dimension()*a_m.dimension();
0140 for(unsigned int i=0;i<D2;i++,pos++) {
0141 *pos = a_conj(*pos); //T = std::complex<>
0142 }
0143 }
0144
0145 template <class MAT>
0146 inline bool is_real(MAT& a_m,typename MAT::elem_t::value_type (*a_imag)(const typename MAT::elem_t&)) {
0147 typedef typename MAT::elem_t T;
0148 T* pos = const_cast<T*>(a_m.data());
0149 unsigned int D2 = a_m.dimension()*a_m.dimension();
0150 for(unsigned int i=0;i<D2;i++,pos++) {if(a_imag(*pos)) return false;}
0151 return true;
0152 }
0153
0154 template <class MAT,class PREC>
0155 inline bool is_real_prec(MAT& a_m,typename MAT::elem_t::value_type (*a_imag)(const typename MAT::elem_t&),
0156 const PREC& a_prec,PREC(*a_fabs)(const typename MAT::elem_t::value_type&)) {\
0157 typedef typename MAT::elem_t T;
0158 T* pos = const_cast<T*>(a_m.data());
0159 unsigned int D2 = a_m.dimension()*a_m.dimension();
0160 for(unsigned int i=0;i<D2;i++,pos++) {if(a_fabs(a_imag(*pos))>=a_prec) return false;}
0161 return true;
0162 }
0163
0164 template <class MAT>
0165 inline bool is_imag(MAT& a_m,typename MAT::elem_t::value_type (*a_real)(const typename MAT::elem_t&)) {
0166 typedef typename MAT::elem_t T;
0167 T* pos = const_cast<T*>(a_m.data());
0168 unsigned int D2 = a_m.dimension()*a_m.dimension();
0169 for(unsigned int i=0;i<D2;i++,pos++) {if(a_real(*pos)) return false;}
0170 return true;
0171 }
0172
0173 template <class CMAT,class RMAT>
0174 inline bool to_real(const CMAT& a_c,RMAT& a_r,typename CMAT::elem_t::value_type (*a_real)(const typename CMAT::elem_t&)) {
0175 if(a_r.dimension()!=a_c.dimension()) return false;
0176 typedef typename CMAT::elem_t CT;
0177 const CT* cpos = a_c.data();
0178 typedef typename RMAT::elem_t RT;
0179 RT* rpos = const_cast<RT*>(a_r.data());
0180 unsigned int D2 = a_c.dimension()*a_c.dimension();
0181 for(unsigned int i=0;i<D2;i++,cpos++,rpos++) *rpos = a_real(*cpos);
0182 return true;
0183 }
0184
0185 template <class RMAT,class CMAT>
0186 inline bool to_complex(const RMAT& a_r,CMAT& a_c) {
0187 if(a_c.dimension()!=a_r.dimension()) return false;
0188 typedef typename RMAT::elem_t RT;
0189 const RT* rpos = a_r.data();
0190 typedef typename CMAT::elem_t CT;
0191 CT* cpos = const_cast<CT*>(a_c.data());
0192 unsigned int D2 = a_r.dimension()*a_r.dimension();
0193 for(unsigned int i=0;i<D2;i++,rpos++,cpos++) *cpos = CT(*rpos);
0194 return true;
0195 }
0196
0197 template <class MAT>
0198 inline void dagger(MAT& a_m,typename MAT::elem_t (*a_conj)(const typename MAT::elem_t&)) {
0199 conjugate<MAT>(a_m,a_conj);
0200 a_m.transpose();
0201 }
0202
0203 template <class CMAT,class RMAT>
0204 inline bool decomplex(const CMAT& a_c,RMAT& a_r,
0205 typename CMAT::elem_t::value_type (*a_real)(const typename CMAT::elem_t&),
0206 typename CMAT::elem_t::value_type (*a_imag)(const typename CMAT::elem_t&)) {
0207 // CMAT = X+iY
0208 // RMAT = | X Y |
0209 // | -Y X |
0210 typedef typename CMAT::elem_t CT; //std::complex<double>
0211 typedef typename RMAT::elem_t RT; //double
0212
0213 unsigned int cdim = a_c.dimension();
0214 if(a_r.dimension()!=2*cdim) {a_r.set_zero();return false;}
0215 RT value;unsigned int r,c;
0216 for(r=0;r<cdim;r++) {
0217 for(c=0;c<cdim;c++) {
0218 const CT& cvalue = a_c.value(r,c);
0219 value = a_real(cvalue);
0220 a_r.set_value(r,c,value);
0221 a_r.set_value(r+cdim,c+cdim,value);
0222 value = a_imag(cvalue);
0223 a_r.set_value(r,c+cdim,value);
0224 a_r.set_value(r+cdim,c,-value);
0225 }
0226 }
0227 return true;
0228 }
0229
0230 template <class VEC_CMAT,class VEC_RMAT>
0231 inline bool decomplex(
0232 const VEC_CMAT& a_vc,VEC_RMAT& a_vr
0233 ,typename VEC_CMAT::value_type::elem_t::value_type (*a_real)(const typename VEC_CMAT::value_type::elem_t&)
0234 ,typename VEC_CMAT::value_type::elem_t::value_type (*a_imag)(const typename VEC_CMAT::value_type::elem_t&)
0235 ) {
0236 // CMAT = X+iY
0237 // RMAT = | X Y |
0238 // | -Y X |
0239 typedef typename VEC_CMAT::size_type sz_t;
0240 sz_t number = a_vc.size();
0241 a_vr.resize(number);
0242 for(sz_t index=0;index<number;index++) {
0243 if(!decomplex(a_vc[index],a_vr[index],a_real,a_imag)) {
0244 a_vr.clear();
0245 return false;
0246 }
0247 }
0248 return true;
0249 }
0250
0251 ///////////////////////////////////////////////////////////////////////////////////////
0252 ///////////////////////////////////////////////////////////////////////////////////////
0253 ///////////////////////////////////////////////////////////////////////////////////////
0254
0255 template <class MAT>
0256 inline MAT commutator(const MAT& a1,const MAT& a2) {
0257 return a1*a2-a2*a1;
0258 }
0259
0260 template <class MAT>
0261 inline MAT anticommutator(const MAT& a1,const MAT& a2) {
0262 return a1*a2+a2*a1;
0263 }
0264
0265 template <class MAT>
0266 inline void commutator(const MAT& a1,const MAT& a2,MAT& a_tmp,MAT& a_res) {
0267 a_res = a1;
0268 a_res *= a2;
0269 a_tmp = a2;
0270 a_tmp *= a1;
0271 a_res -= a_tmp;
0272 }
0273
0274 template <class MAT,class T>
0275 inline void commutator(const MAT& a1,const MAT& a2,MAT& a_tmp,T a_vtmp[],MAT& a_res) {
0276 a_res = a1;
0277 a_res.mul_mtx(a2,a_vtmp);
0278 a_tmp = a2;
0279 a_tmp.mul_mtx(a1,a_vtmp);
0280 a_res -= a_tmp;
0281 }
0282
0283 template <class MAT>
0284 inline void anticommutator(const MAT& a1,const MAT& a2,MAT& a_tmp,MAT& a_res) {
0285 a_res = a1;
0286 a_res *= a2;
0287 a_tmp = a2;
0288 a_tmp *= a1;
0289 a_res += a_tmp;
0290 }
0291
0292 template <class MAT,class T>
0293 inline void anticommutator(const MAT& a1,const MAT& a2,MAT& a_tmp,T a_vtmp[],MAT& a_res) {
0294 a_res = a1;
0295 a_res.mul_mtx(a2,a_vtmp);
0296 a_tmp = a2;
0297 a_tmp.mul_mtx(a1,a_vtmp);
0298 a_res += a_tmp;
0299 }
0300
0301 template <class T,unsigned int D>
0302 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) {
0303 unsigned int order = D;
0304 const T* p1 = a_1.data();
0305 const T* p2 = a_2.data();
0306 const T* pc = a_c.data();
0307 for(unsigned int r=0;r<order;r++) {
0308 for(unsigned int c=0;c<order;c++) {
0309 T _12 = T();
0310 {for(unsigned int i=0;i<order;i++) {
0311 _12 += (*(p1+r+i*order)) * (*(p2+i+c*order));
0312 }}
0313 T _21 = T();
0314 {for(unsigned int i=0;i<order;i++) {
0315 _21 += (*(p2+r+i*order)) * (*(p1+i+c*order));
0316 }}
0317 T diff = (_12-_21) - *(pc+r+c*order);
0318 if(diff<T()) diff *= T(-1);
0319 if(diff>=a_prec) return false;
0320 }
0321 }
0322 return true;
0323 }
0324
0325 template <class T,unsigned int D>
0326 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) {
0327 unsigned int order = D;
0328 const T* p1 = a_1.data();
0329 const T* p2 = a_2.data();
0330 const T* pc = a_c.data();
0331 for(unsigned int r=0;r<order;r++) {
0332 for(unsigned int c=0;c<order;c++) {
0333 T _12 = T();
0334 {for(unsigned int i=0;i<order;i++) {
0335 _12 += (*(p1+r+i*order)) * (*(p2+i+c*order));
0336 }}
0337 T _21 = T();
0338 {for(unsigned int i=0;i<order;i++) {
0339 _21 += (*(p2+r+i*order)) * (*(p1+i+c*order));
0340 }}
0341 T diff = (_12+_21) - *(pc+r+c*order);
0342 if(diff<T()) diff *= T(-1);
0343 if(diff>=a_prec) return false;
0344 }
0345 }
0346 return true;
0347 }
0348
0349 template <class T,unsigned int D>
0350 inline mat<T,D> operator-(const mat<T,D>& a1,const mat<T,D>& a2) {
0351 mat<T,D> res(a1);
0352 res -= a2;
0353 return res;
0354 }
0355 template <class T,unsigned int D>
0356 inline mat<T,D> operator+(const mat<T,D>& a1,const mat<T,D>& a2) {
0357 mat<T,D> res(a1);
0358 res += a2;
0359 return res;
0360 }
0361 template <class T,unsigned int D>
0362 inline mat<T,D> operator*(const mat<T,D>& a1,const mat<T,D>& a2) {
0363 mat<T,D> res(a1);
0364 res *= a2;
0365 return res;
0366 }
0367 template <class T,unsigned int D>
0368 inline mat<T,D> operator*(const T& a_fac,const mat<T,D>& a_m) {
0369 mat<T,D> res(a_m);
0370 res *= a_fac;
0371 return res;
0372 }
0373
0374 template <class T>
0375 inline nmat<T> operator-(const nmat<T>& a1,const nmat<T>& a2) {
0376 nmat<T> res(a1);
0377 res -= a2;
0378 return res;
0379 }
0380 template <class T>
0381 inline nmat<T> operator+(const nmat<T>& a1,const nmat<T>& a2) {
0382 nmat<T> res(a1);
0383 res += a2;
0384 return res;
0385 }
0386 template <class T>
0387 inline nmat<T> operator*(const nmat<T>& a1,const nmat<T>& a2) {
0388 nmat<T> res(a1);
0389 res *= a2;
0390 return res;
0391 }
0392 template <class T>
0393 inline nmat<T> operator*(const T& a_fac,const nmat<T>& a_m) {
0394 nmat<T> res(a_m);
0395 res *= a_fac;
0396 return res;
0397 }
0398
0399 template <class MAT,class REAL>
0400 inline bool mat_fabs(const MAT& a_in,MAT& a_ou,REAL(*a_fabs)(const typename MAT::elem_t&)) {
0401 if(a_in.dimension()!=a_ou.dimension()) {a_ou.set_zero();return false;}
0402 typedef typename MAT::elem_t T;
0403 T* in_pos = const_cast<T*>(a_in.data());
0404 T* ou_pos = const_cast<T*>(a_ou.data());
0405 unsigned int D2 = a_in.dimension()*a_in.dimension();
0406 for(unsigned int i=0;i<D2;i++,in_pos++,ou_pos++) {*ou_pos = a_fabs(*in_pos);}
0407 return true;
0408 }
0409
0410 }
0411
0412 namespace tools {
0413
0414 #define TOOLS_MELEM const typename MAT::elem_t&
0415
0416 ////////////////////////////////////////////////
0417 /// specific D=2 ///////////////////////////////
0418 ////////////////////////////////////////////////
0419 template <class MAT>
0420 inline void matrix_set(MAT& a_m
0421 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01
0422 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11
0423 ){
0424 //a_<R><C>
0425 //vec[R + C * 2];
0426 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0427 vec[0] = a_00;vec[2] = a_01;
0428 vec[1] = a_10;vec[3] = a_11;
0429 }
0430
0431 ////////////////////////////////////////////////
0432 /// specific D=3 ///////////////////////////////
0433 ////////////////////////////////////////////////
0434 template <class MAT>
0435 inline void matrix_set(MAT& a_m
0436 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02 //1 row
0437 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12 //2 row
0438 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22 //3 row
0439 ){
0440 //a_<R><C>
0441 //vec[R + C * 3];
0442 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0443 vec[0] = a_00;vec[3] = a_01;vec[6] = a_02;
0444 vec[1] = a_10;vec[4] = a_11;vec[7] = a_12;
0445 vec[2] = a_20;vec[5] = a_21;vec[8] = a_22;
0446 }
0447
0448 ////////////////////////////////////////////////
0449 /// specific D=4 ///////////////////////////////
0450 ////////////////////////////////////////////////
0451 template <class MAT>
0452 inline void matrix_set(MAT& a_m
0453 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02,TOOLS_MELEM a_03 //1 row
0454 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12,TOOLS_MELEM a_13 //2 row
0455 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22,TOOLS_MELEM a_23 //3 row
0456 ,TOOLS_MELEM a_30,TOOLS_MELEM a_31,TOOLS_MELEM a_32,TOOLS_MELEM a_33 //4 row
0457 ){
0458 //a_<R><C>
0459 //vec[R + C * 4];
0460 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0461 vec[0] = a_00;vec[4] = a_01;vec[ 8] = a_02;vec[12] = a_03;
0462 vec[1] = a_10;vec[5] = a_11;vec[ 9] = a_12;vec[13] = a_13;
0463 vec[2] = a_20;vec[6] = a_21;vec[10] = a_22;vec[14] = a_23;
0464 vec[3] = a_30;vec[7] = a_31;vec[11] = a_32;vec[15] = a_33;
0465 }
0466
0467 ////////////////////////////////////////////////
0468 /// specific D=5 ///////////////////////////////
0469 ////////////////////////////////////////////////
0470 template <class MAT>
0471 inline void matrix_set(MAT& a_m
0472 ,TOOLS_MELEM a_00,TOOLS_MELEM a_01,TOOLS_MELEM a_02,TOOLS_MELEM a_03,TOOLS_MELEM a_04 //1 row
0473 ,TOOLS_MELEM a_10,TOOLS_MELEM a_11,TOOLS_MELEM a_12,TOOLS_MELEM a_13,TOOLS_MELEM a_14 //2 row
0474 ,TOOLS_MELEM a_20,TOOLS_MELEM a_21,TOOLS_MELEM a_22,TOOLS_MELEM a_23,TOOLS_MELEM a_24 //3 row
0475 ,TOOLS_MELEM a_30,TOOLS_MELEM a_31,TOOLS_MELEM a_32,TOOLS_MELEM a_33,TOOLS_MELEM a_34 //4 row
0476 ,TOOLS_MELEM a_40,TOOLS_MELEM a_41,TOOLS_MELEM a_42,TOOLS_MELEM a_43,TOOLS_MELEM a_44 //5 row
0477 ){
0478 //a_<R><C>
0479 //vec[R + C * 5];
0480 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0481 vec[0] = a_00;vec[5] = a_01;vec[10] = a_02;vec[15] = a_03;vec[20] = a_04;
0482 vec[1] = a_10;vec[6] = a_11;vec[11] = a_12;vec[16] = a_13;vec[21] = a_14;
0483 vec[2] = a_20;vec[7] = a_21;vec[12] = a_22;vec[17] = a_23;vec[22] = a_24;
0484 vec[3] = a_30;vec[8] = a_31;vec[13] = a_32;vec[18] = a_33;vec[23] = a_34;
0485 vec[4] = a_40;vec[9] = a_41;vec[14] = a_42;vec[19] = a_43;vec[24] = a_44;
0486 }
0487
0488 ////////////////////////////////////////////////
0489 /// specific D=6 ///////////////////////////////
0490 ////////////////////////////////////////////////
0491 template <class MAT>
0492 inline void matrix_set(MAT& a_m
0493 ,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
0494 ,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
0495 ,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
0496 ,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
0497 ,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
0498 ,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
0499 ){
0500 //a_<R><C>
0501 //vec[R + C * 6];
0502 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0503 vec[0] = a_00;vec[ 6] = a_01;vec[12] = a_02;vec[18] = a_03;vec[24] = a_04;vec[30] = a_05;
0504 vec[1] = a_10;vec[ 7] = a_11;vec[13] = a_12;vec[19] = a_13;vec[25] = a_14;vec[31] = a_15;
0505 vec[2] = a_20;vec[ 8] = a_21;vec[14] = a_22;vec[20] = a_23;vec[26] = a_24;vec[32] = a_25;
0506 vec[3] = a_30;vec[ 9] = a_31;vec[15] = a_32;vec[21] = a_33;vec[27] = a_34;vec[33] = a_35;
0507 vec[4] = a_40;vec[10] = a_41;vec[16] = a_42;vec[22] = a_43;vec[28] = a_44;vec[34] = a_45;
0508 vec[5] = a_50;vec[11] = a_51;vec[17] = a_52;vec[23] = a_53;vec[29] = a_54;vec[35] = a_55;
0509 }
0510
0511 ////////////////////////////////////////////////
0512 /// specific D=10 //////////////////////////////
0513 ////////////////////////////////////////////////
0514 template <class MAT>
0515 inline void matrix_set(MAT& a_m
0516 ,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
0517 ,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
0518 ,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
0519 ,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
0520 ,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
0521 ,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
0522 ,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
0523 ,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
0524 ,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
0525 ,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
0526
0527 ){
0528 //a_<R><C>
0529 //vec[R + C * 10];
0530 typename MAT::elem_t* vec = const_cast<typename MAT::elem_t*>(a_m.data());
0531 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;
0532 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;
0533 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;
0534 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;
0535 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;
0536 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;
0537 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;
0538 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;
0539 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;
0540 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;
0541 }
0542
0543 #undef TOOLS_MELEM
0544
0545 }
0546
0547 ////////////////////////////////////////////////
0548 ////////////////////////////////////////////////
0549 ////////////////////////////////////////////////
0550 #include <ostream>
0551
0552 namespace tools {
0553
0554 //NOTE : print is a Python keyword.
0555 template <class MAT>
0556 inline void dump(std::ostream& a_out,const std::string& aCMT,const MAT& a_matrix) {
0557 if(aCMT.size()) a_out << aCMT << std::endl;
0558 unsigned int D = a_matrix.dimension();
0559 for(unsigned int r=0;r<D;r++) {
0560 for(unsigned int c=0;c<D;c++) {
0561 a_out << " " << a_matrix.value(r,c);
0562 }
0563 a_out << std::endl;
0564 }
0565 }
0566
0567 template <class MAT>
0568 inline bool check_invert(const MAT& a_matrix,std::ostream& a_out) {
0569 MAT I;I.set_identity();
0570 MAT tmp;
0571 if(!a_matrix.invert(tmp)) return false;
0572 tmp.mul_mtx(a_matrix);
0573 if(!tmp.equal(I)) {
0574 dump(a_out,"problem with inv of :",a_matrix);
0575 return false;
0576 }
0577 return true;
0578 }
0579
0580 }
0581
0582 #endif
