Upload files to "ex1D_kahan_summation"

ex1D_kahan_summation/main.cpp

@@ -0,0 +1,33 @@
+#include "mylib.h"
+#include <cmath>
+#include <iostream>
+using namespace std;
+
+int main(int argc, char **argv)
+{
+    for(size_t i = 1; i < 8; ++i)
+    {
+        size_t n = pow(10,i);
+        vector<double> 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;
+}

ex1D_kahan_summation/mylib.cpp

@@ -0,0 +1,34 @@
+#include "mylib.h"
+#include <cassert>       // assert()
+#include <cmath>
+#include <vector>
+using namespace std;
+
+double scalar(vector<double> const &x, vector<double> 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<double> const &x, vector<double> 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;
+}

ex1D_kahan_summation/mylib.h

@@ -0,0 +1,30 @@
+#ifndef FILE_MYLIB
+#define FILE_MYLIB
+#include <vector>
+
+/** 	Inner product
+	@param[in] x	vector
+	@param[in] y	vector
+	@return 	    resulting Euclidian inner product <x,y>
+*/
+double scalar(std::vector<double> const &x, std::vector<double> const &y);
+
+/** 	Inner product using BLAS routines
+	@param[in] x	vector
+	@param[in] y	vector
+	@return 	    resulting Euclidian inner product <x,y>
+*/
+double scalar_cblas(std::vector<double> const &x, std::vector<double> const &y);
+float scalar_cblas(std::vector<float> const &x, std::vector<float> const &y);
+
+
+/** 	L_2 Norm of a vector
+	@param[in] x	vector
+	@return 	    resulting Euclidian norm <x,y>
+*/
+double norm(std::vector<double> const &x);
+
+double Kahan_skalar(std::vector<double> const &x,  std::vector<double> const &y);
+
+
+#endif