///////////////////////////////////////////////////////////////////////////

//

// Copyright (c) 2002-2012, Industrial Light & Magic, a division of Lucas

// Digital Ltd. LLC

// 

// All rights reserved.

// 

// Redistribution and use in source and binary forms, with or without

// modification, are permitted provided that the following conditions are

// met:

// *       Redistributions of source code must retain the above copyright

// notice, this list of conditions and the following disclaimer.

// *       Redistributions in binary form must reproduce the above

// copyright notice, this list of conditions and the following disclaimer

// in the documentation and/or other materials provided with the

// distribution.

// *       Neither the name of Industrial Light & Magic nor the names of

// its contributors may be used to endorse or promote products derived

// from this software without specific prior written permission. 

// 

// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR

// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT

// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,

// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT

// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,

// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY

// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT

// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

//

///////////////////////////////////////////////////////////////////////////

#ifndef INCLUDED_IMATHMATRIX_H

#define INCLUDED_IMATHMATRIX_H

//----------------------------------------------------------------

//

//      2D (3x3) and 3D (4x4) transformation matrix templates.

//

//----------------------------------------------------------------

#include "ImathPlatform.h"

#include "ImathFun.h"

#include "ImathExc.h"

#include "ImathVec.h"

#include "ImathShear.h"

#include "ImathNamespace.h"

#include <cstring>

#include <iostream>

#include <iomanip>

#include <string.h>

#if (defined _WIN32 || defined _WIN64) && defined _MSC_VER

// suppress exception specification warnings

#pragma warning(disable:4290)

#endif

IMATH_INTERNAL_NAMESPACE_HEADER_ENTER

enum Uninitialized {UNINITIALIZED};

template <class T> class Matrix33

{

  public:

    //-------------------

    // Access to elements

    //-------------------

    T           x[3][3];

    T *         operator [] (int i);

    const T *   operator [] (int i) const;

    //-------------

    // Constructors

    //-------------

    Matrix33 (Uninitialized) {}

    Matrix33 ();

                                // 1 0 0

                                // 0 1 0

                                // 0 0 1

    Matrix33 (T a);

                                // a a a

                                // a a a

                                // a a a

    Matrix33 (const T a[3][3]);

                                // a[0][0] a[0][1] a[0][2]

                                // a[1][0] a[1][1] a[1][2]

                                // a[2][0] a[2][1] a[2][2]

    Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i);

                                // a b c

                                // d e f

                                // g h i

    //--------------------------------

    // Copy constructor and assignment

    //--------------------------------

    Matrix33 (const Matrix33 &v);

    template <class S> explicit Matrix33 (const Matrix33<S> &v);

    const Matrix33 &    operator = (const Matrix33 &v);

    const Matrix33 &    operator = (T a);

    //----------------------

    // Compatibility with Sb

    //----------------------

    T *                 getValue ();

    const T *           getValue () const;

    template <class S>

    void                getValue (Matrix33<S> &v) const;

    template <class S>

    Matrix33 &          setValue (const Matrix33<S> &v);

    template <class S>

    Matrix33 &          setTheMatrix (const Matrix33<S> &v);

    //---------

    // Identity

    //---------

    void                makeIdentity();

    //---------

    // Equality

    //---------

    bool                operator == (const Matrix33 &v) const;

    bool                operator != (const Matrix33 &v) const;

    //-----------------------------------------------------------------------

    // Compare two matrices and test if they are "approximately equal":

    //

    // equalWithAbsError (m, e)

    //

    //      Returns true if the coefficients of this and m are the same with

    //      an absolute error of no more than e, i.e., for all i, j

    //

    //      abs (this[i][j] - m[i][j]) <= e

    //

    // equalWithRelError (m, e)

    //

    //      Returns true if the coefficients of this and m are the same with

    //      a relative error of no more than e, i.e., for all i, j

    //

    //      abs (this[i] - v[i][j]) <= e * abs (this[i][j])

    //-----------------------------------------------------------------------

    bool                equalWithAbsError (const Matrix33<T> &v, T e) const;

    bool                equalWithRelError (const Matrix33<T> &v, T e) const;

    //------------------------

    // Component-wise addition

    //------------------------

    const Matrix33 &    operator += (const Matrix33 &v);

    const Matrix33 &    operator += (T a);

    Matrix33            operator + (const Matrix33 &v) const;

    //---------------------------

    // Component-wise subtraction

    //---------------------------

    const Matrix33 &    operator -= (const Matrix33 &v);

    const Matrix33 &    operator -= (T a);

    Matrix33            operator - (const Matrix33 &v) const;

    //------------------------------------

    // Component-wise multiplication by -1

    //------------------------------------

    Matrix33            operator - () const;

    const Matrix33 &    negate ();

    //------------------------------

    // Component-wise multiplication

    //------------------------------

    const Matrix33 &    operator *= (T a);

    Matrix33            operator * (T a) const;

    //-----------------------------------

    // Matrix-times-matrix multiplication

    //-----------------------------------

    const Matrix33 &    operator *= (const Matrix33 &v);

    Matrix33            operator * (const Matrix33 &v) const;

    //-----------------------------------------------------------------

    // Vector-times-matrix multiplication; see also the "operator *"

    // functions defined below.

    //

    // m.multVecMatrix(src,dst) implements a homogeneous transformation

    // by computing Vec3 (src.x, src.y, 1) * m and dividing by the

    // result's third element.

    //

    // m.multDirMatrix(src,dst) multiplies src by the upper left 2x2

    // submatrix, ignoring the rest of matrix m.

    //-----------------------------------------------------------------

    template <class S>

    void                multVecMatrix(const Vec2<S> &src, Vec2<S> &dst) const;

    template <class S>

    void                multDirMatrix(const Vec2<S> &src, Vec2<S> &dst) const;

    //------------------------

    // Component-wise division

    //------------------------

    const Matrix33 &    operator /= (T a);

    Matrix33            operator / (T a) const;

    //------------------

    // Transposed matrix

    //------------------

    const Matrix33 &    transpose ();

    Matrix33            transposed () const;

    //------------------------------------------------------------

    // Inverse matrix: If singExc is false, inverting a singular

    // matrix produces an identity matrix.  If singExc is true,

    // inverting a singular matrix throws a SingMatrixExc.

    //

    // inverse() and invert() invert matrices using determinants;

    // gjInverse() and gjInvert() use the Gauss-Jordan method.

    //

    // inverse() and invert() are significantly faster than

    // gjInverse() and gjInvert(), but the results may be slightly

    // less accurate.

    // 

    //------------------------------------------------------------

    const Matrix33 &    invert (bool singExc = false);

    Matrix33<T>         inverse (bool singExc = false) const;

    const Matrix33 &    gjInvert (bool singExc = false);

    Matrix33<T>         gjInverse (bool singExc = false) const;

    //------------------------------------------------

    // Calculate the matrix minor of the (r,c) element

    //------------------------------------------------

    T                   minorOf (const int r, const int c) const;

    //---------------------------------------------------

    // Build a minor using the specified rows and columns

    //---------------------------------------------------

    T                   fastMinor (const int r0, const int r1, 

                                   const int c0, const int c1) const;

    //------------

    // Determinant

    //------------

    T                   determinant() const;

    //-----------------------------------------

    // Set matrix to rotation by r (in radians)

    //-----------------------------------------

    template <class S>

    const Matrix33 &    setRotation (S r);

    //-----------------------------

    // Rotate the given matrix by r

    //-----------------------------

    template <class S>

    const Matrix33 &    rotate (S r);

    //--------------------------------------------

    // Set matrix to scale by given uniform factor

    //--------------------------------------------

    const Matrix33 &    setScale (T s);

    //------------------------------------

    // Set matrix to scale by given vector

    //------------------------------------

    template <class S>

    const Matrix33 &    setScale (const Vec2<S> &s);

    //----------------------

    // Scale the matrix by s

    //----------------------

    template <class S>

    const Matrix33 &    scale (const Vec2<S> &s);

    //------------------------------------------

    // Set matrix to translation by given vector

    //------------------------------------------

    template <class S>

    const Matrix33 &    setTranslation (const Vec2<S> &t);

    //-----------------------------

    // Return translation component

    //-----------------------------

    Vec2<T>             translation () const;

    //--------------------------

    // Translate the matrix by t

    //--------------------------

    template <class S>

    const Matrix33 &    translate (const Vec2<S> &t);

    //-----------------------------------------------------------

    // Set matrix to shear x for each y coord. by given factor xy

    //-----------------------------------------------------------

    template <class S>

    const Matrix33 &    setShear (const S &h);

    //-------------------------------------------------------------

    // Set matrix to shear x for each y coord. by given factor h[0]

    // and to shear y for each x coord. by given factor h[1]

    //-------------------------------------------------------------

    template <class S>

    const Matrix33 &    setShear (const Vec2<S> &h);

    //-----------------------------------------------------------

    // Shear the matrix in x for each y coord. by given factor xy

    //-----------------------------------------------------------

    template <class S>

    const Matrix33 &    shear (const S &xy);

    //-----------------------------------------------------------

    // Shear the matrix in x for each y coord. by given factor xy

    // and shear y for each x coord. by given factor yx

    //-----------------------------------------------------------

    template <class S>

    const Matrix33 &    shear (const Vec2<S> &h);

    //--------------------------------------------------------

    // Number of the row and column dimensions, since

    // Matrix33 is a square matrix.

    //--------------------------------------------------------

    static unsigned int dimensions() {return 3;}

    //-------------------------------------------------

    // Limitations of type T (see also class limits<T>)

    //-------------------------------------------------

    static T            baseTypeMin()           {return limits<T>::min();}

    static T            baseTypeMax()           {return limits<T>::max();}

    static T            baseTypeSmallest()      {return limits<T>::smallest();}

    static T            baseTypeEpsilon()       {return limits<T>::epsilon();}

    typedef T   BaseType;

    typedef Vec3<T> BaseVecType;

  private:

    template <typename R, typename S>

    struct isSameType

    {

        enum {value = 0};

    };

    template <typename R>

    struct isSameType<R, R>

    {

        enum {value = 1};

    };

};

template <class T> class Matrix44

{

  public:

    //-------------------

    // Access to elements

    //-------------------

    T           x[4][4];

    T *         operator [] (int i);

    const T *   operator [] (int i) const;

    //-------------

    // Constructors

    //-------------

    Matrix44 (Uninitialized) {}

    Matrix44 ();

                                // 1 0 0 0

                                // 0 1 0 0

                                // 0 0 1 0

                                // 0 0 0 1

    Matrix44 (T a);

                                // a a a a

                                // a a a a

                                // a a a a

                                // a a a a

    Matrix44 (const T a[4][4]) ;

                                // a[0][0] a[0][1] a[0][2] a[0][3]

                                // a[1][0] a[1][1] a[1][2] a[1][3]

                                // a[2][0] a[2][1] a[2][2] a[2][3]

                                // a[3][0] a[3][1] a[3][2] a[3][3]

    Matrix44 (T a, T b, T c, T d, T e, T f, T g, T h,

              T i, T j, T k, T l, T m, T n, T o, T p);

                                // a b c d

                                // e f g h

                                // i j k l

                                // m n o p

    Matrix44 (Matrix33<T> r, Vec3<T> t);

                                // r r r 0

                                // r r r 0

                                // r r r 0

                                // t t t 1

    //--------------------------------

    // Copy constructor and assignment

    //--------------------------------

    Matrix44 (const Matrix44 &v);

    template <class S> explicit Matrix44 (const Matrix44<S> &v);

    const Matrix44 &    operator = (const Matrix44 &v);

    const Matrix44 &    operator = (T a);

    //----------------------

    // Compatibility with Sb

    //----------------------

    T *                 getValue ();

    const T *           getValue () const;

    template <class S>

    void                getValue (Matrix44<S> &v) const;

    template <class S>

    Matrix44 &          setValue (const Matrix44<S> &v);

    template <class S>

    Matrix44 &          setTheMatrix (const Matrix44<S> &v);

    //---------

    // Identity

    //---------

    void                makeIdentity();

    //---------

    // Equality

    //---------

    bool                operator == (const Matrix44 &v) const;

    bool                operator != (const Matrix44 &v) const;

    //-----------------------------------------------------------------------

    // Compare two matrices and test if they are "approximately equal":

    //

    // equalWithAbsError (m, e)

    //

    //      Returns true if the coefficients of this and m are the same with

    //      an absolute error of no more than e, i.e., for all i, j

    //

    //      abs (this[i][j] - m[i][j]) <= e

    //

    // equalWithRelError (m, e)

    //

    //      Returns true if the coefficients of this and m are the same with

    //      a relative error of no more than e, i.e., for all i, j

    //

    //      abs (this[i] - v[i][j]) <= e * abs (this[i][j])

    //-----------------------------------------------------------------------

    bool                equalWithAbsError (const Matrix44<T> &v, T e) const;

    bool                equalWithRelError (const Matrix44<T> &v, T e) const;

    //------------------------

    // Component-wise addition

    //------------------------

    const Matrix44 &    operator += (const Matrix44 &v);

    const Matrix44 &    operator += (T a);

    Matrix44            operator + (const Matrix44 &v) const;

    //---------------------------

    // Component-wise subtraction

    //---------------------------

    const Matrix44 &    operator -= (const Matrix44 &v);

    const Matrix44 &    operator -= (T a);

    Matrix44            operator - (const Matrix44 &v) const;

    //------------------------------------

    // Component-wise multiplication by -1

    //------------------------------------

    Matrix44            operator - () const;

    const Matrix44 &    negate ();

    //------------------------------

    // Component-wise multiplication

    //------------------------------

    const Matrix44 &    operator *= (T a);

    Matrix44            operator * (T a) const;

    //-----------------------------------

    // Matrix-times-matrix multiplication

    //-----------------------------------

    const Matrix44 &    operator *= (const Matrix44 &v);

    Matrix44            operator * (const Matrix44 &v) const;

    static void         multiply (const Matrix44 &a,    // assumes that

                                  const Matrix44 &b,    // &a != &c and

                                  Matrix44 &c);         // &b != &c.

    //-----------------------------------------------------------------

    // Vector-times-matrix multiplication; see also the "operator *"

    // functions defined below.

    //

    // m.multVecMatrix(src,dst) implements a homogeneous transformation

    // by computing Vec4 (src.x, src.y, src.z, 1) * m and dividing by

    // the result's third element.

    //

    // m.multDirMatrix(src,dst) multiplies src by the upper left 3x3

    // submatrix, ignoring the rest of matrix m.

    //-----------------------------------------------------------------

    template <class S>

    void                multVecMatrix(const Vec3<S> &src, Vec3<S> &dst) const;

    template <class S>

    void                multDirMatrix(const Vec3<S> &src, Vec3<S> &dst) const;

    //------------------------

    // Component-wise division

    //------------------------

    const Matrix44 &    operator /= (T a);

    Matrix44            operator / (T a) const;

    //------------------

    // Transposed matrix

    //------------------

    const Matrix44 &    transpose ();

    Matrix44            transposed () const;

    //------------------------------------------------------------

    // Inverse matrix: If singExc is false, inverting a singular

    // matrix produces an identity matrix.  If singExc is true,

    // inverting a singular matrix throws a SingMatrixExc.

    //

    // inverse() and invert() invert matrices using determinants;

    // gjInverse() and gjInvert() use the Gauss-Jordan method.

    //

    // inverse() and invert() are significantly faster than

    // gjInverse() and gjInvert(), but the results may be slightly

    // less accurate.

    // 

    //------------------------------------------------------------

    const Matrix44 &    invert (bool singExc = false);

    Matrix44<T>         inverse (bool singExc = false) const;

    const Matrix44 &    gjInvert (bool singExc = false);

    Matrix44<T>         gjInverse (bool singExc = false) const;

    //------------------------------------------------

    // Calculate the matrix minor of the (r,c) element

    //------------------------------------------------

    T                   minorOf (const int r, const int c) const;

    //---------------------------------------------------

    // Build a minor using the specified rows and columns

    //---------------------------------------------------

    T                   fastMinor (const int r0, const int r1, const int r2,

                                   const int c0, const int c1, const int c2) const;

    //------------

    // Determinant

    //------------

    T                   determinant() const;

    //--------------------------------------------------------

    // Set matrix to rotation by XYZ euler angles (in radians)

    //--------------------------------------------------------

    template <class S>

    const Matrix44 &    setEulerAngles (const Vec3<S>& r);

    //--------------------------------------------------------

    // Set matrix to rotation around given axis by given angle

    //--------------------------------------------------------

    template <class S>

    const Matrix44 &    setAxisAngle (const Vec3<S>& ax, S ang);

    //-------------------------------------------

    // Rotate the matrix by XYZ euler angles in r

    //-------------------------------------------

    template <class S>

    const Matrix44 &    rotate (const Vec3<S> &r);

    //--------------------------------------------

    // Set matrix to scale by given uniform factor

    //--------------------------------------------

    const Matrix44 &    setScale (T s);

    //------------------------------------

    // Set matrix to scale by given vector

    //------------------------------------

    template <class S>

    const Matrix44 &    setScale (const Vec3<S> &s);

    //----------------------

    // Scale the matrix by s

    //----------------------

    template <class S>

    const Matrix44 &    scale (const Vec3<S> &s);

    //------------------------------------------

    // Set matrix to translation by given vector

    //------------------------------------------

    template <class S>

    const Matrix44 &    setTranslation (const Vec3<S> &t);

    //-----------------------------

    // Return translation component

    //-----------------------------

    const Vec3<T>       translation () const;

    //--------------------------

    // Translate the matrix by t

    //--------------------------

    template <class S>

    const Matrix44 &    translate (const Vec3<S> &t);

    //-------------------------------------------------------------

    // Set matrix to shear by given vector h.  The resulting matrix

    //    will shear x for each y coord. by a factor of h[0] ;

    //    will shear x for each z coord. by a factor of h[1] ;

    //    will shear y for each z coord. by a factor of h[2] .

    //-------------------------------------------------------------

    template <class S>

    const Matrix44 &    setShear (const Vec3<S> &h);

    //------------------------------------------------------------

    // Set matrix to shear by given factors.  The resulting matrix

    //    will shear x for each y coord. by a factor of h.xy ;

    //    will shear x for each z coord. by a factor of h.xz ;

    //    will shear y for each z coord. by a factor of h.yz ; 

    //    will shear y for each x coord. by a factor of h.yx ;

    //    will shear z for each x coord. by a factor of h.zx ;

    //    will shear z for each y coord. by a factor of h.zy .

    //------------------------------------------------------------

    template <class S>

    const Matrix44 &    setShear (const Shear6<S> &h);

    //--------------------------------------------------------

    // Shear the matrix by given vector.  The composed matrix 

    // will be <shear> * <this>, where the shear matrix ...

    //    will shear x for each y coord. by a factor of h[0] ;

    //    will shear x for each z coord. by a factor of h[1] ;

    //    will shear y for each z coord. by a factor of h[2] .

    //--------------------------------------------------------

    template <class S>

    const Matrix44 &    shear (const Vec3<S> &h);

    //--------------------------------------------------------

    // Number of the row and column dimensions, since

    // Matrix44 is a square matrix.

    //--------------------------------------------------------

    static unsigned int dimensions() {return 4;}

    //------------------------------------------------------------

    // Shear the matrix by the given factors.  The composed matrix 

    // will be <shear> * <this>, where the shear matrix ...

    //    will shear x for each y coord. by a factor of h.xy ;

    //    will shear x for each z coord. by a factor of h.xz ;

    //    will shear y for each z coord. by a factor of h.yz ;

    //    will shear y for each x coord. by a factor of h.yx ;

    //    will shear z for each x coord. by a factor of h.zx ;

    //    will shear z for each y coord. by a factor of h.zy .

    //------------------------------------------------------------

    template <class S>

    const Matrix44 &    shear (const Shear6<S> &h);

    //-------------------------------------------------

    // Limitations of type T (see also class limits<T>)

    //-------------------------------------------------

    static T            baseTypeMin()           {return limits<T>::min();}

    static T            baseTypeMax()           {return limits<T>::max();}

    static T            baseTypeSmallest()      {return limits<T>::smallest();}

    static T            baseTypeEpsilon()       {return limits<T>::epsilon();}

    typedef T   BaseType;

    typedef Vec4<T> BaseVecType;

  private:

    template <typename R, typename S>

    struct isSameType

    {

        enum {value = 0};

    };

    template <typename R>

    struct isSameType<R, R>

    {

        enum {value = 1};

    };

};

//--------------

// Stream output

//--------------

template <class T>

std::ostream &  operator << (std::ostream & s, const Matrix33<T> &m); 

template <class T>

std::ostream &  operator << (std::ostream & s, const Matrix44<T> &m); 

//---------------------------------------------

// Vector-times-matrix multiplication operators

//---------------------------------------------

template <class S, class T>

const Vec2<S> &            operator *= (Vec2<S> &v, const Matrix33<T> &m);

template <class S, class T>

Vec2<S>                    operator * (const Vec2<S> &v, const Matrix33<T> &m);

template <class S, class T>

const Vec3<S> &            operator *= (Vec3<S> &v, const Matrix33<T> &m);

template <class S, class T>

Vec3<S>                    operator * (const Vec3<S> &v, const Matrix33<T> &m);

template <class S, class T>

const Vec3<S> &            operator *= (Vec3<S> &v, const Matrix44<T> &m);

template <class S, class T>

Vec3<S>                    operator * (const Vec3<S> &v, const Matrix44<T> &m);

template <class S, class T>

const Vec4<S> &            operator *= (Vec4<S> &v, const Matrix44<T> &m);

template <class S, class T>

Vec4<S>                    operator * (const Vec4<S> &v, const Matrix44<T> &m);

//-------------------------

// Typedefs for convenience

//-------------------------

typedef Matrix33 <float>  M33f;

typedef Matrix33 <double> M33d;

typedef Matrix44 <float>  M44f;

typedef Matrix44 <double> M44d;

//---------------------------

// Implementation of Matrix33

//---------------------------

template <class T>

inline T *

Matrix33<T>::operator [] (int i)

{

    return x[i];

}

template <class T>

inline const T *

Matrix33<T>::operator [] (int i) const

{

    return x[i];

}

template <class T>

inline

Matrix33<T>::Matrix33 ()

{

    memset (x, 0, sizeof (x));

    x[0][0] = 1;

    x[1][1] = 1;

    x[2][2] = 1;

}

template <class T>

inline

Matrix33<T>::Matrix33 (T a)

{

    x[0][0] = a;

    x[0][1] = a;

    x[0][2] = a;

    x[1][0] = a;

    x[1][1] = a;

    x[1][2] = a;

    x[2][0] = a;

    x[2][1] = a;

    x[2][2] = a;

}

template <class T>

inline

Matrix33<T>::Matrix33 (const T a[3][3]) 

{

    memcpy (x, a, sizeof (x));

}

template <class T>

inline

Matrix33<T>::Matrix33 (T a, T b, T c, T d, T e, T f, T g, T h, T i)

{

    x[0][0] = a;

    x[0][1] = b;

    x[0][2] = c;

    x[1][0] = d;

    x[1][1] = e;

    x[1][2] = f;

    x[2][0] = g;

    x[2][1] = h;

    x[2][2] = i;

}

template <class T>

inline

Matrix33<T>::Matrix33 (const Matrix33 &v)

{

    memcpy (x, v.x, sizeof (x));

}

template <class T>

template <class S>

inline

Matrix33<T>::Matrix33 (const Matrix33<S> &v)

{

    x[0][0] = T (v.x[0][0]);

    x[0][1] = T (v.x[0][1]);

    x[0][2] = T (v.x[0][2]);

    x[1][0] = T (v.x[1][0]);

    x[1][1] = T (v.x[1][1]);

    x[1][2] = T (v.x[1][2]);

    x[2][0] = T (v.x[2][0]);

    x[2][1] = T (v.x[2][1]);

    x[2][2] = T (v.x[2][2]);

}

template <class T>

inline const Matrix33<T> &

Matrix33<T>::operator = (const Matrix33 &v)

{

    memcpy (x, v.x, sizeof (x));

    return *this;

}

template <class T>

inline const Matrix33<T> &

Matrix33<T>::operator = (T a)

{

    x[0][0] = a;

    x[0][1] = a;

    x[0][2] = a;

    x[1][0] = a;

    x[1][1] = a;

    x[1][2] = a;

    x[2][0] = a;

    x[2][1] = a;

    x[2][2] = a;

    return *this;

}

template <class T>

inline T *

Matrix33<T>::getValue ()

{

    return (T *) &x[0][0];

}

template <class T>

inline const T *

Matrix33<T>::getValue () const

{

    return (const T *) &x[0][0];

}

template <class T>

template <class S>

inline void

Matrix33<T>::getValue (Matrix33<S> &v) const

{

    if (isSameType<S,T>::value)

    {

        memcpy (v.x, x, sizeof (x));

    }

    else

    {

        v.x[0][0] = x[0][0];

        v.x[0][1] = x[0][1];

        v.x[0][2] = x[0][2];

        v.x[1][0] = x[1][0];

        v.x[1][1] = x[1][1];

        v.x[1][2] = x[1][2];

        v.x[2][0] = x[2][0];

        v.x[2][1] = x[2][1];

        v.x[2][2] = x[2][2];

    }

}

template <class T>

template <class S>

inline Matrix33<T> &

Matrix33<T>::setValue (const Matrix33<S> &v)

{

    if (isSameType<S,T>::value)

    {

        memcpy (x, v.x, sizeof (x));

    }

    else

    {

        x[0][0] = v.x[0][0];

        x[0][1] = v.x[0][1];

        x[0][2] = v.x[0][2];

        x[1][0] = v.x[1][0];

        x[1][1] = v.x[1][1];

        x[1][2] = v.x[1][2];

        x[2][0] = v.x[2][0];

        x[2][1] = v.x[2][1];

        x[2][2] = v.x[2][2];

    }

    return *this;

}

template <class T>

template <class S>

inline Matrix33<T> &

Matrix33<T>::setTheMatrix (const Matrix33<S> &v)

{

    if (isSameType<S,T>::value)

    {

        memcpy (x, v.x, sizeof (x));

    }

    else

    {

        x[0][0] = v.x[0][0];

        x[0][1] = v.x[0][1];

        x[0][2] = v.x[0][2];

        x[1][0] = v.x[1][0];

        x[1][1] = v.x[1][1];

        x[1][2] = v.x[1][2];

        x[2][0] = v.x[2][0];

        x[2][1] = v.x[2][1];

        x[2][2] = v.x[2][2];

    }

    return *this;

}

template <class T>

inline void

Matrix33<T>::makeIdentity()

{

    memset (x, 0, sizeof (x));

    x[0][0] = 1;

    x[1][1] = 1;

    x[2][2] = 1;

}

template <class T>

bool

Matrix33<T>::operator == (const Matrix33 &v) const

{

    return x[0][0] == v.x[0][0] &&

           x[0][1] == v.x[0][1] &&

           x[0][2] == v.x[0][2] &&

           x[1][0] == v.x[1][0] &&

           x[1][1] == v.x[1][1] &&

           x[1][2] == v.x[1][2] &&

           x[2][0] == v.x[2][0] &&

           x[2][1] == v.x[2][1] &&

           x[2][2] == v.x[2][2];

}

template <class T>

bool

Matrix33<T>::operator != (const Matrix33 &v) const

{

    return x[0][0] != v.x[0][0] ||

           x[0][1] != v.x[0][1] ||

           x[0][2] != v.x[0][2] ||

           x[1][0] != v.x[1][0] ||

           x[1][1] != v.x[1][1] ||

           x[1][2] != v.x[1][2] ||

           x[2][0] != v.x[2][0] ||

           x[2][1] != v.x[2][1] ||

           x[2][2] != v.x[2][2];

}

template <class T>

bool

Matrix33<T>::equalWithAbsError (const Matrix33<T> &m, T e) const

{

    for (int i = 0; i < 3; i++)

        for (int j = 0; j < 3; j++)

            if (!IMATH_INTERNAL_NAMESPACE::equalWithAbsError ((*this)[i][j], m[i][j], e))

                return false;

    return true;

}

template <class T>

bool

Matrix33<T>::equalWithRelError (const Matrix33<T> &m, T e) const

{

    for (int i = 0; i < 3; i++)

        for (int j = 0; j < 3; j++)

            if (!IMATH_INTERNAL_NAMESPACE::equalWithRelError ((*this)[i][j], m[i][j], e))

                return false;

    return true;

}

template <class T>

const Matrix33<T> &

Matrix33<T>::operator += (const Matrix33<T> &v)

{

    x[0][0] += v.x[0][0];

    x[0][1] += v.x[0][1];

    x[0][2] += v.x[0][2];

    x[1][0] += v.x[1][0];

    x[1][1] += v.x[1][1];