Matrix.cpp

#include "Matrix.h"
#include <cassert>

using namespace std;
    int Matrix::rows() const
    {
        return mNumRows;
    }
int Matrix::cols() const
    {
        return mNumCols;
    }
    //constructors
    Matrix::Matrix()
    {
        mNumRows = 0;
        mNumCols = 0;
        mData = nullptr;
    }
    Matrix::Matrix(int numRows, int numCols)
    {
        assert(numRows > 0 && numCols > 0);
        mNumRows = numRows;
        mNumCols = numCols;
        mData = new double*[mNumRows];
        for (int i = 0; i < mNumRows; i++)
        {
            mData[i] = new double[mNumCols];
            for (int j = 0; j < mNumCols; j++)
            {
                mData[i][j] = 0.0;
            }
        }
    }
    //copy constructors
    Matrix::Matrix(const Matrix& other)
    {
        mNumRows = other.mNumRows;
        mNumCols = other.mNumCols;
        mData = new double*[mNumRows];
        for (int i = 0; i < mNumRows; i++)
        {
            mData[i] = new double[mNumCols];
            for (int j = 0; j < mNumCols; j++)
            {
                mData[i][j] = other.mData[i][j];
            }
        }
    }
    //destructors
    Matrix::~Matrix()
    {
        for (int i = 0; i<mNumRows; i++)
        {
            delete[] mData[i];
        }
        delete[] mData;
    }

    //assignment
    Matrix &Matrix::Matrix::operator=(const Matrix& other)
    {
        if (this != &other)
        {
            //clean up old
            for (int i = 0; i < mNumRows; i++)
            {
                delete[] mData[i];
            }
            delete[] mData;

            //copy new
            mNumRows = other.mNumRows;
            mNumCols = other.mNumCols;
            mData = new double*[mNumRows];
            for (int i = 0; i < mNumRows; i++)
            {
                mData[i] = new double[mNumCols];
                for (int j = 0; j < mNumCols; j++)
                {
                    mData[i][j] = other.mData[i][j];
                }
            }
        }
        return *this;
    }
    //Matrix::operator()
    double Matrix::Matrix::operator()(int i, int j) const
    {
        assert(i >= 1 && i <= mNumRows);
        assert( j >= 1 && j <= mNumCols);
        return mData[i - 1][j - 1];
    }
    double& Matrix::operator()(int i, int j)
    {
        assert(i >= 1 && i <= mNumRows);
        assert( j >= 1 && j <= mNumCols);
        return mData[i - 1][j - 1];
    }
    //add
    Matrix Matrix::operator+(const Matrix &other) const
    {
        assert(mNumRows == other.mNumRows && mNumCols == other.mNumCols);
        Matrix result(mNumRows, mNumCols);
        for (int i = 1; i <= mNumRows; i++)  // Change to 1-based indexing
        {
            for (int j = 1; j <= mNumCols; j++)
            {
                result(i, j) = (*this)(i, j) + other(i, j);
            }
        }
        return result;
    }
    //subtract
    Matrix Matrix::operator-(const Matrix &other) const
    {
        assert(mNumRows == other.mNumRows && mNumCols == other.mNumCols);
        Matrix result(mNumRows, mNumCols);
        for (int i = 1; i <= mNumRows; i++)  // Changed to 1-based indexing
        {
            for (int j = 1; j <= mNumCols; j++)
            {
                result(i, j) = (*this)(i, j) - other(i, j);
            }
        }
        return result;
    }
    //scalar multiply
    Matrix Matrix::operator*(double scalar) const
    {
        Matrix result(mNumRows, mNumCols);
        for (int i = 1; i <= mNumRows; i++)  // Changed to 1-based indexing
        {
            for (int j = 1; j <= mNumCols; j++)
            {
                result(i, j) = (*this)(i, j) * scalar;
            }
        }
        return result;
    }
    //multiply
    Vector Matrix::operator*(const Vector& vec) const
    {
        assert(mNumCols == vec.size());
        Vector result(mNumRows);
        for (int i = 1; i <= mNumRows; i++)
        {
            double sum = 0.0;
            for (int j = 1; j <= mNumCols; j++)
            {
                sum += (*this)(i, j) * vec(j);
            }
            result(i) = sum;
        }
        return result;
    }
    //print func
    void Matrix::print() const
    {
        for (int i = 0; i < mNumRows; i++)
        {
            for (int j = 0; j < mNumCols; j++)
            {
                cout << mData[i][j] << " ";
            }
            cout << endl;
        }
    }
    //inverse
    [[nodiscard]] Matrix Matrix::inverse() const
    {
        assert(mNumRows == mNumCols);
        int n = mNumRows;

        Matrix A(*this);
        Matrix I(n, n);
        for (int i = 1; i <= n; i++)
        {
            I(i, i) = 1.0;
        }
        //Gauss-Jordan elimination
        for (int k = 1; k <= n; k++)
        {
            double pivot = A(k, k);
            assert(pivot != 0); //no zero pivot
            for (int j = 1; j <= n; j++)
            {
                A(k, j) /= pivot;
                I(k, j) /= pivot;
            }

            for (int i = 1; i <= n; i++)
            {
                if (i != k)
                {
                    double factor = A(i, k);
                    for (int j = 1; j <= n; j++)
                    {
                        A(i, j) -= factor * A(k, j);
                        I(i,j) -= factor * I(k, j);
                    }
                }
            }
        } return I;
    }
    //for Matrix * Matrix overload
    Matrix Matrix::operator*(const Matrix& other) const
    {
        assert(mNumCols == other.mNumRows);
        Matrix result(mNumRows, other.mNumCols);
        for (int i = 1; i <= mNumRows; i++)
        {
            for (int j = 1; j <= other.mNumCols; j++)
            {
                double sum = 0.0;
                for (int k = 1; k <= mNumCols; k++)
                {
                    sum += (*this)(i, k) * other(k, j);
                }
                result(i, j) = sum;
            }
        }
        return result;
    }
    //pseudo-inverse(Moore-Penrose inverse)
    [[nodiscard]] Matrix Matrix::pseudoInverse() const
    {
        Matrix At(mNumCols, mNumRows);  //tranpose
        for (int i = 1; i <= mNumRows; i++)
        {
            for (int j = 1; j<=mNumCols; j++)
            {
                At(j, i) = (*this)(i, j);
            }
        }
        Matrix AtA = At * (*this);
        Matrix AtA_inv = AtA.inverse();

        return AtA_inv * At;
    }