Finish matrix unit tests

Base/include/Matrix2x2.h

@@ -0,0 +1,47 @@
+/*******************************************************************************
+ * Copyright (c) 2014 Pavel Krajcevski
+ *
+ * This software is provided 'as-is', without any express or implied
+ * warranty. In no event will the authors be held liable for any damages
+ * arising from the use of this software.
+ * 
+ * Permission is granted to anyone to use this software for any purpose,
+ * including commercial applications, and to alter it and redistribute it
+ * freely, subject to the following restrictions:
+ * 
+ *    1. The origin of this software must not be misrepresented; you must not
+ *    claim that you wrote the original software. If you use this software
+ *    in a product, an acknowledgment in the product documentation would be
+ *    appreciated but is not required.
+ *
+ *    2. Altered source versions must be plainly marked as such, and must not be
+ *    misrepresented as being the original software.
+ *
+ *    3. This notice may not be removed or altered from any source
+ *    distribution.
+ *
+ ******************************************************************************/
+
+#ifndef BASE_INCLUDE_MATRIX2X2_H_
+#define BASE_INCLUDE_MATRIX2X2_H_
+
+#include "MatrixSquare.h"
+
+namespace FasTC {
+
+  template <typename T>
+  class Matrix2x2 : public MatrixSquare<T, 2> {
+   public:
+    // Constructors
+    Matrix2x2() { }
+    Matrix2x2(const Matrix2x2<T> &other)
+      : MatrixSquare<T, 2>(other) { }
+    Matrix2x2(const MatrixSquare<T, 2> &other)
+      : MatrixSquare<T, 2>(other) { }
+    Matrix2x2(const MatrixBase<T, 2, 2> &other)
+      : MatrixSquare<T, 2>(other) { }
+  };
+  REGISTER_ONE_TEMPLATE_MATRIX_TYPE(Matrix2x2);
+};
+
+#endif  // BASE_INCLUDE_MATRIX2X2_H_

Base/include/Matrix3x3.h

@@ -31,17 +31,17 @@ namespace FasTC {
 
   template <typename T>
   class Matrix3x3 : public MatrixSquare<T, 3> {
-
    public:
-    
     // Constructors
     Matrix3x3() { }
-    Matrix3x3(const MatrixSquare<T, 3> &other) {
-      for(int i = 0; i < kNumElements; i++) {
-        mat[i] = other[i];
-      }
-    }
+    Matrix3x3(const Matrix3x3<T> &other)
+      : MatrixSquare<T, 3>(other) { }
+    Matrix3x3(const MatrixSquare<T, 3> &other)
+      : MatrixSquare<T, 3>(other) { }
+    Matrix3x3(const MatrixBase<T, 3, 3> &other)
+      : MatrixSquare<T, 3>(other) { }
   };
+  REGISTER_ONE_TEMPLATE_MATRIX_TYPE(Matrix3x3);
 };
 
 #endif  // BASE_INCLUDE_MATRIX3X3_H_

Base/include/Matrix4x4.h

@@ -59,17 +59,17 @@ namespace FasTC {
 
   template <typename T>
   class Matrix4x4 : public MatrixSquare<T, 4> {
-
    public:
-    
     // Constructors
     Matrix4x4() { }
+    Matrix4x4(const Matrix4x4<T> &other)
+      : MatrixSquare<T, 4>(other) { }
     Matrix4x4(const MatrixSquare<T, 4> &other)
       : MatrixSquare<T, 4>(other) { }
-
     Matrix4x4(const MatrixBase<T, 4, 4> &other)
       : MatrixSquare<T, 4>(other) { }
   };
+  REGISTER_ONE_TEMPLATE_MATRIX_TYPE(Matrix4x4);
 };
 
 #endif  // BASE_INCLUDE_MATRIX3X3_H_

Base/include/MatrixBase.h

@@ -158,6 +158,13 @@ namespace FasTC {
     static const EVectorType kVectorType = eVectorType_Matrix; \
   }
 
+  #define REGISTER_ONE_TEMPLATE_MATRIX_SIZED_TYPE(TYPE)        \
+  template<typename T, const int SIZE>                         \
+  class VectorTraits< TYPE <T, SIZE> > {                       \
+  public:                                                      \
+    static const EVectorType kVectorType = eVectorType_Matrix; \
+  }
+
   // Define matrix multiplication for * operator
   template<typename TypeOne, typename TypeTwo>
   class MultSwitch<

Base/include/MatrixSquare.h

@@ -26,6 +26,8 @@
 #define BASE_INCLUDE_MATRIXSQUARE_H_
 
 #include "MatrixBase.h"
+#include <cstdlib>
+#include <ctime>
 
 namespace FasTC {
 
@@ -40,54 +42,82 @@ namespace FasTC {
     MatrixSquare(const MatrixBase<T, N, N> &other)
       : MatrixBase<T, N, N>(other) { }
 
+    MatrixSquare<T, N> Transpose() const {
+      return MatrixBase<T, N, N>::Transpose();
+    }
+
     // Does power iteration to determine the principal eigenvector and eigenvalue.
     // Returns them in eigVec and eigVal after kMaxNumIterations
-    int PowerMethod(VectorBase<T, N> &eigVec, T *eigVal = NULL,
-                    const int kMaxNumIterations = 200) {
+    int PowerMethod(VectorBase<T, N> &eigVec,
+                    T *eigVal = NULL,
+                    const int kMaxNumIterations = 200,
+                    const unsigned int kSeed = time(NULL)) {
+
+      srand(kSeed);
       int numIterations = 0;
 
-      // !SPEED! Find eigenvectors by using the power method. This is good because the
-      // matrix is only 4x4, which allows us to use SIMD...
       VectorBase<T, N> b;
       for(int i = 0; i < N; i++)
-        b[i] = T(1.0);
-      
-      b /= b.Length();
+        b[i] = static_cast<T>(rand());
+      b.Normalize();
 
+      bool badEigenValue = false;
       bool fixed = false;
       numIterations = 0;
       while(!fixed && ++numIterations < kMaxNumIterations) {
 
-        VectorBase<T, N> newB = (*this).operator*(b);
+        VectorBase<T, N> newB = (*this) * b;
 
-        // !HACK! If the principal eigenvector of the covariance matrix
-        // converges to zero, that means that the points lie equally 
-        // spaced on a sphere in this space. In this (extremely rare)
-        // situation, just choose a point and use it as the principal 
-        // direction.
+        // !HACK! If the principal eigenvector of the matrix
+        // converges to zero, that could mean that there is no
+        // principal eigenvector. However, that may be due to
+        // poor initialization of the random vector, so rerandomize
+        // and try again.
         const float newBlen = newB.Length();
         if(newBlen < 1e-10) {
-          eigVec = b;
-          if(eigVal) *eigVal = 0.0;
-          return numIterations;
+          if(badEigenValue) {
+            eigVec = b;
+            if(eigVal) *eigVal = 0.0;
+            return numIterations;
+          }
+
+          VectorBase<T, N> b;
+          for(int i = 0; i < N; i++)
+            b[i] = static_cast<T>(rand());
+
+          b.Normalize();
+          badEigenValue = true;
         }
 
-        T len = newB.Length();
-        newB /= len;
-        if(eigVal)
-          *eigVal = len;
+        // Normalize
+        newB.Normalize();
 
-        if(fabs(1.0f - (b.Dot(newB))) < 1e-5)
+        // If the new eigenvector is close enough to the old one,
+        // then we've converged.
+        if(fabs(1.0f - (b.Dot(newB))) < 1e-8)
           fixed = true;
 
+        // Save and continue.
         b = newB;
       }
 
-      eigVec = b;  
+      // Store the eigenvector in the proper variable.
+      eigVec = b;
+
+      // Store eigenvalue if it was requested
+      if(eigVal) {
+        VectorBase<T, N> result = (*this) * b;
+        *eigVal = result.Length() / b.Length();
+      }
+
       return numIterations;
     }
 
+   private:
+
   };
+  REGISTER_ONE_TEMPLATE_MATRIX_SIZED_TYPE(MatrixSquare);
+
 };
 
 #endif  // BASE_INCLUDE_MATRIXSQUARE_H_

Base/test/TestMatrix.cpp

@@ -123,7 +123,13 @@ TEST(MatrixBase, PointerConversion) {
 }
 
 TEST(MatrixBase, CastVector) {
+  srand(time(NULL));
+
   FasTC::MatrixBase<float, 3, 2> v3f;
+  for(int i = 0; i < 6; i++) {
+    v3f[i] = static_cast<float>(rand());
+  }
+
   FasTC::MatrixBase<double, 3, 2> v3d = v3f;
   FasTC::MatrixBase<int, 3, 2> v3i = v3f;
   for(int i = 0; i < 3*2; i++) {
@@ -198,28 +204,137 @@ TEST(MatrixBase, VectorMultiplication) {
   EXPECT_EQ(v[4], 0 + (3 * 2) + (5 * 3));
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+#include "MatrixSquare.h"
+
 TEST(MatrixSquare, Constructors) {
-  // Stub
-  EXPECT_EQ(0, 1);
+  FasTC::MatrixBase<int, 3, 3> m;
+  m(0, 0) = 1; m(0, 1) = 2; m(0, 2) = 3;
+  m(1, 0) = 2; m(1, 1) = 3; m(1, 2) = 4;
+  m(2, 0) = 3; m(2, 1) = 4; m(2, 2) = 5;
+
+  FasTC::MatrixSquare<int, 3> sqm (m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_EQ(m[i], sqm[i]);
+  }
+
+  FasTC::MatrixSquare<float, 3> fsqm(m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_NEAR(m[i], fsqm[i], kEpsilon);
+  }
+
+  FasTC::MatrixSquare<float, 3> fcsqm(sqm);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_NEAR(fcsqm[i], sqm[i], kEpsilon);
+  }
 }
 
-TEST(MatrixSquare, EigenvalueCalculation) {
-  // Stub
-  EXPECT_EQ(0, 1);
+TEST(MatrixSquare, PowerMethod) {
+  FasTC::MatrixSquare<double, 2> A;
+  A(0, 0) = 0.8f; A(0, 1) = 0.3f;
+  A(1, 0) = 0.2f; A(1, 1) = 0.7f;
+
+  double e;
+  FasTC::VectorBase<double, 2> x;
+  A.PowerMethod(x, &e, 20, 200);
+
+  EXPECT_NEAR(x[0], 0.83205f, 0.0001);
+  EXPECT_NEAR(x[1], 0.5547f, 0.0001);
+
+  EXPECT_NEAR(e, 1.f, 0.0001);
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+#include "Matrix2x2.h"
+
 TEST(Matrix2x2, Constructors) {
-  // Stub
-  EXPECT_EQ(0, 1);
+  FasTC::MatrixBase<int, 2, 2> m;
+  m(0, 0) = 1; m(0, 1) = 2;
+  m(1, 0) = 2; m(1, 1) = 3;
+
+  FasTC::MatrixSquare<int, 2> sqm (m);
+  for(int i = 0; i < 4; i++) {
+    EXPECT_EQ(m[i], sqm[i]);
+  }
+
+  FasTC::Matrix2x2<int> tbtm (m);
+  for(int i = 0; i < 4; i++) {
+    EXPECT_EQ(m[i], tbtm[i]);
+  }
+
+  FasTC::Matrix2x2<float> fsqm(m);
+  for(int i = 0; i < 4; i++) {
+    EXPECT_NEAR(m[i], fsqm[i], kEpsilon);
+  }
+
+  FasTC::Matrix2x2<float> fcsqm(sqm);
+  for(int i = 0; i < 4; i++) {
+    EXPECT_NEAR(fcsqm[i], sqm[i], kEpsilon);
+  }
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+#include "Matrix3x3.h"
+
 TEST(Matrix3x3, Constructors) {
-  // Stub
-  EXPECT_EQ(0, 1);
+  FasTC::MatrixBase<int, 3, 3> m;
+  m(0, 0) = 1; m(0, 1) = 2; m(0, 2) = 3;
+  m(1, 0) = 2; m(1, 1) = 3; m(1, 2) = 4;
+  m(2, 0) = 3; m(2, 1) = 4; m(2, 2) = 5;
+
+  FasTC::MatrixSquare<int, 3> sqm (m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_EQ(m[i], sqm[i]);
+  }
+
+  FasTC::Matrix3x3<int> tbtm (m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_EQ(m[i], tbtm[i]);
+  }
+
+  FasTC::Matrix3x3<float> fsqm(m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_NEAR(m[i], fsqm[i], kEpsilon);
+  }
+
+  FasTC::Matrix3x3<float> fcsqm(sqm);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_NEAR(fcsqm[i], sqm[i], kEpsilon);
+  }
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+#include "Matrix4x4.h"
+
 TEST(Matrix4x4, Constructors) {
-  // Stub
-  EXPECT_EQ(0, 1);
+  FasTC::MatrixBase<int, 4, 4> m;
+  m(0, 0) = 1; m(0, 1) = 2; m(0, 2) = 3; m(0, 3) = 4;
+  m(1, 0) = 2; m(1, 1) = 3; m(1, 2) = 4; m(1, 3) = 5;
+  m(2, 0) = 3; m(2, 1) = 4; m(2, 2) = 5; m(2, 3) = 6;
+  m(3, 0) = 4; m(3, 1) = 5; m(3, 2) = 6; m(3, 3) = 7;
+
+  FasTC::MatrixSquare<int, 4> sqm (m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_EQ(m[i], sqm[i]);
+  }
+
+  FasTC::Matrix4x4<int> tbtm (m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_EQ(m[i], tbtm[i]);
+  }
+
+  FasTC::Matrix4x4<float> fsqm(m);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_NEAR(m[i], fsqm[i], kEpsilon);
+  }
+
+  FasTC::Matrix4x4<float> fcsqm(sqm);
+  for(int i = 0; i < 9; i++) {
+    EXPECT_NEAR(fcsqm[i], sqm[i], kEpsilon);
+  }
 }