File content as of revision 1:085840a3d767:

#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 = 1; i <= _rows; i++) {
        free(_data[i-1]);
    }
    free(_data);
}

Matrix& Matrix::operator=(const Matrix& m)
{
    for (int i = 1; i <= _rows; i++) {
        for (int j = 1; j <= _cols; j++) {
            _data[i-1][j-1] = m._data[i-1][j-1];
        }
    }
    return *this;
}

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 temp(m1);
    return (temp += m2);
    */
    
    //Matrix A(m1);
    //Matrix B(m2);
    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 operator*(const Matrix& A, float scalar)
{
    // 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] += A._data[i-1][j-1]*scalar;
        }
    }
    return C;
}

Matrix operator/(const Matrix& A, float scalar)
{
    // 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] += A._data[i-1][j-1]/scalar;
        }
    }
    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;
}

float norm(const Matrix& A)
{
    // Auxiliary matrix where C = A'
    float n;
    n = 0.0;
    for (int i = 1; i <= A._rows; i++) {
        n += A._data[i-1][0]*A._data[i-1][0];
        }
    return sqrt(n);
}