diff --git a/NIF_Blocks.cpp b/NIF_Blocks.cpp
index b41a938e769d90944988654e1b2f8f7e54bd792d..dba32f1572b447ae9790d543b14a190a6805eebe 100644
--- a/NIF_Blocks.cpp
+++ b/NIF_Blocks.cpp
@@ -73,26 +73,11 @@ void const * ANode::QueryInterface( int id ) const {
 
 Matrix44 ANode::GetLocalTransform() const {
 	//Get transform data from atributes
-	Matrix33 f = GetAttr("Rotation")->asMatrix33();
+	Matrix33 rot = GetAttr("Rotation")->asMatrix33();
 	Float3 tran = GetAttr("Translation")->asFloat3();
 	float scale = GetAttr("Scale")->asFloat();
 
-	//Set up a matrix with rotate and translate information
-	Matrix44 rt;
-	rt[0][0] = f[0][0];	rt[0][1] = f[0][1];	rt[0][2] = f[0][2];	rt[0][3] = 0.0f;
-	rt[1][0] = f[1][0];	rt[1][1] = f[1][1];	rt[1][2] = f[1][2];	rt[1][3] = 0.0f;
-	rt[2][0] = f[2][0];	rt[2][1] = f[2][1];	rt[2][2] = f[2][2];	rt[2][3] = 0.0f;
-	rt[3][0] = tran[0];	rt[3][1] = tran[1];	rt[3][2] = tran[2];	rt[3][3] = 1.0f;
-
-	//Set up another matrix with the scale information
-	Matrix44 s;
-	s[0][0] = scale;	s[0][1] = 0.0f;		s[0][2] = 0.0f;		s[0][3] = 0.0f;
-	s[1][0] = 0.0f;		s[1][1] = scale;	s[1][2] = 0.0f;		s[1][3] = 0.0f;
-	s[2][0] = 0.0f;		s[2][1] = 0.0f;		s[2][2] = scale;	s[2][3] = 0.0f;
-	s[3][0] = 0.0f;		s[3][1] = 0.0f;		s[3][2] = 0.0f;		s[3][3] = 1.0f;
-
-	//Multiply the two for the resulting local transform
-	return MultMatrix44(s, rt);
+	return Matrix44( Vector3( tran[0], tran[1], tran[2] ), rot, scale );
 }
 
 Matrix44 ANode::GetWorldTransform() const {
@@ -107,7 +92,7 @@ Matrix44 ANode::GetWorldTransform() const {
 		Matrix44 par_world = node->GetWorldTransform();
 
 		//Multipy local matrix and parent world matrix for result
-		return MultMatrix44( par_world, local);
+		return par_world * local;
 	}
 	else {
 		//No parent transform, simply return local transform
@@ -132,9 +117,9 @@ Matrix44 ANode::GetLocalBindPos() const {
 		//There is a node parent
 		//multiply its inverse with this block's bind position to get the local bind position
 		Matrix44 par_mat = node->GetWorldBindPos();
-		Matrix44 par_inv = InverseMatrix44( par_mat);
+		Matrix44 par_inv = par_mat.Inverse();
 		
-		return MultMatrix44( bindPosition, par_inv);
+		return bindPosition * par_inv;
 	}
 	else {
 		//No parent transform, simply return local transform
@@ -2492,8 +2477,8 @@ void NiSkinData::StraightenSkeleton() {
 				);
 
 				//Do calculation to get correct bone postion in relation to parent
-				Matrix44 inverse_co = InverseMatrix44(child_offset);
-				Matrix44 child_pos = MultMatrix44( inverse_co, parent_offset);
+				Matrix44 inverse_co = child_offset.Inverse();
+				Matrix44 child_pos = inverse_co * parent_offset;
 
 				//Store result in block's Bind Position Matrix
 				INode * node = (INode*)it2->first->QueryInterface(ID_NODE);
@@ -2566,7 +2551,7 @@ void NiSkinData::RepositionTriShape( blk_ref & tri_shape ) {
 
 		Matrix44 bone_mat = bone_node->GetWorldBindPos();
 
-		Matrix44 result_mat = MultMatrix44( offset_mat, bone_mat);
+		Matrix44 result_mat = offset_mat * bone_mat;
 
 		//GetBuiltUpTransform( bone_blk, end_mat );
 
@@ -2760,8 +2745,8 @@ void NiSkinData::CalculateBoneOffset( INode const * const par_node, IBlock const
 	bone_mat = bone_node->GetWorldBindPos();
 
 	//Inverse bone matrix & multiply with parent node matrix
-	inv_mat = InverseMatrix44(bone_mat);
-	res_mat = MultMatrix44(par_mat, inv_mat);
+	inv_mat = bone_mat.Inverse();
+	res_mat = par_mat * inv_mat;
 
 	//--Extract Scale from first 3 rows--//
 	float scale[3];
@@ -2812,8 +2797,8 @@ void NiSkinData::CalculateOverallOffset( Matrix33 & rot, fVector3 & tr, float &
 	Matrix44 skel_mat = iskel->GetWorldTransform();
 	
 	// Inverse parent node transform & multiply with skeleton matrix
-	Matrix44 inv_mat = InverseMatrix44(par_mat);
-	Matrix44 res_mat = MultMatrix44(inv_mat, skel_mat);
+	Matrix44 inv_mat = par_mat.Inverse();
+	Matrix44 res_mat = inv_mat * skel_mat;
 
 	//--Extract Scale from first 3 rows--//
 	float scale[3];
diff --git a/NIF_Blocks.h b/NIF_Blocks.h
index 5583a82ce70bff802448abdf863ab70c9966b7b9..ae33be99f3ff265509b67330540fb02c3649e50c 100644
--- a/NIF_Blocks.h
+++ b/NIF_Blocks.h
@@ -75,7 +75,7 @@ public:
 	ANode();
 	void Init() { 
 		//Start the bind pose at an identity matrix
-		SetIdentity44(bindPosition);
+		bindPosition = Matrix44::IDENTITY;
 
 		//Start the flags at "Not a skin influence"
 		GetAttr("Flags")->Set(8);
@@ -1887,7 +1887,7 @@ class NiSkinData : public AData, public ISkinData, public ISkinDataInternal {
 
 		NiSkinData() { 
 			//AddAttr( attr_link, "Skin Partition", 0, VER_10_1_0_0 );
-			SetIdentity33(rotation);
+			rotation = Matrix33::IDENTITY;
 			translation[0] = 0.0f;
 			translation[1] = 0.0f;
 			translation[2] = 0.0f;
diff --git a/nif_math.cpp b/nif_math.cpp
index 583aec2c9ac7572c88d832807cd0d44224730fd2..94f6b7b3c404c5371507182b7a6b08c761c637e9 100644
--- a/nif_math.cpp
+++ b/nif_math.cpp
@@ -34,137 +34,19 @@ POSSIBILITY OF SUCH DAMAGE. */
 #include "nif_math.h"
 #include <iomanip>
 
-void PrintMatrix33( Matrix33 const & m, ostream & out ) {
-	out << endl
-		<< "      |" << setw(8) << m[0][0] << "," << setw(8) << m[0][1] << "," << setw(8) << m[0][2] << " |" << endl
-		<< "      |" << setw(8) << m[1][0] << "," << setw(8) << m[1][1] << "," << setw(8) << m[1][2] << " |" << endl
-		<< "      |" << setw(8) << m[2][0] << "," << setw(8) << m[2][1] << "," << setw(8) << m[2][2] << " |" <<endl;
-}
-
-void PrintMatrix44( Matrix44 const & m, ostream & out ) {
-	out << endl
-		<< "      |" << setw(8) << m[0][0] << "," << setw(8) << m[0][1] << "," << setw(8) << m[0][2] << setw(8) << m[0][3] << " |" << endl
-		<< "      |" << setw(8) << m[1][0] << "," << setw(8) << m[1][1] << "," << setw(8) << m[1][2] << setw(8) << m[1][3] << " |" << endl
-		<< "      |" << setw(8) << m[2][0] << "," << setw(8) << m[2][1] << "," << setw(8) << m[2][2] << setw(8) << m[2][3] << " |" << endl
-		<< "      |" << setw(8) << m[3][0] << "," << setw(8) << m[3][1] << "," << setw(8) << m[3][2] << setw(8) << m[3][3] << " |" << endl;
-}
-
-Vector3 MultVector3( Vector3 const & a, Vector3 const & b ) {
-	Vector3 answer;
-	answer.x = a.y * b.z - a.z * b.y;
-	answer.y = a.z * b.x - a.x * b.z;
-	answer.z = a.x * b.y - a.y * b.x;
-
-	return answer;
-}
-
-void SetIdentity33( Matrix33 & m ) {
-	m[0][0] = 1.0f;	m[0][1] = 0.0f;	m[0][2] = 0.0f;
-	m[1][0] = 0.0f;	m[1][1] = 1.0f;	m[1][2] = 0.0f;
-	m[2][0] = 0.0f;	m[2][1] = 0.0f;	m[2][2] = 1.0f;
-}
-
-void SetIdentity44( Matrix44 & m ) {
-	m[0][0] = 1.0f;	m[0][1] = 0.0f;	m[0][2] = 0.0f;	m[0][3] = 0.0f;
-	m[1][0] = 0.0f;	m[1][1] = 1.0f;	m[1][2] = 0.0f;	m[1][3] = 0.0f;
-	m[2][0] = 0.0f;	m[2][1] = 0.0f;	m[2][2] = 1.0f;	m[2][3] = 0.0f;
-	m[3][0] = 0.0f;	m[3][1] = 0.0f;	m[3][2] = 0.0f;	m[3][3] = 1.0f;
-}
-
-Matrix33 MultMatrix33( Matrix33 const & a, Matrix33 const & b ) {
-	Matrix33 result;
-	for (int i = 0; i < 3; i++) {
-		for (int j = 0; j < 3; j++) {
-			float t = 0.0f;
-			for (int k = 0; k < 4; k++) {
-				t += a[i][k] * b[k][j];
-			}
-			result[i][j] = t;
-		}
-	}
-	return result;
-}
-
-Matrix44 MultMatrix44( Matrix44 const & a, Matrix44 const & b ) {
-	Matrix44 result;
-	for (int i = 0; i < 4; i++) {
-		for (int j = 0; j < 4; j++) {
-			float t = 0.0f;
-			for (int k = 0; k < 4; k++) {
-				t += a[i][k] * b[k][j];
-			}
-			result[i][j] = t;
-		}
-	}
-	return result;
-}
-
-float DetMatrix33( Matrix33 const & m ) {
-	return  m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])
-		  - m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])
-		  + m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
-}
-
-float DetMatrix44( Matrix44 const & m ) {
-	Matrix33 sub1(
-		m[1][1], m[1][2], m[1][3],
-		m[2][1], m[2][2], m[2][3],
-		m[3][1], m[3][2], m[3][3]
-	);
-
-	Matrix33 sub2(
-		m[1][0], m[1][2], m[1][3],
-		m[2][0], m[2][2], m[2][3],
-		m[3][0], m[3][2], m[3][3]
-	);
-
-	Matrix33 sub3(
-		m[1][0], m[1][1], m[1][3],
-		m[2][0], m[2][1], m[2][3],
-		m[3][0], m[3][1], m[3][3] 
-	);
-
-	Matrix33 sub4(
-		m[1][0], m[1][1], m[1][2],
-		m[2][0], m[2][1], m[2][2],
-		m[3][0], m[3][1], m[3][2]
-	);
-
-	return  m[0][0] * DetMatrix33( sub1 )
-	      - m[0][1] * DetMatrix33( sub2 )
-	      + m[0][2] * DetMatrix33( sub3 )
-	      - m[0][3] * DetMatrix33( sub4 );
-}
+//Constants
 
-float AdjMatrix44(Matrix44 const & m, int skip_r, int skip_c) {
-	Matrix33 sub;
-	int i = 0, j = 0;
-	for (int r = 0; r < 4; r++) {
-		if (r == skip_c)
-			continue;
-		for (int c = 0; c < 4; c++) {
-			if (c == skip_r)
-				continue;
-			sub[i][j] = m[r][c];
-			j++;
-		}
-		i++;
-		j = 0;
-	}
+const Matrix44 Matrix44::IDENTITY( 1.0f, 0.0f, 0.0f, 0.0f,
+								   0.0f, 1.0f, 0.0f, 0.0f,
+								   0.0f, 0.0f, 1.0f, 0.0f,
+								   0.0f, 0.0f, 0.0f, 1.0f );
 
-	return pow(-1.0f, float(skip_r + skip_c)) * DetMatrix33( sub );
-}
+const Matrix33 Matrix33::IDENTITY( 1.0f, 0.0f, 0.0f,
+								   0.0f, 1.0f, 0.0f,
+								   0.0f, 0.0f, 1.0f );
 
-Matrix44 InverseMatrix44( Matrix44 const & m ) {
-	Matrix44 result;
-	float det = DetMatrix44(m);
-	for (int r = 0; r < 4; r++) {
-		for (int c = 0; c < 4; c++) {
-			result[r][c] = AdjMatrix44(m, r, c) / det;
-		}
-	}
-	return result;
-}
+const Matrix22 Matrix22::IDENTITY( 1.0f, 0.0f,
+								   0.0f, 1.0f );
 
 /*
  * Vector3 Methods
@@ -268,10 +150,22 @@ Vector3 Vector3::CrossProduct( const Vector3 & rh) const {
 //	return *this;
 //}
 
+/*
+ * Matrix22 Methods
+ */
+
+Matrix22::Matrix22() {
+	*this = Matrix22::IDENTITY;
+}
+
 /*
  * Matrix33 Methods
  */
 
+Matrix33::Matrix33() {
+	*this = Matrix33::IDENTITY;
+}
+
 Quaternion Matrix33::AsQuaternion() {
 	Quaternion quat;
 	float tr, s, q[4];
@@ -318,12 +212,145 @@ Quaternion Matrix33::AsQuaternion() {
 	return quat;
 }
 
+float Matrix33::Determinant() const {
+	return  (*this)[0][0] * ( (*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1] )
+		  - (*this)[0][1] * ( (*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0] )
+		  + (*this)[0][2] * ( (*this)[1][0] * (*this)[2][1] - (*this)[1][1] * (*this)[2][0] );
+}
+
+
+
 /*
  * Matrix44 Methods
  */
 
 Matrix44::Matrix44() {
-	*this = IDENTITY44;
+	*this = Matrix44::IDENTITY;
+}
+
+Matrix44::Matrix44( const Vector3 & t, const Matrix33 & r, float scale ) {
+	//Set up a matrix with rotate and translate information
+	Matrix44 rt;
+	rt[0][0] = r[0][0];	rt[0][1] = r[0][1];	rt[0][2] = r[0][2];	rt[0][3] = 0.0f;
+	rt[1][0] = r[1][0];	rt[1][1] = r[1][1];	rt[1][2] = r[1][2];	rt[1][3] = 0.0f;
+	rt[2][0] = r[2][0];	rt[2][1] = r[2][1];	rt[2][2] = r[2][2];	rt[2][3] = 0.0f;
+	rt[3][0] = t.x;     rt[3][1] = t.y;     rt[3][2] = t.z;     rt[3][3] = 1.0f;
+
+	//Set up another matrix with the scale information
+	Matrix44 s;
+	s[0][0] = scale;	s[0][1] = 0.0f;		s[0][2] = 0.0f;		s[0][3] = 0.0f;
+	s[1][0] = 0.0f;		s[1][1] = scale;	s[1][2] = 0.0f;		s[1][3] = 0.0f;
+	s[2][0] = 0.0f;		s[2][1] = 0.0f;		s[2][2] = scale;	s[2][3] = 0.0f;
+	s[3][0] = 0.0f;		s[3][1] = 0.0f;		s[3][2] = 0.0f;		s[3][3] = 1.0f;
+
+	//Multiply the two for the combined transform
+	*this = s * rt;
+}
+
+Matrix44 Matrix44::operator*( const Matrix44 & rh ) const {
+	return Matrix44(*this) *= rh;
+}
+Matrix44 & Matrix44::operator*=( const Matrix44 & rh ) {
+	Matrix44 r;
+	Matrix44 & lh = *this;
+	float t;
+	for (int i = 0; i < 4; i++) {
+		for (int j = 0; j < 4; j++) {
+			t = 0.0f;
+			for (int k = 0; k < 4; k++) {
+				t += lh[i][k] * rh[k][j];
+			}
+			r[i][j] = t;
+		}
+	}
+
+	*this = r;
+	return *this;
+}
+
+Matrix44 Matrix44::operator*( float rh ) const {
+	return Matrix44(*this) *= rh;
+}
+
+Matrix44 & Matrix44::operator*=( float rh ) {
+	for (int i = 0; i < 4; i++) {
+		for (int j = 0; j < 4; j++) {
+			(*this)[i][j] *= rh;
+		}
+	}
+	return *this;
+}
+
+Matrix44 Matrix44::operator+( const Matrix44 & rh ) const {
+	return Matrix44(*this) += rh;
+} 
+
+Matrix44 & Matrix44::operator+=( const Matrix44 & rh ) {
+	for (int i = 0; i < 4; i++) {
+		for (int j = 0; j < 4; j++) {
+			(*this)[i][j] += rh[i][j];
+		}
+	}
+	return *this;
+}
+
+Matrix44 & Matrix44::operator=( const Matrix44 & rh ) {
+	memcpy(rows, rh.rows, sizeof(Float4) * 4);
+	return *this;
+}
+
+bool Matrix44::operator==( const Matrix44 & rh ) const {
+	for (int i = 0; i < 4; i++) {
+		for (int j = 0; j < 4; j++) {
+			if ( (*this)[i][j] != rh[i][j] )
+				return false;
+		}
+	}
+	return true;
+}
+
+Matrix44 Matrix44::Transpose() const {
+	const Matrix44 & t = *this;
+	return Matrix44( t[0][0], t[0][1], t[0][2], t[0][3],
+					 t[1][0], t[1][1], t[1][2], t[1][3],
+					 t[2][0], t[2][1], t[2][2], t[2][3],
+					 t[3][0], t[3][1], t[3][2], t[3][3] );
+}
+
+Matrix33 Matrix44::Submatrix( int skip_r, int skip_c ) const {
+	Matrix33 sub;
+	int i = 0, j = 0;
+	for (int r = 0; r < 4; r++) {
+		if (r == skip_c)
+			continue;
+		for (int c = 0; c < 4; c++) {
+			if (c == skip_r)
+				continue;
+			sub[i][j] = (*this)[r][c];
+			j++;
+		}
+		i++;
+		j = 0;
+	}
+	return sub;
+}
+
+float Matrix44::Adjunct( int skip_r, int skip_c ) const {
+	Matrix33 sub = Submatrix( skip_r, skip_c );
+	return pow(-1.0f, float(skip_r + skip_c)) * sub.Determinant();
+}
+
+Matrix44 Matrix44::Inverse() const {
+	Matrix44 result;
+
+	float det = Determinant();
+	for (int r = 0; r < 4; r++) {
+		for (int c = 0; c < 4; c++) {
+			result[r][c] = Adjunct(r, c) / det;
+		}
+	}
+
+	return result;
 }
 
 /*
diff --git a/nif_math.h b/nif_math.h
index 7642c4e58c5ba269554f28e0954b05f8991a1651..81f5a1c3ce2efefcb02e986d61592298a0814b1e 100644
--- a/nif_math.h
+++ b/nif_math.h
@@ -318,6 +318,9 @@ struct Float2 {
 
 /*! Stores a 2 by 2 matrix used for bump maps. */
 struct Matrix22 {
+	/*! The 2x2 identity matrix */
+	static const Matrix22 IDENTITY;
+
 	Float2 rows[2];  /*!< The two rows of Float2 structures which hold two floating point numbers each. */ 
 	
 	/*! The bracket operator makes it possible to use this structure like a 2x2 C++ array.
@@ -332,7 +335,7 @@ struct Matrix22 {
 	}
 
 	/*! Default Constructor */
-	Matrix22() {}
+	Matrix22();
 
 	/*! This constructor can be used to set all values in this matrix during initialization
 	 * \param m11 The value to set at row 1, column 1.
@@ -427,6 +430,9 @@ struct Float3 {
 
 /*! Stores a 3 by 3 matrix used for rotation. */
 struct Matrix33 {
+	/*! The 3x3 identity matrix */
+	static const Matrix33 IDENTITY;
+
 	Float3 rows[3]; /*!< The three rows of Float3 structures which hold three floating point numbers each. */ 
 	
 	/*! The bracket operator makes it possible to use this structure like a 3x3 C++ array.
@@ -440,8 +446,8 @@ struct Matrix33 {
 		return rows[n];
 	}
 
-	/*! Default constructor. */
-	Matrix33() {}
+	/*! Default constructor.   Initializes matrix to identity.  */
+	Matrix33();
 
 	/*! This constructor can be used to set all values in this matrix during initialization
 	 * \param m11 The value to set at row 1, column 1.
@@ -485,8 +491,17 @@ struct Matrix33 {
 		rows[2][0] = m31; rows[2][1] = m32; rows[2][2] = m33;
 	}
 
+	/*! Returns a quaternion representation of the rotation stored in this matrix. 
+	 * \return A quaternion with an equivalent rotation to the one stored in this matrix.
+	 */
 	Quaternion AsQuaternion();
 
+	/*! Calculates the determinant of this matrix.
+	 * \return The determinant of this matrix.
+	 */
+	float Determinant() const;
+
+	//Undocumented
 	void AsFloatArr( float out[3][3] ) {
 		out[0][0] = rows[0][0]; out[0][1] = rows[0][1]; out[0][2] = rows[0][2];
 		out[1][0] = rows[1][0]; out[1][1] = rows[1][1]; out[1][2] = rows[1][2];
@@ -560,6 +575,9 @@ struct Float4 {
 
 /*! Stores a 4 by 4 matrix used for combined transformations. */
 struct Matrix44 {
+	/*! The 4x4 identity matrix */
+	static const Matrix44 IDENTITY;
+	
 	Float4 rows[4]; /*!< The three rows of Float3 structures which hold three floating point numbers each. */ 
 	
 	/*! The bracket operator makes it possible to use this structure like a 4x4 C++ array.
@@ -611,6 +629,14 @@ struct Matrix44 {
 		rows[3][0] = m41; rows[3][1] = m42; rows[3][2] = m43; rows[3][3] = m44;
 	}
 
+	/*! This constructor allows a 4x4 transform matrix to be initalized from a
+	 * translate vector, a 3x3 rotation matrix, and a scale factor.
+	 * \param translate The translation vector that specifies the new x, y, and z coordinates.
+	 * \param rotate The 3x3 rotation matrix.
+	 * \param scale The scale factor.
+	 */
+	Matrix44( const Vector3 & translate, const Matrix33 & rotation, float scale );
+
 	/*! This function can be used to set all values in this matrix at the same time.
 	 * \param m11 The value to set at row 1, column 1.
 	 * \param m12 The value to set at row 1, column 2.
@@ -641,6 +667,96 @@ struct Matrix44 {
 		rows[3][0] = m41; rows[3][1] = m42; rows[3][2] = m43; rows[3][3] = m44;
 	}
 
+	/* Multiplies this matrix by another.
+	 * \param rh The matrix to multiply this one with.
+	 * \return The result of the multiplication.
+	 */
+	Matrix44 operator*( const Matrix44 & rh ) const;
+
+	/* Multiplies this matrix by another and sets the result to itself.
+	 * \param rh The matrix to multiply this one with.
+	 * \return This matrix is returned.
+	 */
+	Matrix44 & operator*=( const Matrix44 & rh );
+
+	/* Multiplies this matrix by a scalar value.
+	 * \param rh The scalar value to multiply each component of this matrix by.
+	 * \return The result of the multiplication.
+	 */
+	Matrix44 operator*( float rh ) const;
+
+	/* Multiplies this matrix by a scalar value and sets the resutl to itself.
+	 * \param rh The scalar value to multiply each component of this matrix by.
+	 * \return This matrix is returned.
+	 */
+	Matrix44 & operator*=( float rh );
+
+	/* Multiplies this matrix by a vector with x, y, and z components.
+	 * \param rh The vector to multiply this matrix with.
+	 * \return The result of the multiplication.
+	 */
+	Vector3 operator*( const Vector3 & rh ) const;
+
+	/* Adds this matrix to another.
+	 * \param rh The matrix to be added to this one.
+	 * \return The result of the addition.
+	 */
+	Matrix44 operator+( const Matrix44 & rh ) const;
+
+	/* Adds this matrix to another and sets the result to itself.
+	 * \param rh The matrix to be added to this one.
+	 * \return This matrix is returned.
+	 */
+	Matrix44 & operator+=( const Matrix44 & rh );
+
+	/* Sets the values of this matrix to those of the given matrix.
+	 * \param rh The matrix to copy values from.
+	 * \return This matrix is returned.
+	 */
+	Matrix44 & operator=( const Matrix44 & rh );
+
+	/* Compares two 4x4 matricies.  They are considered equal if all components are equal.
+	 * \param rh The matrix to compare this one with.
+	 * \return true if the matricies are equal, false otherwise.
+	 */
+	bool operator==( const Matrix44 & rh ) const;
+
+	/* Allows the contents of the matrix to be printed to an ostream.
+	 * \param lh The ostream to insert the text into.
+	 * \param rh The matrix to insert into the stream.
+	 * \return The given ostream is returned.
+	 */
+	friend ostream & operator<<( ostream & lh, const Matrix44 & rh );
+
+	/*! Calculates the transpose of this matrix.
+	 * \return The transpose of this matrix.
+	 */
+	Matrix44 Transpose() const;
+
+	/*! Calculates the determinant of this matrix.
+	 * \return The determinant of this matrix.
+	 */
+	float Determinant() const;
+
+	/*! Calculates the inverse of this matrix.
+	 * \retun The inverse of this matrix.
+	 */
+	Matrix44 Inverse() const;
+
+	/*! Returns a 3x3 submatrix of this matrix created by skipping the indicated row and column.
+	 * \param skip_r The row to skip.  Must be a value between 0 and 3.
+	 * \param skip_c The colum to skip.  Must be a value between 0 and 3.
+	 * \return The 3x3 submatrix obtained by skipping the indicated row and column.
+	 */
+	Matrix33 Submatrix( int skip_r, int skip_c ) const;
+
+	/*! Calculates the adjunct of this matrix created by skipping the indicated row and column.
+	 * \param skip_r The row to skip.  Must be a value between 0 and 3.
+	 * \param skip_c The colum to skip.  Must be a value between 0 and 3.
+	 * \return The adjunct obtained by skipping the indicated row and column.
+	 */
+	float Adjunct( int skip_r, int skip_c ) const;
+
 	//undocumented
 	void AsFloatArr( float out[4][4] ) {
 		out[0][0] = rows[0][0]; out[0][1] = rows[0][1]; out[0][2] = rows[0][2]; out[0][3] = rows[0][3];
@@ -657,12 +773,6 @@ struct Matrix44 {
     }
 };
 
-/*! The 4x4 identity matrix */
-const Matrix44 IDENTITY44( 1.0f, 0.0f, 0.0f, 0.0f,
-						   0.0f, 1.0f, 0.0f, 0.0f,
-						   0.0f, 0.0f, 1.0f, 0.0f,
-						   0.0f, 0.0f, 0.0f, 1.0f );
-
 /*! Stores a color along with alpha translucency */
 struct Color4 {
 	float r; /*!< The red component of this color.  Should be between 0.0f and 1.0f. */ 
@@ -758,21 +868,8 @@ ostream & operator<<( ostream & out, Matrix22 const & val );
 ostream & operator<<( ostream & out, Float3 const & val );
 ostream & operator<<( ostream & out, Matrix33 const & val );
 ostream & operator<<( ostream & out, Float4 const & val );
-ostream & operator<<( ostream & out, Matrix44 const & val );
 ostream & operator<<( ostream & out, Color4 const & val );
 ostream & operator<<( ostream & out, Quaternion const & val );
 
 
-
-Matrix33 MultMatrix33( Matrix33 const & a, Matrix33 const & b );
-Matrix44 MultMatrix44( Matrix44 const & a, Matrix44 const & b );
-float DetMatrix33( Matrix33 const & m );
-float DetMatrix44( Matrix44 const & m );
-float AdjMatrix44( Matrix44 const & m, int r, int c );
-Matrix44 InverseMatrix44( Matrix44 const & m );
-void SetIdentity33( Matrix33 & m );
-void SetIdentity44( Matrix44 & m );
-void PrintMatrix33( Matrix33 const & m, ostream & out );
-void PrintMatrix44( Matrix44 const & m, ostream & out );
-
 #endif
diff --git a/niflib.cpp b/niflib.cpp
index f683e0a2c67f7111af19a674ee4f4cb8726533a3..f4789c06f9a501347cfc4664db94330ccef475f2 100644
--- a/niflib.cpp
+++ b/niflib.cpp
@@ -553,7 +553,7 @@ void BuildUpBindPositions( blk_ref const & block ) {
 			//Post-multipy the block's bind matrix with the parent's bind matrix
 			Matrix44 par_mat = par_node->GetWorldBindPos();
 			Matrix44 blk_mat = blk_node->GetWorldBindPos();
-			Matrix44 result = MultMatrix44( blk_mat, par_mat);
+			Matrix44 result = blk_mat * par_mat;
 
 			//Store result back to block bind position
 			blk_node->SetWorldBindPos( result );