matrix class

I am trying to maximise the efficiency of my matrix classes. There's a matrix defined at each grid point i. The elements of the matrix at all grid points are are stored inside ptr00, ptr10 and ptr20. Any adivce?

ie for a 3 x 1 matrix

Code:

---

class matrix_col {

 public:

  int i;
  double *ptr00, *ptr10, *ptr20;

  matrix_col(){

    ptr00 = new double[Grids];
    ptr10 = new double[Grids];
    ptr20 = new double[Grids];

    for(i = 0; i < (Grids); i++){

      ptr00[i] = 0;
      ptr10[i] = 0;
      ptr20[i] = 0;
    }
  }

  //copy constructor
  matrix_col(const matrix_col &a){

    int i;
 
    ptr00 = new double[Grids];
    ptr10 = new double[Grids];
    ptr20 = new double[Grids];

    for(i = 0; i < (Grids); i++){

      ptr00[i] = a.ptr00[i];
      ptr10[i] = a.ptr10[i];
      ptr20[i] = a.ptr20[i];
    }
  }

  ~matrix_col(){

    delete [] ptr00;
    delete [] ptr10;
    delete [] ptr20;

  }

  friend matrix_col operator+(const matrix_col &, const matrix_col &);

};


// Overload + for ob1 + ob2.

matrix_col operator+(const matrix_col &left, const matrix_col &right)
{

  int i;
  matrix_col temp;

  for (i = 0; i < Grids; i++){

    temp.ptr00[i] = left.ptr00[i] + right.ptr00[i];
    temp.ptr10[i] = left.ptr10[i] + right.ptr10[i];
    temp.ptr20[i] = left.ptr20[i] + right.ptr20[i];
 
 }

  return temp;
}

---

Have you considered using std::vector<double> instead of double* for the member variables? Your member variable i looks out of place, perhaps it should not be a member variable at all.

I don't really have any advice, except that I think those 3 pointers aren't that neat. You might want to check out the matrix class I wrote a while ago. I can't say it's great, but it's still neat I think. It's far from a final version but I got bored... Maybe I'll finish it someday...

Code:

---

#ifndef MATRIX_H_INCLUDED
#define MATRIX_H_INCLUDED

#include <iostream>
#include <vector>
#include <memory>
#include "vector.h"

template < class T >
class Matrix
{
public:
    Matrix(unsigned int n_cols = 0, unsigned int n_rows = 0);
        Matrix(const Matrix& m);

        void Output( )
        {
                for(int j = 0; j < Rows; j++)
                {
                        for(int i = 0; i < Cols; i++)
                        {
                                std::cout << Data[j][i];
                        }
                        std::cout << std::endl;
                }
        }

        void SetColElements(const std::vector< std::vector< T > >& v, unsigned int col);
        void SetRowElements(const std::vector< std::vector< T > >& v, unsigned int row);
        void SetElement(unsigned int col, unsigned int row, T val);

        void InterchangeRows(unsigned int row_1, unsigned int row_2);
        void InterchangeCols(unsigned int col_1, unsigned int col_2);

        bool LowerTriangle( ) const;
        bool UpperTriangle( ) const;
        bool DiagonalMatrix( ) const;
    bool SquareMatrix( ) const;

        T Determinant( );

        void GaussElim(const Vector< T >&);
        void GaussJordanElim(const Vector< T >&);

        static Matrix< T > IdentityMatrix(unsigned int size);
        static Matrix< T > NullMatrix(unsigned int size);

    Matrix< T > Minor(unsigned int col, unsigned int row);

        Matrix< T >& operator = (const Matrix< T >&);
        Matrix< T >& operator += (const Matrix< T >&);
        Matrix< T >& operator -= (const Matrix< T >&);
        Matrix< T > operator * (const Matrix< T >&) const;
        Matrix< T > operator * (const T&) const;
        Matrix< T > operator *= (const T&);
        Matrix< T > operator - (const Matrix< T >&) const;
        Matrix< T > operator + (const Matrix< T >&) const;
private:
        unsigned int Rows, Cols;
    std::vector< std::vector< T > > Data;
};

template < class T >
Matrix< T >::Matrix< T >(unsigned int n_cols, unsigned int n_rows)
        : Cols(n_cols), Rows(n_rows)
{
        if(Cols > 0 && Rows > 0)
        {
                Data.reserve(Rows);
                for(int i = 0; i < Rows; i++)
                {
                        Data.push_back(std::vector< T >( ));

                        Data[i].reserve(Cols);
                        Data[i].insert(Data[i].begin(), Cols, 0);
                }
        }
        else if(Cols < 0 && Rows < 0)
        {
                Cols = Rows = 0;
        }
}

template < class T >
Matrix< T >::Matrix< T >(const Matrix& m)
{
        *this = m;
}

template < class T >
void Matrix< T >::SetColElements(const std::vector< std::vector< T > >& v, unsigned int col)
{
        if(ls.size() == Rows && Cols != 0 && Rows != 0 && col > 0 && col < Cols)
        {
                for(int j = 0; j < Rows; j++)
                {
                        Data[j][col] = v[j];
                }
        }
}

template < class T >
void Matrix< T >::SetRowElements(const std::vector< std::vector< T > >& v, unsigned int row)
{
        if(ls.size() == Rows && Cols != 0 && Rows != 0 && row > 0 && row < Rows)
        {
                for(int i = 0; i < Cols; i++)
                {
                        Data[row][i] = v[i];
                }
        }
}

template < class T >
void Matrix< T >::SetElement(unsigned int col, unsigned int row, T val)
{
        if(col >= 0 && row >= 0 && col < Cols && row < Rows)
        {
                Data[row][col] = val;
        }
}

template < class T >
void Matrix< T >::InterchangeRows(unsigned int row_1, unsigned int row_2)
{
        if(row_1 != row_2 && (row_1 >= 0 && row_2 >= 0))
        {
                Data[row_1].swap(Data[row_2]);
        }
}

template < class T>
void Matrix< T >::InterchangeCols(unsigned int col_1, unsigned int col_2)
{
        // Parentheses for easier reading only
        // Make sure both columns are different and range in [0, Cols[
        if((col_1 != col_2) && (col_1 >= 0) && (col_2 >= 0) && (col_1 < Cols) && (col_2 < Cols))
        {
                for(int j = 0; j < Rows; j++)
                {
                        T tmp = Data[j][col_1];
                        Data[j][col_1] = Data[j][col_2];
                        Data[j][col_2] = tmp;
                }
        }
}

template < class T >
bool Matrix< T >::LowerTriangle( ) const
{
        // Only square matrixes bigger than 0,0
        if(SquareMatrix() && (Cols != 0 && Rows != 0))
        {
                for(int j = 0; j < Rows; j++)
                {
                        for(int i = 0; i < Cols; i++)
                        {
                                if(i > j && Data[j][i] != 0)
                                {
                                        return false;
                                }
                        }
                }
        }
        else
        {
                return false;
        }
        return true;
}

template < class T >
bool Matrix< T >::UpperTriangle( ) const
{
        // Only square matrixes bigger than 0,0
        if(SquareMatrix() && (Cols != 0 && Rows != 0))
        {
                for(int j = 0; j < Rows; j++)
                {
                        for(int i = 0; i < Cols; i++)
                        {
                                if(i < j && Data[j][i] != 0)
                                {
                                        return false;
                                }
                        }
                }
        }
        else
        {
                return false;
        }
        return true;
}

template < class T >
bool Matrix< T >::DiagonalMatrix( ) const
{
        return (UpperTriangle( ) && LowerTriangle( ));
}

template < class T >
bool Matrix< T >::SquareMatrix( ) const
{
    return (Rows == Cols);
}

template < class T >
T Matrix< T >::Determinant( )
{
        if(SquareMatrix( ))
        {
                if(Cols == 0 && Rows == 0)
                {
                        return static_cast< T >(0);
                }
                if(Cols == 1 && Rows == 1)
                {
                        return Data[0][0];
                }

                if(UpperTriangle( ) || LowerTriangle( ))
                {
                        int total = 1;
                        for(int i = 0; i < Cols; i++)
                        {
                                total *= Data[i][i];
                        }
                        return total;
                }

                if(Cols == 2 && Rows == 2)
                {
                        return (Data[0][0] * Data[1][1] - Data[0][1] * Data[1][0]);
                }
                if(Cols == 3 && Rows == 3)
                {
                        return (Data[0][0] * Data[1][1] * Data[2][2] +
                                        Data[0][2] * Data[1][0] * Data[2][1] +
                                        Data[0][1] * Data[1][2] * Data[2][0] -
                                        Data[0][2] * Data[1][1] * Data[2][0] -
                                        Data[0][0] * Data[1][2] * Data[2][1] -
                                        Data[0][1] * Data[1][0] * Data[2][2]);
                }
                // Let's use the M = LU decomposition method for faster processing
                else
                {
                       
                }
        }
}

template < class T >
void Matrix< T >::GaussElim(const Vector< T >& solution)
{

}

template < class T >
void Matrix< T >::GaussJordanElim(const Vector< T >& solution)
{

}

template < class T >
Matrix< T > Matrix< T >::IdentityMatrix(unsigned int size)
{
        Matrix< T > tmp(size, size);

        for(int j = 0; j < tmp.Rows; j++)
        {
                for(int i = 0; i < tmp.Cols; i++)
                {
                        tmp.Data[j][i] = i == j ? 1 : 0;
                }
        }

        return tmp;
}

template < class T >
Matrix< T > Matrix< T >::NullMatrix(unsigned int size)
{
        return Matrix< T >(size, size);
}

template < class T >
Matrix< T > Matrix< T >::Minor(unsigned int col, unsigned int row)
{
        Matrix< T > tmp(col, row);

        if(SquareMatrix( ) && (Rows != 0 && Cols != 0))
        {
                int y = 0;

                // Move horizontally then down in the matrix
                for(int j = 0; j < Rows; j++)
                {

                        if(j != row)
                        {
                                // We're on a new line...
                                tmp.Data.push_back(std::vector< T >( ));
                                for(int i = 0; i < Cols; i++)
                                {
                                        if(i != col)
                                        {
                                                tmp.Data[y].push_back(Data[j][i]);
                                        }
                                }
                                y++;
                        }
                }
                tmp.Cols = tmp.Data[0].size();
                tmp.Rows = tmp.Data.size();
        }

        return tmp;
}

template < class T >
Matrix< T >& Matrix< T >::operator = (const Matrix< T >& m)
{
        Rows = m.Rows;
        Cols = m.Cols;

        Data.assign(m.Data.begin(), m.Data.end());

        return *this;
}

template < class T >
Matrix< T >& Matrix< T >::operator += (const Matrix< T >& m)
{
        *this = *this + m;
        return *this;
}

template < class T >
Matrix< T >& Matrix< T >::operator -= (const Matrix< T >& m)
{
        *this = *this - m;
        return *this;
}

template < class T >
Matrix< T > Matrix< T >::operator * (const Matrix< T >& m) const
{
        // A Matrix A(a, b) and a matrix B(c, d) yield a matric C(b, c) when multiplied
        Matrix< T > tmp(Rows, m.Cols);

        if(Cols == m.Rows && Rows == m.Cols && Rows > 0 && Cols > 0)
        {
                for(int j_result = 0; j_result < Rows; j_result++)
                {
                        for(int i_result = 0; i_result < m.Cols; i_result++)
                        {
                                int sum = 0;
                                for(int j = 0, i = 0; j < m.Rows && i < Cols; j++, i++)
                                {
                                        sum += Data[j_result][i] * m.Data[j][i_result];
                                }
                                tmp.Data[j_result][i_result] = sum;
                        }
                }
        }

        return tmp;
}

template < class T >
Matrix< T > Matrix< T >::operator * (const T& val) const
{
        Matrix< T > tmp(Cols, Rows);
        tmp = *this;

        if(Rows > 0 && Cols > 0)
        {
                for(int j = 0; j < Rows; j++)
                {
                        for(int i = 0; i < Cols; i++)
                        {
                                tmp.Data[j][i] *= val;
                        }
                }
        }

        return tmp;
}

template < class T >
Matrix< T > Matrix< T >::operator *= (const T& val)
{
        // Slight overhead due to creation of temporary variable
        *this = *this * val;
        return *this;
}

template < class T >
Matrix< T > Matrix< T >::operator - (const Matrix< T >& m) const
{
        Matrix< T > tmp(Cols, Rows);

        if(Rows == m.Rows && Cols == m.Cols && Rows > 0 && Cols > 0)
        {
                for(int j = 0; j < Rows; j++)
                {
                        for(int i = 0; i < Cols; i++)
                        {
                                tmp.Data[j][i] = Data[j][i] - m.Data[j][i];
                        }
                }
        }

        return tmp;
}

template < class T >
Matrix< T > Matrix< T >::operator + (const Matrix< T >& m) const
{
        Matrix< T > tmp(Cols, Rows);

        if(Rows == m.Rows && Cols == m.Cols && Rows > 0 && Cols > 0)
        {
                for(int j = 0; j < Rows; j++)
                {
                        for(int i = 0; i < Cols; i++)
                        {
                                tmp.Data[j][i] = Data[j][i] + m.Data[j][i];
                        }
                }
        }

        return tmp;
}

#endif // MATRIX_H_INCLUDED

---