Homec4science
Diffusion PCSC_project (master) Browse

R7934/Final/io.cpp
073060d11c24master

Last Change
Raw File

io.cpp
View Options

#include <iostream>
#include <cassert>
#include <cmath>
#include "io.hpp"
#include "LinSys.hpp"
using namespace std;
Lin_sys_struct::Lin_sys_struct(int msize) {
A.SetSize(msize,msize);
b.SetSize(msize);
size=msize;
}
Lin_sys_struct Initialize(){
// Initialization
cout<<"\nThis program allows solving linear systems of the type Ax=b.\nStart by giving the size of the linear system : ";
int size;
cin>>size;
Matrix A=Matrix(size);
Vector b=Vector(size);
Lin_sys_struct instance(size);
// Matrix A
cout<<"\nLets fill matrix A."<<endl;
for (int i=0;i<size;i++){
for (int j=0;j<size;j++){
cout<<"Give value of A("<<i+1<<","<<j+1<<") : ";
cin>>A(i+1,j+1);
}
}
TerminalPrint(A,"Matrix A :");
//Vector b
cout<<"\nLets fill vector B."<<endl;
for (int i=0;i<size;i++){
cout<<"Give value of b("<<i+1<<") : ";
cin>>b(i+1);
}
TerminalPrint(b,"Vector b :");
instance.A=A;
instance.b=b;
return instance;
}
string method(){
// Description of methods available
cout<<"Description of available methods:\n----------------------\n------------------------\n";
cout<<"Exact : (suitable for relatively small systems)\n---------------------------------------------\n";
cout<<"Cholesky : Fast but for semipositive definite systems only."<<endl;
cout<<"LU : Works for all definite systems."<<endl<<endl;
cout<<"Approximative : (suitable for large systems)\n----------------------------------------\n";
cout<<"Conjugate_gradient : very fast"<<endl;
cout<<"Gauss_Seidel : Works for diagonally dominant systems."<<endl;
cout<<"Jacobi : Fast but for semipositive definite systems and diagonally dominant system."<<endl;
cout<<"Richardson : Works for all definite systems."<<endl<<endl;
// Choosing of method
string my_method;
cout<<"\nChoose a method by typing a string corresponding to the name of the method, or Exit to abort program\n";
cin>>my_method;
return my_method;
}
//Evaluate if the Matrix A is diagonally dominant
string diagonallyDominant(Matrix A){
string m_diagonallyDominant = "True";
int size = A.NumberOfColumns();
double sum = 0;
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
sum += A(j+1, i+1);
}
sum -= A(i+1, i+1);
if (A(i+1, i+1) < sum) {
m_diagonallyDominant = "False";
}
}
return m_diagonallyDominant;
}
//evaluate the eigenvalues of the matrix A
Vector eigenvalues(Matrix A){
}
void run_algo(Matrix A,Vector b,string my_method){
//string change = "True";
//while (change == "False"){}
if (my_method=="LU") {
LU mySys=LU(A,b);
string command;
cout<<"\nYou chose the LU method. Type a command from the ones available below:\n";
cout<<"U_matrix: outputs the upper triangular decomposition of A.\n";
cout<<"L_matrix: outputs the lower triangular decomposition of A.\n";
cout<<"P_matrix: outputs the permutation matrix of A.\n";
cout<<"Solve: outputs the solution of the system Ax=b\n\n";
cin >> command;
while(command!="Exit") {
if (command == "U_Matrix") { TerminalPrint(mySys.get_U(), "U decomposition of A."); }
else if (command == "L_Matrix") { TerminalPrint(mySys.get_L(), "L decomposition of A."); }
else if (command == "P_Matrix") { TerminalPrint(mySys.get_P(), "Permutation matrix of A."); }
else if (command == "Solve" )
{
Vector x=mySys.Solve(); int counter=0;
for (int i=0;i<b.GetSize();i++){
if (b[i]==0) counter++;
}
if (counter!=b.GetSize()) TerminalPrint(x, "Solution of system Ax=b.");
}
else {cout<<"Command not supported.";}
cout<<"Type Exit to quit the algorithm or type a new command :\n\n";
cin >> command;
}
} else if (my_method=="Cholesky"){
Cholesky mySys=Cholesky(A,b);
string command;
cout<<"\nYou chose the Cholesky method. Type a command from the ones available below:\n";
cout<<"L_matrix: outputs the lower triangular decomposition of A.\n";
cout<<"U_matrix: outputs the upper triangular decomposition of A, which is the transpose of L.\n";
cout<<"P_matrix: outputs the permutation matrix of A.\n";
cout<<"Solve: outputs the solution of the system Ax=b\n\n";
cin >> command;
while(command!="Exit") {
if (command == "U_Matrix") { TerminalPrint(mySys.get_U(), "U decomposition of A."); }
else if (command == "L_Matrix") { TerminalPrint(mySys.get_L(), "L decomposition of A."); }
else if (command == "P_Matrix") { TerminalPrint(mySys.get_P(), "Permutation matrix of A."); }
else if (command == "Solve") { TerminalPrint(mySys.Solve(), "Solution of system Ax=b."); }
else {cout<<"Command not supported.";}
cout<<"Type Exit to quit the algorithm or type a new command :\n\n";
cin >> command;
}
} else if (my_method=="Conjugate_gradient"){
Conjugate_Gradient mySys=Conjugate_Gradient(A,b);
string command;
cout<<"\nYou chose the Conjugate gradient method. Type a command from the ones available below:\n";
cout<<"Set_tol: Sets the convergence criteria to the value given (1e-12 by default). \n";
cout<<"Set_x0 : Sets the initial solution. \n";
cout<<"Set_max_iter : Sets the max number of iterations possible to the value given (5000 by default). \n";
cout<<"Solve: outputs the solution of the system Ax=b. Run this first to unlock following methods:\n\n";
cout<<"iter: outputs the number of iterations needed for convergence. \n";
cout<<"Convergence: outputs a vector with each index corresponding to the error at the iteration. \n";
cin >> command;
while(command!="Exit") {
if (command == "Solve") { TerminalPrint(mySys.Solve(), "Solution of system Ax=b."); }
else if (command == "iter") { TerminalPrint(mySys.get_n_iter(), "Number of iterations needed for convergence."); }
else if (command == "Convergence") {
for (int i=0;i<mySys.get_n_iter();i++){
cout<<mySys.get_e()[i]<<endl;
}
cout<<endl;
}
else if (command == "Set_tol") {
double tol;
cout<<"Type tolerance value you wish to set:\n\n";
cin >> tol;
mySys.set_tol(tol); }
else if (command == "Set_x0") {
Vector x0(b.GetSize());
cout<<"Type x0(1),x0(2),...,x0(n) values you wish to start algorithm with:\n\n";
for (int i=0;i<b.GetSize();i++){
cout<<"Give value of x0("<<i+1<<") : ";
cin>>x0(i+1);
}
TerminalPrint(x0,"Initial solution x0 :"); }
else if (command == "Set_max_iter") {
int max_iter;
cout<<"Type the maximum number of iterations criteria:\n\n";
cin >> max_iter;
mySys.set_iter_max(max_iter); }
else {cout<<"Command not supported.";}
cout<<"Type Exit to quit the algorithm or type a new command :\n\n";
cin >> command;
}
} else if (my_method=="Gauss_Seidel") {
Gauss_Seidel mySys = Gauss_Seidel(A, b);
string command;
string isDD;
isDD = diagonallyDominant(A);
string incorrect = "True";
while (incorrect == "True") {
incorrect = "False";
if (diagonallyDominant == "False") {
cout << "\nYou chose the Gauss_Seidel method.";
cout << "The matrix isn't diagonally dominant therefore the algorithm might not converge.\n";
cout << "You might want to choose another method. Type a command from the ones available below:\n";
cout << "Continue: if you want to procceed \n";
cout << "Change: if you want to change method \n";
cin >> command;
if (command == "Continue") {
cout << "Type a command from the ones available below:\n";
cout<<"Set_tol: Sets the convergence criteria to the value given (1e-12 by default). \n";
cout<<"Set_x0 : Sets the initial solution. \n";
cout<<"Set_max_iter : Sets the max number of iterations possible to the value given (5000 by default). \n";
cout<<"Solve: outputs the solution of the system Ax=b. Run this first to unlock following methods:\n\n";
cout<<"iter: outputs the number of iterations needed for convergence. \n";
cout<<"Convergence: outputs a vector with each index corresponding to the error at the iteration. \n";
cin >> command;
} else if (command == "Change") {
command = "Exit";
} else {
cout << "incorrect command\n";
incorrect = "True";
}
} else {
cout << "\nYou chose the Gauss_Seidel method. Type a command from the ones available below:\n";
cout<<"Set_tol: Sets the convergence criteria to the value given (1e-12 by default). \n";
cout<<"Set_x0 : Sets the initial solution. \n";
cout<<"Set_max_iter : Sets the max number of iterations possible to the value given (5000 by default). \n";
cout<<"Solve: outputs the solution of the system Ax=b. Run this first to unlock following methods:\n\n";
cout<<"iter: outputs the number of iterations needed for convergence. \n";
cout<<"Convergence: outputs a vector with each index corresponding to the error at the iteration. \n";
cin >> command;
}
}
while(command!="Exit") {
if (command == "Solve") { TerminalPrint(mySys.Solve(), "Solution of system Ax=b."); }
else if (command == "iter") { TerminalPrint(mySys.get_n_iter(), "Number of iterations needed for convergence."); }
else if (command == "Convergence") {
for (int i=0;i<mySys.get_n_iter();i++){
cout<<mySys.get_e()[i]<<endl;
}
cout<<endl;
}
else if (command == "Set_tol") {
double tol;
cout<<"Type tolerance value you wish to set:\n\n";
cin >> tol;
mySys.set_tol(tol); }
else if (command == "Set_x0") {
Vector x0(b.GetSize());
cout<<"Type x0(1),x0(2),...,x0(n) values you wish to start algorithm with:\n\n";
for (int i=0;i<b.GetSize();i++){
cout<<"Give value of x0("<<i+1<<") : ";
cin>>x0(i+1);
}
TerminalPrint(x0,"Initial solution x0 :"); }
else if (command == "Set_max_iter") {
int max_iter;
cout<<"Type the maximum number of iterations criteria:\n\n";
cin >> max_iter;
mySys.set_iter_max(max_iter); }
else {cout<<"Command not supported.";}
cout<<"Type Exit to quit the algorithm or type a new command :\n\n";
cin >> command;
}
} else if (my_method=="Jacobi"){
Jacobi mySys=Jacobi(A,b);
string command;
string isDD;
isDD = diagonallyDominant(A);
string incorrect = "True";
while (incorrect == "True") {
incorrect = "False";
if (diagonallyDominant == "False") {
cout << "\nYou chose the Jacobi method.";
cout << "The matrix isn't diagonally dominant therefore the algorithm might not converge.\n";
cout << "You might want to choose another method. Type a command from the ones available below:\n";
cout << "Continue: if you want to procceed \n";
cout << "Change: if you want to change method \n";
cin >> command;
if (command == "Continue") {
cout << "Type a command from the ones available below:\n";
cout << "Set_tol: Sets the convergence criteria to the value given (1e-12 by default). \n";
cout << "Set_x0 : Sets the initial solution. \n";
cout << "Set_max_iter : Sets the max number of iterations possible to the value given (5000 by default). \n";
cout << "Solve: outputs the solution of the system Ax=b. Run this first to unlock following methods:\n\n";
cout << "iter: outputs the number of iterations needed for convergence. \n";
cout << "Convergence: outputs a vector with each index corresponding to the error at the iteration. \n";
cin >> command;
} else if (command == "Change") {
command = "Exit";
} else {
cout << "incorrect command\n";
incorrect = "True";
}
} else {
cout << "\nYou chose the Jacobi method. Type a command from the ones available below:\n";
cout << "Set_tol: Sets the convergence criteria to the value given (1e-12 by default). \n";
cout << "Set_x0 : Sets the initial solution. \n";
cout << "Set_max_iter : Sets the max number of iterations possible to the value given (5000 by default). \n";
cout << "Solve: outputs the solution of the system Ax=b. Run this first to unlock following methods:\n\n";
cout << "iter: outputs the number of iterations needed for convergence. \n";
cout << "Convergence: outputs a vector with each index corresponding to the error at the iteration. \n";
cin >> command;
}
}
while(command!="Exit") {
if (command == "Solve") { TerminalPrint(mySys.Solve(), "Solution of system Ax=b."); }
else if (command == "iter") { TerminalPrint(mySys.get_n_iter(), "Number of iterations needed for convergence."); }
else if (command == "Convergence") {
for (int i=0;i<mySys.get_n_iter();i++){
cout<<mySys.get_e()[i]<<endl;
}
cout<<endl;
}
else if (command == "Set_tol") {
double tol;
cout<<"Type tolerance value you wish to set:\n\n";
cin >> tol;
mySys.set_tol(tol); }
else if (command == "Set_x0") {
Vector x0(b.GetSize());
cout<<"Type x0(1),x0(2),...,x0(n) values you wish to start algorithm with:\n\n";
for (int i=0;i<b.GetSize();i++){
cout<<"Give value of x0("<<i+1<<") : ";
cin>>x0(i+1);
}
TerminalPrint(x0,"Initial solution x0 :"); }
else if (command == "Set_max_iter") {
int max_iter;
cout<<"Type the maximum number of iterations criteria:\n\n";
cin >> max_iter;
mySys.set_iter_max(max_iter); }
else {cout<<"Command not supported.";}
cout<<"Type Exit to quit the algorithm or type a new command :\n\n";
cin >> command;
}
} else if (my_method=="Richardson"){
Richardson mySys=Richardson(A,b);
string command;
cout<<"\nYou chose the Richardson method. Type a command from the ones available below:\n";
cout<<"Set_tol: Sets the convergence criteria to the value given (1e-12 by default). \n";
cout<<"Set_x0 : Sets the initial solution. \n";
cout<<"Set_max_iter : Sets the max number of iterations possible to the value given (5000 by default). \n";
cout<<"Solve: outputs the solution of the system Ax=b. Run this first to unlock following methods:\n\n";
cout<<"iter: outputs the number of iterations needed for convergence. \n";
cout<<"Convergence: outputs a vector with each index corresponding to the error at the iteration. \n";
cin >> command;
while(command!="Exit") {
if (command == "Solve" )
{
Vector x=mySys.Solve(); int counter=0;
for (int i=0;i<b.GetSize();i++){
if (x[i]==0) counter++;
}
if (counter!=b.GetSize()) TerminalPrint(x, "Solution of system Ax=b.");
}
else if (command == "iter") { TerminalPrint(mySys.get_n_iter(), "Number of iterations needed for convergence."); }
else if (command == "Convergence") {
for (int i=0;i<mySys.get_n_iter();i++){
cout<<mySys.get_e()[i]<<endl;
}
cout<<endl;
}
else if (command == "Set_tol") {
double tol;
cout<<"Type tolerance value you wish to set:\n\n";
cin >> tol;
mySys.set_tol(tol); }
else if (command == "Set_x0") {
Vector x0(b.GetSize());
cout<<"Type x0(1),x0(2),...,x0(n) values you wish to start algorithm with:\n\n";
for (int i=0;i<b.GetSize();i++){
cout<<"Give value of x0("<<i+1<<") : ";
cin>>x0(i+1);
}
TerminalPrint(x0,"Initial solution x0 :"); }
else if (command == "Set_max_iter") {
int max_iter;
cout<<"Type the maximum number of iterations criteria:\n\n";
cin >> max_iter;
mySys.set_iter_max(max_iter); }
else {cout<<"Command not supported.";}
cout<<"Type Exit to quit this method and try a new one or type a new command :\n\n";
cin >> command;
}
}else if (my_method=="Exit"){
cout<<"Leaving program .. BYYE\n";
} else {
cout<<"\nMethod not recognized. Select another one.\n\n";
}
}
c4science · Help