AnimaL

Tutorial

Documentation

---

Documentation

• Overview
• Modules
• Namespace List
• Class Hierarchy
• Alphabetical List
• Compound List
• File List
• Namespace Members
• Compound Members
• File Members
• Examples

---

linear.h

Go to the documentation of this file.

#ifndef ____________animal_linear_h________________
#define ____________animal_linear_h________________

#include <animal/container_traits.h>
#include <assert.h>
#include <math.h>


namespace animal {


//=============================================================  

template<class Cont, class Value> inline
void v_assign( Cont& v, const Value& a )
{
    typename container_traits<Cont>::iterator Yi = begin(v), Yend=end(v);
    for( ; Yi!=Yend; ++Yi )
        *(Yi) = a;
}

template<class Cont, class Value> inline
void v_addall( Cont& v, const Value& a )
{
    typename container_traits<Cont>::iterator Yi = begin(v), Yend=end(v);
    for( ; Yi!=Yend; ++Yi )
        *(Yi) += a;
}

template<class Cont,class Real> inline
void v_teq ( Cont& v, const Real& a )
{
    typename container_traits<Cont>::iterator Yi = begin(v), Yend=end(v);
    for( ; Yi!=Yend; ++Yi )
        *(Yi) *= a;
}

template<class Cont> inline
//typename container_traits<Cont>::value_type v_norm( const Cont& v )
typename value_type<Cont>::type v_norm( const Cont& v )
{
    //cout << "v_norm, v_dot(v,v)= "<<v_dot(v,v)<<endl;
    return sqrt(v_dot(v,v));
}

template<class Cont> inline
void v_normalize( Cont& v )
{
    v_teq(v,1/v_norm(v));
}

template<class Cont> inline
typename container_traits<Cont>::value_type v_norm1( const Cont& v )
{
    typename container_traits<Cont const>::iterator xi = begin(v), xend=end(v);
    typename container_traits<Cont>::value_type d = 0;
    for( ; xi!=xend; ++xi )
        d += fabs(*xi);
    return d;
}

template<class Cont> inline
typename container_traits<Cont>::value_type v_normInf( const Cont& v )
{
    typename container_traits<Cont const>::iterator xi = begin(v), xend=end(v);
    typename container_traits<Cont>::value_type d = 0;
    for( ; xi!=xend; ++xi )
        if( fabs(*xi) > d ) d = fabs(*xi);
    return d;
}


template<class Cont> inline
typename container_traits<Cont>::value_type v_sum( const Cont& v )
{
    typename container_traits<Cont const>::iterator xi = begin(v), xend=end(v);
    typename container_traits<Cont>::value_type d = 0;
    for( ; xi!=xend; ++xi )
        d += (*xi);
    return d;
}

template<class Cont> inline
typename container_traits<Cont>::value_type v_sqsum( const Cont& v )
{
    typename container_traits<Cont const>::iterator xi = begin(v), xend=end(v);
    typename container_traits<Cont>::value_type d = 0;
    for( ; xi!=xend; ++xi )
        d += (*xi)*(*xi);
    return d;
}

template<class Cont> inline
typename container_traits<Cont>::value_type v_mean( const Cont& v )
{
    return v_sum(v)/size(v);
}

template<class Cont> inline
typename container_traits<Cont>::value_type v_variance( const Cont& v )
{
    typename container_traits<Cont const>::iterator xi = begin(v), xend=end(v);
    typename container_traits<Cont>::value_type m = v_mean(v);
    typename container_traits<Cont>::value_type d = 0;
    for( ; xi!=xend; ++xi )
    {
        typename container_traits<Cont>::value_type tmp = (*xi)-m;
        tmp*=tmp;
        d+=tmp;
    }
    return d/size(v);
}

//=============================================================  

template<class Cont1, class Cont2> inline
void v_eq ( Cont1& Y, const Cont2& X )
{
  assert( size(Y)==size(X) );
    typename container_traits<Cont1>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont2 const>::iterator Xi = begin(X);
    for( ; Yi!=Yend; ++Yi, ++Xi )
        *(Yi) = (*Xi);
}

template<class Cont1, class Cont2> inline
void v_eq_minus ( Cont1& Y, const Cont2& X )
{
    assert( size(Y)==size(X) );
    typename container_traits<Cont1>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont2 const>::iterator Xi = begin(X);
    for( ; Yi!=Yend; ++Yi, ++Xi )
        *(Yi) = -(*Xi);
}

template<class Cont1, class Cont2, class Real> inline
void v_eq_ab ( Cont1& Y, const Real& a, const Cont2& X )
{
    assert( size(Y)==size(X) );
    typename container_traits<Cont1>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont2 const>::iterator Xi = begin(X);
    for( ; Yi!=Yend; ++Yi, ++Xi )
        *(Yi) = a*(*Xi);//*(Yi) = (*Xi)*a;
}

template<class Cont1, class Cont2, class Cont3> inline
void v_eq_apb ( Cont1& Y, const Cont2& W, const Cont3& X )
{
    assert( size(Y)==size(X) );
    assert( size(Y)==size(W) );
    typename container_traits<Cont1>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont2 const>::iterator Wi=begin(W);
    typename container_traits<Cont3 const>::iterator Xi = begin(X);
    for( ; Yi!=Yend; ++Yi, ++Xi, ++Wi )
        *(Yi) = (*Wi)+(*Xi);
}

template<class Cont1, class Cont2, class Cont3> inline
void v_eq_amb ( Cont1& Y, const Cont2& W, const Cont3& X )
{
    assert( size(Y)==size(X) );
    assert( size(Y)==size(W) );
    typename container_traits<Cont1>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont2 const>::iterator Wi=begin(W);
    typename container_traits<Cont3 const>::iterator Xi = begin(X);
    for( ; Yi!=Yend; ++Yi, ++Xi, ++Wi )
        *(Yi) = (*Wi)-(*Xi);
}

template<class Cont1, class Cont2, class Cont3, class Entry> inline
void v_eq ( Cont1& Y, const Cont2& W, const Entry& a, const Cont3& X )
{
    assert( size(Y)==size(X) );
    assert( size(Y)==size(W) );
    typename container_traits<Cont1>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont3 const>::iterator Xi = begin(X);
    typename container_traits<Cont2 const>::iterator Wi=begin(W);
    for( ; Yi!=Yend; ++Yi, ++Xi, ++Wi )
        *(Yi) = (*Wi)+a*(*Xi);
}

template<class Cont1, class Cont2> inline
void v_peq ( Cont1& Y, const Cont2& X )
{
    assert( size(Y)==size(X) );
    typename container_traits<Cont1>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont2 const>::iterator Xi = begin(X);
    for( ; Yi!=Yend; ++Yi, ++Xi )
        *(Yi) += (*Xi);
}

template<class Cont> inline
void v_meq ( Cont& Y, const Cont& X )
{
    assert( size(Y)==size(X) );
    typename container_traits<Cont>::iterator Yi = begin(Y), Yend=end(Y);
    typename container_traits<Cont const>::iterator Xi = begin(X);
    for( ; Yi!=Yend; ++Yi, ++Xi )
        *(Yi) -= (*Xi);
}

template<class Cont> inline
void v_peq ( Cont& y, typename container_traits<Cont>::value_type a, const Cont& x )
{
    assert( size(y)==size(x) );
    typename container_traits<Cont>::iterator yi = begin(y), yend=end(y);
    typename container_traits<Cont const>::iterator xi = begin(x);
    for( ; yi!=yend; ++xi, ++yi )
        *(yi) += a * (*xi);
}

template<class Cont> inline
void v_meq ( Cont& y, typename container_traits<Cont>::value_type a, const Cont& x )
{
    assert( size(y)==size(x) );
    typename container_traits<Cont>::iterator yi = begin(y), yend=end(y);
    typename container_traits<Cont const>::iterator xi = begin(x);
    for( ; yi!=yend; ++xi, ++yi )
        *(yi) -= a * (*xi);
}

template<class Cont,class Real> inline
void v_peq_ab ( Cont& y, Real a, const Cont& x )
{
    assert( size(y)==size(x) );
    typename container_traits<Cont>::iterator yi = begin(y), yend=end(y);
    typename container_traits<Cont const>::iterator xi = begin(x);
    for( ; yi!=yend; ++xi, ++yi )
        *(yi) += a * (*xi);//*(yi) += (*xi) * a;
}

template<class V1, class V2, class V3> inline
void v_eq_cross ( V1&u, const V2& v, const V3& w )
{
    assert( static_cast<const void*>(&u) != static_cast<const void*>(&v) );
    assert( static_cast<const void*>(&u) != static_cast<const void*>(&w) );
    assert( size(u)==size(v) );
    assert( size(u)==size(w) );
    const int n = size(u);
    for( int i=0, iend=n; i!=iend; ++i )
        u[(i+2)%n] = v[i]*w[(i+1)%n] -v[(i+1)%n]*w[i];
}

template<class V1, class V2, class V3> inline
void v_peq_cross ( V1&u, const V2& v, const V3& w )
{
    assert( static_cast<const void*>(&u) != static_cast<const void*>(&v) );
    assert( static_cast<const void*>(&u) != static_cast<const void*>(&w) );
    assert( size(u)==size(v) );
    assert( size(u)==size(w) );
    const int n = size(u);
    for( int i=0, iend=n; i!=iend; ++i )
        u[(i+2)%n] += v[i]*w[(i+1)%n] -v[(i+1)%n]*w[i];
}

template<class V1, class V2, class V3> inline
void v_meq_cross ( V1&u, const V2& v, const V3& w )
{
    assert( static_cast<const void*>(&u) != static_cast<const void*>(&v) );
    assert( static_cast<const void*>(&u) != static_cast<const void*>(&w) );
    assert( size(u)==size(v) );
    assert( size(u)==size(w) );
    const int n = size(u);
    for( int i=0, iend=n; i!=iend; ++i )
        u[(i+2)%n] -= v[i]*w[(i+1)%n] -v[(i+1)%n]*w[i];
}

template<class Cont> inline
//typename container_traits<Cont>::value_type v_dot( const Cont& x, const Cont& y )
typename value_type<Cont>::type v_dot( const Cont& x, const Cont& y )
{
    assert( size(x)==size(x) );
    typename container_traits<Cont const>::iterator yi = begin(y), yend=end(y);
    typename container_traits<Cont const>::iterator xi = begin(x);
    typename container_traits<Cont>::value_type d = 0;
    for( ; yi!=yend; ++xi, ++yi )
        d += *(yi) * (*xi);
    return d;
}

// vector-vector operations
//======================================================================

template<class M> inline
std::size_t nrows( const M& m )
{
    return container_traits<M const>::size(m);
}

template<class M> inline
std::size_t ncols( const M& m )
{
    return container_traits<typename container_traits<M const>::value_type const>::size(m[0]);
}


template<class M>
struct matrix_traits
{
    typedef typename container_traits<typename container_traits<M>::value_type>::value_type value_type;
};


template<class Cont, class Value> inline
void m_assign( Cont& M, const Value& a )
{
//  container_traits<Cont>::iterator Mi = begin(M), Mend=end(M);
//  for( ; Mi!=Mend; ++Mi )
//      v_assign(*Mi,a);
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] =a;
}


template<class Cont, class Val> inline
//void m_teq ( Cont& M, const typename matrix_traits<Cont>::value_type& a )
void m_teq ( Cont& M, const Val& a )
{
//  container_traits<Cont>::iterator Mi = begin(M), Mend=end(M);
//  for( ; Mi!=Mend; ++Mi )
//      v_teq(*Mi,a);
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] *=a;
}


template<class Cont, class Val> inline
void m_deq ( Cont& M, const Val& a )
{
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] /=a;
}





template<class Cont1, class Cont2> inline
void m_eq_A ( Cont1& M, const Cont2& A )
{
    assert( nrows(M)==nrows(A) );
    assert( ncols(M)==ncols(A) );
    //std::cout<<"r,c: "<<nrows(M)<<","<<ncols(M)<<std::endl;
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] = A[i][j];
    //std::cout<<"m_eq_A OK"<<std::endl;
}

template<class Cont1, class Cont2> inline
//void m_eq_aA ( Cont1& M, typename container_traits<typename container_traits<Cont1>::value_type>::value_type a, const Cont2& A )
void m_eq_aA ( Cont1& M, typename matrix_traits<Cont2>::value_type a, const Cont2& A )
{
    assert( nrows(M)==nrows(A) );
    assert( ncols(M)==ncols(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] = a * A[i][j];
}

template<class Cont1, class Cont2> inline
void m_peq_A ( Cont1& M, const Cont2& A )
{
    assert( nrows(M)==nrows(A) );
    assert( ncols(M)==ncols(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] += A[i][j];
}

template<class Cont1, class Cont2> inline
void m_peq_aA ( Cont1& M, typename matrix_traits<Cont2>::value_type a, const Cont2& A )
{
    assert( nrows(M)==nrows(A) );
    assert( ncols(M)==ncols(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] += a * A[i][j];
}

template<class Cont1, class Cont2> inline
void m_meq_A ( Cont1& M, const Cont2& A )
{
    assert( nrows(M)==nrows(A) );
    assert( ncols(M)==ncols(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] -= A[i][j];
}






/*template<class Cont> inline
typename matrix_traits<Cont>::value_type m_det_A ( const Cont& A, int dim = -1 )
{   



    Cont sous_m;
    typename matrix_traits<Cont>::value_type det = 0;
    if( dim == -1 ) dim = nrows(A);
    if( dim == 1 )  return ( A[0][0] );
    for(int l=0,signe=1 ; l<dim ; ++l)
    {
        for( int i=0, ld=0 ; i<dim ; i++ )
        {
            if( i!=l ) // delete the current line
            {
                for( int c=1 ; c<dim ; c++ ) // begin at 1 to delete the first line
                    sous_m[ld][c-1] = A[i][c];
                ld++;
            }
        det += signe*A[l][0] * m_det_A(sous_m,dim-1);
            signe=-signe;
        }
   }
    return(det);
}*/






template<class Cont1, class Cont2> inline
void m_eq_At ( Cont1& M, const Cont2& A )
{
    assert( nrows(M)==ncols(A) );
    assert( ncols(M)==nrows(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] = A[j][i];
}

template<class Cont1, class Cont2> inline
void m_eq_aAt ( Cont1& M, typename matrix_traits<Cont2>::value_type a, const Cont2& A )
{
    assert( nrows(M)==ncols(A) );
    assert( ncols(M)==nrows(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] = a * A[j][i];
}

template<class Cont1, class Cont2> inline
void m_peq_At ( Cont1& M, const Cont2& A )
{
    assert( nrows(M)==ncols(A) );
    assert( ncols(M)==nrows(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] += A[j][i];
}

template<class Cont1, class Cont2> inline
void m_peq_aAt ( Cont1& M, const typename matrix_traits<Cont2>::value_type a, const Cont2& A )
{
    assert( nrows(M)==ncols(A) );
    assert( ncols(M)==nrows(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] += a * A[j][i];
}

template<class Cont1, class Cont2> inline
void m_meq_At ( Cont1& M, const Cont2& A )
{
    assert( nrows(M)==ncols(A) );
    assert( ncols(M)==nrows(A) );
    for( int i=0, iend=nrows(M); i!=iend; ++i )
        for( int j=0, jend=ncols(M); j!=jend; ++j )
            M[i][j] -= A[j][i];
}




//=============================================================  


template<class M1,class M2, class M3> inline
void m_eq_AB ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j ){
            C[i][j] = 0;
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += A[i][k] * B[k][j];
        }
}

template<class M1,class M2, class M3> inline
void m_eq_aAB ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j ){
            C[i][j] = 0;
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += a * A[i][k] * B[k][j];
        }
}

template<class M1,class M2, class M3> inline
void m_peq_AB ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j )
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += A[i][k] * B[k][j];
}

template<class M1,class M2, class M3> inline
void m_peq_aAB ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j )
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += a * A[i][k] * B[k][j];
}

template<class M1,class M2, class M3> inline
void m_meq_AB ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j )
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] -= A[i][k] * B[k][j];
}




template<class M1,class M2, class M3> inline
void m_eq_ABt ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j )
        {
            C[i][j] = 0;
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += A[i][k] * B[j][k];
        }
}

template<class M1,class M2, class M3> inline
void m_eq_aABt ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j )
        {
            C[i][j] = 0;
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += a * A[i][k] * B[j][k];
        }
}

template<class M1,class M2, class M3> inline
void m_peq_ABt ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j )
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += A[i][k] * B[j][k];
}

template<class M1,class M2, class M3> inline
void m_peq_aABt ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j )
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] += a * A[i][k] * B[j][k];
}

template<class M1,class M2, class M3> inline
void m_meq_ABt ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( nrows(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==ncols(A) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j )
            for( Int k=0, p=ncols(A); k<p; ++k )
                C[i][j] -= A[i][k] * B[j][k];
}



template<class M1,class M2, class M3> inline
void m_eq_AtB ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j ){
            C[i][j] = 0;
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += A[k][i] * B[k][j];
        }
}

template<class M1,class M2, class M3> inline
void m_eq_aAtB ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j ){
            C[i][j] = 0;
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += a * A[k][i] * B[k][j];
        }
}

template<class M1,class M2, class M3> inline
void m_peq_AtB ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j ){
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += A[k][i] * B[k][j];
        }
}

template<class M1,class M2, class M3> inline
void m_peq_aAtB ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j ){
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += a * A[k][i] * B[k][j];
        }
}

template<class M1,class M2, class M3> inline
void m_meq_AtB ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( ncols(B)==ncols(C) );
    assert( nrows(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=ncols(B); j<m; ++j ){
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] -= A[k][i] * B[k][j];
        }
}



template<class M1,class M2, class M3> inline
void m_eq_AtBt ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j ){
            C[i][j] = 0;
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += A[k][i] * B[j][k];
        }
}

template<class M1,class M2, class M3> inline
void m_eq_aAtBt ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j ){
            C[i][j] = 0;
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += a * A[k][i] * B[j][k];
        }
}

template<class M1,class M2, class M3> inline
void m_peq_AtBt ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j ){
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += A[k][i] * B[j][k];
        }
}

template<class M1,class M2, class M3> inline
void m_peq_aAtBt ( M1& C, const typename matrix_traits<M2>::value_type a, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j ){
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] += a * A[k][i] * B[j][k];
        }
}

template<class M1,class M2, class M3> inline
void m_meq_AtBt ( M1& C, const M2& A, const M3& B )
{
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&A) );
    assert( static_cast<const void*>(&C) != static_cast<const void*>(&B) );
    assert( ncols(A)==nrows(C) );
    assert( nrows(B)==ncols(C) );
    assert( ncols(B)==nrows(A) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i )
        for( Int j=0, m=nrows(B); j<m; ++j ){
            for( Int k=0, p=nrows(A); k<p; ++k )
                C[i][j] -= A[k][i] * B[j][k];
        }
}

//=============================================================  






template<class V1, class M, class V2> inline
void v_eq_Ab ( V1& v, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( ncols(A)==size(b) );
    assert( nrows(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i ){
        v[i] = 0;
        for( Int k=0, p=ncols(A); k<p; ++k )
            v[i] += A[i][k] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_eq_aAb ( V1& v, const typename matrix_traits<M>::value_type a, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( ncols(A)==size(b) );
    assert( nrows(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i ){
        v[i] = 0;
        for( Int k=0, p=ncols(A); k<p; ++k )
            v[i] += a * A[i][k] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_peq_aAb ( V1& v, const typename matrix_traits<M>::value_type a, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( ncols(A)==size(b) );
    assert( nrows(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i ){
        for( Int k=0, p=ncols(A); k<p; ++k )
            v[i] += a * A[i][k] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_peq_Ab ( V1& v, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( ncols(A)==size(b) );
    assert( nrows(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i ){
        for( Int k=0, p=ncols(A); k<p; ++k )
            v[i] += A[i][k] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_meq_Ab ( V1& v, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( ncols(A)==size(b) );
    assert( nrows(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=nrows(A); i<n; ++i ){
        for( Int k=0, p=ncols(A); k<p; ++k )
            v[i] -= A[i][k] * b[k];
    }
}




template<class V1, class M, class V2> inline
void v_eq_Atb ( V1& v, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( nrows(A)==size(b) );
    assert( ncols(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i ){
        v[i] = 0;
        for( Int k=0, p=nrows(A); k<p; ++k )
            v[i] += A[k][i] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_eq_aAtb ( V1& v, const typename matrix_traits<M>::value_type a, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( nrows(A)==size(b) );
    assert( ncols(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i ){
        v[i] = 0;
        for( Int k=0, p=nrows(A); k<p; ++k )
            v[i] += a * A[k][i] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_peq_aAtb ( V1& v, const typename matrix_traits<M>::value_type a, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( nrows(A)==size(b) );
    assert( ncols(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i ){
        for( Int k=0, p=nrows(A); k<p; ++k )
            v[i] += a * A[k][i] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_peq_Atb ( V1& v, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( nrows(A)==size(b) );
    assert( ncols(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i ){
        for( Int k=0, p=nrows(A); k<p; ++k )
            v[i] += A[k][i] * b[k];
    }
}

template<class V1, class M, class V2> inline
void v_meq_Atb ( V1& v, const M& A, const V2& b )
{
    assert( static_cast<const void*>(&v) != static_cast<const void*>(&b) );
    assert( nrows(A)==size(b) );
    assert( ncols(A)==size(v) );
    typedef std::size_t Int;
    for( Int i=0, n=ncols(A); i<n; ++i ){
        for( Int k=0, p=nrows(A); k<p; ++k )
            v[i] -= A[k][i] * b[k];
    }
}

// matrix-vector operations
//========================================================================

// LU
template <class M1, class M2>
void m_eq_ludcmp( M1& LU, const M2& A )
{
    assert( static_cast<const void*>(&LU) != static_cast<const void*>(&A) );
    assert( nrows(A) == ncols(A) );
    assert( nrows(A) == ncols(LU) );
    assert( nrows(A) == nrows(LU) );
    
    int i,j,k,n = nrows(A);
    for( i=0; i<n; ++i )
    {   
        for (j=i; j<n; ++j)
        {
            LU[j][i] = A[j][i];                     
            for (k=0; k<i; ++k)
            {
                LU[j][i] -= LU[j][k]*LU[k][i];
            }
        }
        for (j=i+1; j<n; ++j)
        {
            LU[i][j] = A[i][j];                     
            for (k=0; k<i; ++k)
            {
                LU[i][j] -= LU[i][k]*LU[k][j];
            }
            LU[i][j] /= LU[i][i];
        }
    }
    
}






template <class M, class V1, class V2> inline
void v_eq_lusolve( V1& x, const M& LU, const V2& b )
{
    assert( static_cast<const void*>(&x) != static_cast<const void*>(&b) );
    assert( nrows(LU) == ncols(LU) );
    assert( size(x) == ncols(LU) );
    assert( size(b) == nrows(LU) );
    
    int i,j,n=nrows(LU);
    for( i=0; i<n; ++i )
    {
        x[i] = b[i];
        for( j=0; j<i; ++j )
        {
            x[i] -= LU[i][j]*x[j];
        }
        x[i] /= LU[i][i];
    }
    for( i=n-1; i>=0; --i )
    {
        for( j=i+1; j<n; ++j )
        {
            x[i] -= LU[i][j]*x[j];
        }
    }
}

template <class M, class M2, class Vec>
void m_eq_luinv( M& Inv, const M2& LU, Vec& aux )
{
    assert( static_cast<const void*>(&Inv) != static_cast<const void*>(&LU));
    assert( nrows(LU) == ncols(LU) );
    assert( nrows(LU) == nrows(Inv) );
    assert( nrows(LU) == ncols(Inv) );
    assert( nrows(LU) == size(aux) );

    int n=size(aux),i,j,k;
    for( k=0; k<n; ++k )
    {
        v_assign( aux, 0 ); aux[k] = 1;
        for( i=0; i<n; ++i )
        {
            Inv[i][k] = aux[i];
            for( j=0; j<i; ++j )
            {
                Inv[i][k] -= LU[i][j]*Inv[j][k];
            }
            Inv[i][k] /= LU[i][i];
        }
        for( i=n-1; i>=0; --i )
        {
            for( j=i+1; j<n; ++j )
            {
                Inv[i][k] -= LU[i][j]*Inv[j][k];
            }
        }
    }
}




template <class M, class M2>
int ludcmp(M &a, M2 indx )
{ 
double d;
int i, imax = 0, j, k;
  double big, dum, sum, temp;
  int n = nrows(a);
    typedef typename matrix_traits<M>::value_type Real;
  Real vv[n];
  d = 1.0;
  for (i = 0; i < n; i++)
  { big = 0.0;
    for (j = 0; j < n; j++)
      if ((temp = fabs(a[i][j])) > big)
        big = temp;
    if (big == 0.0)
      return -1;
    vv[i] = 1.0/big;
  }
  for (j = 0; j < n; j++)
  { for (i = 0; i < j; i++)
    { sum = a[i][j];
      for (k = 0; k < i; k++)
        sum -= a[i][k]*a[k][j];
      a[i][j] = sum;
    }
    big = 0.0;
    for (i = j; i < n; i++)
    { sum = a[i][j];
      for (k = 0; k < j; k++)
        sum -= a[i][k]*a[k][j];
      a[i][j] = sum;
      if ((dum = vv[i]*fabs(sum)) >= big)
      { big = dum;
        imax = i;
      }
    }
    if (j != imax)
    { for (k = 0; k < n; k++)
      { dum = a[imax][k];
        a[imax][k] = a[j][k];
        a[j][k] = dum;
      }
      d = -d;
      vv[imax] = vv[j];
    }
    indx[j] = imax;
    if (a[j][j] == 0.0)
      a[j][j] = 1.0e-20;
    if (j != n)
    { dum = 1.0/a[j][j];
      for (i = j+1; i < n; i++)
        a[i][j] *= dum;
    }
  }
  for (i = 0; i < n; i++)
  { for (j = 0; j < n; j++)
      if (fabs(a[i][j]) > 1.0e-15)
        break;
    for (k = n-1; k > j; k--)
      if (fabs(a[i][k]) > 1.0e-15)
        break;
  }
  return 0;
}

template <class M, class M2, class M3>
int lubksb(M &a, M2 indx, M3 b)
{ int i, j, k, ip, m, ii = 0;
  double sum;
  int n = nrows(a);
  for (i = 0; i < n; i++)
  { ip = indx[i];
    sum = b[ip];
    b[ip] = b[i];
    if (ii)
    { k = 0;
      if (ii > k)
        k = ii;
      m = n;
      if (i-1 < m)
        m = i-1;
      for (j = k; j < m; j++)
        sum -= a[i][j]*b[j];
    }
    else if (sum)
      ii = i;
    b[i] = sum;
  }
  for (i = n-1; i >= 0; i--)
  { sum = 0;
    k = 0;
    if (i+1 > k)
      k = i+1;
    m = n;
    if (n < m)
      m = n;
    for (j = k; j < m; j++)
      sum -= a[i][j]*b[j];
    b[i] = sum/a[i][i];
  }
  return 0;
}

template <class M>
int ns1__luinv(M &a, M &b)
{ 
    int n = nrows(a);
    typedef typename matrix_traits<M>::value_type Real;
    
    Real col[n];
   int indx[n];
  int i, j, k;
  if (ludcmp(a, indx))
    return 1;
  for (j = 0; j < n; j++)
  { for (i = 0; i < n; i++)
      col[i] = 0.0;
    col[j] = 1.0;
    lubksb(a, indx, col);
    for (i = 0; i < n; i++)
      b[i][j] = col[i];
  }
  for (i = 0; i < n; i++)
  { for (j = 0; j < n; j++)
      if (fabs(b[i][j]) > 1.0e-15)
        break;
    for (k = n-1; k > j; k--)
      if (fabs(b[i][k]) > 1.0e-15)
        break;
    //(*b)[i].resize(j, k);
  }
  return 0;
}






// cholesky

template <class M1, class M2> inline
void m_eq_modchol( M1& L, const M2& A )
{
    assert( static_cast<const void*>(&L) != static_cast<const void*>(&A) );
    assert( nrows(A) == ncols(A) );
    assert( nrows(A) == ncols(L) );
    assert( nrows(A) == nrows(L) );
    
    typedef int Int;
    //typedef container_traits<container_traits<M1>::value_type>::value_type Real;
    typedef typename matrix_traits<M1>::value_type Real;
    
    Int n = nrows(A);
    for( Int i=0; i<n; ++i )
    {
        Real sum = A[i][i];
        for( Int k=i-1; k>=0; --k )
            sum -= L[i][k]*L[i][k];
        if( sum <= 0 ) std::cerr<<"animal::m_eq_modchol nonpositive value " << sum << std::endl;
        L[i][i] = 1.0/sqrt(sum);

        for( Int j=i+1; j<n; ++j )
        {
            Real sum = A[i][j];
            for( Int k=i-1; k>=0; --k )
                sum -= L[i][k]*L[j][k];
            L[j][i] = sum * L[i][i];
        }
    }
}


template <class M, class V1, class V2> inline
void v_eq_cholsolve( V1& x, const M& L, const V2& b )
{
    assert( static_cast<const void*>(&x) != static_cast<const void*>(&b) );
    assert( nrows(L) == ncols(L) );
    assert( size(x) == ncols(L) );
    assert( size(b) == nrows(L) );

    typedef typename container_traits<V1>::value_type Real;
    
    int i,k,n = size(x);
    for( i=0; i<n; ++i )
    {
        Real sum = b[i];
        for( k=i-1; k>=0; --k )
            sum -= L[i][k]*x[k];
        x[i] = sum * L[i][i];
    }
    for( i=n-1; i>=0; --i )
    {
        Real sum = x[i];
        for( k=i+1; k<n; ++k  )
            sum -= L[k][i]*x[k];
        x[i] = sum * L[i][i];
    }
}

template <class M, class M2, class M3> inline
void m_eq_cholinv( M& Inv, const M2& A, M3& Tmp )
{
    assert( static_cast<const void*>(&A) != static_cast<const void*>(&Tmp));
    assert( static_cast<const void*>(&Inv) != static_cast<const void*>(&Tmp));
    m_eq_invsqrt(Tmp,A);
    m_eq_AtB( Inv, Tmp, Tmp );  
}


template <class M1, class M2> inline
void m_eq_invsqrt(M1& B , const M2& A)
{
    assert( static_cast<const void*>(&B) != static_cast<const void*>(&A) );
    assert( nrows(A) == ncols(A) );
    assert( nrows(B) == ncols(A) );
    assert( ncols(B) == ncols(A) );
    
    //typedef container_traits<container_traits<M1>::value_type>::value_type Real;
    typedef typename  matrix_traits<M1>::value_type Real;

    m_eq_modchol(B,A);

    int n=nrows(A);
    for( int i=0; i<n; ++i )
        for( int j=i+1; j<n; ++j )
        {
            Real sum = 0.0;
            for( int k=i; k<j; ++k )
                sum -= B[j][k]*B[k][i];
            B[j][i] = sum * B[j][j];
            B[i][j] = 0.0;
        }   
}  





template <class M1, class M2> inline
void m_eq_cholfactor( M1& L, const M2& A )
{
    assert( static_cast<const void*>(&L) != static_cast<const void*>(&A) );
    assert( nrows(A) == ncols(A) );
    assert( nrows(A) == ncols(L) );
    assert( nrows(A) == nrows(L) );
    
    int n = nrows(A);
    typedef typename matrix_traits<M1>::value_type Real;
    
    for( int i = 0 ; i<n ; ++i)
        for( int j = 0 ; j<n ; ++j)
            L[i][j]=0;  
  
    for (int x = 0;x<n;x++)
    {
        // determine diagonal element
    Real sum = A[x][x];
    for (int y = 0;y<x;y++)
    {
        Real value = L[x][y];
        sum -= value*value;
    }
    L[x][x] = sqrt(sum);

    // determine off-diagonal elements
    for (int y = x+1;y<n;y++)
    {
        sum = A[y][x];
        for (int i = 0;i < x;i++)
            sum -= L[y][i]*L[x][i];
        L[y][x] = sum/L[x][x];
        }
    }
}






template <class M1, class M2> inline
void m_eq_invsqrt2(M1& B , const M2& A)
{
    assert( static_cast<const void*>(&B) != static_cast<const void*>(&A) );
    assert( nrows(A) == ncols(A) );
    assert( nrows(B) == ncols(A) );
    assert( ncols(B) == ncols(A) );
    
    //typedef container_traits<container_traits<M1>::value_type>::value_type Real;
    typedef typename  matrix_traits<M1>::value_type Real;

    m_eq_cholfactor(B,A);

    int n=nrows(A);
    for( int i=0; i<n; ++i )
    {
        B[i][i] = 1.0/A[i][i];
        for( int j=i+1; j<n; ++j )
        {
            Real sum = 0.0;
            for( int k=i; k<j-1; ++k )
                sum -= A[j][k]*A[k][i];
                
            B[j][i] = sum / A[j][j];
            B[i][j] = 0.0;
        }   
    }
} 




template <class M1, class M2> inline
void m_eq_inv(M1& B , const M2& A)
{
    assert( static_cast<const void*>(&B) != static_cast<const void*>(&A) );
    assert( nrows(A) == 3 );
    assert( nrows(A) == 3 );
    assert( nrows(A) == ncols(A) );
    assert( nrows(B) == ncols(A) );
    assert( ncols(B) == ncols(A) );
    
    typedef typename  matrix_traits<M1>::value_type Real;
    
    
    float detval = (A[0][0]*A[1][1]*A[2][2] + A[1][0]*A[2][1]*A[0][2] + A[0][1]*A[1][2]*A[2][0]) - (A[2][0]*A[1][1]*A[0][2]+A[1][0]*A[0][1]*A[2][2]+A[2][1]*A[1][2]*A[0][0]);  // determiant

    if (detval==0)  cerr<<"Cannot inverse a Singular Matrix."<<endl;
    else
    {
    
        /* transpose of the cofactors gives adjoint */
        B[0][0] = (A[1][1]*A[2][2] - A[1][2]*A[2][1]);
        B[0][1] = -1 * (A[0][1]*A[2][2] - A[0][2]*A[2][1]);
        B[0][2] = (A[0][1]*A[1][2] - A[0][2]*A[1][1]);
    
        B[1][0] = -1 * (A[1][0]*A[2][2] - A[1][2]*A[2][0]);
        B[1][1] = (A[0][0]*A[2][2] - A[0][2]*A[2][0]);
        B[1][2] = -1 * (A[0][0]*A[1][2] - A[1][0]*A[0][2]);
    
        B[2][0] = (A[1][0]*A[2][1] - A[1][1]*A[2][0]);
        B[2][1] = -1 * (A[0][0]*A[2][1] -A[0][1]*A[2][0]);
        B[2][2] = (A[0][0]*A[1][1] - A[0][1]*A[1][0]);
    

        B[0][0] = B[0][0] / detval;
        B[0][1] = B[0][1] / detval;
        B[0][2] = B[0][2] / detval;
    
        B[1][0] = B[1][0] / detval;
        B[1][1] = B[1][1] / detval;
        B[1][2] = B[1][2] / detval;
    
        B[2][0] = B[2][0] / detval;
        B[2][1] = B[2][1] / detval;
        B[2][2] = B[2][2] / detval;
    }
}













// linear equation solution
//===========================================================================



}//animal




#endif

---