add bicgstab (not working)

5 files changed, 185 insertions, 18 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/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 1208dad..808ab89 100644
--- a/solver.c
+++ b/solver.c
@@ -35,6 +35,99 @@
 
 // ----------------------------------------------------------------------------
 
+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
  *
@@ -56,29 +149,34 @@ int bicgstab(struct diag* a, double* b, double* x0, struct diag* c, double tol)
 	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 norm_s = 0.0;
 
 	double eps = tol;
 
 	// bootstrap
-	multiply(a, x0, r0);
+	multiply(a, x0, r);
 
 	for (int i = 0; i < len_c; i ++) {
-		r0[i] = b[i] - r0[i];
+		r[i] = b[i] - r[i];
 	}
 
-	memcpy(r, rh, len_c * sizeof(double));
+	memcpy(rh, r, len_c * sizeof(double));
 
 	// loop
 	while (count < 1000) {
 
-		rho1 = dot(r0h, r0, len_c);
+		rho1 = dot(rh, r, len_c);
+		if (rho1 < 1.0e-18) {
+			flag = 2;
+			break;
+		}
 
 		if (count == 1) {
 			memcpy(p, r, len_c * sizeof(double));
@@ -95,33 +193,45 @@ int bicgstab(struct diag* a, double* b, double* x0, struct diag* c, double tol)
 		alph = rho1 / dot(rh, v, len_c);
 
 		for (int i = 0; i < len_c; i ++) {
-			s[i] = r[i] - alpha * v[i];
-			x0[i] = alpha * ph[i];
+			s[i] = r[i] - alph * v[i];
 		}
 
-		norm_s = norm(s, len_c);
+		multiply(c, s, sh);
+		multiply(a, sh, t);
 
-		if (norms < tol) {
-			for (int i = 0; i < len_c; i ++) {
-				tol = norm_s;
-				flag = 0;
-				break;
-			}
+		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(">>> (pcg) :: iteration complete w/ err : %g, flag: %d, steps: %d\n", 
-	eps, flag, count);
+	printf(">>> (bic) :: iteration complete w/ err : %g, \
+		flag: %d, steps: %d\n", eps, flag, count);
 
 	return flag;
 
@@ -209,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/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]);
+	}
+}