enable_if< O<=M &&P<=N >::type *matrix< M, N, T >::get_sub_matrix(matrix< O, P, T > &result, size_t row_offset, size_t col_offset) const {for(size_t row=0;row< O;++row){for(size_t col=0;col< P;++col){result.at(row, col)=at(row_offset+row, col_offset+col);}}return 0;}template< size_t M, size_t N, typename T >template< size_t O, size_t P >typename enable_if< O<=M &&P<=N >::type *matrix< M, N, T >::set_sub_matrix(const matrix< O, P, T > &sub_matrix, size_t row_offset, size_t col_offset){for(size_t row=0;row< O;++row){for(size_t col=0;col< P;++col){at(row_offset+row, col_offset+col)=sub_matrix.at(row, col);}}return 0;}template< size_t M, size_t N, typename T >inline Tmatrix< M, N, T >::det() const {return compute_determinant(*this);}template< size_t M, size_t N, typename T >template< size_t O, size_t P, typename TT >inline bool matrix< M, N, T >::inverse(matrix< O, P, TT > &inverse_, T tolerance, typename enable_if< M==N &&O==P &&O==M &&M >=2 &&M<=4, TT >::type *) const {return compute_inverse(*this, inverse_, tolerance);}template< size_t M, size_t N, typename T >template< size_t O, size_t P >typename enable_if< O==P &&M==N &&O==M &&M >=2 >::type *matrix< M, N, T >::get_adjugate(matrix< O, P, T > &adjugate) const {get_cofactors(adjugate);adjugate=transpose(adjugate);return 0;}template< size_t M, size_t N, typename T >template< size_t O, size_t P >typename enable_if< O==P &&M==N &&O==M &&M >=2 >::type *matrix< M, N, T >::get_cofactors(matrix< O, P, T > &cofactors) const {matrix< M-1, N-1, T > minor_;const size_t _negate=1u;for(size_t row_index=0;row_index< M;++row_index){for(size_t col_index=0;col_index< N;++col_index){if((row_index+col_index)&_negate) cofactors(row_index, col_index)=-get_minor(minor_, row_index, col_index);else cofactors(row_index, col_index)=get_minor(minor_, row_index, col_index);}}return 0;}template< size_t M, size_t N, typename T >template< size_t O, size_t P >Tmatrix< M, N, T >::get_minor(matrix< O, P, T > &minor_, size_t row_to_cut, size_t col_to_cut, typename enable_if< O==M-1 &&P==N-1 &&M==N &&M >=2 >::type *) const {size_t row_offset=0;size_t col_offset=0;for(size_t row_index=0;row_index< M;++row_index){if(row_index==row_to_cut) row_offset=-1;else{for(size_t col_index=0;col_index< M;++col_index){if(col_index==col_to_cut) col_offset=-1;else minor_.at(row_index+row_offset, col_index+col_offset)=at(row_index, col_index);}col_offset=0;}}return compute_determinant(minor_);}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::rotate(const TT angle_, const vector< M-1, T > &axis, typename enable_if< M==N &&M==4, TT >::type *){const T angle=static_cast< T > angle_);const T sine=std::sin(angle);const T cosine=std::cos(angle);const T _zero=0.0;const T one=1.0;const T two=2.0;array[0]=cosine+(one-cosine)*std::pow(axis.array[0], two);array[1]=(one-cosine)*axis.array[0]*axis.array[1]+sine *axis.array[2];array[2]=(one-cosine)*axis.array[0]*axis.array[2]-sine *axis.array[1];array[3]=_zero;array[4]=(one-cosine)*axis.array[0]*axis.array[1]-sine *axis.array[2];array[5]=cosine+(one-cosine)*std::pow(axis.array[1], two);array[6]=(one-cosine)*axis.array[1]*axis.array[2]+sine *axis.array[0];array[7]=_zero;array[8]=(one-cosine)*axis.array[0]*axis.array[2]+sine *axis.array[1];array[9]=(one-cosine)*axis.array[1]*axis.array[2]-sine *axis.array[0];array[10]=cosine+(one-cosine)*std::pow(axis.array[2], two);array[11]=_zero;array[12]=_zero;array[13]=_zero;array[14]=_zero;array[15]=one;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::rotate_x(const TT angle_, typename enable_if< M==N &&M==4, TT >::type *){const T angle=static_cast< T > angle_);const T sine=std::sin(angle);const T cosine=std::cos(angle);T tmp;tmp=array[4]*cosine+array[8]*sine;array[8]=-array[4]*sine+array[8]*cosine;array[4]=tmp;tmp=array[5]*cosine+array[9]*sine;array[9]=-array[5]*sine+array[9]*cosine;array[5]=tmp;tmp=array[6]*cosine+array[10]*sine;array[10]=-array[6]*sine+array[10]*cosine;array[6]=tmp;tmp=array[7]*cosine+array[11]*sine;array[11]=-array[7]*sine+array[11]*cosine;array[7]=tmp;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::rotate_y(const TT angle_, typename enable_if< M==N &&M==4, TT >::type *){const T angle=static_cast< T > angle_);const T sine=std::sin(angle);const T cosine=std::cos(angle);T tmp;tmp=array[0]*cosine-array[8]*sine;array[8]=array[0]*sine+array[8]*cosine;array[0]=tmp;tmp=array[1]*cosine-array[9]*sine;array[9]=array[1]*sine+array[9]*cosine;array[1]=tmp;tmp=array[2]*cosine-array[10]*sine;array[10]=array[2]*sine+array[10]*cosine;array[2]=tmp;tmp=array[3]*cosine-array[11]*sine;array[11]=array[3]*sine+array[11]*cosine;array[3]=tmp;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::rotate_z(const TT angle_, typename enable_if< M==N &&M==4, TT >::type *){const T angle=static_cast< T > angle_);const T sine=std::sin(angle);const T cosine=std::cos(angle);T tmp;tmp=array[0]*cosine+array[4]*sine;array[4]=-array[0]*sine+array[4]*cosine;array[0]=tmp;tmp=array[1]*cosine+array[5]*sine;array[5]=-array[1]*sine+array[5]*cosine;array[1]=tmp;tmp=array[2]*cosine+array[6]*sine;array[6]=-array[2]*sine+array[6]*cosine;array[2]=tmp;tmp=array[3]*cosine+array[7]*sine;array[7]=-array[3]*sine+array[7]*cosine;array[3]=tmp;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::pre_rotate_x(const TT angle_, typename enable_if< M==N &&M==4, TT >::type *){const T angle=static_cast< T > angle_);const T sine=std::sin(angle);const T cosine=std::cos(angle);T tmp;tmp=array[1];array[1]=array[1]*cosine+array[2]*sine;array[2]=tmp *-sine+array[2]*cosine;tmp=array[5];array[5]=array[5]*cosine+array[6]*sine;array[6]=tmp *-sine+array[6]*cosine;tmp=array[9];array[9]=array[9]*cosine+array[10]*sine;array[10]=tmp *-sine+array[10]*cosine;tmp=array[13];array[13]=array[13]*cosine+array[14]*sine;array[14]=tmp *-sine+array[14]*cosine;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::pre_rotate_y(const TT angle_, typename enable_if< M==N &&M==4, TT >::type *){const T angle=static_cast< T > angle_);const T sine=std::sin(angle);const T cosine=std::cos(angle);T tmp;tmp=array[0];array[0]=array[0]*cosine-array[2]*sine;array[2]=tmp *sine+array[2]*cosine;tmp=array[4];array[4]=array[4]*cosine-array[6]*sine;array[6]=tmp *sine+array[6]*cosine;tmp=array[8];array[8]=array[8]*cosine-array[10]*sine;array[10]=tmp *sine+array[10]*cosine;tmp=array[12];array[12]=array[12]*cosine-array[14]*sine;array[14]=tmp *sine+array[14]*cosine;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::pre_rotate_z(const TT angle_, typename enable_if< M==N &&M==4, TT >::type *){const T angle=static_cast< T > angle_);const T sine=std::sin(angle);const T cosine=std::cos(angle);T tmp;tmp=array[0];array[0]=array[0]*cosine+array[1]*sine;array[1]=tmp *-sine+array[1]*cosine;tmp=array[4];array[4]=array[4]*cosine+array[5]*sine;array[5]=tmp *-sine+array[5]*cosine;tmp=array[8];array[8]=array[8]*cosine+array[9]*sine;array[9]=tmp *-sine+array[9]*cosine;tmp=array[12];array[12]=array[12]*cosine+array[13]*sine;array[13]=tmp *-sine+array[13]*cosine;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::scale(const TT _scale[3], typename enable_if< M==N &&M==4, TT >::type *){const T scale0=static_cast< T > _scale[0]);const T scale1=static_cast< T > _scale[1]);const T scale2=static_cast< T > _scale[2]);array[0]*=scale0;array[1]*=scale0;array[2]*=scale0;array[3]*=scale0;array[4]*=scale1;array[5]*=scale1;array[6]*=scale1;array[7]*=scale1;array[8]*=scale2;array[9]*=scale2;array[10]*=scale2;array[11]*=scale2;return *this;}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::scale(const TT x_, const T y_, const T z_, typename enable_if< M==N &&M==4, TT >::type *){const T _x=static_cast< T > x_);array[0]*=_x;array[1]*=_x;array[2]*=_x;array[3]*=_x;array[4]*=y_;array[5]*=y_;array[6]*=y_;array[7]*=y_;array[8]*=z_;array[9]*=z_;array[10]*=z_;array[11]*=z_;return *this;}template< size_t M, size_t N, typename T >template< typename TT >inline matrix< M, N, T > &matrix< M, N, T >::scale(const vector< 3, TT > &scale_, typename enable_if< M==N &&M==4, TT >::type *){return scale(scale_.array);}template< size_t M, size_t N, typename T >template< typename TT >matrix< M, N, T > &matrix< M, N, T >::scale_translation(const TT scale_[3], typename enable_if< M==N &&M==4, TT >::type *){array[12]*=static_cast< T > scale_[0]);array[13]*=static_cast< T > scale_[1]);array[14]*=static_cast< T > scale_[2]);return *this;}template< size_t M, size_t N, typename T >template< typename TT >inline matrix< M, N, T > &matrix< M, N, T >::scale_translation(const vector< 3, TT > &scale_, typename enable_if< M==N &&M==4, TT >::type *){return scale_translation(scale_.array);}template< size_t M, size_t N, typename T >template< typename TT >inline matrix< M, N, T > &matrix< M, N, T >::set_translation(const TT x_, const TT y_, const TT z_, typename enable_if< M==N &&M==4, TT >::type *){array[12]=static_cast< T > x_);array[13]=static_cast< T > y_);array[14]=static_cast< T > z_);return *this;}template< size_t M, size_t N, typename T >template< typename TT >inline matrix< M, N, T > &matrix< M, N, T >::set_translation(const TT trans[3], typename enable_if< M==N &&M==4, TT >::type *){array[12]=static_cast< T > trans[0]);array[13]=static_cast< T > trans[1]);array[14]=static_cast< T > trans[2]);return *this;}template< size_t M, size_t N, typename T >template< typename TT >inline matrix< M, N, T > &matrix< M, N, T >::set_translation(const vector< 3, TT > &translation_, typename enable_if< M==N &&M==4, TT >::type *){return set_translation(translation_.array);}template< size_t M, size_t N, typename T > template< typename TT > inlinevoid matrix< M, N, T >::get_translation(vector< N-1, TT > &translation_) const {for(size_t i=0;i< N-1;++i) translation_.array[i]=array[i+M *(N-1)];}template< size_t M, size_t N, typename T >inline vector< N-1, T > matrix< M, N, T >::get_translation() const {vector< N-1, T > result;get_translation(result);return result;}template< size_t M, size_t N, typename T >template< size_t O >void matrix< M, N, T >::getLookAt(vector< O, T > &eye, vector< O, T > &lookAt, vector< O, T > &up, typename enable_if< M==O+1 &&N==O+1 &&O==3 >::type *) const {matrix< 4, 4, T > inv;inverse(inv);const matrix< 3, 3, T > rotation(transpose(matrix< 3, 3, T >(*this)));eye=inv *vector< 3, T >::ZERO;up=rotation *vector< 3, T >::UP;lookAt=rotation *vector< 3, T >::FORWARD;lookAt.normalize();lookAt=eye+lookAt;}template< size_t M, size_t N, typename T >size_tmatrix< M, N, T >::size() const {return M *N;}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::iteratormatrix< M, N, T >::begin(){return array;}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::iteratormatrix< M, N, T >::end(){return array+size();}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::const_iteratormatrix< M, N, T >::begin() const {return array;}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::const_iteratormatrix< M, N, T >::end() const {return array+size();}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::reverse_iteratormatrix< M, N, T >::rbegin(){return array+size()-1;}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::reverse_iteratormatrix< M, N, T >::rend(){return array-1;}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::const_reverse_iteratormatrix< M, N, T >::rbegin() const {return array+size()-1;}template< size_t M, size_t N, typename T >typename matrix< M, N, T >::const_reverse_iteratormatrix< M, N, T >::rend() const {return array-1;}template< size_t M, size_t N, typename T > template< typename init_functor_t >const matrix< M, N, T > matrix< M, N, T >::get_initialized_matrix(){matrix< M, N, T > matrix_;init_functor_t()(matrix_);return matrix_;}template< size_t M, size_t N, typename T >const matrix< M, N, T >matrix< M, N, T >::IDENTITY(matrix< M, N, T >::get_initialized_matrix< set_to_identity_functor< matrix< M, N, T > > >));template< size_t M, size_t N, typename T >const matrix< M, N, T >matrix< M, N, T >::ZERO(matrix< M, N, T >::get_initialized_matrix< set_to_zero_functor< matrix< M, N, T > > >));template< size_t M, size_t N, typename T >doublematrix< M, N, T >::frobenius_norm() const {double norm=0.0;const_iterator it=begin(), it_end=end();for(;it!=it_end;++it){norm+=*it **it;}return std::sqrt(norm);}template< size_t M, size_t N, typename T >doublematrix< M, N, T >::p_norm(double p) const {double norm=0.0;const_iterator it=begin(), it_end=end();for(;it!=it_end;++it){norm+=std::pow(*it, p);}return std::pow(norm, 1./p);}template< size_t M, size_t N, typename T >template< size_t O >voidmatrix< M, N, T >::khatri_rao_product(const matrix< O, N, T > &right_, matrix< M *O, N, T > &prod_) const {for(size_t col=0;col< N;++col){for(size_t m=0;m< M;++m){for(size_t o=0;o< O;++o){prod_.at(O *m+o, col)=at(m, col)*right_.at(o, col);}}}}template< size_t M, size_t N, typename T >template< size_t O, size_t P >voidmatrix< M, N, T >::kronecker_product(const matrix< O, P, T > &right_, matrix< M *O, N *P, T > &result_) const {for(size_t m=0;m< M;++m){for(size_t n=0;n< N;++n){for(size_t o=0;o< O;++o){for(size_t p=0;p< P;++p){result_.at(O *m+o, P *n+p)=at(m, n)*right_.at(o, p);}}}}}template< size_t M, size_t N, typename T >template< typename TT >voidmatrix< M, N, T >::cast_from(const matrix< M, N, TT > &other){typedef vmml::matrix< M, N, TT > matrix_tt_type;typedef typename matrix_tt_type::const_iterator tt_const_iterator;iterator it=begin(), it_end=end();tt_const_iterator other_it=other.begin();for(;it!=it_end;++it,++other_it){*it=static_cast< T > * | other_it |