diff --git a/src/matrix.py b/src/matrix.py
new file mode 100644
index 0000000..40f93e8
--- /dev/null
+++ b/src/matrix.py
@@ -0,0 +1,315 @@
+# Copyright (C) 2007 Darren Lee Weber
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+# 02111-1307, USA.
+
+__version__ = "$Revision: 1.9 $" # $Date: 2007/06/14 00:24:57 $
+
+class Matrix:
+    """
+    Create and manipulate a matrix object
+ 
+    Matrix(data, dim)
+ 
+    data = list of lists (currently only 2D)
+    dim=(row,col) tuple of int
+ 
+    For example,
+ 
+    #data = [[0.0] * c for i in xrange(r)]
+    data = [[0.0,0.1],[1.0,1.1],[2.0,2.1]]
+    rowN =len(data)
+    colN =len(data[0])
+    m = Matrix(data)
+    m = Matrix(data,dim=(rowN, colN))
+ 
+    d1 = [[0.0, 0.1], [1.0, 1.1], [2.0, 2.1]] # 3x2 matrix
+    d2 = [[0.0, 0.1, 0.2], [1.0, 1.1, 1.2]]   # 2x3 matrix
+    m1 = Matrix(d1)
+    m2 = Matrix(d2)
+    #m3 = m1 + m2                             # dimension error
+    m3 = m1 + m2.transpose()
+    m3 = m1 - m2.transpose()
+    m3 = m1 * m2                        # 3x3
+    m3 = m2 * m1                        # 2x2
+ 
+    m1[2,:]
+    m1[:,2]
+    """
+ 
+    def __init__(self, data=None, dim=None):
+        """
+        create a matrix instance.
+ 
+        m = Matrix([data [, dim]])
+ 
+        <data> is a 2D matrix comprised of a nested list of floats
+        <dim> is a tuple of int values for the row and column size (r,c)
+ 
+        eg:
+        data = [[0.0,0.1],[1.0,1.1],[2.0,2.1]]
+        dim = (3,2) # or (len(data),len(data[0]))
+        """
+ 
+        if data != None:
+            # check data for the cell types (ensure float)?
+            self.data = data
+            r = len(data)
+            c = len(data[0])
+ 
+            # Are all the rows the same length?
+            rowLenCheck = sum([len(data[i]) != c for i in range(r)])
+            if rowLenCheck > 0:
+                raise ValueError
+            else:
+                self.dim = (r,c)
+ 
+            if dim != None:
+                if (dim[0] == r) and (dim[1] == c):
+                    self.dim = (r,c)
+                else:
+                    # over-ride the dim input, do not reshape data!
+                    # print a warning?
+                    self.dim = (r,c)
+        else:
+            if dim != None:
+                if len(dim) == 2:
+                    self.dim = tuple(dim)
+                    r = dim[0]
+                    c = dim[1]
+                else:
+                    # maybe a new exception type?
+                    arg = ("len(dim) != 2: ", dim)
+                    raise ValueError, arg
+ 
+                # BEGIN ALT ----------------------------------------
+                # Does this give unique memory for each element?
+                # self.data = [[0.0] * c for i in xrange(r)]
+ 
+                # It seems that the initialization does not generate
+                # unique memory elements because all list elements
+                # refer to the same number object (0.0), but
+                # modification of any element creates a unique value,
+                # without changing any other values, eg:
+ 
+                ##>>> x = [[0.0] * 3 for i in xrange(2)]
+                ##>>> id(x)
+                # 3079625068L
+                # >>> id(x[0][0])
+                # 136477300
+                # >>> id(x[0][1])
+                # 136477300
+                # >>> id(x[1][1])
+                # 136477300
+                # >>> x[0][0] = 1.0
+                # >>> x
+                # [[1.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
+                # >>>
+                # END ALT ----------------------------------------
+ 
+                # create a zero row vector, with unique memory for each element
+                self.data = [[x * 0.0 for x in range(c)]]
+                for i in range(1,r):
+                    self.data.append([x * 0.0 for x in range(c)])
+            else:
+                self.data = []
+                self.dim = (0,0)
+                #print self.__doc__
+ 
+    def __getitem__(self, i):
+        """
+        matrix[r,c] returns values from matrix.data, eg:
+        data = [[0.0,0.1],[1.0,1.1],[2.0,2.1]]
+        m = Matrix(data)
+        m[2,:]
+        >> [2.0, 2.1000000000000001]
+        """
+        r = i[0]
+        c = i[1]
+        #print "index: (%s, %s)" % (r,c)
+        #print "value: ", self.data[r][c]
+        return self.data[r][c]
+ 
+    def reshape(self, newdim=None):
+        'reshape a matrix object: matrix.reshape(newdim)'
+        print "something to implement later"
+        pass
+ 
+    def transpose(self):
+        'transpose a matrix: m2 = m1.transpose()'
+        m = Matrix(dim=(self.dim[1],self.dim[0]))
+        for r in range(self.dim[0]):
+            for c in range(self.dim[1]):
+                m.data[c][r] = self.data[r][c]
+        return m
+ 
+    def __add__(self, q):
+        '''
+        matrix addition:
+        m3 = matrix1 + matrix2
+        m3 = matrix1 + float
+        m3 = matrix1 + int
+        '''
+        if isinstance(q, Matrix):
+            if self.dim != q.dim:
+                arg = ("p.dim != q.dim", self.dim, q.dim)
+                raise IndexError, arg
+            else:
+                # do the addition
+                m = Matrix(dim=self.dim)
+                for r in range(self.dim[0]): # rows of p and q
+                    m.data[r] = map(lambda x, y: x + y, self.data[r],q.data[r])
+                return m
+        elif isinstance(q, float) or isinstance(q, int):
+            # add a scalar value
+            m = Matrix(dim=self.dim)
+            for r in range(self.dim[0]): # rows
+                m.data[r] = map(lambda x: x + q, self.data[r])
+            return m
+        else:
+            arg = ("q is not a matrix, float or int", q)
+            raise TypeError, arg
+ 
+    def __sub__(self, q):
+        '''
+        matrix subtraction:
+        m3 = matrix1 - matrix2
+        m3 = matrix1 - float
+        m3 = matrix1 - int
+        '''
+        if isinstance(q, Matrix):
+            if self.dim != q.dim:
+                arg = ("p.dim != q.dim", self.dim, q.dim)
+                raise IndexError, arg
+            else:
+                # do the subtraction
+                m = Matrix(dim=self.dim)
+                for r in range(self.dim[0]): # rows of p and q
+                    m.data[r] = map(lambda x, y: x - y, self.data[r], q.data[r])
+                return m
+        elif isinstance(q, float) or isinstance(q, int):
+            # subtract a scalar value
+            m = Matrix(dim=self.dim)
+            for r in range(self.dim[0]): # rows
+                m.data[r] = map(lambda x: x - q, self.data[r])
+            return m
+        else:
+            arg = ("q is not a matrix, float or int", q)
+            raise TypeError, arg
+ 
+    def __mul__(self, q):
+        """
+        multiply two matrices:
+        m = p * q  # p.dim[1] == q.dim[0]
+ 
+        multiply a matrix with a scalar:
+        m = p * q  # where q is a float or int value
+        """
+        if isinstance(q, Matrix):
+            if self.dim[1] != q.dim[0]:
+                arg = ("p.dim[1] != q.dim[0]", self.dim[1], q.dim[0])
+                raise IndexError, arg
+            else:
+                # do the multiplication
+                m = Matrix(dim=(self.dim[0], q.dim[1]))
+                for r in range(self.dim[0]): # rows of p
+                    for c in range(q.dim[1]): # cols of q
+                        # get the dot product of p(r,:) with q(:,c)
+                        pRowVec = self.data[r]
+                        qColVec = [q.data[a][c] for a in xrange(q.dim[0])]
+                        m.data[r][c] = sum(map(lambda x, y: x * y, pRowVec, qColVec))
+                return m
+        elif isinstance(q, float) or isinstance(q, int):
+            # subtract a scalar value
+            m = Matrix(dim=self.dim)
+            for r in range(self.dim[0]): # rows
+                m.data[r] = map(lambda x: x * q, self.data[r])
+            return m
+        else:
+            arg = ("q is not a matrix, float or int", q)
+            raise TypeError, arg
+ 
+    def __div__(self, q):
+        """
+        Divide a matrix with a scalar, eg:
+        m = p / q  # where q is a float or int value
+        This operator will not return a matrix inverse
+        """
+        if isinstance(q, Matrix):
+            # let's not do matrix divide in python, leave the inverse
+            # to a c/c++ library
+            arg = ("q is a matrix: will not calculate inverse", q)
+            raise TypeError, arg
+        elif isinstance(q, float) or isinstance(q, int):
+            # divide a scalar value
+            m = Matrix(dim=self.dim)
+            for r in range(self.dim[0]): # rows
+                m.data[r] = map(lambda x: x / q, self.data[r])
+            return m
+        else:
+            arg = ("q is not a matrix, float or int", q)
+            raise TypeError, arg
+ 
+    def __len__(self):
+        return self.dim[0] * self.dim[1]
+ 
+    def __str__(self):
+        # print the matrix data
+        s = ""
+        for r in range(self.dim[0]):
+            for c in range(self.dim[1]):
+                s += "%f " % (self.data[r][c])
+            s += "\n"
+        return s
+ 
+    def printFormat(self, format):
+        """
+        print the matrix data nicely formatted, eg:
+        matrix.printFormat("%8.4f")
+        """
+        for r in range(self.dim[0]):
+            for c in range(self.dim[1]):
+                print format % (self.data[r][c]),
+            print
+ 
+    def __repr__(self):
+        # return something that will recreate the object
+        return "Matrix(%s, %s)" % (self.data, self.dim)
+
+
+
+# --------------------------------------------------------------------------------
+# Explore the functionality - should be unit testing
+
+testing = 0
+if testing:
+    d1 = [[0.0, 0.1], [1.0, 1.1], [2.0, 2.1]] # 3x2 matrix
+    d2 = [[0.0, 0.1, 0.2], [1.0, 1.1, 1.2]]   # 2x3 matrix
+    m1 = Matrix(d1)
+    m2 = Matrix(d2)
+    #m3 = m1 + m2                             # "dimension" error
+    m3 = m1 + m2.transpose()
+    m3 = m1 - m2.transpose()
+    m3 = m1 * m2                        # 3x3
+    m3 = m2 * m1                        # 2x2
+    m3 += 10.0
+    m3 -= 10.0
+    m3 += 10
+    m3 -= 10
+    m3 /= 10.0
+    m3 *= 10.0
+    m3 /= 10
+    m3 *= 10
+ 
\ No newline at end of file