change to pass-by-reference

12 files changed, 312 insertions, 310 deletions

diff --git a/Makefile b/Makefile
index 50bbef7..3fc7fd7 100644
--- a/Makefile
+++ b/Makefile
@@ -8,7 +8,7 @@ PRG=main
 CC=cc
 CFLAGS=-Wall -Wextra -pedantic -g
 LFLAGS=
-SRC=main.c linalg.c solver.c system_matrix.c
+SRC=test_pcg.c linalg.c solver.c
 INCLUDES=
 LIB=-lm
 
diff --git a/global.h b/global.h
index 97423b0..8100057 100644
--- a/global.h
+++ b/global.h
@@ -26,6 +26,7 @@ struct diag {
 
   int len_c;
   int len_r;
+	int len_d;
   int len_id;
 
 };
diff --git a/linalg.c b/linalg.c
index e96560d..e6ec515 100644
--- a/linalg.c
+++ b/linalg.c
@@ -6,6 +6,7 @@
  *
  * @author ilhan oezgen, wahyd, ilhan.oezgen@wahyd.tu-berlin.de
  * @date 15 apr 17
+ * @date  3 may 17
  **/
 
 // ----------------------------------------------------------------------------
@@ -21,62 +22,54 @@
 /**
  * put a val into the diag matrix dia
  **/
-void put(struct diag dia, int r, int c, double val)
+void put(struct diag* dia, int r, int c, double val)
 {
-	int len_ids =  dia.len_id;
-	int len_d =  dia.len_c;
-	if (dia.len_r > len_d) {
-		len_d = dia.len_r;
-	}
 
-	int* ids    = &dia.ids[0];
-	double* a   = &dia.values[0];
+	int len_ids = dia->len_id;
+	int len_d   = dia->len_d;
 
 	int id = c - r;
 
 	for (int i = 0; i < len_ids; i ++) {
-		if (ids[i] == id) {
-			a[i * len_d + r] = val;
+		if (dia->ids[i] == id) {
+			dia->values[i * len_d + r] = val;
+			break;
 		}
 	}
+
 }
 
 /**
  * get a value (r,c) from diag matrix dia
  **/
-double get(struct diag dia, int r, int c)
+double get(struct diag* dia, int r, int c)
 {
 
-	int len_ids =  dia.len_id;
-	int len_d =  dia.len_c;
-	if (dia.len_r > len_d) {
-		len_d = dia.len_r;
-	}
-
-	int* ids    = &dia.ids[0];
-	double* a   = &dia.values[0];
+	int len_ids =  dia->len_id;
+	int len_d   =  dia->len_d;
 
 	int id = c - r;
 	double val = 0.0;
 
 	for (int i = 0; i < len_ids; i ++) {
-		if (ids[i] == id) {
-			val = a[i * len_d + r];
+		if (dia->ids[i] == id) {
+			val = dia->values[i * len_d + r];
+			break;
 		}
 	}
 
 	return val;
+
 }
 
 /**
  * calculate matrix vector multiplication
  **/
-double* multiply(struct diag dia, double* v)
+void multiply(struct diag* dia, double* v, double* av)
 {
-	int len_c = dia.len_c;
-	int len_r = dia.len_r;
 
-	double* av = calloc(len_c, sizeof(double));
+	int len_c = dia->len_c;
+	int len_r = dia->len_r;
 
 	for (int i = 0; i < len_r; i ++) {
 		for (int j = 0; j < len_c; j ++) {
@@ -84,7 +77,6 @@ double* multiply(struct diag dia, double* v)
 		}
 	}
 
-	return av;
 }
 
 /**
@@ -106,3 +98,4 @@ double norm(double* v, int len)
 {
 	return sqrt(dot(v, v, len));
 }
+
diff --git a/linalg.h b/linalg.h
index 20c2436..3fbce41 100644
--- a/linalg.h
+++ b/linalg.h
@@ -11,11 +11,11 @@
 #ifndef __LINALG_INCLUDED__
 #define __LINALG_INCLUDED__
 
-void put(struct diag, int, int, double);
+void put(struct diag*, int, int, double);
 
-double get(struct diag, int, int);
+double get(struct diag*, int, int);
 
-double* multiply(struct diag, double* v);
+void multiply(struct diag*, double* v, double* av);
 
 double dot(double*, double*, int);
 
diff --git a/main.c b/main.c
index 55135d9..55fb5e2 100644
--- a/main.c
+++ b/main.c
@@ -7,13 +7,15 @@
  * @author ilhan oezgen, wahyd, ilhan.oezgen@wahyd.tu-berlin.de
  * @date 6 apr 17
  *      10 apr 17 : switch to 1d vectors
+ *       2 may 17 : switch to structs
  **/
 
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+
 #include "global.h"
-#include "system_matrix.h"
+#include "sysmat.h"
 #include "linalg.h"
 #include "solver.h"
 
@@ -21,8 +23,8 @@ int main (void)
 {
 
 	// build the domain
-	int nx = 5;
-	int ny = 5;
+	int nx = 3;
+	int ny = 3;
 	int dim = nx * ny;
 
 	double xL = -1.0; // left boundary
@@ -33,8 +35,8 @@ int main (void)
 	double dx = (xR - xL) / nx;
 	double dy = (yU - yD) / ny;
 
-	double kf1 = 1.0e-12;
-	double kf2 = 1.0e-12;
+	double kf1 = 1.0;
+	double kf2 = 1.0;
 
 	// initialize permeability matrix
 
@@ -54,79 +56,87 @@ int main (void)
 
 	// initialize system matrix
 	//
-
-	int ids[5] = { -nx, -1, 0, 1, nx };
-	double* a = system_matrix(nx, ny, kf);
-
-	struct diag dia;
-	dia.values = a;
-	dia.ids = &ids[0];
-	dia.len_c = dim;
-	dia.len_r = dim;
-	dia.len_id = 5;
-
-	// initialize boundary conditions
-	//
-
-	double* amod = (double*) calloc(5 * dim, sizeof(double));
-	amod = memcpy(amod, a, 5 * dim * sizeof(double));
-
-	double* r = (double*) calloc(dim, sizeof(double));
-
-	// dirichlet boundary
-	for (int i = 0; i < nx; i ++) {
-		r[i] = 10.0 * kf[i];
+	// store diagonal ids for each non-zero diagonal
+	int ids[5] = { - nx, -1, 0, 1, nx };
+	// giving a matrix that with shape (5, dim)
+	double* a = (double*) calloc(5 * dim, sizeof(double));
+
+	struct diag m;
+	m.values = a;
+	m.ids    = ids;
+	m.len_c  = dim;
+	m.len_r  = dim;
+	m.len_id = 5;
+
+	m = system_matrix(m, nx, ny, kf);
+	printf("... \n");
+	for (int i = 0; i < dim; i ++) {
 		for (int j = 0; j < dim; j ++) {
-			if (i != j) {
-				put(dia, j, i, 0.0);
-			} else {
-				put(dia, j, i, kf[i]);
-			}
+			printf(" %g ", get(m, i, j));
 		}
-	}
-	
-	dia.values = amod;
-
-	// neumann boundary
-	for (int i = dim - nx; i < dim; i ++) {
-		r[i] = -1.8e-5;
+		printf("\n");
 	}
 
-	// solve the system A*x = b
+	// initialize boundary conditions
 	//
 
-	double* x = (double*) calloc(dim, sizeof(double));
-	for (int i = 0; i < dim; i ++) {
-		x[i] = 1.0;
-	}
-
-	// define preconditioner
-	double* c = (double*) calloc(dim, sizeof(double));
-	for (int i = 0; i < dim; i ++) {
-		c[i] = 1.0 / get(dia, i, i);
-	}
-	int pre_ids[1] = { 0 };
-
-	struct diag pre = { c, pre_ids, dim, dim, 1 };
-
-	int flag = pcg(dia, r, x, pre, 1.0e-6);
-	printf("\n\niteration finished w/ flag : %d\n", flag);
-
-	printf("\nsolution:\n");
-	for (int i = 0; i < ny; i ++) {
-		for (int j = 0; j < nx; j ++) {
-			printf("%g; ", x[i * nx + j]);
-		}
-		printf("\n");
-	}
-
+//	double* amod = (double*) calloc(5 * dim, sizeof(double));
+//	amod = memcpy(amod, m.values, 5 * dim * sizeof(double));
+//
+//	double* r = (double*) calloc(dim, sizeof(double));
+//
+//	// dirichlet boundary
+//	for (int i = 0; i < nx; i ++) {
+//		r[i] = 10.0 * kf[i];
+//		for (int j = 0; j < dim; j ++) {
+//			if (i != j) {
+//				put(m, j, i, 0.0);
+//			} else {
+//				put(m, j, i, kf[i]);
+//			}
+//		}
+//	}
+//
+// m.values = amod;
+//
+//	// neumann boundary
+//	for (int i = dim - nx; i < dim; i ++) {
+//		r[i] = -1.8e-5;
+//	}
+//
+//	// solve the system A*x = b
+//	//
+//
+//	double* x = (double*) calloc(dim, sizeof(double));
+//	for (int i = 0; i < dim; i ++) {
+//		x[i] = 1.0;
+//	}
+//
+
+
+//	// define preconditioner
+//	double* c = (double*) calloc(dim, sizeof(double));
+//	for (int i = 0; i < dim; i ++) {
+//		c[i] = 1.0 / get(m, i, i);
+//	}
+//	int pre_ids[1] = { 0 };
+//
+//	struct mg pre = { c, pre_ids, dim, dim, 1 };
+//
+//	int flag = pcg(m, r, x, pre, 1.0e-6);
+//	printf("\n\niteration finished w/ flag : %d\n", flag);
+//
+//	printf("\nsolution:\n");
+//	for (int i = 0; i < ny; i ++) {
+//		for (int j = 0; j < nx; j ++) {
+//			printf("%g; ", x[i * nx + j]);
+//		}
+//		printf("\n");
+//	}
+//
 	// clean up
-	free(x);
-	free(r);
-	free(amod);
-	free(a);
 	free(kf);
-	free(c);
+	free(a);
 
 	return 0;
 
diff --git a/solver.c b/solver.c
index 99b3ee2..4f5d683 100644
--- a/solver.c
+++ b/solver.c
@@ -42,11 +42,11 @@
  * inout: x0  (in: initial guess, out: final result)
  * inout: tol (in: tolerance crit., out: achieved accuracy)
  */
-int pcg(struct diag a, double* b, double* x0, struct diag c, double tol)
+int pcg(struct diag* a, double* b, double* x0, struct diag* c, double tol)
 {
 
 	// params
-	int len_c = a.len_c;
+	int len_c = a->len_c;
 	int count = 0;
 	int flag  = 1;
 
@@ -54,30 +54,32 @@ int pcg(struct diag a, double* b, double* x0, struct diag c, double tol)
 	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* p  = (double*) calloc(len_c, sizeof(double));
+	double* t  = (double*) calloc(len_c, sizeof(double));
 
 	double alph = 0.0;
 	double beta = 0.0;
 
 	double eps = tol;
 
+
 	// bootstrap
-	r0 = multiply(a, x0);
+	multiply(a, x0, r0);
 
 	for (int i = 0; i < len_c; i ++) {
 		r0[i] = b[i] - r0[i];
 	}
 
-	z0 = multiply(c, r0);
+	multiply(c, r0, z0);
+	
 
 	memcpy(p, z0, len_c * sizeof(double));
 
 	// loop
 	while (count < 100) {
 
-		alph = dot(r0, z0, len_c) / dot(p, multiply(a, p), len_c);
-		t  = multiply(a, p);
+		multiply(a, p, t);
+		alph = dot(r0, z0, len_c) / dot(p, t, len_c);
 
 		for (int i = 0; i < len_c; i ++) {
 			x0[i] = x0[i] + alph * p[i];
@@ -90,7 +92,7 @@ int pcg(struct diag a, double* b, double* x0, struct diag c, double tol)
 			flag = 0;
 			break;
 		} else {
-			z1   = multiply(c, r1);
+			multiply(c, r1, z1);
 			beta = dot(z1, r1, len_c) / dot(z0, r0, len_c);
 
 			for (int i = 0; i < len_c; i ++) {
@@ -102,6 +104,10 @@ int pcg(struct diag a, double* b, double* x0, struct diag c, double tol)
 
 		}
 
+		memset(t,  0.0, len_c * sizeof(double));
+		memset(z1, 0.0, len_c * sizeof(double));
+		memset(r1, 0.0, len_c * sizeof(double));
+
 		count = count + 1;
 	}
 
diff --git a/solver.h b/solver.h
index 0744fe5..8ce4409 100644
--- a/solver.h
+++ b/solver.h
@@ -26,6 +26,6 @@
 #ifndef __SOLVER_INCLUDED__
 #define __SOLVER_INCLUDED__
 
-int pcg(struct diag, double*, double*, struct diag, double);
+int pcg(struct diag*, double*, double*, struct diag*, double);
 
 #endif
diff --git a/sysmat.c b/sysmat.c
new file mode 100644
index 0000000..478113a
--- /dev/null
+++ b/sysmat.c
@@ -0,0 +1,177 @@
+/**
+ * sysmat.c
+ * ---------------
+ *
+ * purpose : return the ``system matrix''
+ *
+ * @author ilhan oezgen, wahyd, ilhan.oezgen@wahyd.tu-berlin.de
+ * @date 6 apr 17
+ *      10 apr 17 : switch to 1d vectors
+ *
+ **/
+
+// ----------------------------------------------------------------------------
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "global.h"
+#include "sysmat.h"
+#include "linalg.h"
+
+// ----------------------------------------------------------------------------
+
+/**
+ * build system matrix
+ *
+ * inout:
+ * struct diag m : system matrix
+ **/
+void system_matrix(struct diag* m, int nx, int ny, double* val)
+{
+
+	// the system matrix is a square matrix with following dimension: 
+	int dim = nx * ny;
+
+	// loop over inner cells with 4 neighbors
+	for (int i = 1; i < ny - 1; i ++) {
+		for (int j = 1; j < nx - 1; j ++) {
+
+			int ik = i * nx + j;
+			int i1 = ik - 1;
+			int i2 = (i - 1) * nx + j;
+			int i3 = ik + 1;
+			int i4 = (i + 1) * nx + j;
+
+			double val1 = harmonic(val[ik], val[i1]);
+			double val2 = harmonic(val[ik], val[i2]);
+			double val3 = harmonic(val[ik], val[i3]);
+			double val4 = harmonic(val[ik], val[i4]);
+
+			put(m, ik, ik,   val1 + val2 + val3 + val4);
+			put(m, ik, i1, - val1);
+			put(m, ik, i2, - val2);
+			put(m, ik, i3, - val3);
+			put(m, ik, i4, - val4);
+
+		}
+	}
+
+	printf("%g\n", get(m, 4, 4));
+	printf("%g\n", get(m, 4, 3));
+	printf("%g\n", get(m, 4, 5));
+	printf("%g\n", get(m, 4, 1));
+	printf("%g\n", get(m, 4, 7));
+
+	return m;
+
+//	// boundary cells with 2 neighbors
+//
+//	// upper left corner : r = 0, c = 0;
+//	double kb1 = harmonic(val[0], val[1]);
+//	double kb2 = harmonic(val[0], val[nx]);
+//
+//	put(m, 0, 0,    kb1 + kb2);
+//	put(m, 0, 1,  - kb1);
+//	put(m, 0, nx, - kb2);
+//
+//	// upper right corner : r = 0, c = nx
+//	kb1 = harmonic(val[nx - 1], val[nx - 2]);
+//	kb2 = harmonic(val[nx - 1], val[2 * nx - 1]);
+//
+//	put(m, nx - 1, nx - 1,       kb1 + kb2);
+//	put(m, nx - 1, nx - 2,     - kb1);
+//	put(m, nx - 1, 2 * nx - 1, - kb2);
+//
+//
+//	// lower left corner : r = ny, c = 0
+//	kb1 = harmonic(val[dim - nx], val[dim - nx + 1]);
+//	kb2 = harmonic(val[dim - nx], val[dim - 2 * nx]);
+//
+//	put(m, dim - nx, dim - nx,       kb1 + kb2);
+//	put(m, dim - nx, dim - nx + 1, - kb1);
+//	put(m, dim - nx, dim - 2 * nx, - kb2);
+//
+//	// lower right corner
+//	kb1 = harmonic(val[dim - 1], val[dim - 2]);
+//	kb2 = harmonic(val[dim - 1], val[dim - 1 - nx]);
+//
+//	put(m, dim - 1, dim - 1,        kb1 + kb2);
+//	put(m, dim - 1, dim - 2,      - kb1);
+//	put(m, dim - 1, dim - 1 - nx, - kb2);
+//
+//
+//	// boundary cells with 3 neighbors
+//
+//	// additional value required
+//	double kb3 = 0.0;
+//
+//	// upper boundary (i = 0)
+//	for (int j = 1; j < nx - 1; j ++) {
+//		kb1 = harmonic(val[j], val[j - 1]);
+//		kb2 = harmonic(val[j], val[j + 1]);
+//		kb3 = harmonic(val[j], val[nx + j]);
+//
+//		put(m, j, j,        kb1 + kb2 + kb3);
+//		put(m, j, j - 1,  - kb1);
+//		put(m, j, j + 1,  - kb2);
+//		put(m, j, nx + j, - kb3);
+//
+//	}
+//
+//	// lower boundary (i = n_max)
+//	for (int j = 1; j < nx - 1; j ++) {
+//
+//		int c = dim - nx;
+//
+//		kb1 = harmonic(val[c + j], val[c + j - 1]);
+//		kb2 = harmonic(val[c + j], val[c + j + 1]);
+//		kb3 = harmonic(val[c + j], val[dim - 2 * nx + j]);
+//
+//		put(m, c + j, c + j,              kb1 + kb2 + kb3);
+//		put(m, c + j, c + j - 1,        - kb1);
+//		put(m, c + j, c + j + 1,        - kb2);
+//		put(m, c + j, dim - 2 * nx + j, - kb3);
+//
+//	}
+//
+//	// left boundary (j = 0)
+//	for (int i = 1; i < ny - 1; i ++) {
+//
+//		int c = i * nx;
+//
+//		kb1 = harmonic(val[c], val[(i - 1) * nx]);
+//		kb2 = harmonic(val[c], val[(i + 1) * nx]);
+//		kb3 = harmonic(val[c], val[i * nx + 1]);
+//
+//		put(m, c, c,              kb1 + kb2 + kb3);
+//		put(m, c, (i - 1) * nx, - kb1);
+//		put(m, c, (i + 1) * nx, - kb2);
+//		put(m, c, i * nx + 1,   - kb3);
+//
+//	}
+//
+//	// right boundary (j = n_max)
+//	for (int i = 1; i < ny - 1; i ++) {
+//
+//		int c = (i + 1) * nx - 1;
+//
+//		kb1 = harmonic(val[c], val[i * nx - 1]);
+//		kb2 = harmonic(val[c], val[(i + 2) * nx - 1]);
+//		kb3 = harmonic(val[c], val[i * nx + nx - 2]);
+//
+//		put(m, c, c,                  kb1 + kb2 + kb3);
+//		put(m, c, i * nx - 1,       - kb1);
+//		put(m, c, (i + 2) * nx - 1, - kb2);
+//		put(m, c, i * nx + nx - 2,  - kb3);
+//
+//	}
+
+	//return m;
+
+}
+
+double harmonic(double ki, double kj)
+{
+	return 2.0 * ki * kj / (ki + kj);
+}
diff --git a/system_matrix.h b/sysmat.h
index 0c034d7..ddfeac7 100644
--- a/system_matrix.h
+++ b/sysmat.h
@@ -11,15 +11,11 @@
  *
  **/
 
-#ifndef __SYSTEM_INCLUDED__
-#define __SYSTEM_INCLUDED__
+#ifndef __SYSMAT_INCLUDED__
+#define __SYSMAT_INCLUDED__
 
-double* system_matrix(int, int, double*);
+struct diag system_matrix(struct diag, int, int, double*);
 
-double avg_kf(double, double);
-
-void put_diag(double*, int*, int, int, int, double);
-
-int sanity_check(double*, int);
+double harmonic(double, double);
 
 #endif
diff --git a/system_matrix.c b/system_matrix.c
deleted file mode 100644
index f84adae..0000000
--- a/system_matrix.c
+++ /dev/null
@@ -1,182 +0,0 @@
-/**
- * sysmat.c
- * ---------------
- *
- * purpose : return the ``system matrix''
- *
- * @author ilhan oezgen, wahyd, ilhan.oezgen@wahyd.tu-berlin.de
- * @date 6 apr 17
- *      10 apr 17 : switch to 1d vectors
- *
- **/
-
-// ----------------------------------------------------------------------------
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#include "sysmat.h"
-#include "linalg.h"
-
-// ----------------------------------------------------------------------------
-
-/**
- * build system matrix
- * 
- * input:
- * nx : number of points in x direction (original system)
- * ny : number of points in y direction (original system)
- **/
-struct diag system_matrix(int nx, int ny, double* val)
-{
-
-	// the system matrix is a square matrix with following dimension:
-	int dim = nx * ny;
-
-	// store diagonal ids for each non-zero diagonal
-	int a_ids[5] = { - nx, -1, 0, 1, nx };
-
-	// giving a matrix that with shape (5, dim)
-	double* a = (double*) calloc(5 * dim, sizeof(double));
-
-	struct diag m;
-	m.values = a;
-	m.ids    = a_ids;
-	m.len_c  = dim;
-	m.len_r  = dim;
-	m.len_id = 5;
-
-	// loop over inner cells with 4 neighbors
-	for (int i = 1; i < ny - 1; i ++) {
-		for (int j = 1; j < nx - 1; j ++) {
-
-			int ik = i * nx + j;
-			int i1 = ik - 1;
-			int i2 = (i - 1) * nx + j;
-			int i3 = ik + 1;
-			int i4 = (i + 1) * nx + j;
-
-			double val1 = harmonic(val[ik], val[i1]);
-			double val2 = harmonic(val[ik], val[i2]);
-			double val3 = harmonic(val[ik], val[i3]);
-			double val4 = harmonic(val[ik], val[i4]);
-
-			put(a, ik, ik,   val1 + val2 + val3 + val4);
-			put(a, ik, i1, - val1);
-			put(a, ik, i2, - val2);
-			put(a, ik, i3, - val3);
-			put(a, ik, i4, - val4);
-
-		}
-	}
-
-	// boundary cells with 2 neighbors
-
-	// upper left corner : r = 0, c = 0;
-	double kb1 = harmonic(val[0], val[1]);
-	double kb2 = harmonic(val[0], val[nx]);
-
-	put(a, 0, 0,    kb1 + kb2);
-	put(a, 0, 1,  - kb1);
-	put(a, 0, nx, - kb2);
-
-	// upper right corner : r = 0, c = nx
-	kb1 = harmonic(val[nx - 1], val[nx - 2]);
-	kb2 = harmonic(val[nx - 1], val[2 * nx - 1]);
-
-	put(a, nx - 1, nx - 1,       kb1 + kb2);
-	put(a, nx - 1, nx - 2,     - kb1);
-	put(a, nx - 1, 2 * nx - 1, - kb2);
-
-
-	// lower left corner : r = ny, c = 0
-	kb1 = harmonic(val[dim - nx], val[dim - nx + 1]);
-	kb2 = harmonic(val[dim - nx], val[dim - 2 * nx]);
-
-	put(a, dim - nx, dim - nx,       kb1 + kb2);
-	put(a, dim - nx, dim - nx + 1, - kb1);
-	put(a, dim - nx, dim - 2 * nx, - kb2);
-
-	// lower right corner
-	kb1 = harmonic(val[dim - 1], val[dim - 2]);
-	kb2 = harmonic(val[dim - 1], val[dim - 1 - nx]);
-
-	put(a, dim - 1, dim - 1,        kb1 + kb2);
-	put(a, dim - 1, dim - 2,      - kb1);
-	put(a, dim - 1, dim - 1 - nx, - kb2);
-
-
-	// boundary cells with 3 neighbors
-
-	// additional value required
-	double kb3 = 0.0;
-
-	// upper boundary (i = 0)
-	for (int j = 1; j < nx - 1; j ++) {
-		kb1 = harmonic(val[j], val[j - 1]);
-		kb2 = harmonic(val[j], val[j + 1]);
-		kb3 = harmonic(val[j], val[nx + j]);
-
-		put(a, j, j,        kb1 + kb2 + kb3);
-		put(a, j, j - 1,  - kb1);
-		put(a, j, j + 1,  - kb2);
-		put(a, j, nx + j, - kb3);
-
-	}
-
-	// lower boundary (i = n_max)
-	for (int j = 1; j < nx - 1; j ++) {
-
-		int c = dim - nx;
-
-		kb1 = harmonic(val[c + j], val[c + j - 1]);
-		kb2 = harmonic(val[c + j], val[c + j + 1]);
-		kb3 = harmonic(val[c + j], val[dim - 2 * nx + j]);
-
-		put(a, c + j, c + j,              kb1 + kb2 + kb3);
-		put(a, c + j, c + j - 1,        - kb1);
-		put(a, c + j, c + j + 1,        - kb2);
-		put(a, c + j, dim - 2 * nx + j, - kb3);
-
-	}
-
-	// left boundary (j = 0)
-	for (int i = 1; i < ny - 1; i ++) {
-
-		int c = i * nx;
-
-		kb1 = harmonic(val[c], val[(i - 1) * nx]);
-		kb2 = harmonic(val[c], val[(i + 1) * nx]);
-		kb3 = harmonic(val[c], val[i * nx + 1]);
-
-		put(a, c, c,              kb1 + kb2 + kb3);
-		put(a, c, (i - 1) * nx, - kb1);
-		put(a, c, (i + 1) * nx, - kb2);
-		put(a, c, i * nx + 1,   - kb3);
-
-	}
-
-	// right boundary (j = n_max)
-	for (int i = 1; i < ny - 1; i ++) {
-
-		int c = (i + 1) * nx - 1;
-
-		kb1 = harmonic(val[c], val[i * nx - 1]);
-		kb2 = harmonic(val[c], val[(i + 2) * nx - 1]);
-		kb3 = harmonic(val[c], val[i * nx + nx - 2]);
-
-		put_a(a, c, c,                  kb1 + kb2 + kb3);
-		put_a(a, c, i * nx - 1,       - kb1);
-		put_a(a, c, (i + 2) * nx - 1, - kb2);
-		put_a(a, c, i * nx + nx - 2,  - kb3);
-
-	}
-
-	return a;
-
-}
-
-double harmonic(double ki, double kj)
-{
-	return 2.0 * ki * kj / (ki + kj);
-}
diff --git a/test_linalg.c b/test_linalg.c
index f7cc48c..d445750 100644
--- a/test_linalg.c
+++ b/test_linalg.c
@@ -26,6 +26,7 @@ int main(void)
 	double* v = (double*) calloc(nd, sizeof(double));
 	ids[0] = -1;
 	ids[1] =  0;
+	double* result = (double*) calloc(nd, sizeof(double));
 
 	a[2] = -12.0;
 
@@ -36,7 +37,7 @@ int main(void)
 		}
 	}
 
-	struct diag m = { a, ids, nd, nd, ni }; // initialize struct
+	struct diag m = { a, ids, nd, nd, nd, ni }; // initialize struct
 
 	double reference_m[5][5] = { { 3.0,   0.0, 0.0, 0.0, 0.0}, 
 															 { 0.0,   3.5, 0.0, 0.0, 0.0},
@@ -54,24 +55,24 @@ int main(void)
 
 	for (int i = 0; i < nd; i ++) {
 		for (int j = 0; j < nd; j ++) {
-			if (fabs(get(m, i, j) - reference_m[i][j]) < 1.0e-12) {
+			if (fabs(get(&m, i, j) - reference_m[i][j]) < 1.0e-12) {
 			printf("\ntest passed\n");
 		  }
 		}
 	}
 
-	put(m, 2, 2, -100.0);
+	put(&m, 2, 2, -100.0);
 	reference_m[2][2] = 100.0;
 
 	for (int i = 0; i < nd; i ++) {
 		for (int j = 0; j < nd; j ++) {
-			if (fabs(get(m, i, j) - reference_m[i][j]) < 1.0e-12) {
+			if (fabs(get(&m, i, j) - reference_m[i][j]) < 1.0e-12) {
 			printf("\ntest passed\n");
 		  }
 		}
 	}
 
-	double* result = multiply(m, v);
+	multiply(&m, v, result);
 
 	double reference[5] = { 15.0, 21.0, -772.0, 36.0, 45.0 };
 
@@ -87,8 +88,8 @@ int main(void)
 
 	printf("\n");
 
+	free(result);
 	free(a);
 	free(ids);
 	free(v);
-	free(result);
 }
diff --git a/test_pcg.c b/test_pcg.c
index aa32e0c..9fe8d3a 100644
--- a/test_pcg.c
+++ b/test_pcg.c
@@ -21,17 +21,17 @@ int main (void)
 	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 };
+	struct diag m  = { a, ids, 2, 2, 2, 3 };
 
 	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
+	struct diag pre = { c, cids, 2, 2, 2, 1 }; // Jacobi preconditioner
 
 	double b[3]    = { 1.0, 2.0 };
-	double x0[2]   = { 10000.0, 1000.0 };
+	double x0[2]   = { 20000.0, 10000.0 };
 
-	pcg(m, b, x0, pre, 1.0e-12);
+	pcg(&m, b, x0, &pre, 1.0e-12);
 
 	printf(">>> x0: \n");
 	for (int i = 0; i < 2; i ++) {