Homec4science
Diffusion PCSC_project (master) Browse

R7934/Final/Vector_Matrix.cpp
a9bd6f84d4a7master

Last Change
Raw File

Vector_Matrix.cpp
View Options

#include <cmath>
#include <iostream>
#include <cassert>
#include <string>
#include "Vector_Matrix.hpp"
#include "LinSys.hpp"
/**************************************
***************************************
***************************************
****** FUNCTION PROTOTYPE *************
***************************************
***************************************
**************************************/
double** AllocateMemory(int rows, int cols);
/**************************************
***************************************
***************************************
************* MATRIX CLASS ************
***************************************
***************************************
***************************************
***************************************
**************************************/
// Overridden copy constructor
// Allocates memory for new matrix, and copies
// entries of other matrix into it
Matrix::Matrix(const Matrix& otherMatrix)
{
/* Copying of Matrix size*/
mNumCols = otherMatrix.NumberOfColumns();
mNumRows = otherMatrix.NumberOfRows();
/*Allocation of memory for copy of matrix*/
mData=AllocateMemory(mNumRows,mNumCols) ;
/* Copying of class data*/
for (int i=0; i<mNumRows; i++){
for (int j=0;j<mNumCols;j++){
mData[i][j]=otherMatrix.mData[i][j];
}
}
}
// Constructor for matrix of a given size
// Allocates memory, and initialises entries
// to zero
Matrix::Matrix(int rows, int cols){
//assert((rows > 0 && cols>0 )&& "Define positive integer values for both arguments defining the size of the matrix.");
mNumCols = cols;
mNumRows=rows;
mData=AllocateMemory(mNumRows,mNumCols) ;
for (int i=0; i<mNumRows; i++) {
for (int j = 0; j < mNumCols; j++) {
mData[i][j] = 0.0;
}
}
}
// Constructor for a square matrix of a given size
// Allocates memory, and initialises entries
// to zero
Matrix::Matrix(int size){
assert(size > 0 && "\nDefine a positive integer value for the square matrix size.");
mNumCols = size;
mNumRows=size;
mData=AllocateMemory(mNumRows,mNumCols) ;
for (int i=0; i<mNumRows; i++) {
for (int j = 0; j < mNumCols; j++) {
mData[i][j] = 0.0;
}
}
}
Matrix::Matrix() {}
/*Matrix::Matrix(){
mNumCols = 2;
mNumRows=2;
mData=AllocateMemory(mNumRows,mNumCols) ;
for (int i=0; i<mNumRows; i++) {
for (int j = 0; j < mNumCols; j++) {
mData[i][j] = 0.0;
}
}
}*/
/*/// Overridden destructor to correctly free memory
Matrix::~Matrix()
{
for (int i=0;i<mNumRows;i++){
delete[] mData[i];
}
delete[] mData;
}*/
// Method to get the number of rows of a matrix
int Matrix::NumberOfRows() const
{
return mNumRows;
}
// Method to get the number of columns of a matrix
int Matrix::NumberOfColumns() const
{
return mNumCols;
}
void Matrix::SetSize(const int rows, const int cols) {
mNumRows=rows;
mNumCols=cols;
mData=AllocateMemory(rows,cols);
}
// Method to check if index is out of bounds
// Method to check if index is out of bounds
int Matrix::CheckIndex(int i_index, int i_val) const{
try {
if (i_val<=-1 || i_val>=mNumRows){
if (i_index==0){
throw (Exception("Index i ","Index i of A(i,j) is out of bounds. (1-based indexing)"));
} else {
throw (Exception("Index j ","Index j of A(i,j) is out of bounds. (1-based indexing)"));
}
} else {
return i_val;
}
} catch (Exception& error){
error.PrintDebug();
std::cout<<"\n\nIndex error.\nSee error report and give alternative index value\n\n";
int j;
std::cin>>j;
while (j<=0 || j>mNumRows){
std::cout<<"\nIndexes are still out of bounds, give alternative values ...\n";
std::cin>>j;
}
return j-1;
}
}
// Read-only variant of [][]
// Note that this uses ‘zero-based’ indexing,
// and a check on the validity of the index
double Matrix::Read(int i, int j) const
{
return mData[CheckIndex(0,i)][CheckIndex(1,j)];
}
// Overloading round brackets
// Note that this uses ‘one-based’ indexing,
// and a check on the validity of the index
double& Matrix::operator()(int i, int j) const
{
return mData[CheckIndex(0,i-1)][CheckIndex(1,j-1)];
}
// Overloading the assignment operator
Matrix& Matrix::operator=(const Matrix& otherVector)
{
assert((mNumRows == otherVector.mNumRows && mNumCols==otherVector.mNumCols) && "Matrices are of different sizes. Cant assign");
for (int i=0; i<mNumRows; i++)
{
for (int j=0;j<mNumCols;j++)
{
mData[i][j] = otherVector.mData[i][j];
}
}
return * this;
}
// Overloading the unary + operator
Matrix Matrix::operator+() const {
Matrix A(mNumRows, mNumCols);
for (int i = 0; i < mNumRows; i++) {
for (int j = 0; j < mNumCols; j++) {
A(i + 1, j + 1) = mData[i][j];
}
}
return A;
}
// Overloading the unary - operator
Matrix Matrix::operator-() const {
Matrix A(mNumRows, mNumCols);
for (int i = 0; i < mNumRows; i++) {
for (int j = 0; j < mNumCols; j++) {
A(i + 1, j + 1) = -mData[i][j];
}
}
return A;
}
// Overloading the binary + operator
Matrix Matrix::operator+(const Matrix& B) const
{
assert((mNumRows == B.mNumRows && mNumCols==B.mNumCols) && "Matrices are of different sizes. Cant perform addition.");
Matrix A(mNumRows, mNumCols);
for (int i = 0; i < mNumRows; i++) {
for (int j = 0; j < mNumCols; j++) {
A.mData[i][j] = mData[i][j]+B.mData[i][j];
}
}
return A;
}
// Overloading the binary - operator
Matrix Matrix::operator-(const Matrix& B) const
{
assert((mNumRows == B.mNumRows && mNumCols==B.mNumCols )&& "Matrices are of different sizes. Cant perform subtraction.");
Matrix A(mNumRows, mNumCols);
for (int i = 0; i < mNumRows; i++) {
for (int j = 0; j < mNumCols; j++) {
A.mData[i][j] = mData[i][j]-B.mData[i][j];
}
}
return A;
}
// Overloading the binary * operator
Matrix Matrix::operator*(const Matrix& B) const
{
assert((mNumRows == B.mNumRows && mNumCols==B.mNumCols )&& "Matrices are of different sizes. Cant perform multiplication.");
Matrix A(mNumRows, mNumCols);
for (int i=0;i<mNumRows;i++){
for (int j=0;j<B.mNumCols;j++){
for (int k=0;k<B.mNumCols;k++){
A.mData[i][j]+=mData[i][k]*B.mData[k][j];
}
}
}
return A;
}
// Overloading the output << operator
std::ostream& operator<<(std::ostream& output, const Matrix& A){
/*Prints to terminal values of a matrix A*/
for (int i=0;i<A.mNumRows;i++){
for (int j=0;j<A.mNumCols;j++){
std::cout<<A.mData[i][j]<<" ";
}
std::cout<<std::endl;
}
std::cout<<std::endl;
}
// Overloading scalar multiplication
Matrix Matrix::operator *(double a) const
{
Matrix A(mNumRows, mNumCols);
for (int i=0;i<mNumRows;i++){
for (int j=0;j<mNumCols;j++){
A.mData[i][j]=mData[i][j]*a;
}
}
return A;
}
// Overloading vector multiplication
Vector Matrix::operator*(const Vector& u) const{
assert(mNumRows==u.GetSize()&& "Matrix and vector are of different sizes. Cant perform multiplication.");
Vector v(u.GetSize());
for (int i=0;i<mNumRows;i++){
for (int j=0;j<mNumCols;j++){
v[i]+=mData[i][j]*u[j];
}
}
return v;
}
/*Computation of determinant : direct formula for 2x2, recursive otherwise*/
double Matrix::Determinant(){
/*Computation of determinant : direct formula for 2x2, recursive otherwise*/
if (mNumRows==2){
return mData[0][0]*mData[1][1]-mData[1][0]*mData[0][1];
} else {
double det=0;
for (int j = 0; j < mNumRows; ++j) {
Matrix A_hat(mNumRows-1);
ReduceMatrix(mData,mNumRows,j,A_hat);
det +=mData[0][j]*std::pow(-1.0,j)*A_hat.Determinant();
A_hat.FreeMemory();
}
return det;
}
}
// Resizing of matrix A for determinant computation
void Matrix::ReduceMatrix(double** A,int size,int col, const Matrix &A_hat){
for (int row = 1; row < size; ++row) {
for (int k = 0; k < col; ++k) {
A_hat(row,k+1) = A[row][k];
}
for (int k = col + 1; k < size; ++k) {
A_hat(row,k) = A[row][k];
}
}
}
void Matrix::FreeMemory() {
/*Frees memory of matrix A*/
for (int i = 0; i < mNumRows; i++) {
delete[] mData[i];
}
delete[] mData;
}
// Allocates memory for matrix A
double** Matrix::AllocateMemory(int rows, int cols) const {
double** A;
A=new double* [rows];
for (int i=0;i<cols;i++){
A[i]=new double[cols];
}
return A;
}
/**************************************
***************************************
***************************************
************* VECTOR CLASS ************
***************************************
***************************************
***************************************
***************************************
**************************************/
// Overridden copy constructor
// Allocates memory for new vector, and copies
// entries of other vector into it
Vector::Vector(const Vector& otherVector)
{
mSize = otherVector.GetSize();
mData = new double [mSize];
for (int i=0; i<mSize; i++)
{
mData[i] = otherVector.mData[i];
}
}
// Constructor for vector of a given size
// Allocates memory, and initialises entries
// to zero
Vector::Vector(int size)
{
assert(size > 0&& "Assign positive integer value for vector size.");
mSize = size;
mData = new double [mSize];
for (int i=0; i<mSize; i++)
{
mData[i] = 0.0;
}
}
Vector::Vector(){}
// Overridden destructor to correctly free memory
Vector::~Vector(){}
// Method to get the size of a vector
int Vector::GetSize() const
{
return mSize;
}
void Vector::SetSize(const int size) {
mSize = size;
mData = new double[mSize];
}
// Method to check if index is out of bounds
int Vector::CheckIndex(int i_index,int i_val) const{
try {
if (i_val<=-1 || i_val>=mSize){
if (i_index==0){
throw (Exception(" i","Index i of vector u[i] is out of bounds. (0-based indexing)"));
}else {
throw (Exception("Index i","Index i of vector u[i] is out of bounds. (1-based indexing)"));
}
} else {
return i_val;
}
} catch (Exception& error){
error.PrintDebug();
std::cout<<"\nIndex error. See error report and give an alternative index value\n\n";
int j;
std::cin>>j;
while (j<=-1 || j>=mSize){
std::cout<<"\nIndex is still out of bounds, give an alternative value ...\n";
std::cin>>j;
}
return j;
}
}
// Overloading square brackets
// Note that this uses ‘zero-based’ indexing,
// and a check on the validity of the index
double& Vector::operator[](int i) const
{
return mData[CheckIndex(0,i)];
}
// Read-only variant of []
// Note that this uses ‘zero-based’ indexing,
// and a check on the validity of the index
double Vector::Read(int i) const
{
return mData[CheckIndex(0,i)];
}
// Overloading round brackets
// Note that this uses ‘one-based’ indexing,
// and a check on the validity of the index
double& Vector::operator()(int i)
{
return mData[CheckIndex(1,i-1)];
}
// Overloading the assignment operator
Vector& Vector::operator=(const Vector& otherVector)
{
//assert(mSize == otherVector.mSize&& "Vectors are of different sizes. Cant assign.");
for (int i=0; i<mSize; i++)
{
mData[i] = otherVector.mData[i];
}
return * this;
}
// Overloading the unary + operator
Vector Vector::operator+() const
{
Vector v(mSize);
for (int i=0; i<mSize; i++)
{
v[i] = mData[i];
}
return v;
}
// Overloading the unary - operator
Vector Vector::operator-() const
{
Vector v(mSize);
for (int i=0; i<mSize; i++)
{
v[i] = -mData[i];
}
return v;
}
// Overloading the binary + operator
Vector Vector::operator+(const Vector& v1) const
{
assert(mSize == v1.mSize&& "Vectors are of different sizes. Cant perform addition.");
Vector v(mSize);
for (int i=0; i<mSize; i++)
{
v[i] = mData[i] + v1.mData[i];
}
return v;
}
// Overloading the binary - operator
Vector Vector::operator-(const Vector& v1) const
{
assert(mSize == v1.mSize&& "Vectors are of different sizes. Cant perform subtraction.");
Vector v(mSize);
for (int i=0; i<mSize; i++)
{
v[i] = mData[i] - v1.mData[i];
}
return v;
}
// Overloading the binary * operator
Vector Vector::operator*(const Vector& v1) const {
assert(mSize == v1.mSize && "Vectors are of different sizes. Cant perform multiplication.");
Vector v(mSize);
for (int i=0; i<mSize; i++)
{
v[i] = mData[i] * v1.mData[i];
}
return v;
}
// Overloading the scalar ^ operator
double Vector::operator^(const Vector& v1) const {
assert(mSize == v1.mSize && "Vectors are of different sizes. Cant perform scalar product.");
double scalar=0;
for (int i=0; i<mSize; i++)
{
scalar += mData[i] * v1.mData[i];
}
return scalar;
}
// Overloading the binary / operator
Vector Vector::operator/(const Vector& v1) const {
assert(mSize == v1.mSize && "Vectors are of different sizes. Cant perform element by element division.");
Vector v(mSize);
for (int i=0; i<mSize; i++)
{
v[i] = mData[i] / v1.mData[i];
}
return v;
}
// Overloading the output << operator
std::ostream& operator<<(std::ostream& output, const Vector& v){
//Prints values of vector
for (int i=0;i<v.GetSize();i++) {
std::cout << v.mData[i]<<std::endl;
}
std::cout<<std::endl;
}
// Overloading scalar multiplication
Vector Vector::operator *(double a) const
{
Vector v(mSize);
for (int i=0; i<mSize; i++)
{
v[i] = a*mData[i];
}
return v;
}
Vector Vector::operator*(const Matrix& A) const{
/*Multiplies vector U[rowsU] by matrix A[rowsA][colsA]*/
assert(mSize==A.NumberOfRows() && "Vector and Matrix are of different sizes. Cant perform multiplication u*A.");
Vector v(A.NumberOfRows());
for (int i=0;i<A.NumberOfRows();i++){
for (int j=0;j<A.NumberOfColumns();j++){
v[i]+=A(j+1,i+1)*mData[j];
}
}
return v;
}
void Vector::Fill(int val) {
for (int i=0;i<mSize;i++){
mData[i]=double(val);
}
}
void Vector::Fill(double val) {
for (int i=0;i<mSize;i++){
mData[i]=val;
}
}
// Method to calculate norm (with default value p=2)
// corresponding to the Euclidean norm
double Vector::CalculateNorm(int p) const
{
double norm_val, sum = 0.0;
for (int i=0; i<mSize; i++)
{
sum += pow(fabs(mData[i]), p);
}
norm_val = pow(sum, 1.0/((double)(p)));
return norm_val;
}
// MATLAB style friend to get the size of a vector
int length(const Vector& v)
{
return v.mSize;
}
/**************************************
***************************************
***************************************
************ MISC. FUNCTIONS **********
***************************************
***************************************
***************************************
***************************************
**************************************/
void TerminalPrint(Vector v,std::string msg) {
/*Prints to terminal values of a vector v and a message */
std::cout<<std::endl<<msg<<std::endl<<v;
}
/*Prints to terminal values of a vector v and a message */
void TerminalPrint(Matrix A,std::string msg) {
std::cout<<std::endl<<msg<<std::endl<<A;
}
/*Prints to terminal values of a scalar and a message */
void TerminalPrint(double scalar,std::string msg) {
std::cout<<std::endl<<msg<<std::endl<<scalar<<std::endl<<std::endl;;
}
/*Prints to terminal values of an int and a message */
void TerminalPrint(int val,std::string msg) {
std::cout<<std::endl<<msg<<std::endl<<val<<std::endl<<std::endl;;
}
/**************************************
***************************************
***************************************
************ ERROR HANDLING************
***************************************
***************************************
***************************************
***************************************
**************************************/
//Constructor
Exception::Exception(std::string tagString,
std::string probString)
{
mTag = tagString;
mProblem = probString;
}
// Terminal print
void Exception::PrintDebug() const
{
std::cerr << "** Error ("<<mTag<<") **\n";
std::cerr << "Problem: " << mProblem << "\n\n";
}
c4science · Help