/Users/barrand/private/dev/softinex/old/inexlib-1.2/inlib/inlib/sqm

Go to the documentation of this file.

// Copyright (C) 2010, Guy Barrand. All rights reserved.
// See the file inlib.license for terms.

#ifndef inlib_sqm
#define inlib_sqm

// sqm is for square matrix

#include "array"

namespace inlib {

template <class T>
class sqm : public array<T> {
public:
  sqm(unsigned a_order):array<T>(2,a_order){
    // a_order is the matrix dimension, which must not be confused
    // with the array dimension which is 2).
  }
  virtual ~sqm() {}
public: 
  sqm(const sqm& a_from):array<T>(a_from){}
  sqm& operator=(const sqm& a_from) {
    array<T>::operator=(a_from);
    return *this;
  }
public:
  bool set_order(unsigned int a_order) {
    std::vector<unsigned int> orders(2);
    for(unsigned int index=0;index<2;index++) orders[index] = a_order;
    return array<T>::configure(orders);
  }

  unsigned int order() const { return array<T>::m_orders[0];}

  T value(unsigned int aR,unsigned int aC) const { 
    //WARNING : range of aR, aC not protected.
    unsigned int ord = order();
    return array<T>::m_vector[aR + aC * ord];
  }

  bool value(unsigned int aR,unsigned int aC,T& a_value) const { 
    unsigned int ord = order();
    if((aR>=ord)||(aC>=ord)) {
      a_value = array<T>::zero();
      return false;
    }
    a_value = array<T>::m_vector[aR + aC * ord];
    return true;
  }

  bool set_value(unsigned int aR,unsigned int aC,const T& a_value) { 
    unsigned int ord = order();
    if((aR>=ord)||(aC>=ord)) return false;
    array<T>::m_vector[aR + aC * ord] = a_value;
    return true;
  }

  void set_value_no_check(unsigned int aR,
                              unsigned int aC,
                              const T& a_value){ 
    array<T>::m_vector[aR + aC * order()] = a_value;
  }

  T trace() const {
    unsigned int ord = order();
    T value = array<T>::zero();
    for(unsigned int c=0;c<ord;c++) value += array<T>::m_vector[c+c*ord];
    return value;
  }

  bool mx_mul(const sqm<T>& a_matrix,sqm<T>& a_result) const {
    // a_result = this * a_matrix
    unsigned int ord = order();
    if(a_matrix.order()!=ord) return false;
    if(a_result.order()!=ord) return false;
    const T* p = &(array<T>::m_vector[0]);            //WARNING : dangerous.
    const T* a_p = &(a_matrix.array<T>::m_vector[0]); //WARNING : dangerous.
    T* r_p = &(a_result.array<T>::m_vector[0]);       //WARNING : dangerous.
    for(unsigned int r=0;r<ord;r++) {
      for(unsigned int c=0;c<ord;c++) {
        T value = array<T>::zero();
        for(unsigned int i=0;i<ord;i++) {
          //value += 
          //  array<T>::m_vector[r+i*ord] * 
          //  a_matrix.array<T>::m_vector[i+c*ord];
          value += (*(p+r+i*ord)) * (*(a_p+i+c*ord)); //optimize.
        }
        //a_result.array<T>::m_vector[r+c*ord] = value;
        *(r_p+r+c*ord) = value;
      }
    }
    return true;
  }

  bool mx_mul(const sqm<T>& a_matrix) {
    sqm<T> res(order());
    bool status = this->mx_mul(a_matrix,res);
    *this = res;
    return status;
  }

  bool mx_left_mul(const sqm<T>& a_matrix,sqm<T>& a_result) const {
    // a_result = a_matrix * this
    unsigned int ord = order();
    if(a_matrix.order()!=ord) return false;
    if(a_result.order()!=ord) return false;
    const T* p = &(array<T>::m_vector[0]);            //WARNING : dangerous.
    const T* a_p = &(a_matrix.array<T>::m_vector[0]); //WARNING : dangerous.
    T* r_p = &(a_result.array<T>::m_vector[0]);       //WARNING : dangerous.
    for(unsigned int r=0;r<ord;r++) {
      for(unsigned int c=0;c<ord;c++) {
        T value = array<T>::zero();
        for(unsigned int i=0;i<ord;i++) {
          //value += 
          //  a_matrix.array<T>::m_vector[r+i*ord] * 
          //  array<T>::m_vector[i+c*ord];
          value += (*(a_p+r+i*ord)) * (*(p+i+c*ord)); //optimize.
        }
        //a_result.array<T>::m_vector[r+c*ord] = value;
        *(r_p+r+c*ord) = value;
      }
    }
    return true;
  }

  bool mx_left_mul(const sqm<T>& a_matrix) {
    sqm<T> res(order());
    bool status = this->mx_left_mul(a_matrix,res);
    *this = res;
    return status;
  }

  virtual T determinant() const {
    unsigned int ord = order();
    if(ord==0) return array<T>::zero();
    else if(ord==1) return value(0,0);
    else if(ord==2) return (value(0,0) * value(1,1) - 
                            value(1,0) * value(0,1)); //Optimize

    unsigned int rord = ord-1;
    std::vector<unsigned int> cs(rord);
    std::vector<unsigned int> rs(rord);

    T det = array<T>::zero();
   {for(unsigned int i=0;i<rord;i++) {cs[i] = i+1;}}
    unsigned int c = 0;
    //if(c>=1) cs[c-1] = c-1;

   {for(unsigned int i=0;i<rord;i++) {rs[i] = i+1;}}
    bool sg = true; //c=0+r=0
    for(unsigned int r=0;r<ord;r++) {
      if(r>=1) rs[r-1] = r-1;
      T subdet = sub_determinant(rs,cs);
      if(sg) 
        det += value(r,c) * subdet;
      else
        det -= value(r,c) * subdet;
      sg = sg?false:true;
    }

    return det;
  }
  T sub_determinant(const std::vector<unsigned int>& aRs,
                          const std::vector<unsigned int>& aCs) const {
    //WARNING : to optimize, we do not check the content of aRs, aCs.
    unsigned int ord = aRs.size();
    if(ord==0) return array<T>::zero();
    else if(ord==1) return value(aRs[0],aCs[0]);
    else if(ord==2) {
      //return (value(aRs[0],aCs[0]) * value(aRs[1],aCs[1]) -
      //        value(aRs[1],aCs[0]) * value(aRs[0],aCs[1])); 
      //Optimize the upper :

      unsigned int r_0 = aRs[0];
      unsigned int r_1 = aRs[1];
      unsigned int c_0 = aCs[0];
      unsigned int c_1 = aCs[1];

      unsigned int ord = order();
      const T* p = &(array<T>::m_vector[0]); //WARNING : dangerous.

      return ( (*(p+r_0+c_0*ord)) * (*(p+r_1+c_1*ord)) -
               (*(p+r_1+c_0*ord)) * (*(p+r_0+c_1*ord)) );
    }

    unsigned int rord = ord-1;
    std::vector<unsigned int> cs(rord);
    std::vector<unsigned int> rs(rord);

    T det = array<T>::zero();
   {for(unsigned int i=0;i<rord;i++) {cs[i] = aCs[i+1];}}
    unsigned int c = 0;
    //if(c>=1) cs[c-1] = c-1;

   {for(unsigned int i=0;i<rord;i++) {rs[i] = aRs[i+1];}}
    bool sg = true; //c=0+r=0
    for(unsigned int r=0;r<ord;r++) {
      if(r>=1) rs[r-1] = aRs[r-1];
      T subdet = sub_determinant(rs,cs);
      if(sg)
        det += value(aRs[r],aCs[c]) * subdet;
      else
        det -= value(aRs[r],aCs[c]) * subdet;
      sg = sg?false:true;
    }

    return det;
  }

  bool invert(sqm<T>& a_result) const {
    //Generic invertion method.
    unsigned int ord = order();
    if(ord==0) return true;

    a_result.set_order(ord);
    if(ord==1) {
      T v;
      if(!value(0,0,v)) return false;
      if(v==array<T>::zero()) return false;
      if(!a_result.set_value(0,0,1./v)) return false;
      return true;
    }
  
    unsigned int rord = ord-1;
    std::vector<unsigned int> cs(rord);
    std::vector<unsigned int> rs(rord);
  
    // Get det with r = 0;
    T det = array<T>::zero();
   {
   {for(unsigned int i=0;i<rord;i++) {rs[i] = i+1;}}
    unsigned int r = 0;
    //if(r>=1) rs[r-1] = r-1;
  
   {for(unsigned int i=0;i<rord;i++) {cs[i] = i+1;}}
    bool sg = true; //r=0+c=0
    for(unsigned int c=0;c<ord;c++) {
      if(c>=1) cs[c-1] = c-1;
      T subdet = sub_determinant(rs,cs);
      T sgn = sg ? array<T>::one() : array<T>::minus_one();
      det += value(r,c) * subdet * sgn;
      T value = subdet * sgn;
      a_result.set_value_no_check(c,r,value);
      sg = sg?false:true;
    }}
  
    if(det==array<T>::zero()) return false;
  
   {for(unsigned int c=0;c<ord;c++) {
      a_result.set_value(c,0,a_result.value(c,0)/det);
    }}
  
   {for(unsigned int i=0;i<rord;i++) {rs[i] = i+1;}}
    bool sgr = false; //r=1+c=0
    for(unsigned int r=1;r<ord;r++) {
      if(r>=1) rs[r-1] = r-1;
     {for(unsigned int i=0;i<rord;i++) {cs[i] = i+1;}}
      bool sg = sgr;
      for(unsigned int c=0;c<ord;c++) {
        if(c>=1) cs[c-1] = c-1;
        T subdet = sub_determinant(rs,cs);
        T sgn = sg ? array<T>::one() : array<T>::minus_one();
        T value = (subdet * sgn)/det;
        a_result.set_value_no_check(c,r,value);
        sg = sg?false:true;
      }
      sgr = sgr?false:true;
    }
  
    return true;
  }

  void transpose() {
    unsigned int ord = order();
    for(unsigned int r=0;r<ord;r++) {
      for(unsigned int c=(r+1);c<ord;c++) {
        T vrc = value(r,c);
        T vcr = value(c,r);
        set_value_no_check(r,c,vcr);
        set_value_no_check(c,r,vrc);
      }
    }
  }

  bool is_symmetric() const {
    unsigned int ord = order();
    for(unsigned int r=0;r<ord;r++) {
      for(unsigned int c=(r+1);c<ord;c++) {
        if(value(r,c)!=value(c,r)) return false;
      }
    }
    return true;
  }
  bool is_antisymmetric() const {
    unsigned int ord = order();
   {for(unsigned int r=0;r<ord;r++) {
      if(value(r,r)!=array<T>::zero()) return false;
    }}
    for(unsigned int r=0;r<ord;r++) {
      for(unsigned int c=(r+1);c<ord;c++) {
        if(value(r,c)!=-value(c,r)) return false;
      }
    }
    return true;
  }

  void set_identity() {
    unsigned int ord = order();
    array<T>::reset();
    for(unsigned int c=0;c<ord;c++) set_value_no_check(c,c,array<T>::one());
  }

  void set_diag(const T& a_value) {
    unsigned int ord = order();
    for(unsigned int c=0;c<ord;c++) set_value_no_check(c,c,a_value);
  }

  //Factories :
  sqm symmetric_part() const {
    sqm res(*this);
    res.transpose();
    res.add(*this);
    res.divide(array<T>::two());
    return res;
  }

  sqm antisymmetric_part() const {
    sqm res(*this);
    res.transpose();
    res.multiply(array<T>::minus_one());
    res.add(*this);
    res.divide(array<T>::two());
    return res;
  }

  sqm& operator*=(const sqm<T>& a_from) {
    sqm<T> res(order());
    /*bool status =*/ this->mx_mul(a_from,res); //FIXME : throw
    *this = res;
    return *this;
  }
  sqm& operator+=(const sqm<T>& a_from) {
    /*bool status =*/ array<T>::add(a_from); //FIXME : throw
    return *this;
  }
  sqm& operator-=(const sqm<T>& a_from) {
    /*bool status =*/ array<T>::subtract(a_from); //FIXME : throw
    return *this;
  }
};

template <class T>
inline sqm<T> operator-(const sqm<T>& a1,const sqm<T>& a2) {
  sqm<T> res(a1);
  res.subtract(a2);
  return res;
}
template <class T>
inline sqm<T> operator+(const sqm<T>& a1,const sqm<T>& a2) {
  sqm<T> res(a1);
  res.add(a2);
  return res;
}
template <class T>
inline sqm<T> operator*(const sqm<T>& a1,const sqm<T>& a2) {
  sqm<T> res(a1.order());
  a1.mx_mul(a2,res);
  return res;
}

template <class T>
inline sqm<T> commutator(const sqm<T>& a1,const sqm<T>& a2) {
  return a1*a2-a2*a1;
}

template <class T>
inline sqm<T> anticommutator(const sqm<T>& a1,const sqm<T>& a2) {
  return a1*a2+a2*a1;
}

template <class T>
inline sqm<T> exp(const sqm<T>& a_matrix,unsigned int aNumber = 100) {
  // result = I + M + M**2/2! + M**3/3! + ....
  sqm<T> result(a_matrix.order());
  result.set_identity();
  sqm<T> m(a_matrix.order());     
  m.set_identity();
  for(unsigned int i=1;i<aNumber;i++) {
    m = m * a_matrix; 
    m.multiply(1./double(i)); 
    result.add(m);
  }
  return result;
}

template <class T>
inline sqm<T> log(const sqm<T>& a_matrix,unsigned int aNumber = 100) {
  // result = (M-I) - (M-I)**2/2 + (M-I)**3/3 +...
  // Warning : cannot converge...
  unsigned int order = a_matrix.order();
  sqm<T> result(order);
  sqm<T> I(order);     
  I.set_identity();
  sqm<T> x(a_matrix);     
  x.subtract(I);
  sqm<T> m(I);
  sqm<T> tmp(order);
  double fact = -1;
  for(unsigned int i=0;i<aNumber;i++) {
    m = m * x; 
    fact *= -1;
    tmp = m;
    tmp.multiply(fact/double(i+1)); 
    result.add(tmp);
  }
  return result;
}

#ifdef __CINT__
// pb with the T in double(*a_func)(T).
#else
template <class T>
inline sqm<double> abs(const sqm<T>& a_matrix,double(*a_func)(T)) {
  const std::vector<T>& avec = a_matrix.vector();
  sqm<double> result(a_matrix.order());
  std::vector<double>& rvec = result.vector();
  unsigned int number = avec.size();
  for(unsigned int index=0;index<number;index++) 
    rvec[index] = a_func(avec[index]);
  return result;
}

template <class T>
inline bool is_log_convergent(const sqm<T>& a_matrix,double(*a_func)(T)) {
  unsigned int order = a_matrix.order();
  sqm<T> I(order);     
  I.set_identity();
  sqm<T> s(a_matrix);
  s.subtract(I);
  sqm<double> a = abs(s,a_func);
  double radius = 1./order;

  std::vector<double>& vec = a.vector();
  typedef typename std::vector<T>::iterator vec_it_t;
  for(vec_it_t it = vec.begin();it!=vec.end();++it) {
    if((*it)>=radius) {
      //::printf("debug : not convergent at : %d : %g rad %g\n",
        //index,a.value(index),radius);
      return false;
    }
  }

  return true;
}
#endif

template <class T>
inline bool determinant(const sqm<T>& a_matrix,T& a_deter) {
  // Brut force determinant by using the kroneckers.
  a_deter = a_matrix.zero();

  unsigned int order = a_matrix.order(); 

  typedef typename std::vector<unsigned int> uints_t;
  typedef typename std::vector<T>::const_iterator const_vec_it_t;

  kronecker<T> k_c(order);
  const std::vector<T>& vec_c = k_c.vector();
  const_vec_it_t it_c;
  uints_t is_c(order);

  //kronecker<T> k_r(order);
  //const std::vector<T>& vec_r = k_r.vector();
  //const_vec_it_t it_r;
  //uints_t is_r(order);

  unsigned int index_c = 0;
  for(it_c = vec_c.begin();it_c!=vec_c.end();++it_c,index_c++) {
    if(*it_c==0.) continue;
    if(!k_c.indices(index_c,is_c)) return false;
    //unsigned int index_r = 0;
    //for(it_r = vec_r.begin();it_r!=vec_r.end();++it_r,index_r++) {
    //  if(*it_r==0.) continue;
    //  if(!k_r.indices(index_r,is_r)) return false;
      T v = array<T>::one();
      for(unsigned iorder=0;iorder<order;iorder++) {
        //v *= a_matrix.value(is_c[iorder],is_r[iorder]); 
        v *= a_matrix.value(is_c[iorder],iorder); 
      }
      a_deter += v;
    //}
  }

  //{T fact = array<T>::one();
  //for(unsigned index=1;index<=order;index++) fact *= index;
  //a_deter /= fact;}

  return true;
}


//
//  En_k with k in [0,n(n-1)/2[
//
//  En_0 :        En_1 :         En_2 :            En_3
//   0  1  0...     0  0 -1...     0  0  0  0...     0  0  0  1... 
//  -1  0  0...     0  0  0...     0  0  1  0...     0  0  0  0...
//   0  0  0...     1  0  0...     0 -1  0  0...     0  0  0  0...
//   0  0  0...     0  0  0...     0  0  0  0...    -1  0  0  0...
//   ..........     ..........     .............     .............
//
//  with r<c r,c in [0,n[ :
//    k(r,c) = c*(c-1)/2 + r
//  Reverse, given k :
//    c is the first such that : c*(c-1)/2 <= k <= c*(c+1)/2 -1
//  then :
//    r = k - c*(c-1)/2
//

// T must accept -1,0,1
template <class T>
class E : public sqm<T> {
public:
  static bool group_constants(unsigned int a_order,inlib::array<T>& a_constants) {
    unsigned int number = a_order * (a_order-1)/2;
  
    std::vector<unsigned int> orders(3,number);
    a_constants.configure(orders);
  
    std::vector< sqm<T> > Es;
    for(unsigned int index=0;index<number;index++)
      Es.push_back(E(a_order,index));

    sqm<T> zero(a_order);
  
    for(unsigned int j=0;j<number;j++) {
      for(unsigned int k=j+1;k<number;k++) {
        sqm<T> cm = commutator(Es[j],Es[k]);
  
        // Find the result in the Es themselves :
        bool found = false;
        for(unsigned int l=0;l<number;l++) {
          if(cm.equal(zero)) {
            found = true; 
            break;
          } else if(cm.equal(Es[l])) {
            std::vector<unsigned int> is(3);
            is[0] = j;
            is[1] = k;
            is[2] = l;
            a_constants.set_value(is,1);
            is[0] = k;
            is[1] = j;
            is[2] = l;
            a_constants.set_value(is,-1);
  
            found = true; 
            break;
          } else {
            sqm<T> tmp(cm);
            tmp.multiply(-1.);
            if(tmp.equal(Es[l])) {
              std::vector<unsigned int> is(3);
              is[0] = j;
              is[1] = k;
              is[2] = l;
              a_constants.set_value(is,-1);
              is[0] = k;
              is[1] = j;
              is[2] = l;
              a_constants.set_value(is,1);
  
              found = true; 
              break;
            }
          }
        }
        if(!found) {
          a_constants.clear();
          return false;
        }
      }
    }
    return true;
  }
  static bool metric(unsigned int a_order,sqm<T>& a_metric) {
    // For order n : -2*I(n*(n-1)/2) for all orders.
    unsigned int number = a_order * (a_order-1)/2;
  
    std::vector< sqm<T> > Es;
    for(unsigned int index=0;index<number;index++)
      Es.push_back(E<T>(a_order,index));
  
    a_metric.set_order(number);
  
    for(unsigned int j=0;j<number;j++) {
      for(unsigned int k=0;k<number;k++) {
        sqm<T> mx = Es[j]*Es[k];
        if(!a_metric.set_value(j,k,mx.trace())) return false;
      }
    }
  
    return true;
  }
  static bool cartan_metric(unsigned int a_order,sqm<T>& a_metric) {
    // For order n : -(2*n-4))*I(n*(n-1)/2)
    unsigned int number = a_order * (a_order-1)/2;
  
    array<T> csts;
    if(!E<T>::group_constants(a_order,csts)) {
      a_metric.clear();
      return false;
    }
  
    // Fill the adjoint representation :
    //      r        r
    //  (Xk)  = cst
    //      c      kc
    std::vector< sqm<T> > xs;
    for(unsigned int index=0;index<number;index++) {
      sqm<T> x(number);
      for(unsigned int j=0;j<number;j++) {
        for(unsigned int k=0;k<number;k++) {
          std::vector<unsigned int> is(3);
          is[0] = index;
          is[1] = k;
          is[2] = j;
          T value;
          if(!csts.value(is,value)) return false;
          if(!x.set_value(j,k,value)) return false;
        }
      }
      xs.push_back(x);
    }
  
    a_metric.set_order(number);
  
    for(unsigned int j=0;j<number;j++) {
      for(unsigned int k=0;k<number;k++) {
        sqm<T> mx = xs[j]*xs[k];
        if(!a_metric.set_value(j,k,mx.trace())) return false;
      }
    }
  
    return true;
  }
  
  static bool casimir(unsigned int a_order,sqm<T>& a_casimir) {
    //  For order n : (n-1)/2 * I(n)
    unsigned int number = a_order * (a_order-1)/2;
  
    sqm<T> metric(number);
    if(!E<T>::metric(a_order,metric)) {
      a_casimir.clear();
      return false;
    }
  
    sqm<T> metric_inv(number);
    if(!metric.invert(metric_inv)) return false;
  
    a_casimir.set_order(a_order);
  
    std::vector< sqm<T> > Es;
    for(unsigned int index=0;index<number;index++)
      Es.push_back(E<T>(a_order,index));
  
    for(unsigned int j=0;j<number;j++) {
      for(unsigned int k=0;k<number;k++) {
        sqm<T> mx = Es[j]*Es[k];
        mx.multiply(metric_inv.value(j,k));
        a_casimir.add(mx);
      }
    }
  
    return true;
  }
public:
  E(unsigned int a_order,unsigned int a_index): sqm<double>(a_order){
    if(!initialize(a_order,a_index)) return; //FIXME : throw
  }
  virtual ~E(){}
public:
  E(const E& a_from): sqm<T>(a_from){}
  E& operator=(const E& a_from) {
    inlib::sqm<T>::operator=(a_from);
    return *this;
  }
private:
  bool initialize(unsigned int a_order,unsigned int a_index) {
    unsigned int number = a_order * (a_order-1)/2;
    if(a_index>number) return false;

    // Find column from index :
    bool found = false;
    unsigned int c = 0;
    unsigned int mn = 0;
    unsigned int mx = 0;
    for(c=1;c<a_order;c++) {
      mn = c*(c-1)/2;
      mx = c*(c+1)/2-1;
      if((mn<=a_index)&&(a_index<=mx)) {
        found = true;
        break;
      }
    }
    if(!found) return false;
    // ok, then compute row :
    unsigned int r = a_index - mn;
    //printf("debug : %d : r %d c %d : mn %d mx %d\n",a_index,r,c,mn,mx);
    // then set value :
    T value = power(-1,r+c+1);
    // If one, then rot 3D vec csts are no more epsilon(i,j,k).
    //T value = 1;
    set_value_no_check(r,c,value);
    set_value_no_check(c,r,-value);
  
    unsigned int icheck = c*(c-1)/2+r;
    if(a_index!=icheck) return false;
    return true;
  }
  static T power(const T& a_A,unsigned int a_B){
    T v = 1;
    for(unsigned int i=0;i<a_B;i++) v *= a_A; 
    return v;
  }
};

// Poincare group in "a_order" dimension with delta metric (identity).
//
// From the operator representation :
//   OpRep_k = 1/2 * J_mu_nu * En_k_mu_nu   with k in [0,n*(n-1)/2[
//   OpRep_alpha = D_alpha                  with alpha in [0,n[
// With :
//   J_mu_nu =   delta_mu_alpha * coord_alpha * D_nu 
//             - delta_nu_alpha * coord_alpha * D_mu
// Not null group constants are :
//   c(j,k,l) = e_n(j,k,l)
//   c(j,alpha,beta) = En_j_beta_alpha
// T must accept -1,0,1
template <class T>
class poincare_cartan_metric : public sqm<T> {
public:
  poincare_cartan_metric(unsigned int a_order)
  :sqm<T>(a_order*(a_order+1)/2){
    unsigned int number = a_order*(a_order+1)/2;

    // Fill the adjoint representation :
    //      c        c
    //  (Xa)  = cst
    //      b      ab
    // a,b,c in [0,n*(n+1)/2[
  
    std::vector< sqm<T> > xs;
  
   {array<T> csts;
    if(!E<T>::group_constants(a_order,csts)) return;
    unsigned int n = a_order*(a_order-1)/2;
    for(unsigned int index=0;index<n;index++) {
      sqm<T> x(number);
  
     {for(unsigned int j=0;j<n;j++) {
        for(unsigned int k=0;k<n;k++) {
          std::vector<unsigned int> is(3);
          is[0] = index;
          is[1] = k;
          is[2] = j;
          T value;
          if(!csts.value(is,value)) return;
          if(!x.set_value(j,k,value)) return;
        }
      }}
  
     {E<T> E(a_order,index);
      for(unsigned int j=0;j<a_order;j++) {
        for(unsigned int k=0;k<a_order;k++) {
  
          T value = E.value(j,k);
  
          if(!x.set_value(n+j,n+k,value)) return;
        }
      }}
  
      xs.push_back(x);
    }}
  
   {for(unsigned int index=0;index<a_order;index++) {
      sqm<T> x(number);
      xs.push_back(x);
    }}
  
    // Set cartan killing metric :
   {for(unsigned int j=0;j<number;j++) {
      for(unsigned int k=0;k<number;k++) {
        inlib::sqm<T> mx = xs[j]*xs[k];
        if(!set_value(j,k,mx.trace())) return;
      }
    }}
  
  }
};

}

#include <ostream>

namespace inlib {

//NOTE : print is a Python keyword.
template <class T>
inline void dump(const sqm<T>& a_matrix,std::ostream& a_out,const std::string& aCMT) {
  if(aCMT.size()) a_out << aCMT << std::endl;
  unsigned int ord = a_matrix.order();
  for(unsigned int r=0;r<ord;r++) {
    for(unsigned int c=0;c<ord;c++) {
      a_out << " " << a_matrix.value(r,c);
    }
    a_out << std::endl;
  }
}

template <class T>
inline bool check_invert(const sqm<T>& a_matrix,std::ostream& a_out) {
  unsigned int ord = a_matrix.order();
  sqm<T> I(ord);
  I.set_identity();
  sqm<T> inv(ord);
  if(!a_matrix.invert(inv)) return false;
  //sqm<T> tmp(ord);
  //a_matrix.mx_mul(inv,tmp);
  sqm<T> tmp = a_matrix * inv;
  if(!tmp.equal(I)) {
    dump(a_matrix,a_out,"problem with inv of :");
    dump(inv,a_out,"found :");
    return false;
  }
  return true;
}

}

#endif