add pcg solver

3 files changed, 114 insertions, 9 deletions

diff --git a/Makefile b/Makefile
index 711785d..3fc7fd7 100644
--- a/Makefile
+++ b/Makefile
@@ -8,7 +8,7 @@ PRG=main
 CC=cc
 CFLAGS=-Wall -Wextra -pedantic -g
 LFLAGS=
-SRC=main.c system_matrix.c linalg.c solver.c
+SRC=test_pcg.c linalg.c solver.c
 INCLUDES=
 LIB=-lm
 
diff --git a/solver.c b/solver.c
index 98c60e9..91fe1d8 100644
--- a/solver.c
+++ b/solver.c
@@ -47,26 +47,76 @@ int pcg(struct diag a, double* b, double* x0, struct diag c, double tol)
 
 	// params
 	int len_c = a.len_c;
-	int len_r = a.len_r;
+	int count = 0;
+	int flag  = 1;
 
 	double* r0 = (double*) calloc(len_c, sizeof(double));
-	double* h0 = (double*) calloc(len_c, sizeof(double));
-	double* d0 = (double*) calloc(len_c, sizeof(double));
+	double* r1 = (double*) calloc(len_c, sizeof(double));
 	double* z0 = (double*) calloc(len_c, sizeof(double));
+	double* z1 = (double*) calloc(len_c, sizeof(double));
+	double* p = (double*) calloc(len_c, sizeof(double));
+	double* t = (double*) calloc(len_c, sizeof(double));
 
-	double la0 = 0.0;
-	double rh0 = 0.0;
+	double alph = 0.0;
+	double beta = 0.0;
 
 	double eps = tol;
 
 	// bootstrap
 	r0 = multiply(a, x0);
 
+	for (int i = 0; i < len_c; i ++) {
+		r0[i] = b[i] - r0[i];
+	}
+
+	z0 = multiply(c, r0);
+	
+	memcpy(p, z0, len_c * sizeof(double));
+
+	// loop
+	while (count < 100) {
+
+		printf(">>> %d\n", count);
+
+		alph = dot(r0, z0, len_c) / dot(p, multiply(a, p), len_c);
+		t  = multiply(a, p);
+
+		for (int i = 0; i < len_c; i ++) {
+			x0[i] = x0[i] + alph * p[i];
+			r1[i] = r0[i] - alph * t[i];
+		}
+
+		eps = norm(r1, len_c);
+
+		if (eps < tol) {
+			flag = 0;
+			break;
+		} else {
+			z1   = multiply(c, r1);
+			beta = dot(z1, r1, len_c) / dot(z0, r0, len_c);
+
+			for (int i = 0; i < len_c; i ++) {
+				p[i] = z1[i] + beta * p[i];
+			}
+
+		memcpy(z0, z1, len_c * sizeof(double));
+		memcpy(r0, r1, len_c * sizeof(double));
+
+		}
+
+		count = count + 1;
+	}
+
 	free(z0);
-	free(d0);
-	free(h0);
+	free(z1);
+	free(p);
 	free(r0);
+	free(r1);
+	free(t);
+
+	printf(">>> iteration complete w/ err : %g\n", eps);
+	printf(">>> iteration complete w/ flag: %d\n", flag);
 
-	return 0;
+	return flag;
 
 }
diff --git a/test_pcg.c b/test_pcg.c
new file mode 100644
index 0000000..4b4eef0
--- /dev/null
+++ b/test_pcg.c
@@ -0,0 +1,55 @@
+/**
+ * test_pcg.c
+ * ----------
+ *
+ * purpose: test solver.pcg()-method
+ * cf. https://en.wikipedia.org/wiki/Conjugate_gradient_method
+ * for the test case used
+ **/
+
+#include <stdio.h>
+#include <stdlib.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, 3 };
+	
+	for (int i = 0; i < 2; i ++) {
+		for (int j = 0; j < 2; j ++) {
+			printf(" %g ", get(m, i, j));
+		}
+		printf("\n");
+	}
+	printf("\n");
+
+	double c[2] = { 1.0 / 4.0, 1.0 / 3.0 };
+	int cids[1] = { 0 };
+
+	struct diag pre = { c, cids, 2, 2, 1 }; // Jacobi preconditioner
+	
+	for (int i = 0; i < 2; i ++) {
+		for (int j = 0; j < 2; j ++) {
+			printf(" %g ", get(pre, i, j));
+		}
+		printf("\n");
+	}
+	printf("\n");
+
+	double b[3]    = { 1.0, 2.0 };
+	double x0[2]   = { 10000.0, 1000.0 };
+
+	printf(">>> iteration finished w/ flag: %d\n", pcg(m, b, x0, pre, 1.0e-3));
+	printf(">>> x0: \n");
+	for (int i = 0; i < 2; i ++) {
+		printf(">>>  %g\n", x0[i]);
+	}
+
+}