Introduction:
Matrix is a simple c++ program including operations and functions for matrices.
The program is a modified version of Matrix class by Jos de Jong.
It supports operators: + - * / and output operator <<
It supports functions: Ones, Zeros, Diag.
About:
matrix program modified by Hamid Soltani. (gmail: hsoltanim)
Code:
Download matrix.cpp
/*
A simple matrix class
c++ code
Author: Jos de Jong, Nov 2007. Updated March 2010
Modified by by Hamid Soltani. (gmail: hsoltanim)
<< opeator added, Inv function improved, some minor changes
https://csvparser.github.io/
Last modified: Aug. 2016.
With this class you can:
- create a 2D matrix with custom size
- get/set the cell values
- use operators +, -, *, /
- use functions Ones(), Zeros(), Diag(), Inv()
- 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);
A = Inv(B);
cols = A.GetCols();
rows = A.GetRows();
you can quick-print the content of a matrix with << operator
*/
#include "stdafx.h"
#include < cstdlib >
#include < cstdio >
#include < math.h >
#include < iostream >
#include < iomanip >
#include < limits >
using namespace std;
// Declarations
class Matrix;
static double lastDet; // Determinant computed by Inv function
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
{
private:
int rows;
int cols;
double** p; // pointer to a matrix with doubles
public:
// constructor
Matrix()
{
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
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
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
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
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 number of rows
int GetRows() const
{
return rows;
}
// returns the number of columns
int GetCols() const
{
return cols;
}
// output operator
friend ostream& operator<<(ostream& os, const Matrix& M)
{
if (M.p != NULL)
for (int r = 0; r < M.rows; r++)
for (int c = 0; c < M.cols; c++)
os << ((c==0) ? ((r == 0) ? "[" : " ") : "")
<< setw(10) << setprecision(4) << M.p[r][c]
<< ((c == M.cols - 1) ? ((r == M.rows - 1) ? "]" : ";\n") : ",");
else
os << "[ ]";
return os;
}
public:
// destructor
~Matrix()
{
// clean up allocated memory
for (int r = 0; r < rows; r++)
{
delete p[r];
}
delete p;
p = NULL;
}
};
/**
* 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;
}
// swap two values
void Swap(double& a, double& b)
{
double temp = a;
a = b;
b = temp;
}
/*
* returns the inverse of Matrix a, stores determinent in DT
*/
Matrix Inv(const Matrix& a, double& DT)
{
Matrix res;
int rows = a.GetRows();
int cols = a.GetRows();
if (rows != cols)
{
throw Exception("Matrix must be square");
return res;
}
// 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
DT = 1;
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 >rows)
{
DT = 0;
throw Exception("Determinant of matrix is zero");
return res;
}
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);
DT *= f;
for (int s = 1; s <= cols; s++)
{
ai(r, s) /= f;
res(r, s) /= f;
}
}
}
}
return res;
}
/*
* returns the inverse of Matrix a, stores determinent in lastDet
*/
Matrix Inv(const Matrix& a)
{
return Inv(a, lastDet);
}
int main(int argc, char *argv[])
{
// below some demonstration of the usage of the Matrix class
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++)
A(r, c) = ++count;
// adjust a value in the matrix (indexes are one-based)
A(2, 1) = 1.23;
// print the whole matrix
cout << "A= \n" << A << "\n";
Matrix B = Ones(rows, cols) + Diag(rows);
cout << "B= \n" << B << "\n";
Matrix B2 = Matrix(rows, 1); // a vector
count = 1;
for (int r = 1; r <= rows; r++)
B2(r, 1) = ++count * 2;
cout << "B2= \n" << B2 << "\n";
cout << "A + B= \n" << A + B << "\n";
cout << "A - B= \n" << A - B << "\n";
cout << "A * B2= \n" << A * B2 << "\n";
cout << "Diag. B2= \n" << Diag(B2) << "\n";
cout << "A= \n" << A << "\n";
cout << "3 - A= \n" << 3 - A << "\n";
Matrix A_inv = Inv(A);
cout << "Inv(A)= \n" << A_inv << "\n";
cout << "Det(A)=" << lastDet << "\n";
cout << "A * Inv(A)= \n" << A * Inv(A) << "\n";
cout << "B / A = \n" << B / A << "\n";
cout << "A / 3 = \n" << A / 3 << "\n";
rows = 2;
cols = 5;
Matrix H = Matrix(rows, cols);
for (int r = 1; r <= rows; r++)
for (int c = 1; c <= cols; c++)
H(r, c) = ++count;
cout << "H= \n" << H << "\n";
int test = 1000;
int band = 20;
cout << "\n\nCreating a " << test << "*" << test
<< " Matrix with " << 2 * band + 1
<< " non-zero middle elements at each row, please wait...";
Matrix M = Matrix(test, test);
for (int i = 1;i <= test;i++)
{
int a = (i <= band) ? 1 : i - band;
int b = (i > test - band) ? test : i + band;
for (int j = a;j <= b;j++)
M(i, j) = 1.0 + rand() % 10;
}
cout << "\n\nComputing Inv M, please wait...";
Matrix MI = Inv(M);
cout << "\n\nComputing M*Inv M, please wait...";
Matrix MMI = M*Inv(M);
double e1 = 1e300;
double d1 = e1;
double e2 = -e1;
double d2 = e2;
for (int i = 1; i <= test;i++)
for (int j = 1; j <= test;j++)
if (i == j)
{
if (d1 > MMI(i, j)) d1 = MMI(i, j);
if (d2 < MMI(i, j)) d2 = MMI(i, j);
} else {
if (e1 > MMI(i, j)) e1 = MMI(i, j);
if (e2 < MMI(i, j)) e2 = MMI(i, j);
}
typedef numeric_limits< double > dbl;
cout << setprecision (dbl::max_digits10)
<< "\n\nDiagonal elements ranging from " << d1 << " to " << d2
<< "\nNon-diagonal elements ranging from " << e1 << " to " << e2;
}
catch (Exception err)
{
cout << "Error: %s\n" << err.msg;
}
catch (...)
{
cout << "An error occured...\n";
}
cout << "\n\nPress Enter to exit...\n";
getchar();
return EXIT_SUCCESS;
}