switch to pcg

2 files changed, 38 insertions, 162 deletions

diff --git a/solver.c b/solver.c
index 58fb90f..98c60e9 100644
--- a/solver.c
+++ b/solver.c
@@ -13,11 +13,13 @@
  * upon successful return, output arguments have the following 
  * values:
  *
- *   x : approximate solution to Ax = b
+ *   x0  : approximate solution to Ax = b
+ *   tol : accuracy of the approximate solution
  *
  * @author ilhan oezgen, wahyd, ilhan.oezgen@wahyd.tu-berlin.de
  * @date 10 apr 17
  * @date 11 apr 17 : add norm function
+ * @date 28 apr 17 : add pcg solver
  *
  **/
 
@@ -26,167 +28,45 @@
 #include <string.h>
 #include <stdlib.h>
 #include <stdio.h>
-#include <math.h>
+
+#include "global.h" // struct diag
 #include "solver.h"
-#include "diagonal.h"
+#include "linalg.h" // multiply, etc.
 
 // ----------------------------------------------------------------------------
 
-
-int pcg(const double* a, const double* c, const double* b, double* x, 
-			  const int* ids, int len_ids, int len_c, int len_r, int max_iter, 
-				double tol)
+/*
+ * preconditioned conjugate gradient solver
+ * symmetric matrices only
+ * c    : preconditioner matrix
+ * 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)
 {
 
-	// bootstrap
-	double* r0 = b - multiply(a, x, ids, len_r, len_c, len_ids);
-	double* h0 = multiply(c, r0, ids, len_r, len_c, len_ids); // to do
-
-	d0 = (double*) calloc(len_c, sizeof(double));
-	memcpy(d0, h0, len_c * sizeof(double));
-
-
-}
-
-
-
-
-
-
-
-
-/**
- * solve linear system Ax = b
- *
- * input:
- *   a        : matrix
- *   b        : right hand-side
- *	 ncol     : number of columns
- *	 nrow     : number of rows
- *	 ndiag    : number of non-zero diagonals
- *   max_iter : maximum number of iterations
- *
- * inout:
- *   x        : solution vector
- *   tol      : convergence criterion / residual error
- **/
-int bicgstab(const double* a, const double* b, double* x, const int* ids, 
-             int len_ids, int len_c, int len_r, int max_iter, double tol)
-{
-
-	int dim = len_c * len_r;
-
-	double* r  = calloc(dim, sizeof(double)); // r = b - A * x
-	double* ax = calloc(dim, sizeof(double)); // A * x
-
-	double* s  = calloc(dim, sizeof(double));
-	double* t  = calloc(dim, sizeof(double));
+	// params
+	int len_c = a.len_c;
+	int len_r = a.len_r;
 
-	ax = multiply(a, x, ids, dim, dim, len_ids);
+	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* z0 = (double*) calloc(len_c, sizeof(double));
 
-	for (int i = 0; i < dim; i ++) {
-			r[i] = b[i] - ax[i];
-	}
+	double la0 = 0.0;
+	double rh0 = 0.0;
 
-	double normb = norm(b, dim);
-	if (fabs(normb) < 1.0e-15) {
-		normb = 1.0;
-	}
+	double eps = tol;
 
-	double err = norm(r, dim) / normb;
-
-	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));
-
-	for (int i = 1; i <= max_iter; i ++) {
-		
-		rho1 = dot(rtilde, r, dim);
-
-		if (fabs(rho1) < 1.0e-15) {
-			tol = err;
-			free_mem(t, s, v, p, rtilde, r, ax);
-			return 2;
-		}
-
-		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, ids, dim, dim, len_ids);
-
-		alph = rho1 / dot(rtilde, v, dim);
-
-		for (int j = 0; j < dim; 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] + alph * p[j];
-			}
-			tol = err;
-			max_iter = i;
-			free_mem(t, s, v, p, rtilde, r, ax);
-			return 0;
-		}
-
-		t = multiply(a, s, ids, dim, dim, len_ids);
-
-		ome1 = dot(t, s, dim) / dot(t, t, dim);
-
-		for (int j = 0; j < dim; j ++) {
-			x[j] = x[j] + ome1 * s[j];
-		}
-
-		for (int j = 0; j < dim; j ++) {
-			r[j] = s[j] - ome1 * t[j];
-		}
-
-		rho0 = rho1;
-		ome0 = ome1;
-
-		err = norm(r, dim) / normb;
-		if (err < tol) {
-			tol = err;
-			max_iter = i;
-			free_mem(t, s, v, p, rtilde, r, ax);
-			return 0;
-		}
-
-		printf(" >>> %d\n", i);
-		printf(" >>> err: %g\n", err);
+	// bootstrap
+	r0 = multiply(a, x0);
 
-	}
+	free(z0);
+	free(d0);
+	free(h0);
+	free(r0);
 
-	tol = err;
-	free_mem(t, s, v, p, rtilde, r, ax);
-	return 1;
-}
+	return 0;
 
-/**
- * free memory
- **/
-void free_mem(double* t, double* s, double* v, double* p, 
-		          double* rtilde, double* r, double* ax) 
-{
-	free(t);
-	free(s);
-	free(v);
-	free(p);
-	free(rtilde);
-	free(r);
-	free(ax);
 }
diff --git a/solver.h b/solver.h
index f842c40..be4fb6e 100644
--- a/solver.h
+++ b/solver.h
@@ -1,5 +1,5 @@
 /**
- * solver.h
+ * solver.c
  * ----------
  *
  * gradient-based iterative solver for linear systems of equations
@@ -13,23 +13,19 @@
  * upon successful return, output arguments have the following 
  * values:
  *
- *   x : approximate solution to Ax = b
+ *   x0  : approximate solution to Ax = b
+ *   tol : accuracy of the approximate solution
  *
  * @author ilhan oezgen, wahyd, ilhan.oezgen@wahyd.tu-berlin.de
  * @date 10 apr 17
  * @date 11 apr 17 : add norm function
- * @date 27 apr 17 : rename and add pcg solver
+ * @date 28 apr 17 : add pcg solver
  *
  **/
 
-#ifndef __BICGST_INCLUDED__
-#define __BICGST_INCLUDED__
+#ifndef __SOLVER_INCLUDED__
+#define __SOLVER_INCLUDED__
 
-int pcg(const double*);
-
-int bicgstab(const double*, const double*, double*, const int*, int, int, 
-             int, int, double);
-
-void free_mem(double*, double*, double*, double*, double*, double*, double*);
+int pcg(struct diag, double*, double*, struct diag, double);
 
 #endif