fix bicgstab

3 files changed, 95 insertions, 158 deletions

diff --git a/solver.c b/solver.c
index 808ab89..e5c2899 100644
--- a/solver.c
+++ b/solver.c
@@ -20,6 +20,7 @@
  * @date 10 apr 17
  * @date 11 apr 17 : add norm function
  * @date 28 apr 17 : add pcg solver
+ * @date  8 jun 17 : add bicgstab solver
  *
  **/
 
@@ -28,6 +29,7 @@
 #include <string.h>
 #include <stdlib.h>
 #include <stdio.h>
+#include <math.h>
 
 #include "global.h" // struct diag
 #include "solver.h"
@@ -35,207 +37,143 @@
 
 // ----------------------------------------------------------------------------
 
-int bic(struct diag* a, double* b, double* x0, double tol)
+int bic(struct diag* a, double* b, double* x0, struct diag* m, 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* r      = (double*) calloc(len_c, sizeof(double));
+	double* r_tld  = (double*) calloc(len_c, sizeof(double));
+	double* p      = (double*) calloc(len_c, sizeof(double));
+	double* v      = (double*) calloc(len_c, sizeof(double));
+	double* p_hat  = (double*) calloc(len_c, sizeof(double));
+	double* s      = (double*) calloc(len_c, sizeof(double));
+	double* s_hat  = (double*) calloc(len_c, sizeof(double));
+	double* t      = (double*) calloc(len_c, sizeof(double));
+
+	double eps   = tol;
+	double omega = 1.0;
 	double rho   = 1.0;
-	double alp   = 1.0;
-	double ome   = 1.0;
+	double rho_1 = 1.0;
+	double beta  = 1.0;
+	double alpha = 1.0;
 
-	double rh1   = 0.0;
-	double bet   = 0.0;
-	double eps   = 0.0;
 
-	multiply(a, x0, r);
+	double bnrm = norm( b, len_c);
+
+	if (fabs(bnrm) < tol) {
+		bnrm = 1.0;
+	}
 
+	// r = b - A x
+	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));
+	eps = norm(r, len_c) / bnrm;
 
-	memset(v, 0.0, len_c * sizeof(double));
-	memset(p, 0.0, len_c * sizeof(double));
+	if (eps < tol) {
+		flag = 0;
+		printf(">>> (bic) :: iteration complete w/ err : %g, \
+			flag: %d, steps: %d\n", eps, flag, count);
+		return flag;
+	}
 
-	for (int i = 1; i < 1000; i ++) {
+	memcpy(r_tld, r, len_c * sizeof(double));
 
-		eps = norm(r, len_c);
-		if (eps < tol) {
-			flag = 0;
-			break;
-		}
+	for (count = 1; count < 1000; count ++) {
 
-		rh1 = dot(rhat, r, len_c);
-		bet = (rh1 / rho) * (alp / ome);
+		rho = dot(r_tld, r, len_c);
 
-		for (int k = 0; k < len_c; k ++) {
-			p[k] = r[k] + bet * (p[k] - ome * v[k]);
+		if (fabs(rho) < tol) {
+			break;
 		}
 
-		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];
+		if (count > 1) {
+			beta = (rho / rho_1) * (alpha * omega);
+			for (int j = 0; j < len_c; j ++) {
+				p[j] = r[j] + beta * (p[j] - omega * v[j]);
+			}
+		} else {
+			memcpy(p, r, len_c * sizeof(double));
 		}
 
-		multiply(a, s, t);
-
-		ome = dot(t, s, len_c) / dot(t, t, len_c);
+		multiply(m, p, p_hat);
+		multiply(a, p_hat, v);
 
-		for (int k = 0; k < len_c; k ++) {
-			x0[k] = h[k] + ome * s[k];
-			r[k]  = s[k] - ome * t[k];
+		alpha = rho / dot(r_tld, v, len_c);
+		for (int j = 0; j < len_c; j ++) {
+			s[j] = r[j] - alpha * v[j];
 		}
 
-		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];
-	}
+		double snrm = norm(s, len_c);
 
-	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]);
+		if (snrm < tol) {
+			for (int j = 0; j < len_c; j ++) {
+				x0[j] = x0[j] + alpha + p_hat[j];
 			}
+			eps = snrm / bnrm;
+			break;
 		}
 
-		multiply(c, p, ph);
-		multiply(a, ph, v);
+		multiply(m, s, s_hat);
+		multiply(a, s_hat, t);
 
-		alph = rho1 / dot(rh, v, len_c);
+		omega = dot(t, s, len_c) / dot(t, t, len_c);
 
-		for (int i = 0; i < len_c; i ++) {
-			s[i] = r[i] - alph * v[i];
+		for (int j = 0; j < len_c; j ++) {
+			x0[j] = x0[j] + alpha * p_hat[j] + omega * s_hat[j];
+			r[j] = s[j] - omega * t[j];
 		}
 
-		multiply(c, s, sh);
-		multiply(a, sh, t);
+		eps = norm(r, len_c) / bnrm;
 
-		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];
+		if (eps < tol) {
+			break;
 		}
 
-		//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;
+		if (fabs(omega) < tol) {
 			break;
 		}
 
-		count = count + 1;
+		rho_1 = rho;
+		memset(p_hat, 0.0, len_c * sizeof(double));
+		memset(v, 0.0, len_c * sizeof(double));
+		memset(s_hat, 0.0, len_c * sizeof(double));
+		memset(t, 0.0, len_c * sizeof(double));
 
 	}
 
-	free(sh);
-	free(s);
-	free(v);
-	free(ph);
-	free(p);
-	free(rh);
-	free(r);
-	free(t);
+	if (eps < tol) {
+		flag = 0;
+	} else if (fabs(omega) < tol) {
+		flag = -2;
+	} else if (fabs(rho) < tol) {
+		flag = -1;
+	} else {
+		flag = 1;
+	}
 
 	printf(">>> (bic) :: iteration complete w/ err : %g, \
 		flag: %d, steps: %d\n", eps, flag, count);
 
+	free(r);
+	free(r_tld);
+	free(p);
+	free(v);
+	free(p_hat);
+	free(s);
+	free(s_hat);
+	free(t);
+
 	return flag;
 
 }
+
+
 /*
  * preconditioned conjugate gradient solver
  * symmetric matrices only
@@ -320,8 +258,7 @@ int pcg(struct diag* a, double* b, double* x0, struct diag* c, double tol)
 	free(t);
 
 	printf(">>> (pcg) :: iteration complete w/ err : %g, \
-		flag: %d, steps: %d\n", 
-	eps, flag, count);
+		flag: %d, steps: %d\n", eps, flag, count);
 
 	return flag;
 
diff --git a/solver.h b/solver.h
index 9fc5bb5..2404897 100644
--- a/solver.h
+++ b/solver.h
@@ -30,6 +30,6 @@ 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);
+int bic(struct diag*, double*, double*, struct diag*, double);
 
 #endif
diff --git a/test_bic.c b/test_bic.c
index b90334c..bd5399b 100644
--- a/test_bic.c
+++ b/test_bic.c
@@ -23,15 +23,15 @@ int main (void)
 
 	struct diag m  = { a, ids, 2, 2, 2, 3 };
 
-	double c[2] = { 1.0 / sqrt(4.0), 1.0 / sqrt(3.0) };
+	double c[2] = { 1.0, 1.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 };
+	double b[2]    = { 1.0, 2.0 };
+	double x0[2]   = { 20000.0, 1000.0 };
 
-	bicgstab(&m, b, x0, &pre, 1.0e-12);
+	bic(&m, b, x0, &pre, 1.0e-12);
 
 	printf(">>> x0: \n");
 	for (int i = 0; i < 2; i ++) {