Mat class

[afdave14]

I have never seen this before. I get no compiling errors but when I run this, the values in the Matrix are "nan". Here's the code (I've omitted some code to make it short enough. Let me know if you need some code):

Mat.h

#ifndef MAT_H
#define MAT_H

class IndexError{};
class SizeError{};

template<class D>
class Mat
{
public:
    Mat(const Mat< D > & num)
        : _rowsize(num.rowsize()), _colsize(num.colsize())
    {
        allocate_arrays();

        for (int i = 0; i < _rowsize; i++)
        {
            for (int j = 0; j < _colsize; j++)
            {
                _p[i][j] = num(i, j);
            }
        }
    }
    Mat(int rowsize=1, int colsize=1)
        :_rowsize(rowsize), _colsize(colsize)
    {
        allocate_arrays();

        for (int i = 0; i < _rowsize; i++)
        {
            for (int j = 0; j < _colsize; j++)
            {
                _p[i][j] = 0;
            }
        }
    }
    Mat(int rowsize, int colsize, D x[])
    {
        _rowsize = rowsize;
        _colsize = colsize;

        allocate_arrays();

        int i = 0, j = 0;
        while (i < _rowsize)
        {
            for (int k = 0; k < _colsize; k++, j++)
            {
                _p[i][k] = x[j];
            }
            i++;
        }
    }
    ~Mat()
    {
        for (int i = 0; i < _rowsize; i++)
            delete [] _p[i];

        delete [] _p;
    }
    Mat & operator=(const Mat< D > &);
    D & operator()(const int, const int) const;
    Mat & operator+=(const Mat< D > &);
    Mat operator+(const Mat< D > &) const;
    Mat & operator-=(const Mat< D > &);
    Mat operator-(const Mat< D > &) const;
    Mat operator*(const Mat< D > &) const;
    Mat operator*(const D &) const;
    void insert(int, int, D) const;
    int rowsize() const { return _rowsize; }
    int colsize() const { return _colsize; }
    //D get_value(const int, const int) const;
    Mat inv() const;
private:
    D ** _p;
    int _rowsize;
    int _colsize;
    void allocate_arrays()
    {
        _p = new D*[_rowsize];
        for (int i = 0; i < _rowsize; i++)
            _p[i] = new D[_colsize];
    }
};




template< class D >
std::ostream & operator<<(std::ostream & cout, const Mat< D > & matrix)
{
    ...
}

template< class D >
Mat< D > operator*(int a, const Mat< D > & b)
{
    Mat< D > prod(b);

    for (int i = 0; i < prod.rowsize(); i++)
    {
        for (int j = 0; j < b.colsize(); j++)
        {
            prod(i, j) = (prod(i, j) * a);
        }
    }

    return prod;
}

template< class D >
int len(const Mat< D > & num)
{
    ...
}

template< class D >
D & Mat<D>::operator()(const int row_index, const int col_index) const
{
    if (row_index >= _rowsize || col_index >= _colsize) { std::cout << "I FOUND IT" << std::endl; throw IndexError(); }

    return (_p[row_index][col_index]);
}

template< class D >
Mat< D > & Mat<D>::operator=(const Mat< D > & num)
{
    ...
}

template< class D >
Mat< D > & Mat<D>::operator+=(const Mat< D > & num)
{
    if (_rowsize != num.rowsize()) throw SizeError();
    for (int i = 0; i < _rowsize; i++)
    {
        for (int j = 0; j < _colsize; j++)
            this->_p[i][j] += num(i, j);
    }

    return *this;
}

template< class D >
Mat< D > Mat<D>::operator+(const Mat< D > & num) const
{
    if (_rowsize != num.rowsize()) throw SizeError();
    Mat sum(_rowsize, _colsize);
    
    return(sum += num);
}

template< class D >
Mat< D > & Mat<D>::operator-=(const Mat< D > & num)
{
    if (_rowsize != num.rowsize()) throw SizeError();
    for (int i = 0; i < _rowsize; i++)
    {
        for (int j = 0; j < _colsize; j++)
            this->_p[i][j] = this->_p[i][j] - num(i, j);
    }

    return *this;
}

template< class D >
Mat< D > Mat<D>::operator-(const Mat< D > & num) const
{
    if (_rowsize != num.rowsize()) throw SizeError();
    Mat sum(_rowsize, _colsize);
    
    return (sum -= num);
}

template< class D >
Mat< D > Mat<D>::operator*(const Mat< D > & num) const
{
    if (_colsize != num.rowsize()) throw SizeError();
    
    Mat prod(_rowsize, num.colsize());

    for (int i = 0; i < prod.rowsize(); i++)
    {
        for (int j = 0; j < prod.colsize(); j++)
        {
            for (int k = 0; k < _colsize; k++)
            {
                prod(i, j) += _p[i][k] * num(k, j);
            }
        }
    }

    return prod;
}

template< class D >
Mat< D > Mat<D>::operator*(const D & num) const
{
    for (int i = 0; i < _rowsize; i++)
    {
        for (int j = 0; j < _colsize; j++)
            _p[i][j] *= num;
    }

    return *this;
}

template< class D >
void Mat< D >::insert(int row_index, int col_index, D value) const
{
    ...
}

template< class D >
Mat< D > Mat<D>::inv() const
{
    Mat< D > inverse(_rowsize, _colsize);
    int order = this->rowsize();
    double det = 1.0 / CalcDeterminant(*this, order);
    
    Mat< D > minor(order - 1, order - 1);
    
    for (int j = 0; j < order; j++)
    {
        for (int i = 0; i < order; i++)
        {
            GetMinor(*this, minor, j, i, order);
            inverse(i, j) = det * CalcDeterminant(minor, order - 1);
            if ((i + j) % 2 == 1)
            {
                inverse(i, j) = -inverse(i, j);
            }
        }
    }
    std::cout << inverse << '\n' << std::endl;
    return inverse;
}

template< class D >
double det(const Mat< D > & mat)
{
    // order must be >= 0
    // stop the recursion when matrix is a single element
    int order = mat.rowsize();
    
    if( order == 1 )
        return mat(0, 0);
    
    // the determinant value
    double detr = 0;
    
    // allocate the cofactor matrix
    Mat< D > minor(order - 1, order - 1);
    
    for(int i = 0; i < order; i++ )
    {
        // get minor of element (0,i)
        GetMinor(mat, minor, 0, i, order);
        
        // the recusion is here!
        detr += (i % 2 == 1 ? -1.0 : 1.0) * mat(0, i) * CalcDeterminant(minor, order-1);
    }
    
    return detr;
}


// matrix inversioon
// the result is put in Y
template< class D >
void MatrixInversion(Mat< D > A, int order, Mat< D > Y)
{
    // get the determinant of a
    double det = 1.0 / CalcDeterminant(A,order);
    
    // memory allocation
    int * temp = new int[(order - 1) * (order - 1)];
    Mat< D > minor(order - 1, order - 1);
    
    for(int j = 0; j < order; j++)
    {
        for(int i = 0; i < order; i++)
        {
            // get the co-factor (matrix) of A(j,i)
            GetMinor(A, minor, j, i, order);
            Y(i, j) = det * CalcDeterminant(minor, order - 1);
            
            if((i + j) % 2 == 1)
                Y(i, j) = -Y(i, j);
        }
    }
    
    // release memory

    delete [] temp;
}


// calculate the cofactor of element (row,col)
template< class D >
int GetMinor(Mat< D > src, Mat< D > dest, int row, int col, int order)
{
    // indicate which col and row is being copied to dest
    int colCount = 0, rowCount = 0;
    
    for(int i = 0; i < order; i++ )
    {
        if( i != row )
        {
            colCount = 0;
            
            for(int j = 0; j < order; j++ )
            {
                // when j is not the element
                if( j != col )
                {
                    dest(rowCount, colCount) = src(i, j);
                    colCount++;
                }
            }
            rowCount++;
        }
    }
    return 1;
}


// Calculate the determinant recursively.
template< class D >
double CalcDeterminant(Mat< D > mat, int order)
{
    // order must be >= 0
    // stop the recursion when matrix is a single element
    if( order == 1 )
        return mat(0, 0);
    
    // the determinant value
    double det = 0;
    
    // allocate the cofactor matrix
    Mat< D > minor(order - 1, order - 1);

    for(int i = 0; i < order; i++ )
    {
        // get minor of element (0,i)
        GetMinor(mat, minor, 0, i, order);
        
        // the recusion is here!
        det += (i % 2 == 1 ? -1.0 : 1.0) * mat(0, i) * CalcDeterminant(minor, order-1);
    }
    std::cout << det << std::endl;
    return det;
}

#endif 

[afdave14]

And here's main.cpp:

#include <iostream>
#include "Mat.h"

int main()
{
    double x[] = {1.23, 2.34, 3.45, 4.56};//, 5.67, 6.78, 7.89, 8.91, 9.11};
    double y[] = {6.78, 5.67, 4.56, 3.45};//, 2.34, 1.23};
    Mat< double > a(2, 2, x);

    Mat< double > b(2, 2);
    b = a.inv();
    std::cout << b << std::endl;
    
    Mat< double >c = a * b;
    std::cout << '\n' << c << std::endl;

    return 0;
}

[afdave14]

Well, I think I figured it out. I forgot to make the dest Matrix in the GetMinor function a reference, and so when I changed it in the GetMinor function it didn't do jack. (I was wondering why I kept getting 0). Values aren't exact, though, unfortunately. Any comments would be nice. I'm sure somebody would like to give some constructive criticism and I'm all ears

[tition]

Some comments on it:
1) initializing a matrix with zero rows and zero columns is much better. No memory allocation in matrix creation = much more efficient.
2) Your way of computing determinants is the most inefficient standard method. On each recursion step, you do n recursive calls, where n is the size of the matrix. This means for an n by n matrix you will be doing ~n! operations - your code will die on a 20 by 20 matrix.

Instead, use gaussian elimination:

http://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode

In order to calculate the determinant using Gaussian elimination, initialize

D Determinant=1;
and do
Determinant*=A[i,j];
right before the step "divide each entry in row i by A[i,j]" in the algorithm.

3) Matrix inverses are also calculated using Gaussian elimination. I can explain that too, but gotta go right now: you can look it up yourself too:
http://en.wikipedia.org/wiki/Gaussian_elimination#Finding_the_inverse_of_a_matrix
If you have trouble with the implementation, I will gladly try to help.



I also have an educational question to fill in my c++ gaps: can you explain the syntax

    Mat(const Mat< D > & num)
        : _rowsize(num.rowsize()), _colsize(num.colsize())

, also, what is its purpose?