_Matrix3.cc

Go to the documentation of this file.

#include <vecmath/GMatrix.h>
#include <vecmath/Matrix3SVD.h>

template<> const double _Matrix3<double>::EPS = 1e-16;
template<> const double _Matrix3<double>::ERR_EPS = 1e-08;
template<> const float _Matrix3<float>::EPS = 1e-16;
template<> const float _Matrix3<float>::ERR_EPS = 1e-08;

template <class Type>
void _Matrix3<Type>::add(Type d)
{
    m00 += d;
    m01 += d;
    m02 += d;
    m10 += d;
    m11 += d;
    m12 += d;
    m20 += d;
    m21 += d;
    m22 += d;
}

template <class Type>
void _Matrix3<Type>::add(Type d, const _Matrix3<Type>& m)
{
    set(m.m00 + d, m.m01 + d, m.m02 + d,
        m.m10 + d, m.m11 + d, m.m12 + d,
        m.m20 + d, m.m21 + d, m.m22 + d);
}

template <class Type>
void _Matrix3<Type>::add(const _Matrix3<Type>& m)
{
    m00 += m.m00;
    m01 += m.m01;
    m02 += m.m02;
    m10 += m.m10;
    m11 += m.m11;
    m12 += m.m12;
    m20 += m.m20;
    m21 += m.m21;
    m22 += m.m22;
}

template <class Type>
void _Matrix3<Type>::add(const _Matrix3<Type>& m, const _Matrix3<Type>& m1)
{
    set(m.m00 + m1.m00, m.m01 + m1.m01, m.m02 + m1.m02,
        m.m10 + m1.m10, m.m11 + m1.m11, m.m12 + m1.m12,
        m.m20 + m1.m20, m.m21 + m1.m21, m.m22 + m1.m22);
}

template <class Type>
void _Matrix3<Type>::compute_svd(double m[9], double scale[3], double rot[9], bool flag)
{
#if 0
    GMatrix gad(3, 3, ad, true);
    GMatrix gad1(3, 3, ad1, true);
    GMatrix gad2(3, 3, ad2, true);
    
    gad.computeSVD(gad, gad, gad1, gad2, flag);
#endif
    
    double scale_a[9];
    Matrix3SVD::compute_svd(m, scale_a, rot, false);
    scale[0] = scale_a[0];
    scale[1] = scale_a[1];
    scale[2] = scale_a[2];
//  cerr << "Matrix3x3 compute_svd not properly implemented - check!" << endl;
}

template <class Type>
bool _Matrix3<Type>::epsilonEquals(const _Matrix3<Type>& m, Type epsilon) const
{
    return (Math::epsilonEquals(m00, m.m00, epsilon) &&
            Math::epsilonEquals(m01, m.m01, epsilon) &&
            Math::epsilonEquals(m02, m.m02, epsilon) &&
            Math::epsilonEquals(m10, m.m10, epsilon) &&
            Math::epsilonEquals(m11, m.m11, epsilon) &&
            Math::epsilonEquals(m12, m.m12, epsilon) &&
            Math::epsilonEquals(m20, m.m20, epsilon) &&
            Math::epsilonEquals(m21, m.m21, epsilon) &&
            Math::epsilonEquals(m22, m.m22, epsilon));
}

template <class Type>
bool _Matrix3<Type>::equals(const _Matrix3<Type>& m) const
{
    return (m00 == m.m00 && m01 == m.m01 && m02 == m.m02
            && m10 == m.m10 && m11 == m.m11 && m12 == m.m12
            && m20 == m.m20 && m21 == m.m21 && m22 == m.m22);
}

template <class Type>
void _Matrix3<Type>::getColumn(int column, _Vector3<Type>& v) const
{
    switch (column) {
    case 0:
        v.set(m00, m10, m20);
        break;
    case 1:
        v.set(m01, m11, m21);
        break;
    case 2:
        v.set(m02, m12, m22);
        break;
    default:
        break;
    }
}

template <class Type>
void _Matrix3<Type>::getColumn(int column, Type ad[]) const
{
    switch (column) {
    case 0:
        ad[0] = m00;
        ad[1] = m10;
        ad[2] = m20;
        break;
    case 1:
        ad[0] = m01;
        ad[1] = m11;
        ad[2] = m21;
        break;
    case 2:
        ad[0] = m02;
        ad[1] = m12;
        ad[2] = m22;
        break;
    default:
        break;
    }
}

template <class Type>
Type _Matrix3<Type>::getElement(int row, int column) const
{
    switch (row) {
    case 0:
        switch (column) {
        case 0: return m00;
        case 1: return m01;
        case 2: return m02;
        }
        break;
    case 1:
        switch (column) {
        case 0: return m10;
        case 1: return m11;
        case 2: return m12;
        }
        break;
    case 2:
        switch (column) {
        case 0: return m20;
        case 1: return m21;
        case 2: return m22;
        }
        break;
    }
    return 0;
}

template <class Type>
void _Matrix3<Type>::getRow(int row, _Vector3<Type>& v) const
{
    switch (row) {
    case 0:
        v.set(m00, m01, m02);
        break;
    case 1:
        v.set(m10, m11, m12);
        break;
    case 2:
        v.set(m20, m21, m22);
        break;
    default:
        break;
    }
}

template <class Type>
void _Matrix3<Type>::getRow(int row, Type ad[3]) const
{
    switch (row) {
    case 0:
        ad[0] = m00;
        ad[1] = m01;
        ad[2] = m02;
        break;
    case 1:
        ad[0] = m10;
        ad[1] = m11;
        ad[2] = m12;
        break;
    case 2:
        ad[0] = m20;
        ad[1] = m21;
        ad[2] = m22;
        break;
    default:
        break;
    }
}

template <class Type>
void _Matrix3<Type>::invert(const _Matrix3<Type>& m)
{
    Type d = m.determinant();

    if (this != &m) {
        m00 = (m.m11 * m.m22 - m.m21 * m.m12) / d;
        m01 = (m.m02 * m.m21 - m.m22 * m.m01) / d;
        m02 = (m.m01 * m.m12 - m.m11 * m.m02) / d;

        m10 = (m.m12 * m.m20 - m.m22 * m.m10) / d;
        m11 = (m.m00 * m.m22 - m.m20 * m.m02) / d;
        m12 = (m.m02 * m.m10 - m.m12 * m.m00) / d;

        m20 = (m.m10 * m.m21 - m.m20 * m.m11) / d;
        m21 = (m.m01 * m.m20 - m.m21 * m.m00) / d;
        m22 = (m.m00 * m.m11 - m.m10 * m.m01) / d;
    } else {
        Type ad[9];

        ad[0] = (m.m11 * m.m22 - m.m21 * m.m12) / d;
        ad[1] = (m.m02 * m.m21 - m.m22 * m.m01) / d;
        ad[2] = (m.m01 * m.m12 - m.m11 * m.m02) / d;

        ad[3] = (m.m12 * m.m20 - m.m22 * m.m10) / d;
        ad[4] = (m.m00 * m.m22 - m.m20 * m.m02) / d;
        ad[5] = (m.m02 * m.m10 - m.m12 * m.m00) / d;

        ad[6] = (m.m10 * m.m21 - m.m20 * m.m11) / d;
        ad[7] = (m.m01 * m.m20 - m.m21 * m.m00) / d;
        ad[8] = (m.m00 * m.m11 - m.m10 * m.m01) / d;
        set(ad);
    }
}

template <class Type>
void _Matrix3<Type>::mul(Type d)
{
    m00 *= d;
    m01 *= d;
    m02 *= d;
    m10 *= d;
    m11 *= d;
    m12 *= d;
    m20 *= d;
    m21 *= d;
    m22 *= d;
}

template <class Type>
void _Matrix3<Type>::mul(Type d, const _Matrix3<Type>& m)
{
    set(m.m00 * d, m.m01 * d, m.m02 * d,
        m.m10 * d, m.m11 * d, m.m12 * d,
        m.m20 * d, m.m21 * d, m.m22 * d);
}


template <class Type>
void _Matrix3<Type>::mul(const _Matrix3<Type>& m, const _Matrix3<Type>& m1)
{
    if(this != &m && this != &m1) {
        m00 = m.m00 * m1.m00 + m.m01 * m1.m10 + m.m02 * m1.m20;
        m01 = m.m00 * m1.m01 + m.m01 * m1.m11 + m.m02 * m1.m21;
        m02 = m.m00 * m1.m02 + m.m01 * m1.m12 + m.m02 * m1.m22;
        m10 = m.m10 * m1.m00 + m.m11 * m1.m10 + m.m12 * m1.m20;
        m11 = m.m10 * m1.m01 + m.m11 * m1.m11 + m.m12 * m1.m21;
        m12 = m.m10 * m1.m02 + m.m11 * m1.m12 + m.m12 * m1.m22;
        m20 = m.m20 * m1.m00 + m.m21 * m1.m10 + m.m22 * m1.m20;
        m21 = m.m20 * m1.m01 + m.m21 * m1.m11 + m.m22 * m1.m21;
        m22 = m.m20 * m1.m02 + m.m21 * m1.m12 + m.m22 * m1.m22;
    } else {
        Type d0 = m.m00 * m1.m00 + m.m01 * m1.m10 + m.m02 * m1.m20;
        Type d1 = m.m00 * m1.m01 + m.m01 * m1.m11 + m.m02 * m1.m21;
        Type d2 = m.m00 * m1.m02 + m.m01 * m1.m12 + m.m02 * m1.m22;
        Type d3 = m.m10 * m1.m00 + m.m11 * m1.m10 + m.m12 * m1.m20;
        Type d4 = m.m10 * m1.m01 + m.m11 * m1.m11 + m.m12 * m1.m21;
        Type d5 = m.m10 * m1.m02 + m.m11 * m1.m12 + m.m12 * m1.m22;
        Type d6 = m.m20 * m1.m00 + m.m21 * m1.m10 + m.m22 * m1.m20;
        Type d7 = m.m20 * m1.m01 + m.m21 * m1.m11 + m.m22 * m1.m21;
        Type d8 = m.m20 * m1.m02 + m.m21 * m1.m12 + m.m22 * m1.m22;
        set(d0, d1, d2,
            d3, d4, d5,
            d6, d7, d8);
    }
}

template <class Type>
void _Matrix3<Type>::mulNormalize(const _Matrix3<Type>& m, const _Matrix3<Type>& m1)
{
    double tmp[9];
    double tmp_scale[9];
    double tmp_rot[9];
    tmp[0] = m.m00 * m1.m00 + m.m01 * m1.m10 + m.m02 * m1.m20;
    tmp[1] = m.m00 * m1.m01 + m.m01 * m1.m11 + m.m02 * m1.m21;
    tmp[2] = m.m00 * m1.m02 + m.m01 * m1.m12 + m.m02 * m1.m22;
    tmp[3] = m.m10 * m1.m00 + m.m11 * m1.m10 + m.m12 * m1.m20;
    tmp[4] = m.m10 * m1.m01 + m.m11 * m1.m11 + m.m12 * m1.m21;
    tmp[5] = m.m10 * m1.m02 + m.m11 * m1.m12 + m.m12 * m1.m22;
    tmp[6] = m.m20 * m1.m00 + m.m21 * m1.m10 + m.m22 * m1.m20;
    tmp[7] = m.m20 * m1.m01 + m.m21 * m1.m11 + m.m22 * m1.m21;
    tmp[8] = m.m20 * m1.m02 + m.m21 * m1.m12 + m.m22 * m1.m22;
    compute_svd(tmp, tmp_scale, tmp_rot, false);
    set(tmp_rot);
}

template <class Type>
void _Matrix3<Type>::mulTransposeBoth(const _Matrix3<Type> m, const _Matrix3<Type>& m1)
{
    if(this != &m && this != &m1) {
        m00 = m.m00 * m1.m00 + m.m10 * m1.m01 + m.m20 * m1.m02;
        m01 = m.m00 * m1.m10 + m.m10 * m1.m11 + m.m20 * m1.m12;
        m02 = m.m00 * m1.m20 + m.m10 * m1.m21 + m.m20 * m1.m22;
        m10 = m.m01 * m1.m00 + m.m11 * m1.m01 + m.m21 * m1.m02;
        m11 = m.m01 * m1.m10 + m.m11 * m1.m11 + m.m21 * m1.m12;
        m12 = m.m01 * m1.m20 + m.m11 * m1.m21 + m.m21 * m1.m22;
        m20 = m.m02 * m1.m00 + m.m12 * m1.m01 + m.m22 * m1.m02;
        m21 = m.m02 * m1.m10 + m.m12 * m1.m11 + m.m22 * m1.m12;
        m22 = m.m02 * m1.m20 + m.m12 * m1.m21 + m.m22 * m1.m22;
    } else {
        Type d0 = m.m00 * m1.m00 + m.m10 * m1.m01 + m.m20 * m1.m02;
        Type d1 = m.m00 * m1.m10 + m.m10 * m1.m11 + m.m20 * m1.m12;
        Type d2 = m.m00 * m1.m20 + m.m10 * m1.m21 + m.m20 * m1.m22;
        Type d3 = m.m01 * m1.m00 + m.m11 * m1.m01 + m.m21 * m1.m02;
        Type d4 = m.m01 * m1.m10 + m.m11 * m1.m11 + m.m21 * m1.m12;
        Type d5 = m.m01 * m1.m20 + m.m11 * m1.m21 + m.m21 * m1.m22;
        Type d6 = m.m02 * m1.m00 + m.m12 * m1.m01 + m.m22 * m1.m02;
        Type d7 = m.m02 * m1.m10 + m.m12 * m1.m11 + m.m22 * m1.m12;
        Type d8 = m.m02 * m1.m20 + m.m12 * m1.m21 + m.m22 * m1.m22;
        set(d0, d1, d2,
            d3, d4, d5,
            d6, d7, d8);
    }
}

template <class Type>
void _Matrix3<Type>::mulTransposeLeft(const _Matrix3<Type>& m, const _Matrix3<Type>& m1)
{
    if(this != &m && this != &m1) {
        m00 = m.m00 * m1.m00 + m.m10 * m1.m10 + m.m20 * m1.m20;
        m01 = m.m00 * m1.m01 + m.m10 * m1.m11 + m.m20 * m1.m21;
        m02 = m.m00 * m1.m02 + m.m10 * m1.m12 + m.m20 * m1.m22;
        m10 = m.m01 * m1.m00 + m.m11 * m1.m10 + m.m21 * m1.m20;
        m11 = m.m01 * m1.m01 + m.m11 * m1.m11 + m.m21 * m1.m21;
        m12 = m.m01 * m1.m02 + m.m11 * m1.m12 + m.m21 * m1.m22;
        m20 = m.m02 * m1.m00 + m.m12 * m1.m10 + m.m22 * m1.m20;
        m21 = m.m02 * m1.m01 + m.m12 * m1.m11 + m.m22 * m1.m21;
        m22 = m.m02 * m1.m02 + m.m12 * m1.m12 + m.m22 * m1.m22;
    } else {
        Type d0 = m.m00 * m1.m00 + m.m10 * m1.m10 + m.m20 * m1.m20;
        Type d1 = m.m00 * m1.m01 + m.m10 * m1.m11 + m.m20 * m1.m21;
        Type d2 = m.m00 * m1.m02 + m.m10 * m1.m12 + m.m20 * m1.m22;
        Type d3 = m.m01 * m1.m00 + m.m11 * m1.m10 + m.m21 * m1.m20;
        Type d4 = m.m01 * m1.m01 + m.m11 * m1.m11 + m.m21 * m1.m21;
        Type d5 = m.m01 * m1.m02 + m.m11 * m1.m12 + m.m21 * m1.m22;
        Type d6 = m.m02 * m1.m00 + m.m12 * m1.m10 + m.m22 * m1.m20;
        Type d7 = m.m02 * m1.m01 + m.m12 * m1.m11 + m.m22 * m1.m21;
        Type d8 = m.m02 * m1.m02 + m.m12 * m1.m12 + m.m22 * m1.m22;
        set(d0, d1, d2,
            d3, d4, d5,
            d6, d7, d8);
    }
}

template <class Type>
void _Matrix3<Type>::mulTransposeRight(const _Matrix3<Type>& m, const _Matrix3<Type>& m1)
{
    if(this != &m && this != &m1) {
        m00 = m.m00 * m1.m00 + m.m01 * m1.m01 + m.m02 * m1.m02;
        m01 = m.m00 * m1.m10 + m.m01 * m1.m11 + m.m02 * m1.m12;
        m02 = m.m00 * m1.m20 + m.m01 * m1.m21 + m.m02 * m1.m22;
        m10 = m.m10 * m1.m00 + m.m11 * m1.m01 + m.m12 * m1.m02;
        m11 = m.m10 * m1.m10 + m.m11 * m1.m11 + m.m12 * m1.m12;
        m12 = m.m10 * m1.m20 + m.m11 * m1.m21 + m.m12 * m1.m22;
        m20 = m.m20 * m1.m00 + m.m21 * m1.m01 + m.m22 * m1.m02;
        m21 = m.m20 * m1.m10 + m.m21 * m1.m11 + m.m22 * m1.m12;
        m22 = m.m20 * m1.m20 + m.m21 * m1.m21 + m.m22 * m1.m22;
    } else {
        Type d0 = m.m00 * m1.m00 + m.m01 * m1.m01 + m.m02 * m1.m02;
        Type d1 = m.m00 * m1.m10 + m.m01 * m1.m11 + m.m02 * m1.m12;
        Type d2 = m.m00 * m1.m20 + m.m01 * m1.m21 + m.m02 * m1.m22;
        Type d3 = m.m10 * m1.m00 + m.m11 * m1.m01 + m.m12 * m1.m02;
        Type d4 = m.m10 * m1.m10 + m.m11 * m1.m11 + m.m12 * m1.m12;
        Type d5 = m.m10 * m1.m20 + m.m11 * m1.m21 + m.m12 * m1.m22;
        Type d6 = m.m20 * m1.m00 + m.m21 * m1.m01 + m.m22 * m1.m02;
        Type d7 = m.m20 * m1.m10 + m.m21 * m1.m11 + m.m22 * m1.m12;
        Type d8 = m.m20 * m1.m20 + m.m21 * m1.m21 + m.m22 * m1.m22;
        set(d0, d1, d2,
            d3, d4, d5,
            d6, d7, d8);
    }
}

template <class Type>
void _Matrix3<Type>::negate(const _Matrix3<Type>& m)
{
    set(-m.m00, -m.m01, -m.m02,
        -m.m10, -m.m11, -m.m12,
        -m.m20, -m.m21, -m.m22);
}

template <class Type>
void _Matrix3<Type>::normalize()
{
    double tmp_scale[9];
    double tmp_rot[9];
    getScaleRotate(tmp_scale, tmp_rot);
    set(tmp_rot);
}

template <class Type>
void _Matrix3<Type>::normalize(const _Matrix3<Type>& m)
{
    double tmp[9];
    double tmp_scale[9];
    double tmp_rot[9];
    m.get(tmp);
    compute_svd(tmp, tmp_scale, tmp_rot, false);
    set(tmp_rot);
}

template <class Type>
void _Matrix3<Type>::normalizeCP(const _Matrix3<Type>& m)
{
    Type d = sqrt(m.m00 * m.m00 + m.m10 * m.m10 + m.m20 * m.m20);
    m00 = m.m00 / d;
    m10 = m.m10 / d;
    m20 = m.m20 / d;
    d = sqrt(m.m01 * m.m01 + m.m11 * m.m11 + m.m21 * m.m21);
    m01 = m.m01 / d;
    m11 = m.m11 / d;
    m21 = m.m21 / d;
    m02 = m10 * m21 - m11 * m20;
    m12 = m01 * m20 - m00 * m21;
    m22 = m00 * m11 - m01 * m10;
}

template <class Type>
void _Matrix3<Type>::rotX(Type angle)
{
    Type s = sin(angle);
    Type c = cos(angle);

    set(1.0, 0.0, 0.0,
        0.0,   c,  -s,
        0.0,   s,   c);
}

template <class Type>
void _Matrix3<Type>::rotY(Type angle)
{
    Type s = sin(angle);
    Type c = cos(angle);

    set(c  , 0.0,   s,
        0.0, 1.0, 0.0,
        -s , 0.0,   c);
}

template <class Type>
void _Matrix3<Type>::rotZ(Type angle)
{
    Type s = sin(angle);
    Type c = cos(angle);

    set(c  ,  -s, 0.0,
        s  ,   c, 0.0,
        0.0, 0.0, 1.0);
}

template <class Type>
void _Matrix3<Type>::setAxisAngle(Type x, Type y, Type z, Type angle)
{
    Type d = sqrt(x * x + y * y + z * z);

    if (Math::epsilonEquals(d)) {
        setIdentity();
    } else {
        Type d1 = x / d;
        Type d2 = y / d;
        Type d3 = z / d;
        Type d4 = sin(angle);
        Type d5 = cos(angle);
        Type d6 = 1.0 - d5;
        Type d7 = x * z;
        Type d8 = x * y;
        Type d9 = y * z;
        set(d6 * d1 * d1 + d5, d6 * d8 - d4 * d3, d6 * d7 + d4 * d2,
            d6 * d8 + d4 * d3, d6 * d2 * d2 + d5, d6 * d9 - d4 * d1,
            d6 * d7 - d4 * d2, d6 * d9 + d4 * d1, d6 * d3 * d3 + d5);
    }
}

template <class Type>
void _Matrix3<Type>::setQuat(Type x, Type y, Type z, Type w)
{
    set(1.0 - 2.0 * y * y - 2.0 * z * z,
        2.0 * (x * y + w * z),
        2.0 * (x * z - w * y),
        2.0 * (x * y - w * z),
        1.0 - 2.0 * x * x - 2.0 * z * z,
        2.0 * (y * z + w * x),
        2.0 * (x * z + w * y),
        2.0 * (y * z - w * x),
        1.0 - 2.0 * x * x - 2.0 * y * y);
}

template <class Type>
void _Matrix3<Type>::setRow(int row, Type x, Type y, Type z)
{
    switch(row) {
    case 0:
        m00 = x;
        m01 = y;
        m02 = z;
        break;
    case 1:
        m10 = x;
        m11 = y;
        m12 = z;
        break;
    case 2:
        m20 = x;
        m21 = y;
        m22 = z;
        break;
    default:
        break;
    }
}

template <class Type>
void _Matrix3<Type>::setRow(int row, const Type v[3])
{
    switch(row) {
    case 0:
        m00 = v[0];
        m01 = v[1];
        m02 = v[2];
        break;
    case 1:
        m10 = v[0];
        m11 = v[1];
        m12 = v[2];
        break;
    case 2:
        m20 = v[0];
        m21 = v[1];
        m22 = v[2];
        break;
    default:
        break;
    }
}

template <class Type>
void _Matrix3<Type>::setRow(int row, const _Vector3<Type>& v)
{
    switch(row) {
    case 0:
        m00 = v.x;
        m01 = v.y;
        m02 = v.z;
        break;
    case 1:
        m10 = v.x;
        m11 = v.y;
        m12 = v.z;
        break;
    case 2:
        m20 = v.x;
        m21 = v.y;
        m22 = v.z;
        break;
    default:
        break;
    }
}

template <class Type>
void _Matrix3<Type>::sub(const _Matrix3<Type>& m)
{
    m00 -= m.m00;
    m01 -= m.m01;
    m02 -= m.m02;
    m10 -= m.m10;
    m11 -= m.m11;
    m12 -= m.m12;
    m20 -= m.m20;
    m21 -= m.m21;
    m22 -= m.m22;
}

template <class Type>
void _Matrix3<Type>::sub(const _Matrix3<Type>& m, const _Matrix3<Type>& m1)
{
    set(m.m00 - m1.m00, m.m01 - m1.m01, m.m02 - m1.m02,
        m.m10 - m1.m10, m.m11 - m1.m11, m.m12 - m1.m12,
        m.m20 - m1.m20, m.m21 - m1.m21, m.m22 - m1.m22);
}

template <class Type>
void _Matrix3<Type>::transform(const _Tuple3<Type>& t, _Tuple3<Type>& t1) const
{
    Type x = t.x, y = t.y;
    t1.x = m00 * x + m01 * y + m02 * t.z;
    t1.y = m10 * x + m11 * y + m12 * t.z;
    t1.z = m20 * x + m21 * y + m22 * t.z;
}

template <class Type>
void _Matrix3<Type>::transpose(const _Matrix3<Type>& m)
{
    if (this != &m) {
        set(m.m00, m.m10, m.m20,
            m.m01, m.m11, m.m21,
            m.m02, m.m12, m.m22);
    } else
        transpose();
}

template <class Type>
std::ostream& operator <<(std::ostream& os, const _Matrix3<Type>& m)
{
    os << "{ ";
    for (int i = 0; i < 9;) {
        double d = m.getElement(i/3, i%3);
        if (Math::abs(d) < m.EPS)
            d = 0.0;
        os << d;
        i++;
        if (i % 3)
            os << ", ";
        else {
            if (i < 9)
                os << std::endl << "  ";
            else
                os << " }" << std::endl;
        }
    }
    return os;
}

template <class Type>
std::ostream& operator <<(std::ostream& os, const _Matrix3<Type>* m)
{
    return os << *m;
}

template class _Matrix3<double>;
template class _Matrix3<float>;
template std::ostream& operator <<(std::ostream&, const _Matrix3<double>&);
template std::ostream& operator <<(std::ostream&, const _Matrix3<float>&);
template std::ostream& operator <<(std::ostream&, const _Matrix3<double>*);
template std::ostream& operator <<(std::ostream&, const _Matrix3<float>*);

---