diff --git a/src/Makefile b/src/Makefile index 1e668e6..29d1bfa 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,17 +1,19 @@ CC=g++ +OBJS= generalAlgorithm.o specialAlgorithm.o + .PHONY: clean -all: simpleFile analyticPlot +all: main analyticPlot -simpleFile: simpleFile.o +main: main.o $(OBJS) $(CC) -o $@ $^ analyticPlot: analyticPlot.o $(CC) -o $@ $^ %.o: %.cpp - $(CC) -c $< -o $@ + $(CC) -c -o $@ $^ clean: rm *.o diff --git a/src/generalAlgorithm.cpp b/src/generalAlgorithm.cpp new file mode 100644 index 0000000..4ff167d --- /dev/null +++ b/src/generalAlgorithm.cpp @@ -0,0 +1,61 @@ +#include "generalAlgorithm.hpp" +#include +#include +#include +#include +#include + +arma::vec* general_algorithm( + arma::vec* sub_diag, + arma::vec* main_diag, + arma::vec* sup_diag, + arma::vec* g_vec +) +{ + int n = g_vec->n_elem; + double d; + + for (int i = 1; i < n; i++) { + d = (*sub_diag)(i-1) / (*main_diag)(i-1); + (*main_diag)(i) -= d*(*sup_diag)(i-1); + (*g_vec)(i) -= d*(*g_vec)(i-1); + } + + (*g_vec)(n-1) /= (*main_diag)(n-1); + + for (int i = n-2; i >= 0; i--) { + (*g_vec)(i) = ((*g_vec)(i) - (*sup_diag)(i) * (*g_vec)(i+1)) / (*main_diag)(i); + } + return g_vec; +} + +void general_algorithm_main() +{ + arma::vec sub_diag, main_diag, sup_diag, g_vec, *v_vec; + std::ofstream ofile; + int steps; + double step_size; + + for (int i = 0; i < 6; i++) { + steps = std::pow(10, i+1); + step_size = 1./(double) steps; + main_diag.resize(steps - 1); + g_vec.resize(steps - 1); + sub_diag.resize(steps - 2); + sup_diag.resize(steps - 2); + + sub_diag.fill(-1.); + sup_diag.fill(-1.); + main_diag.fill(2.); + + v_vec = general_algorithm(&sub_diag, &main_diag, &sup_diag, &g_vec); + + ofile.open("general_algorithm_" + std::to_string(steps) + ".txt"); + + for (int j=0; j < v_vec->n_elem; j++) { + ofile << std::setprecision(4) << std::scientific << step_size*(i+1) << "," + << std::setprecision(4) << std::scientific << (*v_vec)(i) << std::endl; + } + ofile.close(); + } +} diff --git a/src/generalAlgorithm.hpp b/src/generalAlgorithm.hpp new file mode 100644 index 0000000..cf9fd35 --- /dev/null +++ b/src/generalAlgorithm.hpp @@ -0,0 +1,15 @@ +#ifndef __GENERAL_ALG__ +#define __GENERAL_ALG__ + +#include + +arma::vec* general_algorithm( + arma::vec* sub_diag, + arma::vec* main_diag, + arma::vec* sup_diag, + arma::vec* g_vec +); + +void general_algorithm_main(); + +#endif diff --git a/src/main.cpp b/src/main.cpp index 65bbf4a..27fbeeb 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,24 +1,124 @@ -#include "GeneralAlgorithm.hpp" #include -#include +#include +#include +#include +#include +#include +#include + +#include "generalAlgorithm.hpp" +#include "specialAlgorithm.hpp" + +#define TIMING_ITERATIONS 5 + +void error( + std::string filename, + arma::vec* x_vec, + arma::vec* v_vec, + arma::vec* a_vec +) +{ + std::ofstream ofile; + ofile.open(filename); + + if (!ofile.is_open()) { + exit(1); + } + + for (int i=0; i < a_vec->n_elem; i++) { + double sub = (*a_vec)(i) - (*v_vec)(i); + ofile << std::setprecision(8) << std::scientific << (*x_vec)(i) + << std::setprecision(8) << std::scientific << std::log10(std::abs(sub)) + << std::setprecision(8) << std::scientific << std::log10(std::abs(sub/(*a_vec)(i))) + << std::endl; + } + + ofile.close(); +} double f(double x) { - return 100. * std::exp(-10.*x); + return 100*std::exp(-10*x); } -double a_sol(double x) { - return 1. - (1. - std::exp(-10)) * x - std::exp(-10*x); +void build_g_vec(int n_steps, arma::vec* g_vec) { + g_vec->resize(n_steps-1); + + double step_size = 1./ (double) n_steps; + for (int i=0; i < n_steps-1; i++) { + (*g_vec)(i) = f((i+1)*step_size); + } } -int main() { - arma::mat A = arma::eye(3,3); +void build_array( + int n_steps, + arma::vec* sub_diag, + arma::vec* main_diag, + arma::vec* sup_diag, + arma::vec* g_vec +) +{ + sub_diag->resize(n_steps-2); + main_diag->resize(n_steps-1); + sup_diag->resize(n_steps-2); - GeneralAlgorithm ga(3, &A, f, a_sol, 0., 1.); + sub_diag->fill(-1); + main_diag->fill(2); + sup_diag->fill(-1); - ga.solve(); - std::cout << "Time: " << ga.time(5) << std::endl; - ga.error(); + build_g_vec(n_steps, g_vec); +} - return 0; + +void timing() { + arma::vec sub_diag, main_diag, sup_diag, g_vec; + int n_steps; + + std::ofstream ofile; + ofile.open("timing.txt"); + + // Timing + for (int i=1; i <= 8; i++) { + n_steps = std::pow(10, i); + clock_t g_1, g_2, s_1, s_2; + double g_res = 0, s_res = 0; + + for (int j=0; j < TIMING_ITERATIONS; j++) { + build_array(n_steps, &sub_diag, &main_diag, &sup_diag, &g_vec); + + g_1 = clock(); + + general_algorithm(&sub_diag, &main_diag, &sup_diag, &g_vec); + + g_2 = clock(); + + g_res += (double) (g_2 - g_1) / CLOCKS_PER_SEC; + build_g_vec(n_steps, &g_vec); + + s_1 = clock(); + + special_algorithm(-1., 2., -1., &g_vec); + + s_2 = clock(); + + s_res += (double) (s_2 - s_1) / CLOCKS_PER_SEC; + + } + ofile + << n_steps << "," + << g_res / (double) TIMING_ITERATIONS << "," + << s_res / (double) TIMING_ITERATIONS << std::endl; + } + + ofile.close(); +} + +void error_file() { } + +int main() +{ + timing(); + general_algorithm_main(); + special_algorithm_main(); +} diff --git a/src/simpleFile.cpp b/src/simpleFile.cpp deleted file mode 100644 index 344e7a6..0000000 --- a/src/simpleFile.cpp +++ /dev/null @@ -1,162 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#define TIMING_ITERATIONS 5 - -arma::vec* general_algorithm( - arma::vec* sub_diag, - arma::vec* main_diag, - arma::vec* sup_diag, - arma::vec* g_vec -) -{ - int n = g_vec->n_elem; - double d; - - for (int i = 1; i < n; i++) { - d = (*sub_diag)(i-1) / (*main_diag)(i-1); - (*main_diag)(i) -= d*(*sup_diag)(i-1); - (*g_vec)(i) -= d*(*g_vec)(i-1); - } - - (*g_vec)(n-1) /= (*main_diag)(n-1); - - for (int i = n-2; i >= 0; i--) { - (*g_vec)(i) = ((*g_vec)(i) - (*sup_diag)(i) * (*g_vec)(i+1)) / (*main_diag)(i); - } - return g_vec; -} - - -arma::vec* special_algorithm( - double sub_sig, - double main_sig, - double sup_sig, - arma::vec* g_vec -) -{ - int n = g_vec->n_elem; - arma::vec diag = arma::vec(n); - - for (int i = 1; i < n; i++) { - // Calculate values for main diagonal based on indices - diag(i-1) = (double)(i+1) / i; - (*g_vec)(i) += (*g_vec)(i-1) / diag(i-1); - } - // The last element in main diagonal has value (i+1)/i = (n+1)/n - (*g_vec)(n-1) /= (double)(n+1) / (n); - - for (int i = n-2; i >= 0; i--) { - (*g_vec)(i) = ((*g_vec)(i) + (*g_vec)(i+1))/ diag(i); - } - - return g_vec; -} - -void error( - std::string filename, - arma::vec* x_vec, - arma::vec* v_vec, - arma::vec* a_vec -) -{ - std::ofstream ofile; - ofile.open(filename); - - if (!ofile.is_open()) { - exit(1); - } - - for (int i=0; i < a_vec->n_elem; i++) { - double sub = (*a_vec)(i) - (*v_vec)(i); - ofile << std::setprecision(8) << std::scientific << (*x_vec)(i) - << std::setprecision(8) << std::scientific << std::log10(std::abs(sub)) - << std::setprecision(8) << std::scientific << std::log10(std::abs(sub/(*a_vec)(i))) - << std::endl; - } - - ofile.close(); - -} - -double f(double x) { - return 100*std::exp(-10*x); -} - -void build_array( - int n_steps, - arma::vec* sub_diag, - arma::vec* main_diag, - arma::vec* sup_diag, - arma::vec* g_vec -) -{ - sub_diag->resize(n_steps-2); - main_diag->resize(n_steps-1); - sup_diag->resize(n_steps-2); - - sub_diag->fill(-1); - main_diag->fill(2); - sup_diag->fill(-1); - - g_vec->resize(n_steps-1); - - double step_size = 1./ (double) n_steps; - for (int i=0; i < n_steps-1; i++) { - (*g_vec)(i) = f((i+1)*step_size); - } - -} - -void timing() { - arma::vec sub_diag, main_diag, sup_diag, g_vec; - int n_steps; - - std::ofstream ofile; - ofile.open("timing.txt"); - - // Timing - for (int i=1; i <= 8; i++) { - n_steps = std::pow(10, i); - clock_t g_1, g_2, s_1, s_2; - double g_res = 0, s_res = 0; - - for (int j=0; j < TIMING_ITERATIONS; j++) { - build_array(n_steps, &sub_diag, &main_diag, &sup_diag, &g_vec); - - g_1 = clock(); - - general_algorithm(&sub_diag, &main_diag, &sup_diag, &g_vec); - - g_2 = clock(); - - g_res += (double) (g_2 - g_1) / CLOCKS_PER_SEC; - build_array(n_steps, &sub_diag, &main_diag, &sup_diag, &g_vec); - - s_1 = clock(); - - special_algorithm(-1., 2., -1., &g_vec); - - s_2 = clock(); - - s_res += (double) (s_2 - s_1) / CLOCKS_PER_SEC; - - } - ofile - << n_steps << "," - << g_res / (double) TIMING_ITERATIONS << "," - << s_res / (double) TIMING_ITERATIONS << std::endl; - } - - ofile.close(); -} - -int main() -{ - timing(); -} diff --git a/src/specialAlgorithm.cpp b/src/specialAlgorithm.cpp new file mode 100644 index 0000000..ad88699 --- /dev/null +++ b/src/specialAlgorithm.cpp @@ -0,0 +1,52 @@ +#include "specialAlgorithm.hpp" +#include + +arma::vec* special_algorithm( + double sub_sig, + double main_sig, + double sup_sig, + arma::vec* g_vec +) +{ + int n = g_vec->n_elem; + arma::vec diag = arma::vec(n); + + for (int i = 1; i < n; i++) { + // Calculate values for main diagonal based on indices + diag(i-1) = (double)(i+1) / i; + (*g_vec)(i) += (*g_vec)(i-1) / diag(i-1); + } + // The last element in main diagonal has value (i+1)/i = (n+1)/n + (*g_vec)(n-1) /= (double)(n+1) / (n); + + for (int i = n-2; i >= 0; i--) { + (*g_vec)(i) = ((*g_vec)(i) + (*g_vec)(i+1))/ diag(i); + } + + return g_vec; +} + +void special_algorithm_main() +{ + arma::vec g_vec, *v_vec; + std::ofstream ofile; + int steps; + double sub_sig, main_sig, sup_sig, step_size; + + for (int i = 0; i < 6; i++) { + steps = std::pow(10, i+1); + step_size = 1./(double) steps; + g_vec.resize(steps - 1); + + v_vec = special_algorithm(sub_sig, main_sig, sup_sig, &g_vec); + + ofile.open("special_algorithm_" + std::to_string(steps) + ".txt"); + + for (int j=0; j < v_vec->n_elem; j++) { + ofile << std::setprecision(4) << std::scientific << step_size*(i+1) << "," + << std::setprecision(4) << std::scientific << (*v_vec)(i) << std::endl; + } + ofile.close(); + } +} + diff --git a/src/specialAlgorithm.hpp b/src/specialAlgorithm.hpp new file mode 100644 index 0000000..1578e64 --- /dev/null +++ b/src/specialAlgorithm.hpp @@ -0,0 +1,15 @@ +#ifndef __SPECIAL_ALG__ +#define __SPECIAL_ALG__ + +#include + +arma::vec* special_algorithm( + double sub_sig, + double main_sig, + double sup_sig, + arma::vec* g_vec +); + +void special_algorithm_main(); + +#endif