7 Replies - 56654 Views - Last Post: 06 February 2012 - 05:36 PM Rate Topic: -----

#1 krish_iitd  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 2
  • Joined: 24-June 08

determinant of matrix in c++

Posted 24 June 2008 - 11:37 AM

:^: hello , i am interested to find the simple function to find the determinant of square matrix in c++ language as possible as simple ,can any one Please help me .
Is This A Good Question/Topic? 0
  • +

Replies To: determinant of matrix in c++

#2 joske  Icon User is offline

  • D.I.C Regular
  • member icon

Reputation: 43
  • View blog
  • Posts: 297
  • Joined: 04-September 07

Re: determinant of matrix in c++

Posted 24 June 2008 - 12:40 PM

Here is some code i wrote some time ago as an exercise. It contains a simple matrix class, including function det() to calculate the determinant. For 2x2 matrices the determinant is very easy, for larger matrices it is a little bit more complex, amongst others you have to calculate all minor of the matrix.

It is not a neat C++ program, as I haven't split the code in a header and source file.
/*
A simple matrix class
c++ code
Author: Jos de Jong, Nov 2007

With this class you can:
    - create a matrix
    - get/set the values
    - use operators +, -, *
    - use functions ones(), zeros(), diag(), det(), inv(), size()
    - print the content of the matrix
    
Usage:
    you can create a matrix by:
        Matrix A;
        Matrix A = Matrix(rows, cols);
        Matrix A = B;

    you can get and set matrix elements by:
        A(2,3) = 5.6;        // set an element of Matix A
        value = A(3,1);      // get an element of Matrix A
        A = B;               // copy content of Matrix B to Matrix A

    you can apply operations with matrices and doubles:
        A = B + C;
        A = B - C;
        A = B * C;
        A = -B;

    the following functions are available:
        A = ones(rows, cols);
        A = zeros(rows, cols);
        A = diag(n);
        A = diag(B)/>;
        d = det(A);
        A = inv(B)/>;
        cols = A.get_cols();
        rows = A.get_rows();
        cols = size(A, 1);
        rows = size(A, 2);
        
    you can quick-print the content of a matrix in the console with:
        A.print();

Todo:
    inverse of a matrix. see [url="http://mathworld.wolfram.com/MatrixInverse.html"]http://mathworld.wolfram.com/MatrixInverse.html[/url]
    division with a matrix
    multiplying a matrix with a double
*/

#include <cstdlib>
#include <cstdio>
#include <math.h>

using namespace std;


/* 
 * a simple exception class
 * you can create an exeption by entering
 *      throw Exception("...Error description...");
 * and get the error message from the data msg for displaying:
 *      err.msg
 */
class Exception {
   public:
      const char* msg;
      Exception(const char* arg) : msg(arg) { };
};


class Matrix
{
public:
    // constructors
    Matrix()
    {
        //printf("Executing constructor Matrix() ...\n");
        // create a Matrix object without content
        p = NULL;
        rows = 0;
        cols = 0;
    }
    Matrix(const int row_count, const int column_count)
    {
        // create a Matrix object with given number of rows and columns
        //printf("Executing constructor Matrix(const int column_count, const int row_count) ...\n");
        p = NULL;

        if (row_count > 0 && column_count > 0)
        {
            rows = row_count;
            cols = column_count;

            p = new double*[rows];
            for (int r = 0; r < rows; r++)
            {
                p[r] = new double[rows];

                // initially fill in zeros for all values in the matrix;
                for (int c = 0; c < cols; c++)
                {
                    p[r][c] = 0;
                }
            }
        }
    }
    // assignment operator
    Matrix(const Matrix& a)
    {
        //printf("Executing constructor Matrix(const Matrix& a) ...\n");
        rows = a.rows;
        cols = a.cols;
        p = new double*[a.rows];
        for (int r = 0; r < a.rows; r++)
        {
            p[r] = new double[a.cols];

            // copy the values from the matrix a
            for (int c = 0; c < a.cols; c++)
            {
                p[r][c] = a.p[r][c];
            }
        }
    }

    // index operator. You can use this class like myMatrix(col, row)
    // the indexes are one-based, not zero based.
	double& operator()(const int r, const int c)
	{
        //printf("Executing Matrix operator() ...\n");
      	if (p != NULL && r > 0 && r <= rows && c > 0 && c <= cols)
      	{
          	return p[r-1][c-1];
        }
      	else
      	{
            throw Exception("Subscript out of range");
        }
    }

    // assignment operator
    Matrix& operator= (const Matrix& a)
    {
        //printf("Executing Matrix operator= ...\n");
        rows = a.rows;
        cols = a.cols;
        p = new double*[a.rows];
        for (int r = 0; r < a.rows; r++)
        {
            p[r] = new double[a.cols];

            // copy the values from the matrix a
            for (int c = 0; c < a.cols; c++)
            {
                p[r][c] = a.p[r][c];
            }
        }
        return *this;
    }

    // operator addition
    friend Matrix operator+(const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.rows == b.rows && a.cols == b.cols)
        {
         	Matrix res(a.rows, a.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c = 0; c < a.cols; c++)
                {
                    res.p[r][c] = a.p[r][c] + b.p[r][c];
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    // add a double value (elements wise)
    Matrix& add(double v)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                p[r][c] += v;
            }
        }
     	return *this;
    }

    // subtract a double value (elements wise)
    Matrix& subtract(double v)
    {
        return add(-v);
    }

    // multiply a double value (elements wise)
    Matrix& multiply(double v)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                p[r][c] *= v;
            }
        }
     	return *this;
    }

    // operator subtraction
    friend Matrix operator- (const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.rows == b.rows && a.cols == b.cols)
        {
         	Matrix res(a.rows, a.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c = 0; c < a.cols; c++)
                {
                    res.p[r][c] = a.p[r][c] - b.p[r][c];
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }
    // operator unary minus
    friend Matrix operator- (const Matrix& a)
    {
     	Matrix res(a.rows, a.cols);

        for (int r = 0; r < a.rows; r++)
        {
            for (int c = 0; c < a.cols; c++)
            {
                res.p[r][c] = -a.p[r][c];
            }
        }
    
        return res;
    }
    
    // operator multiplication
    friend Matrix operator* (const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.cols == b.rows)
        {
         	Matrix res(a.rows, b.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c_res = 0; c_res < b.cols; c_res++)
                {
                    for (int c = 0; c < a.cols; c++)
                    {
                        res.p[r][c_res] += a.p[r][c] * b.p[c][c_res];
                    }
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    /**
     * returns the minor from the given matrix where
     * the selected row and column are removed
     */
    Matrix minor(const int row, const int col)
    {
        //printf("minor(row=%i, col=%i), rows=%i, cols=%i\n", row, col, rows, cols);
        Matrix res;
      	if (row > 0 && row <= rows && col > 0 && col <= cols)
      	{
            res = Matrix(rows - 1, cols - 1);
            
            // copy the content of the matrix to the minor, except the selected
            for (int r = 1; r <= (rows - (row >= rows)); r++)
            {
                for (int c = 1; c <= (cols - (col >= cols)); c++)
                {
                    //printf("r=%i, c=%i, value=%f, rr=%i, cc=%i \n", r, c, p[r-1][c-1], r - (r > row), c - (c > col));
                    res(r - (r > row), c - (c > col)) = p[r-1][c-1];
                }
            }        
        }
        else
        {
            throw Exception("Index for minor out of range");
        }
        
        return res;
    }
    
    /* 
     * returns the size of the i-th dimension of the matrix.
     * i.e. for i=1 the function returns the number of rows, 
     * and for i=2 the function returns the number of columns
     * else the function returns 0
     */
    int size(const int i)
    {
        if (i == 1)
        {
            return rows;
        }
        else if (i == 2)
        {
            return cols;
        }
        return 0;
    }

    // returns the number of rows
    int get_rows()
    {
        return rows;
    }
    
    // returns the number of columns
    int get_cols()
    {
        return cols;
    }

    // print the contents of the matrix
    void print()
    {
        if (p != NULL)
        {
            printf("[");
            for (int r = 0; r < rows; r++)
            {
                if (r > 0)
                {
                    printf(" ");
                }
                for (int c = 0; c < cols-1; c++)
                {
                    printf("%.2f, ", p[r][c]);
                }
                if (r < rows-1)
                {
                    printf("%.2f;\n", p[r][cols-1]);
                }
                else
                {
                    printf("%.2f]\n", p[r][cols-1]);
                }
            }
        }
        else
        {
            // matrix is empty
            printf("[ ]\n");
        }
    }

public:
    // destructor
    ~Matrix()
    {
        // clean up allocated memory
        for (int r = 0; r < rows; r++)
        {
            delete p[r];
        }
        delete p;
        p = NULL;
    }

private:
    int rows;
    int cols;
    double** p;         // pointer to a matrix with doubles
};

/* 
 * i.e. for i=1 the function returns the number of rows,
 * and for i=2 the function returns the number of columns
 * else the function returns 0
 */
int size(Matrix& a, const int i)
{
    return a.size(i);
}


// addition of Matrix with double
Matrix operator+ (const Matrix& a, const double B)/>
{
    Matrix res = a;
    res.add(B)/>;
    return res;
}
// addition of double with Matrix
Matrix operator+ (const double b, const Matrix& a)
{
    Matrix res = a;
    res.add(B)/>;
    return res;
}

// subtraction of Matrix with double
Matrix operator- (const Matrix& a, const double B)/>
{
    Matrix res = a;
    res.subtract(B)/>;
    return res;
}
// subtraction of double with Matrix
Matrix operator- (const double b, const Matrix& a)
{
    Matrix res = -a;
    res.add(B)/>;
    return res;
}

// multiplication of Matrix with double
Matrix operator* (const Matrix& a, const double B)/>
{
    Matrix res = a;
    res.multiply(B)/>;
    return res;
}
// multiplication of double with Matrix
Matrix operator* (const double b, const Matrix& a)
{
    Matrix res = a;
    res.multiply(B)/>;
    return res;
}


/**
 * returns a matrix with size cols x rows with ones as values
 */
Matrix ones(const int rows, const int cols)
{
    Matrix res = Matrix(rows, cols);

    for (int r = 1; r <= rows; r++)
    {
        for (int c = 1; c <= cols; c++)
        {
            res(r, c) = 1;
        }
    }
    return res;
}

/**
 * returns a matrix with size cols x rows with zeros as values
 */
Matrix zeros(const int rows, const int cols)
{
    return Matrix(rows, cols);
}


/**
 * returns a diagonal matrix with size n x n with ones at the diagonal
 * @param   v a vector
 * @return  a diagonal matrix with ones on the diagonal
 */
Matrix diag(const int n)
{
    Matrix res = Matrix(n, n);
    for (int i = 1; i <= n; i++)
    {
        res(i, i) = 1;
    }
    return res;
}
/**
 * returns a diagonal matrix
 * @param   v a vector
 * @return  a diagonal matrix with the given vector v on the diagonal
 */
Matrix diag(Matrix& v)
{
    Matrix res;
    if (v.get_cols() == 1)
    {
        // the given matrix is a vector n x 1
        int rows = v.get_rows();
        res = Matrix(rows, rows);
        
        // copy the values of the vector to the matrix
        for (int r=1; r <= rows; r++)
        {
            res(r, r) = v(r, 1);
        }
    }
    else if (v.get_rows() == 1)
    {
        // the given matrix is a vector 1 x n
        int cols = v.get_cols();
        res = Matrix(cols, cols);

        // copy the values of the vector to the matrix
        for (int c=1; c <= cols; c++)
        {
            res(c, c) = v(1, c);
        }        
    }
    else
    {
        throw Exception("Parameter for diag must be a vector");
    }
    return res;
}

/*
 * returns the determinant of Matrix a
 */
double det(Matrix& a)
{
    double d = 0;       // value of the determinant
    int rows = a.get_rows();
    int cols = a.get_cols();
    
    if (rows == cols)
    {
        // this is a square matrix
        if (rows == 1)
        {
            // this is a 1 x 1 matrix
            d = a(1, 1);
        }
        else if (rows == 2)
        {
            // this is a 2 x 2 matrix
            // the determinant of [a11,a12;a21,a22] is det = a11*a22-a21*a12
            d = a(1, 1) * a(2, 2) - a(2, 1) * a(1, 2);
        }
        else
        {
            // this is a matrix of 3 x 3 or larger
            for (int c = 1; c <= cols; c++)
            {
                Matrix M = a.minor(1, c);
                //d += pow(-1, 1+c)  * a(1, c) * det(M);
                d += (c%2 + c%2 - 1) * a(1, c) * det(M);  // faster than with pow()
            }
        }
    }
    else
    {
        throw Exception("Matrix must be square");
    }
    return d;
}

// swap two values
void swap(double& a, double& B)/>
{
    double temp = a;
    a = b;
    b = temp;
}

/*
 * returns the inverse of Matrix a
 */
Matrix inv(Matrix& a)
{
    Matrix res;
    double d = 0;       // value of the determinant
    int rows = a.get_rows();
    int cols = a.get_cols();
    
    d = det(a);
    if (rows == cols && d != 0)
    {
        // this is a square matrix
        if (rows == 1)
        {
            // this is a 1 x 1 matrix
            res = Matrix(rows, cols);
            res(1, 1) = 1 / a(1, 1);
        }
        else if (rows == 2)
        {
            // this is a 2 x 2 matrix
            res = Matrix(rows, cols);
            res(1, 1) = a(2, 2);
            res(1, 2) = -a(1, 2);
            res(2, 1) = -a(2, 1);
            res(2, 2) = a(1, 1);
            res = (1/d) * res;
        }
        else
        {
            // this is a matrix of 3 x 3 or larger
            // calculate inverse using gauss-jordan elimination
            //      [url="http://mathworld.wolfram.com/MatrixInverse.html"]http://mathworld.wolfram.com/MatrixInverse.html[/url]
            //      [url="http://math.uww.edu/~mcfarlat/inverse.htm"]http://math.uww.edu/~mcfarlat/inverse.htm[/url]
            res = diag(rows);     // a diagonal matrix with ones at the diagonal
            Matrix ai = a;        // make a copy of Matrix a

            for (int c = 1; c <= cols; c++)
            {
                // element (c, c) should be non zero. if not, swap content 
                // of lower rows
                int r;
                for (r = c; r <= rows && ai(r, c) == 0; r++)
                {
                }                
                if (r != c)
                {
                    // swap rows
                    for (int s = 1; s <= cols; s++)
                    {
                        swap(ai(c, s), ai(r, s));
                        swap(res(c, s), res(r, s));
                    }
                }

                // eliminate non-zero values on the other rows at column c
                for (int r = 1; r <= rows; r++)
                {
                    if(r != c)
                    {
                        // eleminate value at column c and row r
                        if (ai(r, c) != 0)
                        {
                            double f = - ai(r, c) / ai(c, c);
                            
                            // add (f * row c) to row r to eleminate the value
                            // at column c
                            for (int s = 1; s <= cols; s++)
                            {
                                ai(r, s) += f * ai(c, s);
                                res(r, s) += f * res(c, s);
                            }
                        }                    
                    }
                    else
                    {
                        // make value at (c, c) one,
                        // divide each value on row r with the value at ai(c,c)
                        double f = ai(c, c);
                        for (int s = 1; s <= cols; s++)
                        {
                            ai(r, s) /= f;
                            res(r, s) /= f;
                        }                        
                    }
                }                
            }
        }
    }
    else
    {
        if (rows == cols)
        {
            throw Exception("Matrix must be square");
        }
        else
        {
            throw Exception("Determinant of matrix is zero");
        }
    }
    return res;
}


int main(int argc, char *argv[])
{
    try
    {
        // create an empty matrix of 3x3 (will initially contain zeros)
        int cols = 3;
        int rows = 3;
        Matrix A = Matrix(cols, rows);

        // fill in some values in matrix a
        int count = 0;
        for (int r = 1; r <= rows; r++)
        {
            for (int c = 1; c <= cols; c++)
            {
                count ++;
                A(r, c) = count;
            }
        }

        // adjust a value in the matrix  (indexes are one-based)
        A(2,1) = 1.23;

        // read a value from the matrix  (indexes are one-based)
        double centervalue = A(2,2);
        printf("centervalue = %f \n", centervalue);
        printf("\n");

        // print the whole matrix
        printf("A = \n");
        A.print();
        printf("\n");

        Matrix B = ones(rows, cols) + diag(rows);
        printf("B = \n");
        B.print();
        printf("\n");

        Matrix B2 = Matrix(rows, 1);        // a vector
        count = 1;
        for (int r = 1; r <= rows; r++)
        {
            count ++;
            B2(r, 1) = count * 2;
        }
        printf("B2 = \n");
        B2.print();
        printf("\n");
        
        Matrix C;
        C = A + B;
        printf("A + B = \n");
        C.print();
        printf("\n");

        C = A - B;
        printf("A - B = \n");
        C.print();
        printf("\n");

        C = A * B2;
        printf("A * B2 = \n");
        C.print();
        printf("\n");

        // create a diagonal matrix
        Matrix E = diag(B2);
        printf("E = \n");
        E.print();
        printf("\n");

        // calculate determinant
        Matrix D = Matrix(2, 2);
        D(1,1) = 2;
        D(1,2) = 4;
        D(2,1) = 1;
        D(2,2) = -2;
        printf("D = \n");
        D.print();
        printf("det(D) = %f\n\n", det(D));

        printf("A = \n");
        A.print();
        printf("\n");
        printf("det(A) = %f\n\n", det(A));
        
        Matrix F;
        F = 3 - A ;
        printf("3 - A = \n");
        F.print();
        printf("\n");
        
        // test inverse
        Matrix G = Matrix(2, 2);
        G(1, 1) = 1;
        G(1, 2) = 2;
        G(2, 1) = 3;
        G(2, 2) = 4;
        printf("G = \n");
        G.print();
        printf("\n");
        Matrix G_inv = inv(G);
        printf("inv(G) = \n");
        G_inv.print();
        printf("\n");
        
        Matrix A_inv = inv(A);
        printf("inv(A) = \n");
        A_inv.print();
        printf("\n");
    }
    catch (Exception err)
    {
        printf("Error: %s\n", err.msg);
    }
    catch (...)
    {
        printf("An error occured...\n");
    }

    printf("\n");
    system("PAUSE");
    return EXIT_SUCCESS;
}


Was This Post Helpful? 0
  • +
  • -

#3 joske  Icon User is offline

  • D.I.C Regular
  • member icon

Reputation: 43
  • View blog
  • Posts: 297
  • Joined: 04-September 07

Re: determinant of matrix in c++

Posted 20 March 2010 - 04:17 AM

Here some updated code.

Pete pointed me to an error in the contstructor (there was a line p[r] = new double[rows]; instead of p[r] = new double[cols];. Thanks!

Furthermore, some slight changes to make the code cross platform.

/*
A simple matrix class
c++ code
Author: Jos de Jong, Nov 2007. Updated March 2010

With this class you can:
    - create a matrix
    - get/set the values
    - use operators +, -, *
    - use functions ones(), zeros(), diag(), det(), inv(), size()
    - print the content of the matrix

Usage:
    you can create a matrix by:
        Matrix A;
        Matrix A = Matrix(rows, cols);
        Matrix A = B;

    you can get and set matrix elements by:
        A(2,3) = 5.6;        // set an element of Matix A
        value = A(3,1);      // get an element of Matrix A
        A = B;               // copy content of Matrix B to Matrix A

    you can apply operations with matrices and doubles:
        A = B + C;
        A = B - C;
        A = B * C;
        A = -B;

    the following functions are available:
        A = Ones(rows, cols);
        A = Zeros(rows, cols);
        A = Diag(n);
        A = Diag(B)/>;
        d = Det(A);
        A = Inv(B)/>;
        cols = A.GetCols();
        rows = A.GetRows();
        cols = Size(A, 1);
        rows = Size(A, 2);

    you can quick-print the content of a matrix in the console with:
        A.Print();

Todo:
    inverse of a matrix. see http://mathworld.wolfram.com/MatrixInverse.html
    division with a matrix
    multiplying a matrix with a double
*/

#include <cstdlib>
#include <cstdio>
#include <math.h>

#define PAUSE {printf("Press \"Enter\" to continue\n"); fflush(stdin); getchar(); fflush(stdin);}


/*
 * a simple exception class
 * you can create an exeption by entering
 *      throw Exception("...Error description...");
 * and get the error message from the data msg for displaying:
 *      err.msg
 */
class Exception {
   public:
      const char* msg;
      Exception(const char* arg) : msg(arg) { };
};


class Matrix
{
public:
    // constructors
    Matrix()
    {
        //printf("Executing constructor Matrix() ...\n");
        // create a Matrix object without content
        p = NULL;
        rows = 0;
        cols = 0;
    }
    Matrix(const int row_count, const int column_count)
    {
        // create a Matrix object with given number of rows and columns
        //printf("Executing constructor Matrix(const int column_count, const int row_count) ...\n");
        p = NULL;

        if (row_count > 0 && column_count > 0)
        {
            rows = row_count;
            cols = column_count;

            p = new double*[rows];
            for (int r = 0; r < rows; r++)
            {
                p[r] = new double[cols];

                // initially fill in zeros for all values in the matrix;
                for (int c = 0; c < cols; c++)
                {
                    p[r][c] = 0;
                }
            }
        }
    }
    // assignment operator
    Matrix(const Matrix& a)
    {
        //printf("Executing constructor Matrix(const Matrix& a) ...\n");
        rows = a.rows;
        cols = a.cols;
        p = new double*[a.rows];
        for (int r = 0; r < a.rows; r++)
        {
            p[r] = new double[a.cols];

            // copy the values from the matrix a
            for (int c = 0; c < a.cols; c++)
            {
                p[r][c] = a.p[r][c];
            }
        }
    }

    // index operator. You can use this class like myMatrix(col, row)
    // the indexes are one-based, not zero based.
	double& operator()(const int r, const int c)
	{
        //printf("Executing Matrix operator() ...\n");
      	if (p != NULL && r > 0 && r <= rows && c > 0 && c <= cols)
      	{
          	return p[r-1][c-1];
        }
      	else
      	{
            throw Exception("Subscript out of range");
        }
    }

    // assignment operator
    Matrix& operator= (const Matrix& a)
    {
        //printf("Executing Matrix operator= ...\n");
        rows = a.rows;
        cols = a.cols;
        p = new double*[a.rows];
        for (int r = 0; r < a.rows; r++)
        {
            p[r] = new double[a.cols];

            // copy the values from the matrix a
            for (int c = 0; c < a.cols; c++)
            {
                p[r][c] = a.p[r][c];
            }
        }
        return *this;
    }

    // operator addition
    friend Matrix operator+(const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.rows == b.rows && a.cols == b.cols)
        {
         	Matrix res(a.rows, a.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c = 0; c < a.cols; c++)
                {
                    res.p[r][c] = a.p[r][c] + b.p[r][c];
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    // add a double value (elements wise)
    Matrix& Add(double v)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                p[r][c] += v;
            }
        }
     	return *this;
    }

    // subtract a double value (elements wise)
    Matrix& Subtract(double v)
    {
        return Add(-v);
    }

    // multiply a double value (elements wise)
    Matrix& Multiply(double v)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                p[r][c] *= v;
            }
        }
     	return *this;
    }

    // operator subtraction
    friend Matrix operator- (const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.rows == b.rows && a.cols == b.cols)
        {
         	Matrix res(a.rows, a.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c = 0; c < a.cols; c++)
                {
                    res.p[r][c] = a.p[r][c] - b.p[r][c];
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }
    // operator unary minus
    friend Matrix operator- (const Matrix& a)
    {
     	Matrix res(a.rows, a.cols);

        for (int r = 0; r < a.rows; r++)
        {
            for (int c = 0; c < a.cols; c++)
            {
                res.p[r][c] = -a.p[r][c];
            }
        }

        return res;
    }

    // operator multiplication
    friend Matrix operator* (const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.cols == b.rows)
        {
         	Matrix res(a.rows, b.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c_res = 0; c_res < b.cols; c_res++)
                {
                    for (int c = 0; c < a.cols; c++)
                    {
                        res.p[r][c_res] += a.p[r][c] * b.p[c][c_res];
                    }
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    /**
     * returns the minor from the given matrix where
     * the selected row and column are removed
     */
    Matrix Minor(const int row, const int col)
    {
        //printf("minor(row=%i, col=%i), rows=%i, cols=%i\n", row, col, rows, cols);
        Matrix res;
      	if (row > 0 && row <= rows && col > 0 && col <= cols)
      	{
            res = Matrix(rows - 1, cols - 1);

            // copy the content of the matrix to the minor, except the selected
            for (int r = 1; r <= (rows - (row >= rows)); r++)
            {
                for (int c = 1; c <= (cols - (col >= cols)); c++)
                {
                    //printf("r=%i, c=%i, value=%f, rr=%i, cc=%i \n", r, c, p[r-1][c-1], r - (r > row), c - (c > col));
                    res(r - (r > row), c - (c > col)) = p[r-1][c-1];
                }
            }
        }
        else
        {
            throw Exception("Index for minor out of range");
        }

        return res;
    }

    /*
     * returns the size of the i-th dimension of the matrix.
     * i.e. for i=1 the function returns the number of rows,
     * and for i=2 the function returns the number of columns
     * else the function returns 0
     */
    int Size(const int i)
    {
        if (i == 1)
        {
            return rows;
        }
        else if (i == 2)
        {
            return cols;
        }
        return 0;
    }

    // returns the number of rows
    int GetRows()
    {
        return rows;
    }

    // returns the number of columns
    int GetCols()
    {
        return cols;
    }

    // print the contents of the matrix
    void Print()
    {
        if (p != NULL)
        {
            printf("[");
            for (int r = 0; r < rows; r++)
            {
                if (r > 0)
                {
                    printf(" ");
                }
                for (int c = 0; c < cols-1; c++)
                {
                    printf("%.2f, ", p[r][c]);
                }
                if (r < rows-1)
                {
                    printf("%.2f;\n", p[r][cols-1]);
                }
                else
                {
                    printf("%.2f]\n", p[r][cols-1]);
                }
            }
        }
        else
        {
            // matrix is empty
            printf("[ ]\n");
        }
    }

public:
    // destructor
    ~Matrix()
    {
        // clean up allocated memory
        for (int r = 0; r < rows; r++)
        {
            delete p[r];
        }
        delete p;
        p = NULL;
    }

private:
    int rows;
    int cols;
    double** p;         // pointer to a matrix with doubles
};

/*
 * i.e. for i=1 the function returns the number of rows,
 * and for i=2 the function returns the number of columns
 * else the function returns 0
 */
int Size(Matrix& a, const int i)
{
    return a.Size(i);
}


// addition of Matrix with double
Matrix operator+ (const Matrix& a, const double B)/>
{
    Matrix res = a;
    res.Add(B)/>;
    return res;
}
// addition of double with Matrix
Matrix operator+ (const double b, const Matrix& a)
{
    Matrix res = a;
    res.Add(B)/>;
    return res;
}

// subtraction of Matrix with double
Matrix operator- (const Matrix& a, const double B)/>
{
    Matrix res = a;
    res.Subtract(B)/>;
    return res;
}
// subtraction of double with Matrix
Matrix operator- (const double b, const Matrix& a)
{
    Matrix res = -a;
    res.Add(B)/>;
    return res;
}

// multiplication of Matrix with double
Matrix operator* (const Matrix& a, const double B)/>
{
    Matrix res = a;
    res.Multiply(B)/>;
    return res;
}
// multiplication of double with Matrix
Matrix operator* (const double b, const Matrix& a)
{
    Matrix res = a;
    res.Multiply(B)/>;
    return res;
}


/**
 * returns a matrix with size cols x rows with ones as values
 */
Matrix Ones(const int rows, const int cols)
{
    Matrix res = Matrix(rows, cols);

    for (int r = 1; r <= rows; r++)
    {
        for (int c = 1; c <= cols; c++)
        {
            res(r, c) = 1;
        }
    }
    return res;
}

/**
 * returns a matrix with size cols x rows with zeros as values
 */
Matrix Zeros(const int rows, const int cols)
{
    return Matrix(rows, cols);
}


/**
 * returns a diagonal matrix with size n x n with ones at the diagonal
 * @param   v a vector
 * @return  a diagonal matrix with ones on the diagonal
 */
Matrix Diag(const int n)
{
    Matrix res = Matrix(n, n);
    for (int i = 1; i <= n; i++)
    {
        res(i, i) = 1;
    }
    return res;
}
/**
 * returns a diagonal matrix
 * @param   v a vector
 * @return  a diagonal matrix with the given vector v on the diagonal
 */
Matrix Diag(Matrix& v)
{
    Matrix res;
    if (v.GetCols() == 1)
    {
        // the given matrix is a vector n x 1
        int rows = v.GetRows();
        res = Matrix(rows, rows);

        // copy the values of the vector to the matrix
        for (int r=1; r <= rows; r++)
        {
            res(r, r) = v(r, 1);
        }
    }
    else if (v.GetRows() == 1)
    {
        // the given matrix is a vector 1 x n
        int cols = v.GetCols();
        res = Matrix(cols, cols);

        // copy the values of the vector to the matrix
        for (int c=1; c <= cols; c++)
        {
            res(c, c) = v(1, c);
        }
    }
    else
    {
        throw Exception("Parameter for diag must be a vector");
    }
    return res;
}

/*
 * returns the determinant of Matrix a
 */
double Det(Matrix& a)
{
    double d = 0;       // value of the determinant
    int rows = a.GetRows();
    int cols = a.GetRows();

    if (rows == cols)
    {
        // this is a square matrix
        if (rows == 1)
        {
            // this is a 1 x 1 matrix
            d = a(1, 1);
        }
        else if (rows == 2)
        {
            // this is a 2 x 2 matrix
            // the determinant of [a11,a12;a21,a22] is det = a11*a22-a21*a12
            d = a(1, 1) * a(2, 2) - a(2, 1) * a(1, 2);
        }
        else
        {
            // this is a matrix of 3 x 3 or larger
            for (int c = 1; c <= cols; c++)
            {
                Matrix M = a.Minor(1, c);
                //d += pow(-1, 1+c)  * a(1, c) * Det(M);
                d += (c%2 + c%2 - 1) * a(1, c) * Det(M);  // faster than with pow()
            }
        }
    }
    else
    {
        throw Exception("Matrix must be square");
    }
    return d;
}

// swap two values
void Swap(double& a, double& B)/>
{
    double temp = a;
    a = b;
    b = temp;
}

/*
 * returns the inverse of Matrix a
 */
Matrix Inv(Matrix& a)
{
    Matrix res;
    double d = 0;       // value of the determinant
    int rows = a.GetRows();
    int cols = a.GetRows();

    d = Det(a);
    if (rows == cols && d != 0)
    {
        // this is a square matrix
        if (rows == 1)
        {
            // this is a 1 x 1 matrix
            res = Matrix(rows, cols);
            res(1, 1) = 1 / a(1, 1);
        }
        else if (rows == 2)
        {
            // this is a 2 x 2 matrix
            res = Matrix(rows, cols);
            res(1, 1) = a(2, 2);
            res(1, 2) = -a(1, 2);
            res(2, 1) = -a(2, 1);
            res(2, 2) = a(1, 1);
            res = (1/d) * res;
        }
        else
        {
            // this is a matrix of 3 x 3 or larger
            // calculate inverse using gauss-jordan elimination
            //      http://mathworld.wolfram.com/MatrixInverse.html
            //      http://math.uww.edu/~mcfarlat/inverse.htm
            res = Diag(rows);     // a diagonal matrix with ones at the diagonal
            Matrix ai = a;        // make a copy of Matrix a

            for (int c = 1; c <= cols; c++)
            {
                // element (c, c) should be non zero. if not, swap content
                // of lower rows
                int r;
                for (r = c; r <= rows && ai(r, c) == 0; r++)
                {
                }
                if (r != c)
                {
                    // swap rows
                    for (int s = 1; s <= cols; s++)
                    {
                        Swap(ai(c, s), ai(r, s));
                        Swap(res(c, s), res(r, s));
                    }
                }

                // eliminate non-zero values on the other rows at column c
                for (int r = 1; r <= rows; r++)
                {
                    if(r != c)
                    {
                        // eleminate value at column c and row r
                        if (ai(r, c) != 0)
                        {
                            double f = - ai(r, c) / ai(c, c);

                            // add (f * row c) to row r to eleminate the value
                            // at column c
                            for (int s = 1; s <= cols; s++)
                            {
                                ai(r, s) += f * ai(c, s);
                                res(r, s) += f * res(c, s);
                            }
                        }
                    }
                    else
                    {
                        // make value at (c, c) one,
                        // divide each value on row r with the value at ai(c,c)
                        double f = ai(c, c);
                        for (int s = 1; s <= cols; s++)
                        {
                            ai(r, s) /= f;
                            res(r, s) /= f;
                        }
                    }
                }
            }
        }
    }
    else
    {
        if (rows == cols)
        {
            throw Exception("Matrix must be square");
        }
        else
        {
            throw Exception("Determinant of matrix is zero");
        }
    }
    return res;
}


int main(int argc, char *argv[])
{
    try
    {
        // create an empty matrix of 3x3 (will initially contain zeros)
        int cols = 3;
        int rows = 3;
        Matrix A = Matrix(cols, rows);

        // fill in some values in matrix a
        int count = 0;
        for (int r = 1; r <= rows; r++)
        {
            for (int c = 1; c <= cols; c++)
            {
                count ++;
                A(r, c) = count;
            }
        }

        // adjust a value in the matrix  (indexes are one-based)
        A(2,1) = 1.23;

        // read a value from the matrix  (indexes are one-based)
        double centervalue = A(2,2);
        printf("centervalue = %f \n", centervalue);
        printf("\n");

        // print the whole matrix
        printf("A = \n");
        A.Print();
        printf("\n");

        Matrix B = Ones(rows, cols) + Diag(rows);
        printf("B = \n");
        B.Print();
        printf("\n");

        Matrix B2 = Matrix(rows, 1);        // a vector
        count = 1;
        for (int r = 1; r <= rows; r++)
        {
            count ++;
            B2(r, 1) = count * 2;
        }
        printf("B2 = \n");
        B2.Print();
        printf("\n");

        Matrix C;
        C = A + B;
        printf("A + B = \n");
        C.Print();
        printf("\n");

        C = A - B;
        printf("A - B = \n");
        C.Print();
        printf("\n");

        C = A * B2;
        printf("A * B2 = \n");
        C.Print();
        printf("\n");

        // create a diagonal matrix
        Matrix E = Diag(B2);
        printf("E = \n");
        E.Print();
        printf("\n");

        // calculate determinant
        Matrix D = Matrix(2, 2);
        D(1,1) = 2;
        D(1,2) = 4;
        D(2,1) = 1;
        D(2,2) = -2;
        printf("D = \n");
        D.Print();
        printf("Det(D) = %f\n\n", Det(D));

        printf("A = \n");
        A.Print();
        printf("\n");
        printf("Det(A) = %f\n\n", Det(A));

        Matrix F;
        F = 3 - A ;
        printf("3 - A = \n");
        F.Print();
        printf("\n");

        // test inverse
        Matrix G = Matrix(2, 2);
        G(1, 1) = 1;
        G(1, 2) = 2;
        G(2, 1) = 3;
        G(2, 2) = 4;
        printf("G = \n");
        G.Print();
        printf("\n");
        Matrix G_inv = Inv(G);
        printf("Inv(G) = \n");
        G_inv.Print();
        printf("\n");

        Matrix A_inv = Inv(A);
        printf("Inv(A) = \n");
        A_inv.Print();
        printf("\n");

        rows = 2;
        cols = 5;
        Matrix H = Matrix(rows, cols);
        for (int r = 1; r <= rows; r++)
        {
            for (int c = 1; c <= cols; c++)
            {
                count ++;
                H(r, c) = count;
            }
        }
        printf("H = \n");
        H.Print();
        printf("\n");
    }
    catch (Exception err)
    {
        printf("Error: %s\n", err.msg);
    }
    catch (...)
    {
        printf("An error occured...\n");
    }

    printf("\n");
    PAUSE;

    return EXIT_SUCCESS;
}


Was This Post Helpful? 0
  • +
  • -

#4 joske  Icon User is offline

  • D.I.C Regular
  • member icon

Reputation: 43
  • View blog
  • Posts: 297
  • Joined: 04-September 07

Re: determinant of matrix in c++

Posted 21 March 2010 - 09:28 AM

I cleaned up the code a little bit, and added matrix division.

/*
A simple matrix class
c++ code
Author: Jos de Jong, Nov 2007. Updated March 2010

With this class you can:
    - create a matrix
    - get/set the cell values
    - use operators +, -, *, /
    - use functions Ones(), Zeros(), Diag(), Det(), Inv(), Size()
    - print the content of the matrix

Usage:
    you can create a matrix by:
        Matrix A;
        Matrix A = Matrix(rows, cols);
        Matrix A = B;

    you can get and set matrix elements by:
        A(2,3) = 5.6;        // set an element of Matix A
        value = A(3,1);      // get an element of Matrix A
        value = A.get(3,1);  // get an element of a constant Matrix A
        A = B;               // copy content of Matrix B to Matrix A

    you can apply operations with matrices and doubles:
        A = B + C;
        A = B - C;
        A = -B;
        A = B * C;
        A = B / C;

    the following functions are available:
        A = Ones(rows, cols);
        A = Zeros(rows, cols);
        A = Diag(n);
        A = Diag(B)/>;
        d = Det(A);
        A = Inv(B)/>;
        cols = A.GetCols();
        rows = A.GetRows();
        cols = Size(A, 1);
        rows = Size(A, 2);

    you can quick-print the content of a matrix in the console with:
        A.Print();
*/

#include <cstdlib>
#include <cstdio>
#include <math.h>

#define PAUSE {printf("Press \"Enter\" to continue\n"); fflush(stdin); getchar(); fflush(stdin);}

// Declarations
class Matrix;
double Det(const Matrix& a);
Matrix Diag(const int n);
Matrix Diag(const Matrix& v);
Matrix Inv(const Matrix& a);
Matrix Ones(const int rows, const int cols);
int Size(const Matrix& a, const int i);
Matrix Zeros(const int rows, const int cols);


/*
 * a simple exception class
 * you can create an exeption by entering
 *      throw Exception("...Error description...");
 * and get the error message from the data msg for displaying:
 *      err.msg
 */
class Exception
{
   public:
      const char* msg;
      Exception(const char* arg)
       : msg(arg)
      {
      }
};



class Matrix
{
public:
    // constructor
    Matrix()
    {
        //printf("Executing constructor Matrix() ...\n");
        // create a Matrix object without content
        p = NULL;
        rows = 0;
        cols = 0;
    }

    // constructor
    Matrix(const int row_count, const int column_count)
    {
        // create a Matrix object with given number of rows and columns
        p = NULL;

        if (row_count > 0 && column_count > 0)
        {
            rows = row_count;
            cols = column_count;

            p = new double*[rows];
            for (int r = 0; r < rows; r++)
            {
                p[r] = new double[cols];

                // initially fill in zeros for all values in the matrix;
                for (int c = 0; c < cols; c++)
                {
                    p[r][c] = 0;
                }
            }
        }
    }

    // assignment operator
    Matrix(const Matrix& a)
    {
        rows = a.rows;
        cols = a.cols;
        p = new double*[a.rows];
        for (int r = 0; r < a.rows; r++)
        {
            p[r] = new double[a.cols];

            // copy the values from the matrix a
            for (int c = 0; c < a.cols; c++)
            {
                p[r][c] = a.p[r][c];
            }
        }
    }

    // index operator. You can use this class like myMatrix(col, row)
    // the indexes are one-based, not zero based.
	  double& operator()(const int r, const int c)
	  {
      	if (p != NULL && r > 0 && r <= rows && c > 0 && c <= cols)
      	{
          	return p[r-1][c-1];
        }
      	else
      	{
            throw Exception("Subscript out of range");
        }
    }

    // index operator. You can use this class like myMatrix.get(col, row)
    // the indexes are one-based, not zero based.
    // use this function get if you want to read from a const Matrix
	  double get(const int r, const int c) const
	  {
      	if (p != NULL && r > 0 && r <= rows && c > 0 && c <= cols)
      	{
          	return p[r-1][c-1];
        }
      	else
      	{
            throw Exception("Subscript out of range");
        }
    }

    // assignment operator
    Matrix& operator= (const Matrix& a)
    {
        rows = a.rows;
        cols = a.cols;
        p = new double*[a.rows];
        for (int r = 0; r < a.rows; r++)
        {
            p[r] = new double[a.cols];

            // copy the values from the matrix a
            for (int c = 0; c < a.cols; c++)
            {
                p[r][c] = a.p[r][c];
            }
        }
        return *this;
    }

    // add a double value (elements wise)
    Matrix& Add(const double v)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                p[r][c] += v;
            }
        }
     	  return *this;
    }

    // subtract a double value (elements wise)
    Matrix& Subtract(const double v)
    {
        return Add(-v);
    }

    // multiply a double value (elements wise)
    Matrix& Multiply(const double v)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                p[r][c] *= v;
            }
        }
     	  return *this;
    }

    // divide a double value (elements wise)
    Matrix& Divide(const double v)
    {
     	  return Multiply(1/v);
    }

    // addition of Matrix with Matrix
    friend Matrix operator+(const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.rows == b.rows && a.cols == b.cols)
        {
         	Matrix res(a.rows, a.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c = 0; c < a.cols; c++)
                {
                    res.p[r][c] = a.p[r][c] + b.p[r][c];
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    // addition of Matrix with double
    friend Matrix operator+ (const Matrix& a, const double B)/>
    {
        Matrix res = a;
        res.Add(B)/>;
        return res;
    }
    // addition of double with Matrix
    friend Matrix operator+ (const double b, const Matrix& a)
    {
        Matrix res = a;
        res.Add(B)/>;
        return res;
    }

    // subtraction of Matrix with Matrix
    friend Matrix operator- (const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.rows == b.rows && a.cols == b.cols)
        {
         	Matrix res(a.rows, a.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c = 0; c < a.cols; c++)
                {
                    res.p[r][c] = a.p[r][c] - b.p[r][c];
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    // subtraction of Matrix with double
    friend Matrix operator- (const Matrix& a, const double B)/>
    {
        Matrix res = a;
        res.Subtract(B)/>;
        return res;
    }
    // subtraction of double with Matrix
    friend Matrix operator- (const double b, const Matrix& a)
    {
        Matrix res = -a;
        res.Add(B)/>;
        return res;
    }

    // operator unary minus
    friend Matrix operator- (const Matrix& a)
    {
     	Matrix res(a.rows, a.cols);

        for (int r = 0; r < a.rows; r++)
        {
            for (int c = 0; c < a.cols; c++)
            {
                res.p[r][c] = -a.p[r][c];
            }
        }

        return res;
    }

    // operator multiplication
    friend Matrix operator* (const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match
        if (a.cols == b.rows)
        {
         	Matrix res(a.rows, b.cols);

            for (int r = 0; r < a.rows; r++)
            {
                for (int c_res = 0; c_res < b.cols; c_res++)
                {
                    for (int c = 0; c < a.cols; c++)
                    {
                        res.p[r][c_res] += a.p[r][c] * b.p[c][c_res];
                    }
                }
            }
         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    // multiplication of Matrix with double
    friend Matrix operator* (const Matrix& a, const double B)/>
    {
        Matrix res = a;
        res.Multiply(B)/>;
        return res;
    }
    // multiplication of double with Matrix
    friend Matrix operator* (const double b, const Matrix& a)
    {
        Matrix res = a;
        res.Multiply(B)/>;
        return res;
    }

    // division of Matrix with Matrix
    friend Matrix operator/ (const Matrix& a, const Matrix& B)/>
    {
        // check if the dimensions match: must be square and equal sizes
        if (a.rows == a.cols && a.rows == a.cols && b.rows == b.cols)
        {
         	Matrix res(a.rows, a.cols);

          res = a * Inv(B)/>;

         	return res;
        }
        else
        {
            // give an error
            throw Exception("Dimensions does not match");
        }

        // return an empty matrix (this never happens but just for safety)
        return Matrix();
    }

    // division of Matrix with double
    friend Matrix operator/ (const Matrix& a, const double B)/>
    {
        Matrix res = a;
        res.Divide(B)/>;
        return res;
    }

    // division of double with Matrix
    friend Matrix operator/ (const double b, const Matrix& a)
    {
        Matrix b_matrix(1, 1);
        b_matrix(1,1) = b;

        Matrix res = b_matrix / a;
        return res;
    }

    /**
     * returns the minor from the given matrix where
     * the selected row and column are removed
     */
    Matrix Minor(const int row, const int col) const
    {
        Matrix res;
      	if (row > 0 && row <= rows && col > 0 && col <= cols)
      	{
            res = Matrix(rows - 1, cols - 1);

            // copy the content of the matrix to the minor, except the selected
            for (int r = 1; r <= (rows - (row >= rows)); r++)
            {
                for (int c = 1; c <= (cols - (col >= cols)); c++)
                {
                    res(r - (r > row), c - (c > col)) = p[r-1][c-1];
                }
            }
        }
        else
        {
            throw Exception("Index for minor out of range");
        }

        return res;
    }

    /*
     * returns the size of the i-th dimension of the matrix.
     * i.e. for i=1 the function returns the number of rows,
     * and for i=2 the function returns the number of columns
     * else the function returns 0
     */
    int Size(const int i) const
    {
        if (i == 1)
        {
            return rows;
        }
        else if (i == 2)
        {
            return cols;
        }
        return 0;
    }

    // returns the number of rows
    int GetRows() const
    {
        return rows;
    }

    // returns the number of columns
    int GetCols() const
    {
        return cols;
    }

    // print the contents of the matrix
    void Print() const
    {
        if (p != NULL)
        {
            printf("[");
            for (int r = 0; r < rows; r++)
            {
                if (r > 0)
                {
                    printf(" ");
                }
                for (int c = 0; c < cols-1; c++)
                {
                    printf("%.2f, ", p[r][c]);
                }
                if (r < rows-1)
                {
                    printf("%.2f;\n", p[r][cols-1]);
                }
                else
                {
                    printf("%.2f]\n", p[r][cols-1]);
                }
            }
        }
        else
        {
            // matrix is empty
            printf("[ ]\n");
        }
    }

public:
    // destructor
    ~Matrix()
    {
        // clean up allocated memory
        for (int r = 0; r < rows; r++)
        {
            delete p[r];
        }
        delete p;
        p = NULL;
    }

private:
    int rows;
    int cols;
    double** p;         // pointer to a matrix with doubles
};

/*
 * i.e. for i=1 the function returns the number of rows,
 * and for i=2 the function returns the number of columns
 * else the function returns 0
 */
int Size(const Matrix& a, const int i)
{
    return a.Size(i);
}


/**
 * returns a matrix with size cols x rows with ones as values
 */
Matrix Ones(const int rows, const int cols)
{
    Matrix res = Matrix(rows, cols);

    for (int r = 1; r <= rows; r++)
    {
        for (int c = 1; c <= cols; c++)
        {
            res(r, c) = 1;
        }
    }
    return res;
}

/**
 * returns a matrix with size cols x rows with zeros as values
 */
Matrix Zeros(const int rows, const int cols)
{
    return Matrix(rows, cols);
}


/**
 * returns a diagonal matrix with size n x n with ones at the diagonal
 * @param   v a vector
 * @return  a diagonal matrix with ones on the diagonal
 */
Matrix Diag(const int n)
{
    Matrix res = Matrix(n, n);
    for (int i = 1; i <= n; i++)
    {
        res(i, i) = 1;
    }
    return res;
}
/**
 * returns a diagonal matrix
 * @param   v a vector
 * @return  a diagonal matrix with the given vector v on the diagonal
 */
Matrix Diag(const Matrix& v)
{
    Matrix res;
    if (v.GetCols() == 1)
    {
        // the given matrix is a vector n x 1
        int rows = v.GetRows();
        res = Matrix(rows, rows);

        // copy the values of the vector to the matrix
        for (int r=1; r <= rows; r++)
        {
            res(r, r) = v.get(r, 1);
        }
    }
    else if (v.GetRows() == 1)
    {
        // the given matrix is a vector 1 x n
        int cols = v.GetCols();
        res = Matrix(cols, cols);

        // copy the values of the vector to the matrix
        for (int c=1; c <= cols; c++)
        {
            res(c, c) = v.get(1, c);
        }
    }
    else
    {
        throw Exception("Parameter for diag must be a vector");
    }
    return res;
}

/*
 * returns the determinant of Matrix a
 */
double Det(const Matrix& a)
{
    double d = 0;       // value of the determinant
    int rows = a.GetRows();
    int cols = a.GetRows();

    if (rows == cols)
    {
        // this is a square matrix
        if (rows == 1)
        {
            // this is a 1 x 1 matrix
            d = a.get(1, 1);
        }
        else if (rows == 2)
        {
            // this is a 2 x 2 matrix
            // the determinant of [a11,a12;a21,a22] is det = a11*a22-a21*a12
            d = a.get(1, 1) * a.get(2, 2) - a.get(2, 1) * a.get(1, 2);
        }
        else
        {
            // this is a matrix of 3 x 3 or larger
            for (int c = 1; c <= cols; c++)
            {
                Matrix M = a.Minor(1, c);
                //d += pow(-1, 1+c)  * a(1, c) * Det(M);
                d += (c%2 + c%2 - 1) * a.get(1, c) * Det(M);  // faster than with pow()
            }
        }
    }
    else
    {
        throw Exception("Matrix must be square");
    }
    return d;
}

// swap two values
void Swap(double& a, double& B)/>
{
    double temp = a;
    a = b;
    b = temp;
}

/*
 * returns the inverse of Matrix a
 */
Matrix Inv(const Matrix& a)
{
    Matrix res;
    double d = 0;       // value of the determinant
    int rows = a.GetRows();
    int cols = a.GetRows();

    d = Det(a);
    if (rows == cols && d != 0)
    {
        // this is a square matrix
        if (rows == 1)
        {
            // this is a 1 x 1 matrix
            res = Matrix(rows, cols);
            res(1, 1) = 1 / a.get(1, 1);
        }
        else if (rows == 2)
        {
            // this is a 2 x 2 matrix
            res = Matrix(rows, cols);
            res(1, 1) = a.get(2, 2);
            res(1, 2) = -a.get(1, 2);
            res(2, 1) = -a.get(2, 1);
            res(2, 2) = a.get(1, 1);
            res = (1/d) * res;
        }
        else
        {
            // this is a matrix of 3 x 3 or larger
            // calculate inverse using gauss-jordan elimination
            //      http://mathworld.wolfram.com/MatrixInverse.html
            //      http://math.uww.edu/~mcfarlat/inverse.htm
            res = Diag(rows);     // a diagonal matrix with ones at the diagonal
            Matrix ai = a;        // make a copy of Matrix a

            for (int c = 1; c <= cols; c++)
            {
                // element (c, c) should be non zero. if not, swap content
                // of lower rows
                int r;
                for (r = c; r <= rows && ai(r, c) == 0; r++)
                {
                }
                if (r != c)
                {
                    // swap rows
                    for (int s = 1; s <= cols; s++)
                    {
                        Swap(ai(c, s), ai(r, s));
                        Swap(res(c, s), res(r, s));
                    }
                }

                // eliminate non-zero values on the other rows at column c
                for (int r = 1; r <= rows; r++)
                {
                    if(r != c)
                    {
                        // eleminate value at column c and row r
                        if (ai(r, c) != 0)
                        {
                            double f = - ai(r, c) / ai(c, c);

                            // add (f * row c) to row r to eleminate the value
                            // at column c
                            for (int s = 1; s <= cols; s++)
                            {
                                ai(r, s) += f * ai(c, s);
                                res(r, s) += f * res(c, s);
                            }
                        }
                    }
                    else
                    {
                        // make value at (c, c) one,
                        // divide each value on row r with the value at ai(c,c)
                        double f = ai(c, c);
                        for (int s = 1; s <= cols; s++)
                        {
                            ai(r, s) /= f;
                            res(r, s) /= f;
                        }
                    }
                }
            }
        }
    }
    else
    {
        if (rows == cols)
        {
            throw Exception("Matrix must be square");
        }
        else
        {
            throw Exception("Determinant of matrix is zero");
        }
    }
    return res;
}


int main(int argc, char *argv[])
{
    try
    {
        // create an empty matrix of 3x3 (will initially contain zeros)
        int cols = 3;
        int rows = 3;
        Matrix A = Matrix(cols, rows);

        // fill in some values in matrix a
        int count = 0;
        for (int r = 1; r <= rows; r++)
        {
            for (int c = 1; c <= cols; c++)
            {
                count ++;
                A(r, c) = count;
            }
        }

        // adjust a value in the matrix  (indexes are one-based)
        A(2,1) = 1.23;

        // read a value from the matrix  (indexes are one-based)
        double centervalue = A(2,2);
        printf("centervalue = %f \n", centervalue);
        printf("\n");

        // print the whole matrix
        printf("A = \n");
        A.Print();
        printf("\n");

        Matrix B = Ones(rows, cols) + Diag(rows);
        printf("B = \n");
        B.Print();
        printf("\n");

        Matrix B2 = Matrix(rows, 1);        // a vector
        count = 1;
        for (int r = 1; r <= rows; r++)
        {
            count ++;
            B2(r, 1) = count * 2;
        }
        printf("B2 = \n");
        B2.Print();
        printf("\n");

        Matrix C;
        C = A + B;
        printf("A + B = \n");
        C.Print();
        printf("\n");

        C = A - B;
        printf("A - B = \n");
        C.Print();
        printf("\n");

        C = A * B2;
        printf("A * B2 = \n");
        C.Print();
        printf("\n");

        // create a diagonal matrix
        Matrix E = Diag(B2);
        printf("E = \n");
        E.Print();
        printf("\n");

        // calculate determinant
        Matrix D = Matrix(2, 2);
        D(1,1) = 2;
        D(1,2) = 4;
        D(2,1) = 1;
        D(2,2) = -2;
        printf("D = \n");
        D.Print();
        printf("Det(D) = %f\n\n", Det(D));

        printf("A = \n");
        A.Print();
        printf("\n");
        printf("Det(A) = %f\n\n", Det(A));

        Matrix F;
        F = 3 - A ;
        printf("3 - A = \n");
        F.Print();
        printf("\n");

        // test inverse
        Matrix G = Matrix(2, 2);
        G(1, 1) = 1;
        G(1, 2) = 2;
        G(2, 1) = 3;
        G(2, 2) = 4;
        printf("G = \n");
        G.Print();
        printf("\n");
        Matrix G_inv = Inv(G);
        printf("Inv(G) = \n");
        G_inv.Print();
        printf("\n");

        Matrix A_inv = Inv(A);
        printf("Inv(A) = \n");
        A_inv.Print();
        printf("\n");

        Matrix A_A_inv = A * Inv(A);
        printf("A * Inv(A) = \n");
        A_A_inv.Print();
        printf("\n");

        Matrix B_A = B / A;
        printf("B / A = \n");
        B_A.Print();
        printf("\n");

        Matrix A_3 = A / 3;
        printf("A / 3 = \n");
        A_3.Print();
        printf("\n");

        rows = 2;
        cols = 5;
        Matrix H = Matrix(rows, cols);
        for (int r = 1; r <= rows; r++)
        {
            for (int c = 1; c <= cols; c++)
            {
                count ++;
                H(r, c) = count;
            }
        }
        printf("H = \n");
        H.Print();
        printf("\n");
    }
    catch (Exception err)
    {
        printf("Error: %s\n", err.msg);
    }
    catch (...)
    {
        printf("An error occured...\n");
    }

    printf("\n");
    PAUSE;

    return EXIT_SUCCESS;
}


Was This Post Helpful? 0
  • +
  • -

#5 Guest_Chris*


Reputation:

Re: determinant of matrix in c++

Posted 19 November 2010 - 07:23 AM

Yikes!
After finding this thread, I got scared so wrote my own:

/* const int dim  is the dimensionality of the matrix */
	double det = 0.0;
	for (int i = 0; i < dim; i++)
	{
	    double a = 1.0, b = 1.0;
	    for (int row = 0; row < dim; row++)
		{
		    a *= m[row][(i+row)%dim];
			b *= m[row][(dim-1) - (i+row)%dim];
		}
		det += a - b;
	}
	return det;


MOD EDIT: When posting code...USE CODE TAGS!!!

:code:

This post has been edited by JackOfAllTrades: 19 November 2010 - 08:22 AM

Was This Post Helpful? 0

#6 Yiannis_Spyridakis  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 2
  • Joined: 19-June 11

Re: determinant of matrix in c++

Posted 19 June 2011 - 05:11 AM

View Postjoske, on 21 March 2010 - 09:28 AM, said:

//...
public:
    // destructor
    ~Matrix()
    {
        // clean up allocated memory
        for (int r = 0; r < rows; r++)
        {
            delete p[r];
        }
        delete p;
        p = NULL;
    }




Thanks Joske, this is very useful.

The destructor needs one correction, I think: You need to call delete[] instead of delete when deallocating both the whole array and individual rows (otherwise you leak memory)

So the destructor becomes:
//...
public:
    // destructor
    ~Matrix()
    {
        // clean up allocated memory
        for (int r = 0; r < rows; r++)
        {
            delete[] p[r];
        }
        delete[] p;
        p = NULL;
    }




Cheers,
Yiannis
Was This Post Helpful? 0
  • +
  • -

#7 Yiannis_Spyridakis  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 2
  • Joined: 19-June 11

Re: determinant of matrix in c++

Posted 19 June 2011 - 05:57 AM

View PostYiannis_Spyridakis, on 19 June 2011 - 05:11 AM, said:

...
The destructor needs one correction, I think: You need to call delete[] instead of delete when deallocating both the whole array and individual rows (otherwise you leak memory)

So the destructor becomes:
...
Cheers,
Yiannis


Similarily, the assignment operator needs to check for previously allocated memory and deallocate it if necessary. That is, if you 1) allocate a Matrix and then 2) assign it to another matrix the code as it now is will cause a memory leak.

It's often convenient to introduce a 'destroy' private function that does all the deallocation and implement the copy constructor by delegating to the assignment operator.
Here is how I would implement it (I had to rename the class to MyMatrix to avoid name clashes. Also I am currently refreshing my C++, so I wouldn't blindly trust this code :whistling: )

// Copy constructor
MyMatrix::MyMatrix(const MyMatrix& RHS)
	: p(nullptr), cols(0), rows(0)
{
	this->operator=(RHS);
}

// Assignment operator
MyMatrix& MyMatrix::operator=(const MyMatrix& RHS){
	
	if (p){
		this->destroy();
	}

	rows = RHS.rows;
	cols = RHS.cols;
	p = new double*[rows];
	
	int r, c;
	for (r = 0; r < rows; ++r){
		// Allocate a new row of elements
		p[r] = new double[cols];
		// Copy the values from matrix RHS
		for (c = 0; c < cols; ++c){
			p[r][c] = RHS.p[r][c];
		}
	}
	return *this;
}

// Private helper (used for deallocation)
void MyMatrix::destroy(){
	// Clean up allocated memory
	for (int r = 0; r < rows; ++r){
		delete[] p[r];
	}
	delete[] p; p = nullptr;
	cols = rows = 0;
}

// Destructor
MyMatrix::~MyMatrix()
{
	this->destroy();
}



Was This Post Helpful? 0
  • +
  • -

#8 rnarozniak  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 1
  • Joined: 06-February 12

Re: determinant of matrix in c++

Posted 06 February 2012 - 05:36 PM

this thread is kind of old, but I would like to make a couple of suggestions re your matrix class:

1. the recursive, laplacian expansion algorithm you are using scales at O(n!). The preferred approach is to use the fact that the determinant of an n^2 matrix is simply the product of the elements in the diag of the upper triangular matrix. This reduces the order to ~O(2n^2/3). Your recursive algorithm will be out of control for, say, n>12.

2. I would also suggest using LU decomposition for your inversion algorithm. This reduces the number of flops from about O(n^3) to O(n^2).

I hope this helps.

Regards,

This post has been edited by modi123_1: 06 February 2012 - 06:09 PM
Reason for edit:: removed excessively long quote

Was This Post Helpful? 0
  • +
  • -

Page 1 of 1