File content as of revision 0:431ee55036ca:

#include "Matrix.h"

Matrix::Matrix() : _rows(1), _cols(1)
{
  // Allocate memory
  _data = (float **)malloc(sizeof(float *));
  _data[0] = (float *)malloc(sizeof(float));
  // Initialize values
  _data[0][0] = 0.0;
}

Matrix::Matrix(int rows, int cols) : _rows(rows), _cols(cols)
{

  // Allocate memory
  _data = (float **)malloc(_rows * sizeof(float *)); // = new float*[_rows];
  for (int i = 1; i <= _rows; i++)
  {
    _data[i-1] = (float *)malloc(_cols * sizeof(float)); // = new float[_cols];
  }
  // Initialize values
  for (int i = 1; i <= _rows; i++)
  {
    for (int j = 1; j <= _cols; j++)
    {
      _data[i-1][j-1] = 0.0;
    }
  }
}

Matrix::~Matrix()
{
    // Free memory
    for (int i = 0; i < _rows; ++i) {
        free(_data[i]);
    }
    free(_data);
}

float& Matrix::operator()(int row, int col)
{
  // Return cell data
  return _data[row-1][col-1];
}

Matrix operator+(const Matrix& A, const Matrix& B)
{
    // Auxiliary matrix where C = A+B
    Matrix C(A._rows,A._cols);
    for (int i = 1; i <= A._rows; i++)
    {
      for (int j = 1; j <= A._cols; j++)
      {
          C._data[i-1][j-1] += A._data[i-1][j-1]+B._data[i-1][j-1];
      }
    }
    return C;
}

Matrix operator-(const Matrix& A, const Matrix& B)
{
    // Auxiliary matrix where C = A-B
    Matrix C(A._rows,A._cols);
    for (int i = 1; i <= A._rows; i++)
    {
      for (int j = 1; j <= A._cols; j++)
      {
          C._data[i-1][j-1] += A._data[i-1][j-1]-B._data[i-1][j-1];
      }
    }
    return C;
}

Matrix operator*(const Matrix& A, const Matrix& B)
{
    // Auxiliary matrix where C = A*B
    Matrix C(A._rows,B._cols);
    for (int i = 1; i <= A._rows; i++)
    {
      for (int j = 1; j <= B._cols; j++)
      {
        for (int k = 1; k <= A._cols; k++)
        {
          C._data[i-1][j-1] += A._data[i-1][k-1]*B._data[k-1][j-1];
        }
      }
    }
    return C;
}

Matrix operator*(float scalar, const Matrix& A)
{
    // Auxiliary matrix where C = scalar*A
    Matrix C(A._rows,A._cols);
    for (int i = 1; i <= A._rows; i++)
    {
      for (int j = 1; j <= A._cols; j++)
      {
        C._data[i-1][j-1] += scalar*A._data[i-1][j-1];
      }
    }
    return C;
}

Matrix transpose(const Matrix& A)
{
    // Auxiliary matrix where C = A'
    Matrix C(A._cols, A._rows);
    for (int i = 1; i <= C._rows; i++)
    {
        for (int j = 1; j <= C._cols; j++)
        {
            C._data[i-1][j-1] = A._data[j-1][i-1];
        }
    }
    return C;
}

Matrix inverse(const Matrix& A)
{
    // Apply A = LDL' factorization where L is a lower triangular matrix and D 
    // is a block diagonal matrix
    Matrix L(A._cols, A._rows);
    Matrix D(A._cols, A._rows);
    float L_sum;
    float D_sum;
    for (int i = 1; i <= A._rows; i++)
    {
        for (int j = 1; j <= A._rows; j++)
        {
            if (i >= j)
            {
                if (i == j)
                {
                    L._data[i-1][j-1] = 1.0;
                    D_sum = 0.0;
                    for (int k = 1; k <= (i-1); k++)
                    {
                        D_sum += L._data[i-1][k-1]*L._data[i-1][k-1]*D._data[k-1][k-1];
                    }
                    D._data[i-1][j-1] = A._data[i-1][j-1] - D_sum;
                }
                else
                {
                    L_sum = 0.0;
                    for (int k = 1; k <= (j-1); k++)
                    {
                        L_sum += L._data[i-1][k-1]*L._data[j-1][k-1]*D._data[k-1][k-1];
                    }
                    L._data[i-1][j-1] = (1.0/D._data[j-1][j-1])*(A._data[i-1][j-1]-L_sum);
                }
            }
        }
    }
    // Compute the inverse of L and D matrices 
    Matrix L_inv(A._cols, A._rows);
    Matrix D_inv(A._cols, A._rows);
    float L_inv_sum;
    for (int i = 1; i <= A._rows; i++)
    {
        for (int j = 1; j <= A._rows; j++)
        {
            if (i >= j)
            {
                if (i == j)
                {
                    L_inv._data[i-1][j-1] = 1.0/L._data[i-1][j-1];
                    D_inv._data[i-1][j-1] = 1.0/D._data[i-1][j-1];
                }
                else
                {
                    L_inv_sum = 0.0;
                    for (int k = j; k <= (i-1); k++)
                    {
                        L_inv_sum += L._data[i-1][k-1]*L_inv._data[k-1][j-1];
                    }
                    L_inv._data[i-1][j-1] = -L_inv_sum;
                }
            }
        }
    }
    // Compute the inverse of A matrix in terms of the inverse of L and D 
    // matrices 
    return transpose(L_inv)*D_inv*L_inv;
}