#include "mymath.h"
/**-------PRIVATE------------------------------------------------*/
// Determines the distance between two values in the vector
usint Matrix::distance()
{
std::ostringstream out;
std::string str;
usint abst = 2;
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
out << round(_arr[i][j], _precision);
// cast value in a string
str = out.str();
usint k = str.length() + 1;
if (k > abst)
{
abst = k;
}
out.str(""); // clear stream
out.clear(); // reset flags
}
}
return abst;
}
// Set the matrix type
void Matrix::setType()
{
// set default value
this->_type = DEFAULT_MATRIX;
// check if matrix is a ROW_MATRIX
if (_row == 1 and _col > 1)
{
this->_type = ROW_MATRIX;
}
// check if matrix is a COL_MATRIX
else if (_col == 1 and _row > 1)
{
this->_type = COL_MATRIX;
}
// check if matrix is a UNIT-,DIAGONAL- or TRIANGULAR_MATRIX
else if (_square)
{
bool top = true;
bool down = true;
bool unit = true;
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
// check the top-right triangle
if (j > i)
{
if (_arr[i][j] != 0)
{
top = false;
}
}
// check the bottom-left triangle
if (j < i)
{
if (_arr[i][j] != 0)
{
down = false;
}
}
// check the diagonal
if (i == j)
{
if (_arr[i][j] != 1)
{
unit = false;
}
}
}
}
// Evaluation
if (top and down and unit)
{
this->_type = UNIT_MATRIX;
}
else if (top and down and !unit)
{
this->_type = DIAGONAL_MATRIX;
}
else if (top xor down)
{
this->_type = TRIANGULAR_MATRIX;
}
}
}
// Swapping rows
void Matrix::swapRows(usint aRow1, usint aRow2)
{
for (usint i = 0; i < _col; i++)
{
double tmp = _arr[aRow1][i];
_arr[aRow1][i] = _arr[aRow2][i];
_arr[aRow2][i] = tmp;
}
}
/**-------PUBLIC-------------------------------------------------*/
// Standard constructor
Matrix::Matrix(const usint aRow, const usint aCol)
: _row(aRow), _col(aCol), _precision(2), _square(aRow == aCol), _type(DEFAULT_MATRIX) // fast assignments
{
if (_row == 0 or _col == 0)
{
throw std::length_error("Matrix::Matrix() - rows and colums must be greater than zero !");
}
// allocate memory
_arr = new double*[_row];
for (usint i = 0; i < _row; i++) // rows gets columns
{
_arr[i] = new double[_col];
}
// set default value => Zero matrix
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] = 0.0;
}
}
}
// Copy constructor
Matrix::Matrix(const Matrix& aCopy)
: _row(aCopy._row), _col(aCopy._col), _precision(aCopy._precision),
_square(aCopy._square), _type(aCopy._type)
{
// allocate memory
_arr = new double*[_row];
for (usint i = 0; i < _row; i++)
{
_arr[i] = new double[_col];
}
// fill vector with values
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] = aCopy._arr[i][j];
}
}
}
// Destructor
Matrix::~Matrix()
{
// deallocate memory
for (usint i = 0; i < _row; i++)
{
delete[] _arr[_row];
}
delete[] _arr;
}
/**------- Public methods----------------------------------------*/
// Fill vector with values
void Matrix::fillIn()
{
double value(0);
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
std::cout << "[" <<i << "][" << j <<"]: ";
// try to get a valid floating number
while (!(std::cin >> value))
{
std::cout << "---> Wrong input. Please type only floating numbers ! \n";
std::cin.clear();
std::cin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
std::cout << "[" <<i << "][" << j <<"]: ";
}
// assign the value
_arr[i][j] = value;
}
}
std::cout << std::endl;
// set type
this->setType();
}
void Matrix::setPrecision(usint aVal)
{
_precision = aVal;
}
/**-------OPERATORS----------------------------------------------*/
//------------------- Assignment
const Matrix& Matrix::operator= (const Matrix& aRef)
{
if (this != &aRef) // Verhindert Selbstkopie
{
// call the destructor
this->~Matrix();
// allocate memory
_arr = new double*[_row];
for (usint i = 0; i < _row; i++)
{
_arr[i] = new double[_col];
}
// copy values
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] = aRef._arr[i][j];
}
}
}
// return reference
return *this;
}
//------------------- Negation
Matrix Matrix::operator- () const
{
// call copy constructor
Matrix temp(*this);
// negate each value
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
temp._arr[i][j] = -temp._arr[i][j];
}
}
// return copy
return temp;
}
//------------------- Operator to add two vectors together
const Matrix& Matrix::operator+= (const Matrix& aRef)
{
// vectors must have same structure
if (_row != aRef._row or _col != aRef._col)
throw std::length_error("Matrix::operator+=");
// add
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] += aRef._arr[i][j];
}
}
// return reference
return *this;
}
Matrix Matrix::operator+ (const Matrix& aRef) const
{
// vectors must have same structure
if (_row != aRef._row or _col != aRef._col)
throw std::length_error("Matrix::operator+");
Matrix temp(*this);
temp += aRef;
return temp;
}
const Matrix& Matrix::operator+= (double aVal)
{
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] += aVal;
}
}
return *this;
}
Matrix Matrix::operator+ (double aVal) const
{
Matrix temp(*this);
temp += aVal;
return temp;
}
//------------------- Operator for subtract two vectors
const Matrix& Matrix::operator-= (const Matrix& aRef)
{
// vectors must have same structure
if (_row != aRef._row or _col != aRef._col)
throw std::length_error("Matrix::operator-=");
// subtract
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] -= aRef._arr[i][j];
}
}
// return reference
return *this;
}
Matrix Matrix::operator- (const Matrix& aRef) const
{
// vectors must have same structure
if (_row != aRef._row or _col != aRef._col)
throw std::length_error("Matrix::operator-");
Matrix temp(*this);
temp -= aRef;
return temp;
}
const Matrix& Matrix::operator-= (double aVal)
{
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] -= aVal;
}
}
return *this;
}
Matrix Matrix::operator- (double aVal) const
{
Matrix temp(*this);
temp -= aVal;
return temp;
}
//--------------- Multiplication with "Falkschen Schema"
Matrix Matrix::operator* (const Matrix& aRef) const
{
// vectors must have special structure
if (_col != aRef._row)
throw std::length_error("Matrix::operator*");
// create temporary object
Matrix temp(_row, aRef._col);
for (usint i = 0; i < aRef._col; i++)
{
for (usint j = 0; j < _row; j++)
{
for (usint k = 0; k < _col; k++)
{
temp._arr[j][i] += (_arr[j][k] * aRef._arr[k][i]);
}
}
}
// return object
return temp;
}
const Matrix& Matrix::operator*= (double aVal)
{
for (usint i = 0; i < _row; i++)
{
for (usint j = 0; j < _col; j++)
{
_arr[i][j] *= aVal;
}
}
return *this;
}
Matrix Matrix::operator* (double aVal) const
{
Matrix temp(*this);
temp *= aVal;
return temp;
}
// Overloading the output operator
std::ostream& operator<< (std::ostream& aOut, Matrix& aRef)
{
// Abstand ermitteln [abst]
usint abst = aRef.distance();
for (usint i = 0; i < aRef._row; i++)
{
aOut << "|";
for (usint j = 0; j < aRef._col; j++)
{
aOut << std::setw(abst) << round(aRef._arr[i][j], aRef._precision);
}
aOut << " |\n";
}
aOut << std::endl;
return aOut;
}
/**-------Special functions--------------------------------------------*/
// try to create a form of triangle
Matrix Matrix::toTriangle()
{
// transforming only by squarish matrix
if (!this->_square)
{
throw std::logic_error("This matrix is not squarish, thatswhy it can't be transform to triangle.");
}
if (this->_type == TRIANGULAR_MATRIX)
{
return *this;
}
Matrix tmp(*this);
for (usint i = 0; i < _row -1; i++)
{
for (usint j = i+1; j < _row; j++)
{
if (tmp._arr[i][i] == 0)
{
tmp.swapRows(i,j);
}
double z = -(tmp._arr[j][i]/tmp._arr[i][i]);
for (usint k = i; k < _col; k++)
{
tmp._arr[j][k] += (tmp._arr[i][k] * z);
}
}
}
tmp._type = TRIANGULAR_MATRIX;
return tmp;
}
// Try to compute the determinant
const double Matrix::determinant()
{
// Only squarish matrix have a determinant
if (!_square)
{
throw std::logic_error("This matrix is not squarish, thatswhy they haven't a determinant.");
}
double deter;
if (this->_type == TRIANGULAR_MATRIX)
{
deter = _arr[0][0];
for (usint i = 1; i < _row; i++)
{
deter *= _arr[i][i];
}
}
else
{
Matrix tmp = this->toTriangle();
deter = tmp._arr[0][0];
for (usint i = 1; i < tmp._row; i++)
{
deter *= tmp._arr[i][i];
}
}
return deter;
}
// Get the rang of a matrix
const usint Matrix::rang()
{
if (!_square)
{
throw std::logic_error("This program can only compute the rang of a squarish matrix.");
}
Matrix tmp = this->toTriangle();
usint rg = 0;
for (usint i = 0; i < tmp._row; i++)
{
double count = 0;
for (usint j = 0; j < tmp._col; j++)
{
count += tmp._arr[i][j];
}
if (count != 0.0)
{
rg++;
}
}
return rg;
}
// Try to transform to inverse matrix with "adjunkte"
Matrix Matrix::toInverse()
{
// if matrix have no square format then dont exits the inverse
if (!_square)
{
throw std::logic_error("This matrix is not squarish, thatswhy it haven't an inverse matrix.");
}
const double det = this->determinant();
if (!det)
{
throw std::logic_error("The determinant is zero so it dont't exists an inverse matrix.");
}
// create inverse matrix
Matrix inverseMatrix(_row, _col);
// create temporary matrix for computing the complementary matrix
usint size = _row - 1;
Matrix tmp(size, size);
// loop for every row
for (usint i = 0; i < _row; i++)
{
// loop for every column
for (usint j = 0; j < _col; j++)
{
// position counter for collected values (which are not striked)
usint x = 0, y = 0;
// interior loop for striking row and column
for (usint k = 0; k < _row; k++)
{
for (usint l = 0; l < _col; l++)
{
if (i != k and j != l)
{
tmp._arr[x][y] = this->_arr[k][l];
y++;
if (y == size)
{
y = 0;
x++;
}
}
}
}
// complementary matrix value
inverseMatrix._arr[j][i] = (std::pow(-1.0, j+i) * tmp.determinant()) / det;
}
}
return inverseMatrix;
}
// Get the matrix type as string
std::string Matrix::getType()
{
switch(_type)
{
case ROW_MATRIX: return "ROW-Matrix";
case COL_MATRIX: return "COLUMN-Matrix";
case UNIT_MATRIX: return "UNIT-Matrix";
case DIAGONAL_MATRIX: return "DIAGONAL-Marix";
case TRIANGULAR_MATRIX: return "TRIANGULAR-Matrix";
default: return "DEFAULT-Matrix";
}
}
///--------------------- extern function---------------------------
// Round the value with precision
const double round(const double value, const int precision)
{
const double multiplier(std::pow(10.0, precision));
return std::floor(value * multiplier + 0.5) / multiplier;
}
