From 9098bcc99a1136cab40c0c7406d26852a7270923 Mon Sep 17 00:00:00 2001 From: "jakob.schratter" Date: Wed, 22 Oct 2025 15:22:37 +0200 Subject: [PATCH] Upload files to "ex1D_kahan_summation" --- ex1D_kahan_summation/main.cpp | 33 +++++++++++++++++++++++++++++++++ ex1D_kahan_summation/mylib.cpp | 34 ++++++++++++++++++++++++++++++++++ ex1D_kahan_summation/mylib.h | 30 ++++++++++++++++++++++++++++++ 3 files changed, 97 insertions(+) create mode 100644 ex1D_kahan_summation/main.cpp create mode 100644 ex1D_kahan_summation/mylib.cpp create mode 100644 ex1D_kahan_summation/mylib.h diff --git a/ex1D_kahan_summation/main.cpp b/ex1D_kahan_summation/main.cpp new file mode 100644 index 0000000..8d93e96 --- /dev/null +++ b/ex1D_kahan_summation/main.cpp @@ -0,0 +1,33 @@ +#include "mylib.h" +#include +#include +using namespace std; + +int main(int argc, char **argv) +{ + for(size_t i = 1; i < 8; ++i) + { + size_t n = pow(10,i); + vector x(n); + for (size_t k = 0; k < n; ++k) + x[k] = 1.0/(k + 1); + + + // compute scalar products + double sum_1 = scalar(x, x); + double sum_2 = Kahan_skalar(x, x); + + // compute error + double err_1 = abs(sum_1 - pow(M_PI,2)/6); + double err_2 = abs(sum_2 - pow(M_PI,2)/6); + + cout << "n = " << n << endl; + cout << "Normal scalar product: " << sum_1 << "\terror: " << err_1 << endl; + cout << "Kahan scalar product: " << sum_2 << "\terror: " << err_2 << endl; + cout << endl; + } + + + + return 0; +} diff --git a/ex1D_kahan_summation/mylib.cpp b/ex1D_kahan_summation/mylib.cpp new file mode 100644 index 0000000..11fc394 --- /dev/null +++ b/ex1D_kahan_summation/mylib.cpp @@ -0,0 +1,34 @@ +#include "mylib.h" +#include // assert() +#include +#include +using namespace std; + +double scalar(vector const &x, vector const &y) +{ + assert(x.size() == y.size()); + size_t const N = x.size(); + double sum = 0.0; + for (size_t i = 0; i < N; ++i) + { + sum += x[i] * y[i]; + } + return sum; +} + + + +double Kahan_skalar(vector const &x, vector const &y) +{ + double sum = 0; + double c = 0; + size_t n = x.size(); + for (size_t i = 0; i < n; ++i) + { + double z = x[i]*y[i] - c; // c is the part that got lost in the last iteration + double t = sum + z; // when adding sum + z, the lower digits are lost if sum is large + c = (t - sum) - z; // now we recover the lower digits to add in the next iteration + sum = t; + } + return sum; +} \ No newline at end of file diff --git a/ex1D_kahan_summation/mylib.h b/ex1D_kahan_summation/mylib.h new file mode 100644 index 0000000..7914250 --- /dev/null +++ b/ex1D_kahan_summation/mylib.h @@ -0,0 +1,30 @@ +#ifndef FILE_MYLIB +#define FILE_MYLIB +#include + +/** Inner product + @param[in] x vector + @param[in] y vector + @return resulting Euclidian inner product +*/ +double scalar(std::vector const &x, std::vector const &y); + +/** Inner product using BLAS routines + @param[in] x vector + @param[in] y vector + @return resulting Euclidian inner product +*/ +double scalar_cblas(std::vector const &x, std::vector const &y); +float scalar_cblas(std::vector const &x, std::vector const &y); + + +/** L_2 Norm of a vector + @param[in] x vector + @return resulting Euclidian norm +*/ +double norm(std::vector const &x); + +double Kahan_skalar(std::vector const &x, std::vector const &y); + + +#endif