prepare for bicgstab

1 files changed, 91 insertions, 0 deletions

diff --git a/solver.c b/solver.c
index 4427217..1208dad 100644
--- a/solver.c
+++ b/solver.c
@@ -36,6 +36,97 @@
 // ----------------------------------------------------------------------------
 
 /*
+ * 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 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);
+
+	for (int i = 0; i < len_c; i ++) {
+		r0[i] = b[i] - r0[i];
+	}
+
+	memcpy(r, rh, len_c * sizeof(double));
+
+	// loop
+	while (count < 1000) {
+
+		rho1 = dot(r0h, r0, len_c);
+
+		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] - alpha * v[i];
+			x0[i] = alpha * ph[i];
+		}
+
+		norm_s = norm(s, len_c);
+
+		if (norms < tol) {
+			for (int i = 0; i < len_c; i ++) {
+				tol = norm_s;
+				flag = 0;
+				break;
+			}
+		}
+
+		count = count + 1;
+
+	}
+
+	free(s);
+	free(v);
+	free(ph);
+	free(p);
+	free(rh);
+	free(r);
+
+	printf(">>> (pcg) :: iteration complete w/ err : %g, flag: %d, steps: %d\n", 
+	eps, flag, count);
+
+	return flag;
+
+}
+/*
  * preconditioned conjugate gradient solver
  * symmetric matrices only
  * c    : preconditioner matrix