From 55b8290f383a1b1aaf07e1976e97655e41614466 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 20 Oct 2011 18:02:35 +0200
Subject: [PATCH] Play around with the Gauss-Seidel method
---
src/Makefile.am | 6 +++++
src/gauss-seidel-sample.cc | 50 ++++++++++++++++++++++++++++++++++++
src/gaussseidel.hh | 52 ++++++++++++++++++++++++++++++++++++++
3 files changed, 108 insertions(+)
create mode 100644 src/gauss-seidel-sample.cc
create mode 100644 src/gaussseidel.hh
diff --git a/src/Makefile.am b/src/Makefile.am
index 4adfdcb1..5d7f6d03 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -5,6 +5,12 @@ check_PROGRAMS = \
test-python \
test-gradient-method
+bin_PROGRAMS = \
+ gauss-seidel-sample
+
+gauss_seidel_sample_SOURCES = \
+ gauss-seidel-sample.cc gaussseidel.hh
+
test_python_SOURCES = \
test-python.cc
diff --git a/src/gauss-seidel-sample.cc b/src/gauss-seidel-sample.cc
new file mode 100644
index 00000000..cd60bc76
--- /dev/null
+++ b/src/gauss-seidel-sample.cc
@@ -0,0 +1,50 @@
+/* -*- mode:c++; mode:semantic -*- */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include "nicefunction.hh"
+#include "mynonlinearity.hh"
+
+#include "gaussseidel.hh"
+
+int main() {
+ size_t const dimension = 3;
+ typedef Dune::FieldMatrix<double, dimension, dimension> LocalMatrix;
+ typedef Dune::FieldVector<double, dimension> LocalVector;
+ size_t const nodes = 10;
+ typedef Dune::FieldMatrix<LocalMatrix, nodes, nodes> GlobalMatrix;
+ typedef Dune::FieldVector<LocalVector, nodes> GlobalVector;
+
+ GlobalMatrix A;
+ // semi-random data
+ for (size_t i = 0; i < nodes; ++i)
+ for (size_t k = 0; k < nodes; ++k)
+ for (size_t j = 0; j < dimension; ++j)
+ for (size_t l = 0; l < dimension; ++l)
+ A[i][k][j][l] = ((i == k) && (j == l) ? 1 : 0);
+
+ GlobalVector b;
+ // semi-random data
+ for (size_t i = 0; i < nodes; ++i)
+ for (size_t j = 0; j < dimension; ++j)
+ b[i][j] = 1;
+
+ Dune::SampleFunction const f = Dune::SampleFunction(); // semi-random function
+ Dune::MyNonlinearity<dimension> const phi(f);
+
+ GlobalVector unu(LocalVector(0));
+ // semi-random data
+ for (size_t m = 0; m < nodes; ++m)
+ for (size_t j = 0; j < dimension; ++j)
+ unu[m][j] = m + j;
+
+ GaussSeidel(unu, A, b, phi);
+ std::cout << std::endl << std::endl;
+ GaussSeidel(unu, A, b, phi);
+ return 0;
+}
diff --git a/src/gaussseidel.hh b/src/gaussseidel.hh
new file mode 100644
index 00000000..46db7f1c
--- /dev/null
+++ b/src/gaussseidel.hh
@@ -0,0 +1,52 @@
+/* -*- mode:c++; mode:semantic -*- */
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include "nicefunction.hh"
+#include "mynonlinearity.hh"
+#include "samplefunctional.hh"
+
+// Just for debugging
+template <int dim, int nodes>
+double J_eval(
+ Dune::FieldMatrix<Dune::FieldMatrix<double, dim, dim>, nodes, nodes> const &
+ A,
+ Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &b,
+ Dune::MyNonlinearity<dim> phi,
+ Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &v) {
+ double ret = 0;
+ for (size_t i = 0; i < nodes; ++i) {
+ for (size_t k = 0; k < nodes; ++k) {
+ Dune::FieldMatrix<double, dim, dim> const localA = A[i][k];
+ Dune::FieldVector<double, dim> tmp;
+ localA.mv(v[k], tmp);
+ ret += 1 / 2 * (tmp * v[i]);
+ }
+ ret -= b[i] * v[i];
+ ret +=
+ phi(v[i]); // an example of phi^i(v[i]) (phi is independent of i here)
+ }
+ return ret;
+}
+
+template <int dim, int nodes>
+void GaussSeidel(
+ Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> &unu,
+ Dune::FieldMatrix<Dune::FieldMatrix<double, dim, dim>, nodes, nodes> const &
+ A,
+ Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &b,
+ Dune::MyNonlinearity<dim> const &phi) {
+ for (size_t m = 0; m < nodes; ++m) {
+ typename Dune::SampleFunctional<dim>::SmallVector const localb = b[m];
+ typename Dune::SampleFunctional<dim>::SmallMatrix const localA = A[m][m];
+ Dune::SampleFunctional<dim> J(localA, localb, phi);
+ typename Dune::SampleFunctional<dim>::SmallVector correction;
+ for (size_t local_iteration = 0; local_iteration < 2;
+ ++local_iteration) { // 2 is random here
+ Dune::minimise(J, unu[m], correction);
+ unu[m] += correction;
+ }
+ std::cout << J_eval(A, b, phi, unu) << std::endl; // debugging
+ }
+}
--
GitLab