Merge remote-tracking branch 'origin/master'

8 files changed, 369 insertions, 36 deletions

diff --git a/Makefile b/Makefile
index c2007b9..15fbc38 100644
--- a/Makefile
+++ b/Makefile
@@ -7,7 +7,7 @@ PRG=main
 CC=cc
 CFLAGS=-Wall -Wextra -pedantic
 LFLAGS=-fopenmp
-SRC=main.c linalg.c solver.c sysmat.c
+SRC=test_bic.c linalg.c solver.c sysmat.c
 RUNFLAGS=OMP_NUM_THREADS=4
 INCLUDES=
 LIB=-lm
diff --git a/main.c b/main.c
index 75f6446..275fea3 100644
--- a/main.c
+++ b/main.c
@@ -32,6 +32,17 @@ int main (void)
 	double yD = -1.0; // lower boundary
 	double yU =  1.0; // upper boundary
 
+	// boundary type
+	int north = 0;
+	int south = 1;
+	int east  = 2;
+	int west  = 2;
+
+	double val_north = 40.0;
+	double val_south = 1.8e-5;
+	double val_east  = 0.0;
+	double val_west  = 0.0;
+
 	double dx = (xR - xL) / nx;
 	double dy = (yU - yD) / ny;
 
@@ -46,7 +57,6 @@ int main (void)
 	double* a  = (double*) calloc(5 * dim, sizeof(double));
 	double* c  = (double*) calloc(dim, sizeof(double));
 	double* r  = (double*) calloc(dim, sizeof(double));
-	double* ar = (double*) calloc(dim, sizeof(double));
 	double* x  = (double*) calloc(dim, sizeof(double));
 	double* vx  = (double*) calloc(dim, sizeof(double));
 	double* vy  = (double*) calloc(dim, sizeof(double));
@@ -78,39 +88,8 @@ int main (void)
 
 	// initialize boundary conditions
 	//
-
-
-	// dirichlet boundary
-	for (int i = 0; i < nx; i ++) {
-		x[i] = 40.0;
-	}
-
-	multiply(&m, x, ar);
-
-	for (int i = 0; i < dim; i ++) {
-		r[i] = r[i] - ar[i];
-		if (i < nx) {
-			r[i] = x[i];
-		}
-	}
-
-	free(ar);
-
-	for (int i = 0; i < nx; i ++) {
-		for (int j = 0; j < dim; j ++) {
-			if (i != j) {
-				put(&m, j, i, 0.0);
-				put(&m, i, j, 0.0);
-			} else {
-				put(&m, j, i, 1.0);
-			}
-		}
-	}
-
-	// neumann boundary
-	for (int i = dim - nx; i < dim; i ++) {
-		r[i] = 1.8e-5;
-	}
+	boundary(dim, nx, ny, north, south, east, west, val_north, val_south, 
+		val_east, val_west, &m, r, x);
 
 	// solve the system A*x = b
 	//
diff --git a/result/log b/result/log
deleted file mode 100644
index e69de29..0000000
--- a/result/log
+++ /dev/null
diff --git a/solver.c b/solver.c
index 4427217..808ab89 100644
--- a/solver.c
+++ b/solver.c
@@ -35,6 +35,207 @@
 
 // ----------------------------------------------------------------------------
 
+int bic(struct diag* a, double* b, double* x0, double tol)
+{
+
+	int len_c = a->len_c;
+	int count = 1;
+	int flag  = 1;
+
+	double* r    = (double*) calloc(len_c, sizeof(double));
+	double* rhat = (double*) calloc(len_c, sizeof(double));
+	double* h    = (double*) calloc(len_c, sizeof(double));
+	double* p    = (double*) calloc(len_c, sizeof(double));
+	double* v    = (double*) calloc(len_c, sizeof(double));
+	double* s    = (double*) calloc(len_c, sizeof(double));
+	double* t    = (double*) calloc(len_c, sizeof(double));
+
+	double rho   = 1.0;
+	double alp   = 1.0;
+	double ome   = 1.0;
+
+	double rh1   = 0.0;
+	double bet   = 0.0;
+	double eps   = 0.0;
+
+	multiply(a, x0, r);
+
+	for (int i = 0; i < len_c; i ++) {
+		r[i] = b[i] - r[i];
+	}
+
+	memcpy(rhat, r, len_c * sizeof(double));
+
+	memset(v, 0.0, len_c * sizeof(double));
+	memset(p, 0.0, len_c * sizeof(double));
+
+	for (int i = 1; i < 1000; i ++) {
+
+		eps = norm(r, len_c);
+		if (eps < tol) {
+			flag = 0;
+			break;
+		}
+
+		rh1 = dot(rhat, r, len_c);
+		bet = (rh1 / rho) * (alp / ome);
+
+		for (int k = 0; k < len_c; k ++) {
+			p[k] = r[k] + bet * (p[k] - ome * v[k]);
+		}
+
+		multiply(a, p, v);
+		alp = rh1 / dot(rhat, v, len_c);
+
+		for (int k = 0; k < len_c; k ++) {
+			h[k] = x0[k] + alp * p[k];
+			s[k] = r[k] - alp * v[k];
+		}
+
+		multiply(a, s, t);
+
+		ome = dot(t, s, len_c) / dot(t, t, len_c);
+
+		for (int k = 0; k < len_c; k ++) {
+			x0[k] = h[k] + ome * s[k];
+			r[k]  = s[k] - ome * t[k];
+		}
+
+		rho = rh1;
+
+		count = count + 1;
+
+		memset(s, 0.0, len_c * sizeof(double));
+		memset(t, 0.0, len_c * sizeof(double));
+		memset(h, 0.0, len_c * sizeof(double));
+
+	}
+
+	printf(">>> (bic) :: iteration complete w/ err : %g, \
+		flag: %d, steps: %d\n", eps, flag, count);
+
+	free(r);
+	free(rhat);
+	free(h);
+	free(p);
+	free(v);
+	free(s);
+	free(t);
+
+	return flag;
+
+}
+
+
+
+/*
+ * preconditioned biconjugate gradient solver
+ *
+ * c    : preconditioner matrix
+ * inout: x0  (in: initial guess, out: final result)
+ * inout: tol (in: tolerance crit., out: achieved accuracy)
+ */
+int bicgstab(struct diag* a, double* b, double* x0, struct diag* c, double tol)
+{
+
+	// params
+	int len_c = a->len_c;
+	int count = 1;
+	int flag  = 1;
+
+	double* r  = (double*) calloc(len_c, sizeof(double));
+	double* rh = (double*) calloc(len_c, sizeof(double));
+	double* p  = (double*) calloc(len_c, sizeof(double));
+	double* ph = (double*) calloc(len_c, sizeof(double));
+	double* v  = (double*) calloc(len_c, sizeof(double));
+	double* s  = (double*) calloc(len_c, sizeof(double));
+	double* sh = (double*) calloc(len_c, sizeof(double));
+	double* t  = (double*) calloc(len_c, sizeof(double));
+
+	double rho1   = 0.0;
+	double rho2   = 0.0;
+	double beta   = 0.0;
+	double alph   = 0.0;
+	double omeg   = 0.0;
+
+	double eps = tol;
+
+	// bootstrap
+	multiply(a, x0, r);
+
+	for (int i = 0; i < len_c; i ++) {
+		r[i] = b[i] - r[i];
+	}
+
+	memcpy(rh, r, len_c * sizeof(double));
+
+	// loop
+	while (count < 1000) {
+
+		rho1 = dot(rh, r, len_c);
+		if (rho1 < 1.0e-18) {
+			flag = 2;
+			break;
+		}
+
+		if (count == 1) {
+			memcpy(p, r, len_c * sizeof(double));
+		} else {
+			beta = (rho1 / rho2) * (alph / omeg);
+			for (int i = 0; i < len_c; i ++) {
+				p[i] = r[i] + beta * (p[i] - omeg * v[i]);
+			}
+		}
+
+		multiply(c, p, ph);
+		multiply(a, ph, v);
+
+		alph = rho1 / dot(rh, v, len_c);
+
+		for (int i = 0; i < len_c; i ++) {
+			s[i] = r[i] - alph * v[i];
+		}
+
+		multiply(c, s, sh);
+		multiply(a, sh, t);
+
+		omeg = dot(t, s, len_c) / dot(t, t, len_c);
+
+		for (int i = 0; i < len_c; i ++) {
+			x0[i] = x0[i] + alph * ph[i] + omeg * sh[i];
+			r[i]  = s[i] - omeg * t[i];
+		}
+
+		//memset(s, 0.0, len_c * sizeof(double));
+		//memset(t, 0.0, len_c * sizeof(double));
+
+		rho2 = rho1;
+		eps  = norm(r, len_c);
+
+		if (eps < tol) {
+			flag = 0;
+			break;
+		}
+
+		count = count + 1;
+
+	}
+
+	free(sh);
+	free(s);
+	free(v);
+	free(ph);
+	free(p);
+	free(rh);
+	free(r);
+	free(t);
+
+	printf(">>> (bic) :: iteration complete w/ err : %g, \
+		flag: %d, steps: %d\n", eps, flag, count);
+
+	return flag;
+
+}
 /*
  * preconditioned conjugate gradient solver
  * symmetric matrices only
@@ -118,7 +319,8 @@ int pcg(struct diag* a, double* b, double* x0, struct diag* c, double tol)
 	free(r1);
 	free(t);
 
-	printf(">>> (pcg) :: iteration complete w/ err : %g, flag: %d, steps: %d\n", 
+	printf(">>> (pcg) :: iteration complete w/ err : %g, \
+		flag: %d, steps: %d\n", 
 	eps, flag, count);
 
 	return flag;
diff --git a/solver.h b/solver.h
index 8ce4409..9fc5bb5 100644
--- a/solver.h
+++ b/solver.h
@@ -28,4 +28,8 @@
 
 int pcg(struct diag*, double*, double*, struct diag*, double);
 
+int bicgstab(struct diag*, double*, double*, struct diag*, double);
+
+int bic(struct diag*, double*, double*, double);
+
 #endif
diff --git a/sysmat.c b/sysmat.c
index d0c172b..d390c37 100644
--- a/sysmat.c
+++ b/sysmat.c
@@ -161,6 +161,99 @@ void system_matrix(struct diag* m, int nx, int ny, double* val)
 
 }
 
+
+/**
+ * set boundary conditions
+ *
+ * inout: m
+ **/
+void boundary(int dim, int nx, int ny, int north, int south, int east, 
+	int west, double val_north, double val_south, double val_east, 
+	double val_west, struct diag* m, double* r, double* x)
+{
+
+	// 0: dirichlet | 1: neumann | 2: closed wall
+
+
+
+	if (north == 0) {
+
+		for (int i = 0; i < nx;  i ++) {
+			x[i] = val_north;
+		}
+
+		double* ar = (double*) calloc(dim, sizeof(double));
+
+		multiply(m, x, ar);
+
+		for (int i = 0; i < dim; i ++) {
+			r[i] = r[i] - ar[i];
+			if (i < nx) {
+				r[i] = x[i];
+			}
+		}
+
+		free(ar);
+
+		for (int i = 0; i < nx; i ++) {
+			for (int j = 0; j < dim; j ++) {
+				if (i != j) {
+					put(m, j, i, 0.0);
+					put(m, i, j, 0.0);
+				} else {
+					put(m, j, i, 1.0);
+				}
+			}
+		}
+
+	} else if (north == 1) {
+
+		for (int i = 0; i < nx; i ++) {
+			r[i] = val_north;
+		}
+
+	}
+
+	if (south == 0) {
+
+		for (int i = dim - nx; i < dim; i ++) {
+			x[i] = val_south;
+		}
+
+		double* ar = (double*) calloc(dim, sizeof(double));
+
+		multiply(m, x, ar);
+
+		for (int i = 0; i < dim; i ++) {
+			r[i] = r[i] - ar[i];
+			if ((i >= dim - nx) && (i < dim)) {
+				r[i] = x[i];
+			}
+		}
+
+		free(ar);
+
+		for (int i = dim - nx; i < dim; i ++) {
+			for (int j = 0; j < dim; j ++) {
+				if (i != j) {
+					put(m, j, i, 0.0);
+					put(m, i, j, 0.0);
+				} else {
+					put(m, j, i, 1.0);
+				}
+			}
+		}
+
+	} else if (south == 1) {
+
+		for (int i = dim - nx; i < dim; i ++) {
+			r[i] = val_south;
+		}
+
+	}
+
+}
+
 /**
  * compute darcy velocities
  *
diff --git a/sysmat.h b/sysmat.h
index 3837d38..bcc2514 100644
--- a/sysmat.h
+++ b/sysmat.h
@@ -16,6 +16,9 @@
 
 void system_matrix(struct diag*, int, int, double*);
 
+void boundary(int, int, int, int, int, int, int, double, double, double, 
+	double, struct diag*, double*, double*);
+
 void darcy(double*, double*, int, int, double, double, double*, double*);
 
 double harmonic(double, double);
diff --git a/test_bic.c b/test_bic.c
new file mode 100644
index 0000000..b90334c
--- /dev/null
+++ b/test_bic.c
@@ -0,0 +1,52 @@
+/**
+ * test_bic.c
+ * ----------
+ *
+ * purpose: test solver.bicgstab()-method
+ * cf. https://en.wikipedia.org/wiki/Conjugate_gradient_method
+ * for the test case used
+ **/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "global.h"
+#include "linalg.h"
+#include "solver.h"
+
+int main (void)
+{
+
+	// init system
+	double a[6] = { 0.0, 1.0, 4.0, 3.0, 1.0, 0.0 };
+	int  ids[3] = { -1, 0, 1 };
+
+	struct diag m  = { a, ids, 2, 2, 2, 3 };
+
+	double c[2] = { 1.0 / sqrt(4.0), 1.0 / sqrt(3.0) };
+	int cids[1] = { 0 };
+
+	struct diag pre = { c, cids, 2, 2, 2, 1 }; // Jacobi preconditioner
+
+	double b[3]    = { 1.0, 2.0 };
+	double x0[2]   = { 20000.0, 10000.0 };
+
+	bicgstab(&m, b, x0, &pre, 1.0e-12);
+
+	printf(">>> x0: \n");
+	for (int i = 0; i < 2; i ++) {
+		printf(">>>  %g\n", x0[i]);
+	}
+
+	if (fabs(x0[0] - (1.0/11.0)) <= 1.0e-10) {
+		printf(">>> test passed.\n");
+	} else {
+		printf(">>> test failed. expected: %g, found: %g\n", (1.0/11.0), x0[0]);
+	}
+
+	if (fabs(x0[1] - (7.0/11.0)) <= 1.0e-10) {
+		printf(">>> test passed.\n");
+	} else {
+		printf(">>> test failed. expected: %g, found: %g\n", (7.0/11.0), x0[1]);
+	}
+}