debug and extract methods

8 files changed, 283 insertions, 160 deletions

diff --git a/Makefile b/Makefile
index 4e12936..1055d3d 100644
--- a/Makefile
+++ b/Makefile
@@ -4,11 +4,11 @@
 # make clean     : delete main program and all object files
 # make mrpropper : delete main program, all object files and results
 
-PRG=main
+PRG=test
 CC=cc
 CFLAGS=-Wall -Wextra -pedantic -g
 LFLAGS=
-SRC=main.c system_matrix.c bicgstab.c
+SRC=test.c system_matrix.c diagonal.c bicgstab.c
 INCLUDES=
 LIB=-lm
 
diff --git a/bicgstab.c b/bicgstab.c
index 733c0d4..f2265ea 100644
--- a/bicgstab.c
+++ b/bicgstab.c
@@ -28,6 +28,7 @@
 #include <stdio.h>
 #include <math.h>
 #include "bicgstab.h"
+#include "diagonal.h"
 
 // ----------------------------------------------------------------------------
 
@@ -52,17 +53,9 @@ int bicgstab(const double* a, const double* b, double* x, const int* diag_ids,
 
 	int dim = ncol * nrow;
 
-	double err   = 0.0;
-	double rho1  = 1.0;
-	double rho2  = 1.0;
-	double alpha = 1.0;
-	double beta  = 1.0;
-	double omega = 1.0;
-
 	double* r  = calloc(dim, sizeof(double)); // r = b - A * x
 	double* ax = calloc(dim, sizeof(double)); // A * x
-	double* p  = calloc(dim, sizeof(double));
-	double* v  = calloc(dim, sizeof(double));
+
 	double* s  = calloc(dim, sizeof(double));
 	double* t  = calloc(dim, sizeof(double));
 
@@ -72,59 +65,54 @@ int bicgstab(const double* a, const double* b, double* x, const int* diag_ids,
 			r[i] = b[i] - ax[i];
 	}
 
-	double* rtilde = (double*) calloc(dim, sizeof(double));
-	memcpy(rtilde, r, dim * sizeof(double));
-
 	double normb = norm(b, dim);
-
-	if (normb == 0.0) {
+	if (fabs(normb) < 1.0e-15) {
 		normb = 1.0;
 	}
 
+	double err = norm(r, dim) / normb;
 
-	// if the error is already less then our tolerance, 
-	// we are done here.
-	err = norm(r, dim) / normb;
-
-	if (err <= tol) {
-		tol = err;
-		max_iter = 0;
-		free_mem(t, s, v, p, rtilde, r, ax);
-		return 0;
-	}
+	double* rtilde = (double*) calloc(dim, sizeof(double));
+	memcpy(rtilde, r, dim * sizeof(double));
+	
+	double rho0 = 1.0;
+	double rho1 = 0.0;
+	double alph = 1.0;
+	double ome0 = 1.0;
+	double ome1 = 0.0;
+	
+	double* p  = calloc(dim, sizeof(double));
+	double* v  = calloc(dim, sizeof(double));
 
-	// else, our initial guess was not good enough.
 	for (int i = 1; i <= max_iter; i ++) {
-
+		
 		rho1 = dot(rtilde, r, dim);
 
-		printf(" >>> %d\n", i);
-		printf("rho1 : %f\n", rho1);
-
-		if (sqrt(rho1*rho1) < 1.0e-18) {
+		if (fabs(rho1) < 1.0e-15) {
 			tol = err;
 			free_mem(t, s, v, p, rtilde, r, ax);
 			return 2;
 		}
 
-		if (i == 1) {
-			memcpy(p, r, dim * sizeof(double));
-		} else {
-			beta = (rho1 / rho2) * (alpha / omega);
-			for (int j = 0; j < dim; j ++) {
-				p[j] = r[j] + beta * (p[j] - omega * v[j]);
-			}
+		double beta = (rho1/rho0) * (alph/ome0);
+
+		for (int j = 0; j < dim; j ++) {
+			p[j] = r[j] + beta * (p[j] - ome0 * v[j]);
 		}
 
 		v = multiply(a, p, diag_ids, ncol, nrow, ndiag);
-		alpha = rho1 / dot(rtilde, v, dim);
+
+		alph = rho1 / dot(rtilde, v, dim);
+
 		for (int j = 0; j < dim; j ++) {
-			s[j] = r[j] - alpha * v[j];
+			s[j] = r[j] - alph * v[j];
 		}
+
 		err = norm(s, dim) / normb;
+
 		if (err < tol) {
 			for (int j = 0; j < dim; j ++) {
-				x[j] = x[j] + alpha * p[j];
+				x[j] = x[j] + alph * p[j];
 			}
 			tol = err;
 			max_iter = i;
@@ -133,15 +121,19 @@ int bicgstab(const double* a, const double* b, double* x, const int* diag_ids,
 		}
 
 		t = multiply(a, s, diag_ids, ncol, nrow, ndiag);
-		omega = dot(t,s, dim) / dot(t,t, nrow);
+
+		ome1 = dot(t, s, dim) / dot(t, t, dim);
+
 		for (int j = 0; j < dim; j ++) {
-			x[j] = x[j] + omega * s[j];
+			x[j] = x[j] + ome1 * s[j];
 		}
+
 		for (int j = 0; j < dim; j ++) {
-			r[j] = s[j] - omega * t[j];
+			r[j] = s[j] - ome1 * t[j];
 		}
 
-		rho2 = rho1;
+		rho0 = rho1;
+		ome0 = ome1;
 
 		err = norm(r, dim) / normb;
 		if (err < tol) {
@@ -151,6 +143,9 @@ int bicgstab(const double* a, const double* b, double* x, const int* diag_ids,
 			return 0;
 		}
 
+		printf(" >>> %d\n", i);
+		printf(" >>> err: %g\n", err);
+
 	}
 
 	tol = err;
@@ -158,104 +153,6 @@ int bicgstab(const double* a, const double* b, double* x, const int* diag_ids,
 	return 1;
 }
 
-
-/**
- * return the L2-norm of the vector
- *
- * input:
- *	v : vector
- *	n : dimension of v
- *
- * output:
- *	L2-norm (magnitude) of v
- **/
-double norm(const double* v, double n)
-{
-	double normv = 0.0;
-	for (int i = 0; i < n; i ++) {
-		normv = normv + (v[i] * v[i]);
-	}
-
-	return sqrt(normv);
-}
-
-
-/**
- * calculate matrix vector multiplication
- *
- * input:
- *	a       : matrix in diagonal form
- *	v       : vector of dimension nrow
- *	diag_ids: diagonal ids
- *	ncol    : number of columns
- *	nrow    : number of rows
- *	ndiag   : number of non-zero diagonals
- *
- * output:
- *	product
- **/
-double* multiply(const double* a, const double* v, const int* ids, 
-                 int ncol, int nrow, int ndia)
-{
-	double* av = calloc(ncol, sizeof(double));
-	double dim = ncol * nrow;
-
-	for (int i = 0; i < nrow; i ++) {
-		for (int j = 0; j < ncol; j ++) {
-			av[i] = av[i] + get_dia(a, ids, ndia, dim, i, j) * v[j];
-		}
-	}
-
-	return av;
-}
-
-/**
- * put a val into the matrix a
- **/
-void put_dia(const double* a, const int* ids, int ndia, int dim, 
-             int r, int c, double val)
-{
-	int id = c - r;
-
-	for (int i = 0; i < ndia; i ++) {
-		if (ids[i] == id) {
-			dia[i * dim + r] = val;
-		}
-	}
-}
-
-/**
- * get a value (r,c) from matrix a
- **/
-double get_dia(const double* a, const int* ids, int ndia, int dim,
-                int r, int c)
-{
-	int id = c - r;
-	double val = 0.0;
-
-	for (int i = 0; i < ndia; i ++) {
-		if (ids[i] == id) {
-			val = a[i * dim + r];
-		}
-	}
-
-	return val;
-}
-
-/**
- * return the dot product of vectors v1 and v2
- **/
-double dot(const double* v1, const double* v2, int n)
-{
-	double dot_product = 0.0;
-	for (int i = 0; i < n; i ++) {
-		dot_product = dot_product + v1[i] * v2[i];
-	}
-	return dot_product;
-}
-
-
-
 /**
  * free memory
  **/
diff --git a/bicgstab.h b/bicgstab.h
index 6616154..a83983f 100644
--- a/bicgstab.h
+++ b/bicgstab.h
@@ -27,12 +27,6 @@
 int bicgstab(const double*, const double*, double*, const int*, int, int, 
              int, int, double);
 
-double norm(const double*, double);
-
 void free_mem(double*, double*, double*, double*, double*, double*, double*);
 
-double dot(const double*, const double*, int);
-
-double* multiply(const double*, const double*, const int*, int, int, int);
-
 #endif
diff --git a/diagonal.c b/diagonal.c
new file mode 100644
index 0000000..9d2c642
--- /dev/null
+++ b/diagonal.c
@@ -0,0 +1,90 @@
+/**
+ * diagonal.c
+ * ----------
+ *
+ * provide functionality of a sparse matrix in diagonal format
+ *
+ * @author ilhan oezgen, wahyd, ilhan.ozgen@wahyd.tu-berlin.de
+ * @date 15 apr 17
+ **/
+
+// ----------------------------------------------------------------------------
+
+#include<stdlib.h>
+#include<math.h>
+
+// ----------------------------------------------------------------------------
+
+/**
+ * put a val into the matrix a
+ **/
+void put_dia(double* a, const int* ids, int n, int len_dia, 
+             int r, int c, double val)
+{
+	int id = c - r;
+	
+	for (int i = 0; i < n; i ++) {
+		if (ids[i] == id) {
+			a[i * len_dia + r] = val;
+		}
+	}
+}
+
+/**
+ * get a value (r,c) from matrix a
+ **/
+double get_dia(const double* a, const int* ids, int n, int len_dia,
+               int r, int c)
+{
+	int id = c - r;
+	double val = 0.0;
+
+	for (int i = 0; i < n; i ++) {
+		if (ids[i] == id) {
+			val = a[i * len_dia + r];
+		}
+	}
+
+	return val;
+}
+
+/**
+ * calculate matrix vector multiplication
+ **/
+double* multiply(const double* a, const double* v, const int* ids, 
+                 int len_r, int len_c, int n)
+{
+	double* av = calloc(len_c, sizeof(double));
+	double len_dia = len_c; // length of main diagonal
+	if (len_c < len_r) {
+		len_dia = len_r;
+	}
+
+	for (int i = 0; i < len_r; i ++) {
+		for (int j = 0; j < len_c; j ++) {
+			av[i] = av[i] + get_dia(a, ids, n, len_dia, i, j) * v[j];
+		}
+	}
+
+	return av;
+}
+
+/**
+ * return the dot product of vectors v1 and v2
+ **/
+double dot(const double* v1, const double* v2, int len)
+{
+	double dot_product = 0.0;
+	for (int i = 0; i < len; i ++) {
+		dot_product = dot_product + v1[i] * v2[i];
+	}
+	return dot_product;
+}
+
+/**
+ * return the L2-norm of the vector
+ **/
+double norm(const double* v, int len)
+{
+	return sqrt(dot(v, v, len));
+}
diff --git a/diagonal.h b/diagonal.h
new file mode 100644
index 0000000..70165df
--- /dev/null
+++ b/diagonal.h
@@ -0,0 +1,24 @@
+/**
+ * diagonal.h
+ * ----------
+ *
+ * provide functionality of a sparse matrix in diagonal format
+ *
+ * @author ilhan oezgen, wahyd, ilhan.ozgen@wahyd.tu-berlin.de
+ * @date 15 apr 17
+ **/
+
+#ifndef __DIAGON_INCLUDED__
+#define __DIAGON_INCLUDED__
+
+void put_dia(double*, const int*, int, int, int, int, double);
+ 
+double get_dia(double*, const int*, int, int, int, int);
+
+double* multiply(const double*, const double*, const int*, int, int, int);
+
+double dot(const double*, const double*, int);
+
+double norm(const double*, int);
+
+#endif
diff --git a/main.c b/main.c
index 4a97634..c552806 100644
--- a/main.c
+++ b/main.c
@@ -13,6 +13,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include "system_matrix.h"
+#include "diagonal.h"
 #include "bicgstab.h"
 
 int main (void)
@@ -31,8 +32,8 @@ int main (void)
 	double dx = (xR - xL) / nx;
 	double dy = (yU - yD) / ny;
 
-	double kf1 = 1.0;
-	double kf2 = 10.0;
+	double kf1 = 1.0e-12;
+	double kf2 = 1.0e-11;
 
 	// initialize permeability matrix
 
@@ -82,14 +83,9 @@ int main (void)
 		r[i] = 10.0 * kf[i];
 		for (int j = 0; j < dim; j ++) {
 			if (i != j) {
-				int id = j - i;
-				for (int k = 0; k < 5; k ++) {
-					if ( diag_ids[k] == id ) {
-						amod[k * dim + i] = 0.0;
-					}
-				}
+				put_dia(amod, diag_ids, 5, dim, j, i, 0.0);
 			} else {
-				amod[2 * dim + i] = kf[i]; // main diagonal
+				put_dia(amod, diag_ids, 5, dim, j, i, kf[i]);
 			}
 		}
 	}
@@ -125,7 +121,7 @@ int main (void)
 		x[i] = 1.0;
 	}
 
-	int flag = bicgstab(a, r, x, diag_ids, nx, ny, 5, 1000, 1.0e-3);
+	int flag = bicgstab(a, r, x, diag_ids, nx, ny, 5, 15, 1.0e-3);
 	printf("\n\niteration finished w/ flag : %d\n", flag);
 
 	printf("\nsolution:\n");
diff --git a/system_matrix.c b/system_matrix.c
index b7dc3aa..abb9436 100644
--- a/system_matrix.c
+++ b/system_matrix.c
@@ -16,6 +16,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include "system_matrix.h"
+#include "diagonal.h"
 
 // ----------------------------------------------------------------------------
 //
diff --git a/test.c b/test.c
new file mode 100644
index 0000000..f878417
--- /dev/null
+++ b/test.c
@@ -0,0 +1,121 @@
+/**
+ * test class
+ **/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "bicgstab.h"
+#include "diagonal.h"
+
+int main (void)
+{
+
+	double** a = (double**) calloc(5, sizeof(double*));
+	for (int i = 0; i < 5; i++) {
+		a[i] = (double*) calloc(5, sizeof(double));
+	}
+
+  a[0][0] = 2.0;
+	a[1][1] = 3.0;
+	a[2][2] = 4.0;
+	a[3][3] = 5.0;
+	a[4][4] = 6.0;
+
+	a[0][1] = 1.0;
+	a[1][2] = 1.1;
+	a[2][3] = 1.2;
+	a[3][4] = 1.3;
+
+	a[1][0] = -1.0;
+	a[2][1] = -1.1;
+	a[3][2] = -1.2;
+	a[4][3] = -1.3;
+
+	for (int i = 0; i < 5; i ++) {
+		for (int j = 0; j < 5; j ++) {
+			if (a[i][j] >= 0.0) {
+			printf(" +%1.2f ", a[i][j]);
+			} else {
+			printf(" %1.2f ", a[i][j]);
+			}
+		}
+		printf("\n");
+	}
+
+	double* dia = (double*) calloc(3*5, sizeof(double*));
+	int ids[3] = { -1, 0, 1 };
+
+	for (int i = 0; i < 5; i ++) {
+		for (int j = 0; j < 5; j ++) {
+			put_dia(dia, ids, 3, 5, i, j, a[i][j]);
+		}
+	}
+
+	printf("\n");
+
+	for (int i = 0; i < 5; i ++) {
+		for (int j = 0; j < 5; j ++) {
+			double val = get_dia(dia, ids, 3, 5, i, j);
+			if (fabs(val - a[i][j]) > 1.0e-12) {
+				printf("test failed!");
+			}
+			if (val >= 0.0) {
+				printf(" +%1.2f ", val);
+			} else {
+				printf(" %1.2f ", val);
+			}
+		}
+		printf("\n");
+	}
+
+
+	for (int i = 0; i < 5; i ++) {
+		free(a[i]);
+	}
+
+	double* v = (double*) calloc(5, sizeof(double*));
+	
+	v[0] = 1.0;
+	v[1] = 2.0;
+	v[2] = 0.0;
+	v[3] = 4.0;
+	v[4] = 5.0;
+
+	printf("\n");
+	for (int i = 0; i < 5; i ++) {
+		printf(" %f ", v[i]);
+	}
+	printf("\n");
+
+	printf("\n");
+	if (fabs(norm(v,5) - 6.7823299831252681) < 1.0e-12) {
+		printf("norm test passed\n");
+	} else {
+		printf("norm: %g \n", norm(v, 5));
+	}
+
+	if (fabs(dot(v,v,5) - 46.0) < 1.0e-12) {
+		printf("dot test passed\n");
+	} else {
+		printf("dot : %g \n", dot(v, v, 5));
+	}
+
+	double* dia_v = multiply(dia, v, ids, 5, 5, 3);
+
+	double ref[5] = { 4.0, 5.0, 2.6, 26.5, 24.8 };
+
+	printf("\n");
+
+	for (int i = 0; i < 5; i ++) {
+		if (fabs(ref[i] - dia_v[i]) > 1.0e-12) {
+			printf("test failed in multiply");
+		}
+		printf(" %f ", dia_v[i]);
+	}
+	printf("\n");
+
+	free(v);
+	free(dia);
+
+}