ex 8

BSP_3_8_jacobi_seq/Makefile

@@ -0,0 +1,54 @@
+#
+# use GNU-Compiler tools
+COMPILER=GCC_
+# COMPILER=GCC_SEQ_
+# alternatively from the shell
+# export COMPILER=GCC_
+# or, alternatively from the shell
+# make COMPILER=GCC_
+
+MAIN = main
+SOURCES = ${MAIN}.cpp vdop.cpp geom.cpp\
+	getmatrix.cpp jacsolve.cpp userset.cpp
+# 	dexx.cpp debugd.cpp skalar.cpp  vecaccu.cpp accudiag.cpp
+
+OBJECTS = $(SOURCES:.cpp=.o)
+
+PROGRAM	= ${MAIN}.${COMPILER}
+
+# uncomment the next to lines for debugging and detailed performance analysis
+CXXFLAGS += -g
+# -pg slows down the code on my laptop when using CLANG_
+#LINKFLAGS += -pg
+#CXXFLAGS += -Q --help=optimizers
+#CXXFLAGS += -fopt-info
+
+include ../${COMPILER}default.mk
+
+#############################################################################
+# additional specific cleaning in this directory
+clean_all::
+	@rm -f t.dat*
+
+
+#############################################################################
+# special testing
+# NPROCS	= 4
+#
+TFILE	= t.dat
+# TTMP	= t.tmp
+#
+graph: $(PROGRAM)
+# 	@rm -f $(TFILE).*
+	# next two lines only sequentially
+	./$(PROGRAM)
+	@mv  $(TFILE).000 $(TFILE)
+# 	$(MPIRUN) $(MPIFLAGS) -np $(NPROCS) $(PROGRAM)
+# 	@echo " "; echo "Manipulate data for graphics."; echo " "
+# 	@cat $(TFILE).* > $(TTMP)
+# 	@sort -b -k 2    $(TTMP)   -o $(TTMP).1
+# 	@sort -b -k 1    $(TTMP).1 -o $(TTMP).2
+# 	@awk  -f nl.awk  $(TTMP).2  > $(TFILE)
+# 	@rm -f $(TTMP).* $(TTMP) $(TFILE).*
+#
+	-gnuplot jac.dem

BSP_3_8_jacobi_seq/jacsolve.h

@@ -0,0 +1,18 @@
+#ifndef JACSOLVE_FILE
+#define JACSOLVE_FILE
+#include "getmatrix.h"
+#include <vector>
+
+
+/**
+ * Solves linear system of equations  K @p u = @p f  via the Jacobi iteration.
+ * We use a distributed symmetric  CSR matrix @p SK and initial guess of the
+ * solution is set to 0.
+ * @param[in] SK	CSR matrix
+ * @param[in] f		distributed local vector storing the right hand side
+ * @param[out] u	accumulated local vector storing the solution.
+*/
+void JacobiSolve(CRS_Matrix const &SK, std::vector<double> const &f, std::vector<double> &u);
+
+
+#endif

BSP_3_8_jacobi_seq/main.cpp

@@ -0,0 +1,129 @@
+//		MPI code in C++.
+//		See [Gropp/Lusk/Skjellum, "Using MPI", p.33/41 etc.]
+//		and  /opt/mpich/include/mpi2c++/comm.h  for details
+
+#include "geom.h"
+#include "getmatrix.h"
+#include "jacsolve.h"
+#include "userset.h"
+#include "vdop.h"
+
+#include <chrono>             // timing
+#include <cmath>
+#include <iostream>
+using namespace std;
+using namespace std::chrono;  // timing
+
+
+int main(int, char ** )
+{
+    const int numprocs = 1;
+    const int myrank   = 0;
+
+    if (myrank == 0)
+    {
+        cout << "\n There are " << numprocs << " processes running.\n \n";
+    }
+
+    const auto procx = static_cast<int>(sqrt(numprocs + 0.0));
+    const int  procy = procx;
+
+    if (procy * procx != numprocs)
+    {
+        cout << "\n Wrong number of processors !\n \n";
+    }
+    else
+    {
+// #####################################################################
+//      Here starts the real code
+// #####################################################################
+        //bool ScaleUp = !true;
+        int nx, ny, NXglob, NYglob; /* number of local intervals on (xl,xr)=:nx, (yb,yt)=:ny */
+        //nx = 1024;
+        //ny = 1024;
+        nx = 100;
+        ny = 100;
+        NXglob = nx * procx;
+        NYglob = ny * procy;
+        cout << "Intervalls: " << NXglob << " x " << NYglob << endl;
+
+// ##################### STL ###########################################
+        {
+            Mesh_2d_3_square const mesh(nx, ny);
+            //mesh.Debug();
+
+            CRS_Matrix SK(mesh);                   // CRS matrix
+            //SK.Debug();
+            vector<double> uv(SK.Nrows(), 0.0);    // temperature
+            vector<double> fv(SK.Nrows(), 0.0);    // r.h.s.
+
+            SK.CalculateLaplace(fv);
+            //SK.Debug();
+
+            //mesh.SetU(uv);         // deprecated
+            //mesh.SetF(fv);         // deprecated
+            // Two ways to initialize the vector
+            //mesh.SetValues(uv,f_zero);             // functional
+            mesh.SetValues(uv, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
+
+            SK.ApplyDirichletBC(uv, fv);
+            //SK.Compare2Old(nnode, id, ik, sk);
+            //SK.Debug();
+
+            auto tstart = system_clock::now();       // start timer
+
+            JacobiSolve(SK, fv, uv );          // solve the system of equations
+
+            auto tend = system_clock::now();         // end timer
+            auto duration = duration_cast<microseconds>(tend - tstart);
+            auto t1 = static_cast<double>(duration.count()) / 1e6 ;   // t1 in seconds
+            cout << "JacobiSolve: timing in sec. : " << t1 << endl;
+
+            //CompareVectors(uv, nnode, u, 1e-6);    // Check correctness
+
+            //mesh.SaveVectorP("t.dat", uv);
+            //mesh.Visualize(uv);
+
+        }
+// ##################### STL ###########################################
+        {
+            //Mesh_2d_3_matlab const mesh("square_tiny.txt");
+            Mesh_2d_3_matlab const mesh("square_100.txt");
+            //Mesh_2d_3_matlab const mesh("L_shape.txt");
+            //mesh.Debug();
+
+            CRS_Matrix SK(mesh);                   // CRS matrix
+            //SK.Debug();
+
+            vector<double> uv(SK.Nrows(), 0.0);    // temperature
+            vector<double> fv(SK.Nrows(), 0.0);    // r.h.s.
+
+            SK.CalculateLaplace(fv);
+            //SK.Debug();
+
+            //mesh.SetU(uv);         // deprecated
+            // Two ways to initialize the vector
+            //mesh.SetValues(uv,f_zero);             // user function
+            mesh.SetValues(uv, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
+
+            SK.ApplyDirichletBC(uv, fv);
+            //SK.Compare2Old(nnode, id, ik, sk);
+            //SK.Debug();
+
+            auto tstart = system_clock::now();       // start timer
+
+            JacobiSolve(SK, fv, uv );          // solve the system of equations
+
+            auto tend = system_clock::now();         // end timer
+            auto duration = duration_cast<microseconds>(tend - tstart);
+            auto t1 = static_cast<double>(duration.count()) / 1e6 ;   // t1 in seconds
+            cout << "JacobiSolve: timing in sec. : " << t1 << endl;
+
+            //mesh.Write_ascii_matlab("uv.txt", uv);
+            //mesh.Visualize(uv);
+        }
+        return 0;
+    }
+}
+
+