Matrix.cpp

#include <iostream>
#include <math.h>
#include <iomanip>
#include <assert.h>
#include <ctime>
#include <cstdlib>
#include "Matrix.hpp"

using namespace std;


/**
 * Matrix default constructor
 */
Matrix::Matrix()
{
    n = new int;
    m = new int;
    *n = 3;
    *m = 3;
    
    allocEmptyMatrix(*n, *m);
}

/**
 * @param   int n - number of rows
 * @param   int m - number of columns
 */
Matrix::Matrix(int a_n, int a_m)
{
    this->n = new int;
    this->m = new int;
    *n = a_n;
    *m = a_m;
    
    allocEmptyMatrix(*n, *m);
}

/**
 * @param   double** tdarray - two dimensional array
 * @param   int n - number of rows
 * @param   int m - number of columns
 */
Matrix::Matrix(double** tdarray, int a_n, int a_m)
{
    this->n = new int;
    this->m = new int;
    *n = a_n;
    *m = a_m;
    
    allocEmptyMatrix(*n, *m);
    
    for (int i = 0; i < *n; i++) {
        for (int j = 0; j < *m; j++) {
            this->matrix[i][j] = tdarray[i][j];
        }
    }
}

/**
 * @param   const Matrix - Matrix to be copied
 */
Matrix::Matrix(const Matrix & ptr)
{
    n = new int;
    m = new int;
    *n = *(ptr.n);
    *m = *(ptr.m);
    
    allocEmptyMatrix(*n, *m);
    
    for (int i = 0; i < *n; i++) {
        for (int j = 0; j < *m; j++) {
            this->matrix[i][j] = ptr.matrix[i][j];
        }
    }
}

/**
 * Matrix destructor
 */
Matrix::~Matrix()
{
    for (int i = 0; i < *n; i++) {
        delete [] matrix[i];
    }
    delete [] matrix;

    delete n;
    delete m;
    n = m = 0;
}

/**
 * @param   int n - number of rows
 * @param   int m - number of columns
 * @return  void
 */
void Matrix::allocEmptyMatrix(int n, int m)
{
    matrix = new double*[n];
    for (int i = 0; i < n; i++) {
        matrix[i] = new double[m];
    }
}

/**
 * Primitive output method
 * @return  void
 */
void Matrix::print() const
{
    for (int i = 0; i < *n; i++) {
        cout << "|";
        for (int j = 0; j < *m; j++) {
            cout << setw(8) << matrix[i][j];
        }
        cout << "|" << endl;
    }
}

/**
 * @param   int n - number of rows
 * @param   int m - number of columns
 * @return  Matrix
 */
Matrix Matrix::ones(int n, int m)
{
    Matrix A(n, m);
    
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            A.matrix[i][j] = 1;
        }
    }
    
    return A;
}

/**
 * @param   int n - number of rows
 * @param   int m - number of columns
 * @return  Matrix
 */
Matrix Matrix::unity(int n, int m)
{
    Matrix A(n, m);
    
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            if (i == j) {
                A.matrix[i][j] = 1.0;
            }
        }
    }
    return A;
}

/**
 * @param   int n - number of rows
 * @param   int m - number of columns
 * @return  Matrix
 */
Matrix Matrix::random(int n, int m)
{   
    Matrix A(n, m);
	
    srand(time(NULL));   
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            A.matrix[i][j] = rand() % 2;
        }
    }
    
    return A;
}

/**
 * @return Matrix
 */
Matrix Matrix::transposition()
{
    Matrix A(*m, *n);
    
    for (int i = 0; i < *n; i++) {
        for (int j = 0; j < *m; j++) {
            A.matrix[j][i] = this->matrix[i][j];
        }
    }
    
    return A;
}

/**
 * @param   Matrix
 * @return  Matrix
 */
Matrix Matrix::add(const Matrix A) const
{
    if (!isSameType(A)) {
        cout << "You cannot add two Matrixes with diferent type (n, m)" << endl;
        return A;
    }
    
    Matrix C(*n, *m);

    for (int i = 0; i < *n; i++) {
        for (int j = 0; j < *m; j++) {
            C.matrix[i][j] = this->matrix[i][j] + A.matrix[i][j];
        } 
    }
    
    return C;
}

/**
 * @param   Matrix
 * @return  Matrix
 */
Matrix Matrix::subtract(const Matrix A) const
{
    if (!isSameType(A)) {
        cout << "You cannot subtract two Matrixes with diferent type (n, m)" << endl;
        return A;
    }
    
    Matrix C(*n, *m);

    for (int i = 0; i < *n; i++) {
        for (int j = 0; j < *m; j++) {
            C.matrix[i][j] = this->matrix[i][j] - A.matrix[i][j];
        }
    }
    
    return C;
}

/**
 * @param   Matrix
 * @return  Matrix
 */
Matrix Matrix::multiply(const Matrix A) const
{    
    Matrix C(*n, *(A.m));
    double temp = .0;

    for (int i = 0; i < *(C.n); i++) {
        for (int j = 0; j < *(C.m); j++) {
            for (int k = 0; k < *m && k < *(A.n); k++) {
                temp  = temp + this->matrix[i][k] * A.matrix[k][j];
            }
            C.matrix[i][j] = temp;
            temp = .0;
        }
    }
    
    return C;
}


/**
 * @param   Matrix
 * @return  bool
 */
bool Matrix::equals(const Matrix A) const
{
    if (!isSameType(A)) {
        return false;
    }
    
    for (int i = 0; i < *n; i++) {
        for (int j = 0; j < *m; j++) {
            if (this->matrix[i][j] != A.matrix[i][j]) {
                return false;
            }
        }
    }
    return true;
}

/**
 * @return  double
 * @throws  InvalidStateException
 */
double Matrix::determinant()
{
    assert(*n == *m);

    double det = 1;
    int sign = 1;
    Matrix A = gaussElimination();
    sign = A.formatToUpperTriangleForm();
    det *= (double) sign;
    A.print();
    for (int i = 0; i < *(A.n); i++) {
        det *= A.matrix[i][i];
    } 
    
    if (det == -0) { // formating zero hack
        return .0;
    }
    return det;
}

/**
 * @throws  ChangedRows()
 * @return  int - sign (+-1)
 */
int Matrix::formatToUpperTriangleForm()
{
    int indexRow = 0;
    int row = 0;
    int sign = 1;
    for (int j = 0; j < *m; j++) {
        if (isEmptyColumn(j)) {
            row++;
            continue;   
        }
        indexRow = getHighestRowIndex(j);
        if (indexRow > row) {
            switchRows(indexRow, row);    
            sign *= -1;
        }
        row++;
    }
    
    return sign;
}


/**
 * @return  Matrix
 */
Matrix Matrix::gaussElimination()
{
    Matrix A = Matrix(this->matrix, *n, *m);
    // highest number index, second highest number index
    int hni, shni = 0;
    double alpha = .0;
    for (int j = 0; j < *(A.m); j++) {
        while(!A.isSingle(j) && !A.isEmptyColumn(j)) {
            hni = A.getHighestRowIndex(j);
            for (int i = 0; i < *(A.n); i++) {
                if (A.matrix[i][j] != 0) {
                    if (hni == shni) continue;
                    if (A.matrix[i][j] >= A.matrix[shni][j]) {
                        shni = i;
                    }
                }
            }
            if (A.matrix[shni][j] == 0) {
                break;
            }

            alpha = A.matrix[hni][j] / A.matrix[shni][j] * -1;
            A.addRows(hni, alpha, shni);
            shni = 0;
        }        
    }
    return A;
}

/**
 * @throws  InvalidStateException
 * @param   int - column
 * @return  int
 */
int Matrix::getHighestRowIndex(int column)
{
    assert(column < *m);
    
    int highestRow = fabs(matrix[0][column]);
    int index = 0;
    for (int i = 0; i < *n; i++) {
        if (fabs(matrix[i][column]) > highestRow) {
            highestRow = fabs(matrix[i][column]);
            index = i;
        }
    }
    return index;
}

/**
 * @param   int - row to switch
 * @param   int - row to switch
 * @throws  InvalidStateException
 * @return  void
 */
void Matrix::switchRows(int r1, int r2)
{
    assert(r1 < *n && r2 < *n);

    double temp;
    for (int j = 0; j < *m; j++) {
        temp = matrix[r1][j];
        matrix[r1][j] = matrix[r2][j];
        matrix[r2][j] = temp;
    }
}

/**
 * @param   int - row to be added
 * @param   int - row which will add
 * @throws  InvalidStateException
 * @return  void
 */
void Matrix::addRows(int r1, int r2) 
{
    assert(r1 < *n && r2 < *n);
    
    for (int j = 0; j < *m; j++) {
        matrix[r1][j] += matrix[r2][j];
    }
}

/**
 * @param   int - row to be added into
 * @param   double - constant - must differ from zero
 * @param   int - row which will add
 * @throws  InvalidStateException
 * @return  void
 */
void Matrix::addRows(int r1, double alpha, int r2)
{
    assert(r1 < *n && r2 < *n);
    
    assert(alpha);

    for (int j = 0; j < *m; j++) {
        matrix[r1][j] += alpha * matrix[r2][j];
    }
    
}
/**
 * @param   Matrix
 * @return  bool
 */
bool Matrix::isSameType(const Matrix A) const
{
    if (*(A.n) == *n && *(A.m) == *m) {
        return true;
    }
    return false;
}

/**
 * @param   int - indef of the column
 * @return  bool
 */
bool Matrix::isSingle(int column) const
{
    int num = 0;
    for (int i = 0; i < *n; i++) {
        if (num > 1) return false;
        if (this->matrix[i][column] != 0) {
            num++;
        }
    }
    return (num == 1);
}

/**
 * @param   int - index of the column
 * @return  bool
 */
bool Matrix::isEmptyColumn(int column) const
{
    for (int i = 0; i < *n; i++) {
        if (matrix[i][column]) {
            return false;
        }
    }
    return true;
}

/**
 * Implements Doolitle's algorithm - LU decomposition
 * @see     http://www.google.com/url?sa=t&source=web&ct=res&cd=4&ved=0CCYQFjAD&url=http%3A%2F%2Fwww.engr.colostate.edu%2F~thompson%2FhPage%2FCourseMat%2FTutorials%2FCompMethods%2Fdoolittle.pdf&ei=GgPwS_7vHo_7OdblqKQI&usg=AFQjCNHvPEvfvPaSYith1lUAlj0YEaEeEQ&sig2=R_uj-74fYn6wjGHNVuDS_A
 * @param   Matrix A
 * @return  void
 */
void Matrix::LUDecomposition(Matrix A)
{
    int n = *(A.n);

    Matrix U(n, n);
    Matrix L(n, n);
 
    for (int i = 0; i < n; i++) {
        for (int j = i; j < n; j++) {
            U.matrix[i][j] = A.matrix[i][j];
            for (int k = 0; k < i; k++) {
                U.matrix[i][j] = U.matrix[i][j] - L.matrix[i][k] * U.matrix[k][j];
            }
        }
        for (int j = i; j < n; j++) {
            L.matrix[j][i] = A.matrix[j][i];
            for (int k = 0; k < i; k++) {
                L.matrix[j][i] = L.matrix[j][i] - L.matrix[j][k] * U.matrix[k][i];
            }
            L.matrix[j][i] = L.matrix[j][i] / U.matrix[i][i];
        }
    }

    Matrix Result = L.multiply(U); 
    
    L.print();
    cout << " * "<< endl;
    U.print();
    cout << endl;
    cout << " = " << endl;
    Result.print();
}

Matrix.hpp

#ifndef MATRIX_HPP
#define MATRIX_HPP

class Matrix {

public:
    Matrix();
    Matrix(int n, int m);
    Matrix(double** matrix, int n, int m);
    Matrix(const Matrix &);
    ~Matrix();
    
    void allocEmptyMatrix(int n, int m);
    void print() const;
    
    static Matrix unity(int n = 3, int m = 3);
    static Matrix ones(int n = 3, int m = 3);
    static Matrix random(int n = 3, int m = 3);
    static void LUDecomposition(Matrix A);
    
    Matrix transposition();
    Matrix add(const Matrix A) const;
    Matrix subtract(const Matrix A) const;
    Matrix multiply(const Matrix A) const;
    
    bool equals(const Matrix A) const;
    double determinant();
    
protected:
    bool isSameType(const Matrix A) const;

    int getHighestRowIndex(int column);
    
    bool isSingle(int column) const;
    bool isEmptyColumn(int column) const;
    
    void switchRows(int r1, int r2);
    void addRows(int r1, int r2);
    void addRows(int r1, double alpha, int r2);

    int formatToUpperTriangleForm();
    
    Matrix gaussElimination();
    
    /** @var double matrix  **/
    double** matrix;
    /** @var int number of rows **/
    int* n;
    /** @var int number of columns **/
    int* m;

};

#endif