ignore dat files, add num6

2020-06-11 23:22:05 +02:00
parent b522ade33c
commit d9f0815f84
8 changed files with 446 additions and 8 deletions
@@ -0,0 +1,152 @@
+#include <iostream>    // notwendig zur Ausgabe
+#include <vector>
+#include "hdnum.hh"    // hdnum header
+
+
+namespace hdnum {
+
+    template<typename REAL>
+    class SparseMatrix {
+        struct MatrixEntry {
+            int i;
+            int j;
+            REAL value;
+        };
+
+        public:
+            void AddEntry (int i, int j, REAL value) {
+                if (i < 0 || j < 0)
+                    HDNUM_ERROR("Zeilen und Spalten dürfen nicht negativ sein!");
+                // create new matrix entry
+                MatrixEntry entry;
+                entry.i = i;
+                entry.j = j;
+                entry.value = value;
+                // add entry to entries vector
+                entries.insert(entries.begin(), entry);
+            }
+
+            template<typename V>
+            void mv (Vector<V>& y, const Vector<V>& x) {
+                // iterate through all entries
+                for(typename std::vector<MatrixEntry>::iterator it = entries.begin(); it != entries.end(); ++it) {
+                    if ((*it).i > y.size() || (*it).j > x.size())
+                        HDNUM_ERROR("Vektor(en) hat/haben falsche Dimension!");
+                    // modify corresponding vector entry
+                    y[(*it).i] += x[(*it).j] * (*it).value;
+                }
+            }
+
+        private:
+            // vector storing non-zero elements
+            std::vector<MatrixEntry> entries;
+    };
+}
+
+// zum Test, fluss matrizen aus rohrleitungsnetzwerk
+hdnum::SparseMatrix<double> flussMatrixSparse(int N) {
+    hdnum::SparseMatrix<double> A;
+
+    // berechnung der kantenlaenge
+    int n = floor(sqrt(N+1));
+
+    for(int i = 0; i < N; i++) {
+        int edges = 0;
+        if ((i+1)%n != 0) { // nicht linker rand
+            edges++;
+            if (i-1 >= 0) A.AddEntry(i, i-1, -1);
+        }
+        if ((i+2)%n != 0) { // nicht rechter rand
+            edges++;
+            A.AddEntry(i, i+1, -1);
+        }
+        if (i+n < N) { // nicht unterer rand
+            edges++;
+            A.AddEntry(i, i+n, -1);
+        }
+        if (i-n >= -1) { // nicht oberer rand
+            edges++;
+            if (i-n >= 0) A.AddEntry(i, i-n, -1);
+        }
+        A.AddEntry(i, i, edges);
+    }
+    return A;
+}
+
+// Funktion zum Aufstellen der Matrix (aus rohrleitungsnetzwerk)
+hdnum::DenseMatrix<double> flussMatrix(int N) {
+    hdnum::DenseMatrix<double> A(N,N);
+
+    // berechnung der kantenlaenge
+    int n = floor(sqrt(N+1));
+
+    for(int i = 0; i < N; i++) {
+        int edges = 0;
+        if ((i+1)%n != 0) { // nicht linker rand
+            edges++;
+            if (i-1 >= 0) A(i, i-1) = -1; // falls nicht der referenzknoten
+        }
+        if ((i+2)%n != 0) { // nicht rechter rand
+            edges++;
+            A(i, i+1) = -1;
+        }
+        if (i+n < N) { // nicht unterer rand
+            edges++;
+            A(i, i+n) = -1;
+        }
+        if (i-n >= -1) { // nicht oberer rand
+            edges++;
+            if (i-n >= 0) A(i, i-n) = -1; // falls nicht der referenzknoten
+        }
+        A(i, i) = edges;
+    }
+    return A;
+}
+
+// funktion zum testen
+void testMatrixVectorProduct(int N) {
+    hdnum::SparseMatrix<double> A = flussMatrixSparse(N);
+    hdnum::DenseMatrix<double> B = flussMatrix(N);
+
+    // test vektor
+    hdnum::Vector<double> x(N);
+    for(int i=0; i<N; i++) {
+        x[i] = i+1;
+    }
+    hdnum::Vector<double> y1(N);
+    hdnum::Vector<double> y2(N);
+
+    // Pretty-printing fuer Vektoren
+    x.scientific(false);
+    x.width(15);
+
+    // test ob Ax = Bx
+    A.mv(y1, x);
+    B.mv(y2, x);
+    std::cout << "Mit Sparse Matrix: " << y1 << std::endl;
+    std::cout << "Mit Dense Matrix: " << y2;
+}
+
+int main () {
+    // testMatrixVectorProduct(15);
+
+    // messe zeiten fuer n=4,...,14 und N = 2^n
+    for(int n=4; n<=14; n++) {
+        int N = pow(2, n);
+
+        // test vektor
+        hdnum::Vector<double> x(N);
+        for(int i=0; i<N; i++) {
+            x[i] = i+1;
+        }
+        hdnum::Vector<double> y(N);
+
+        hdnum::SparseMatrix<double> A = flussMatrixSparse(N);
+        hdnum::DenseMatrix<double> B = flussMatrix(N);
+
+        // timer hier fuer die SparseMatrix implementierung
+        hdnum::Timer myTimer;
+        A.mv(y, x);
+        std::cout << N << " " << myTimer.elapsed() << std::endl;
+    }
+}