Homec4science
Diffusion PCSC_project (master) Browse

R7934/Final_v1.4/LinSys.cpp
d46187a89a17master

Last Change
Raw File

LinSys.cpp
View Options

/*!
* \file LinSys.cpp
* \brief Functions associated to the different methods of resolution
*
* \author Nabil & Xavier
* \version 1.0
* \date 14.12.2018
*
* */
/*********************************************************/
#include <cmath>
#include <iostream>
#include <cassert>
#include <string>
#include "LinSys.hpp"
using namespace std;
/*********************************************************/
/**************************************
************* LINSYS CLASS ************
**************************************/
//Default constructor
Lin_Sys::Lin_Sys() {}
//Constructor
Lin_Sys::Lin_Sys(Matrix A ,Vector b)
{
mSize=A.NumberOfRows();
mpb.SetSize(mSize);
mpA.SetSize(mSize,mSize);
mpb=b;
mpA=A;
}
//Another constructor
Lin_Sys::Lin_Sys(Matrix A ,Vector b,Vector x_0)
{
mSize=A.NumberOfRows();
mpb.SetSize(mSize);
mpA.SetSize(mSize,mSize);
mpb=b;
mpA=A;
}
//Another constructor
Lin_Sys::Lin_Sys(Matrix A ,Vector b,Vector x_0,double tol)
{
mSize=A.NumberOfRows();
mpb.SetSize(mSize);
mpA.SetSize(mSize,mSize);
mpb=b;
mpA=A;
}
//Another constructor
Lin_Sys::Lin_Sys(Matrix A ,Vector b,Vector x_0,double tol,int max_iter)
{
mSize=A.NumberOfRows();
mpb.SetSize(mSize);
mpA.SetSize(mSize,mSize);
mpb=b;
mpA=A;
}
//Default destructor
Lin_Sys::~Lin_Sys() {}
//Function to get the left hand side
Matrix Lin_Sys::GetMatrix() { return mpA;}
//Function to get the right hand side
Vector Lin_Sys::GetRightHandSide() { return mpb;}
//Function to get the size of the system
int Lin_Sys::GetSize(){ return mSize;}
//Overloaded solve function
Vector Lin_Sys::Solve () {}
/**************************************
********* EXACT METHOD CLASS **********
**************************************/
//Constructor
Exact::Exact(Matrix A,Vector b) : Lin_Sys(A,b)
{
int size = A.NumberOfRows();
U.SetSize(size, size);
L.SetSize(size, size);
P.SetSize(size,size);
b_star.SetSize(size);
z.SetSize(size);
// Initialization of matrix L and P
for (int i = 0; i < L.NumberOfRows(); i++) {
L(i + 1, i + 1) = 1;
P(i+1,i+1)=1;
}
// Initialization of matrix U
U = A;
// Initialization of vector b*
b_star = b;
// Partial Pivoting of U
double temp1,temp2;
for (int i = size; i > 1; i--)
{
if (U(i-1,1) < U(i,1)) {
for (int j = 1; j <= size; j++)
{
temp1 = U(i, j);
U(i, j) = U(i - 1, j);
U(i - 1, j) = temp1;
temp1 = P(i, j);
P(i, j) = P(i - 1, j);
P(i - 1, j) = temp1;
}
temp2 = b_star(i);
b_star(i) = b_star(i - 1);
b_star(i - 1) = temp2;
}
}
}
//Default destructor
Exact::~Exact(){}
//Function to get the upper triangular decomposition
Matrix Exact::get_U() {return U;}
//Function to get the lower triangular decomposition
Matrix Exact::get_L() {return L;}
//Function to get the permutation matrix
Matrix Exact::get_P() {return P;}
//Function to get the pivoted right hand side
Vector Exact::get_b_star(){return b_star;}
//Overloaded solve function
Vector Exact::Solve() {return z;}
/**************************************
************* LU CLASS ************
**************************************/
//Default constructor
LU::LU(Matrix A,Vector b) : Exact(A,b)
{
// Gauss Jordan Pivoting
for (int i = 0; i < GetSize() - 1; i++)
{
for (int j = i + 1; j < GetSize(); j++)
{
L(j + 1, i + 1) = U(j + 1, i + 1) / U(i + 1, i + 1);
for (int k = 0; k < GetSize(); k++)
{
U(j+1,k+1) -= U(i+1,k+1) * L(j + 1, i + 1);
}
b_star[j]-=b_star[i]*L(j + 1, i + 1);
}
}
}
//Default Destructor
LU::~LU(){}
//Overloaded solve function
Vector LU::Solve()
{
z.Fill(0);
for (int i=GetSize()-1;i>-1;i--)
{
z[i]+=b_star[i];
for (int j=i+1;j<GetSize();j++)
{
z[i]-=U(i+1,j+1)*z[j];
}
z[i]/=U(i+1,i+1);
}
return z;
}
/**************************************
************* CHOLESKY CLASS **********
**************************************/
//Constructor
Cholesky::Cholesky(Matrix A,Vector b) : Exact(A,b)
{
// Construction of L matrix
for (int i = 0; i < GetSize(); i++)
{
for (int k = 0; k < (i + 1); k++)
{
double sum = 0;
for (int j = 0; j < k; j++) {
sum += (L(i + 1, j + 1) * L(k + 1, j + 1));
}
if (i == k)
{
L(i + 1, k + 1) = sqrt(A(i + 1, i + 1) - sum);
} else
{
L(i + 1, k + 1) = (1.0 / L(k + 1, k + 1)) * (A(i + 1, k + 1) - sum);
}
}
}
//Construction of L transpose matrix
for (int i = 0; i < GetSize(); i++)
{
for (int j = 0; j < GetSize(); j++)
{
U(i + 1, j + 1) = L(j + 1, i + 1);
}
}
//Construction of b* vector
b_star.Fill(0);
for (int i=0;i<GetSize();i++)
{
b_star[i]+=GetRightHandSide()[i];
for (int j=i-1;j>=0;j--)
{
b_star[i]-=L(i+1,j+1)*b_star[j];
}
b_star[i]/=L(i+1,i+1);
}
}
//Destructor
Cholesky::~Cholesky() {}
//Overloaded solve function
Vector Cholesky::Solve() {
z.Fill(0);
for (int i=b_star.GetSize()-1;i>-1;i--)
{
z[i]+=b_star[i];
for (int j=i+1;j<b_star.GetSize();j++)
{
z[i]-=U(i+1,j+1)*z[j];
}
z[i]/=U(i+1,i+1);
}
return z;
}
/**************************************
********* APPROX METHOD CLASS *********
**************************************/
//Default constructor
Approx::Approx(Matrix A,Vector b) : Lin_Sys(A,b)
{
int size = A.NumberOfRows();
x0.SetSize(size);x_old.SetSize(size);
x_new.SetSize(size);e.SetSize(size);
tol = 1e-12;
n_iter_max = 5000;
convergence.SetSize(n_iter_max);
}
//Another constructor
Approx::Approx(Matrix A,Vector b,Vector x_0) : Lin_Sys(A,b,x_0)
{
int size = A.NumberOfRows();
x0.SetSize(size);x_old.SetSize(size);
x_new.SetSize(size);e.SetSize(size);
x0=x_0;
tol = 1e-12;
n_iter_max = 5000;
convergence.SetSize(n_iter_max);
}
//Another constructor
Approx::Approx(Matrix A,Vector b,Vector x_0,double tolerance) : Lin_Sys(A,b,x_0,tolerance)
{
int size = A.NumberOfRows();
x0.SetSize(size);x_old.SetSize(size);
x_new.SetSize(size);e.SetSize(size);
x0=x_0;
tol = tolerance;
n_iter_max = 5000;
convergence.SetSize(n_iter_max);
}
//Another constructor
Approx::Approx(Matrix A,Vector b,Vector x_0,double tolerance,int max_iter) : Lin_Sys(A,b,x_0,tolerance,max_iter)
{
int size = A.NumberOfRows();
x0.SetSize(size);x_old.SetSize(size);
x_new.SetSize(size);e.SetSize(size);
x0=x_0;
tol = tolerance;
n_iter_max = max_iter;
convergence.SetSize(n_iter_max);
}
//Destructor
Approx::~Approx(){}
//Function to get the number of iteration for convergence
int Approx::get_n_iter() {return n_iter;}
//Function to get the maximum number of iteration for convergence
int Approx::get_n_iter_max() {return n_iter_max;}
//Function to get a vector containing the evolution of the error
Vector Approx::get_convergence() {return convergence;}
//Function to set the tolerance to a user_defined value
void Approx::set_tol(double my_tol){tol=my_tol;}
//Function to set the maximum number of iterations allowed to a user_defined value
void Approx::set_iter_max(int my_iter_max){n_iter_max=my_iter_max;}
//Overloaded solve function
Vector Approx::Solve() {return x_new;}
/**************************************
***** Conjugate_gradient CLASS ********
**************************************/
//Constructor
Conjugate_Gradient::Conjugate_Gradient(Matrix A,Vector b): Approx(A,b)
{
int size = A.NumberOfRows();
r_old.SetSize(size);
p_new.SetSize(size);
p_old.SetSize(size);
}
//Constructor
Conjugate_Gradient::Conjugate_Gradient(Matrix A,Vector b,Vector x_0): Approx(A,b,x_0)
{
int size = A.NumberOfRows();
r_old.SetSize(size);
p_new.SetSize(size);
p_old.SetSize(size);
x0=x_0;
}
//Constructor
Conjugate_Gradient::Conjugate_Gradient(Matrix A,Vector b,Vector x_0,double tolerance): Approx(A,b,x_0,tolerance)
{
int size = A.NumberOfRows();
r_old.SetSize(size);
p_new.SetSize(size);
p_old.SetSize(size);
x0=x_0;
}
//Constructor
Conjugate_Gradient::Conjugate_Gradient(Matrix A,Vector b,Vector x_0,double tolerance,int max_iter): Approx(A,b,x_0,tolerance,max_iter)
{
int size = A.NumberOfRows();
r_old.SetSize(size);
p_new.SetSize(size);
p_old.SetSize(size);
x0=x_0;
}
//Destructor
Conjugate_Gradient::~Conjugate_Gradient() {}
//Overloaded solve function
Vector Conjugate_Gradient::Solve()
{
n_iter=0;
x_new=x0;
double div_counter=0;
r_old= GetRightHandSide()-GetMatrix()*x0;
p_old= r_old;
e.Fill(1);
while((e^e)>tol && n_iter<n_iter_max)
{
double temp = r_old^r_old;
double alpha = temp /(p_old^(GetMatrix()*p_old));
x_new =x_new+p_old*alpha;
e=r_old-GetMatrix()*p_old*alpha;
double beta = (e^e)/temp;
p_new=e+p_old*beta;
r_old=e;p_old=p_new;
convergence[n_iter]=e^e;
if (n_iter>0 && div_counter>n_iter/2) {div_counter++;}
n_iter++;
}
if (div_counter>n_iter/10)
{
std::cout<<"Conjugate_Gradient method failed to converge. You can use the Convergence command to analyze divergence\n"
"Try with another set of parameters or change method.\n";
x_new.Fill(0);
} else
{
std::cout<<"Conjugate_Gradient method succeeded in converging"<<std::endl;
}
return x_new;
}
/**************************************
************* Gauss Seidel CLASS ******
**************************************/
//Constructor
Gauss_Seidel::Gauss_Seidel(Matrix A,Vector b): Approx(A,b) {}
//Constructor
Gauss_Seidel::Gauss_Seidel(Matrix A,Vector b,Vector x_0): Approx(A,b,x_0) {x0=x_0;}
//Constructor
Gauss_Seidel::Gauss_Seidel(Matrix A,Vector b,Vector x_0,double tolerance): Approx(A,b,x_0,tolerance) {
x0=x_0;
tol=tolerance;
}
//Constructor
Gauss_Seidel::Gauss_Seidel(Matrix A,Vector b,Vector x_0,double tolerance,int max_iter): Approx(A,b,x_0,tolerance,max_iter) {
x0=x_0;
tol=tolerance;
n_iter_max=max_iter;
}
//Destructor
Gauss_Seidel::~Gauss_Seidel() {}
//Overloaded solve function
Vector Gauss_Seidel::Solve()
{
//Initialization
int size=GetSize();
n_iter=0;
e=GetRightHandSide()-GetMatrix()*x0;
double div_counter=0;
while((e^e)>tol && n_iter<n_iter_max)
{
for (int i = 1; i <= size; ++i)
{
x_new(i) = GetRightHandSide()(i) / GetMatrix()(i,i);
for (int j = 1; j <= size; ++j)
{
if (i == j)
{
} else {
x_new(i)-= (GetMatrix()(i,j) * x_new(j))/GetMatrix()(i,i);
}
}
}
x_old=x_new;
e=GetRightHandSide()-GetMatrix()*x_new;
convergence[n_iter]=e^e;
if (n_iter>0 && convergence[n_iter]-convergence[n_iter-1]>0) {div_counter++;}
n_iter++;
}
if (div_counter>n_iter/10)
{
std::cout<<"Gauss Seidel method failed to converge. You can use the Convergence command to analyze divergence\n"
"Try with another set of parameters or change method.\n";
x_new.Fill(0);
} else
{
std::cout<<"Gauss Seidel method succeeded in converging"<<std::endl;
}
return x_new;
}
/**************************************
************* Jacobi CLASS ************
**************************************/
//Constructors
Jacobi::Jacobi(Matrix A,Vector b): Approx(A,b) {}
//Constructor
Jacobi::Jacobi(Matrix A,Vector b,Vector x_0): Approx(A,b,x_0) {x0=x_0;}
//Constructor
Jacobi::Jacobi(Matrix A,Vector b,Vector x_0,double tolerance): Approx(A,b,x_0,tolerance) {
x0=x_0;
tol=tolerance;
}
//Constructor
Jacobi::Jacobi(Matrix A,Vector b,Vector x_0,double tolerance,int max_iter): Approx(A,b,x_0,tolerance,max_iter) {
x0=x_0;
tol=tolerance;
n_iter_max=max_iter;
}
//Destructor
Jacobi::~Jacobi() {}
//Overloaded solve function
Vector Jacobi::Solve()
{
//Initialization
int size=GetSize();
n_iter=0;
e=GetRightHandSide()-GetMatrix()*x0;
double div_counter=0;
while((e^e)>tol && n_iter<n_iter_max)
{
for (int i = 1; i <= size; ++i)
{
x_new(i) = GetRightHandSide()(i) / GetMatrix()(i,i);
for (int j = 1; j <= size; ++j)
{
if (i == j)
{
} else {
x_new(i)-= (GetMatrix()(i,j) * x_old(j))/GetMatrix()(i,i);
}
}
}
x_old=x_new;
e=GetRightHandSide()-GetMatrix()*x_new;
convergence[n_iter]=e^e;
if (n_iter>0 && convergence[n_iter]-convergence[n_iter-1]>0) {div_counter++;}
n_iter++;
}
if (div_counter>n_iter/10)
{
std::cout<<"Jacobi method failed to converge. You can use the Convergence command to analyze divergence\n"
"Try with another set of parameters or change method.\n";
x_new.Fill(0);
} else
{
std::cout<<"Jacobi method succeeded in converging"<<std::endl;
}
return x_new;
}
/**************************************
************* Richardson CLASS ********
**************************************/
//Constructor
Richardson::Richardson(Matrix A,Vector b): Approx(A,b) {}
//Constructor
Richardson::Richardson(Matrix A,Vector b,Vector x_0): Approx(A,b) {x0=x_0;}
//Constructor
Richardson::Richardson(Matrix A,Vector b,Vector x_0,double tolerance): Approx(A,b,x_0,tolerance) {
x0=x_0;
tol=tolerance;
}
//Constructor
Richardson::Richardson(Matrix A,Vector b,Vector x_0,double tolerance,int max_iter): Approx(A,b,x_0,tolerance,max_iter) {
x0=x_0;
tol=tolerance;
n_iter_max=max_iter;
}
//Destructor
Richardson::~Richardson() {}
//Overloaded solve function
Vector Richardson::Solve()
{
//Initialization
int size=GetSize();
n_iter=0;
e=GetRightHandSide()-GetMatrix()*x0;
double div_counter=0;
while((e^e)>tol && n_iter<n_iter_max && (e^e)<1e10)
{
for (int i = 1; i <= size; ++i)
{
e(i)=GetRightHandSide()(i)/GetMatrix()(i,i);
x_new(i) = x_old(i)+e(i);
for (int j = 1; j <= size; ++j)
{
e(i)-=(GetMatrix()(i,j) * x_old(j))/GetMatrix()(i,i);
x_new(i)-= (GetMatrix()(i,j) * x_old(j))/GetMatrix()(i,i);
}
}
x_old=x_new;
e=GetRightHandSide()-GetMatrix()*x_new;
convergence[n_iter]=e^e;
if (n_iter>0 && convergence[n_iter]-convergence[n_iter-1]>0) {
div_counter++;
}
n_iter++;
}
if (div_counter>n_iter/10)
{
std::cout<<"Richardson method failed to converge. You can use the Convergence command to analyze divergence\n"
"Try with another set of parameters or change method.\n";
x_new.Fill(0);
} else
{
std::cout<<"Richardson method succeeded in converging"<<std::endl;
}
return x_new;
}
c4science · Help