Release v0.0.12

This next release fixes a number of problems with Campy for a more complex example: steepest descent. This example encompasses a number of advanced capabilities of Campy and C#, which is best explained with the implementation shown below. In this example, you will note use of value types, reference types, generics, and multiple Parallel.For() calls.

using System;
using System.Text;
using System.Collections.Generic;
using System.Collections.ObjectModel;

namespace test
{
    class Program
    {
        static void Main(string[] args)
        {
            var A = new SquareMatrix(new Collection<double>() { 3, 2, 2, 6 });
            var b = new Vector(new Collection<double>() {2, -8});
            var x = new Vector(new Collection<double>() {-2, -2});
	    var r = SD.SteepestDescent(A, b, x);
            System.Console.WriteLine(r.ToString());
        }
    }

    class SquareMatrix
    {
        public int N { get; private set; }
        private List<double> data;
        public SquareMatrix(int n)
        {
            N = n;
            data = new List<double>();
            for (int i = 0; i < n*n; ++i) data.Add(0);
        }

        public SquareMatrix(Collection<double> c)
        {
            data = new List<double>(c);
            var s = Math.Sqrt(c.Count);
            N = (int)Math.Floor(s);
            if (s != (double)N)
            {
                throw new Exception("Need to provide square matrix sized initializer.");
            }
        }

        public static Vector operator *(SquareMatrix a, Vector b)
        {
            Vector result = new Vector(a.N);
            Campy.Parallel.For(result.N, i =>
            {
                for (int j = 0; j < result.N; ++j)
                    result[i] += a.data[i * result.N + j] * b[j];
            });
            return result;
        }
    }

    class Vector
    {
        public int N { get; private set; }
        private List<double> data;

        public Vector(int n)
        {
            N = n;
            data = new List<double>();
            for (int i = 0; i < n; ++i) data.Add(0);
        }

        public double this[int i]
        {
            get
            {
                return data[i];
            }
            set
            {
                data[i] = value;
            }
        }

        public Vector(Collection<double> c)
        {
            data = new List<double>(c);
            N = c.Count;
        }

        public static double operator *(Vector a, Vector b)
        {
            double result = 0;
            for (int i = 0; i < a.N; ++i) result += a[i] * b[i]; return result; } public static Vector operator *(double a, Vector b) { Vector result = new Vector(b.N); Campy.Parallel.For(b.N, i => { result[i] = a * b[i]; });
            return result;
        }

        public static Vector operator -(Vector a, Vector b)
        {
            Vector result = new Vector(a.N);
            Campy.Parallel.For(a.N, i => { result[i] = a[i] - b[i]; });
            return result;
        }

        public static Vector operator +(Vector a, Vector b)
        {
            Vector result = new Vector(a.N);
            Campy.Parallel.For(a.N, i => { result[i] = a[i] + b[i]; });
            return result;
        }

        public override string ToString()
        {
            StringBuilder sb = new StringBuilder();
            for (int i = 0; i < data.Count; ++i)
            {
                sb.Append(data[i] + " ");
            }
            return sb.ToString();
        }
    }

    class SD
    {
        public static Vector SteepestDescent(SquareMatrix A, Vector b, Vector x)
        {
            // Similar to http://ta.twi.tudelft.nl/nw/users/mmbaumann/projects/Projekte/MPI2_slides.pdf
            for (;;)
            {
                Vector r = b - A * x;
                double rr = r * r;
                double rAr = r * (A * r);
		if (Math.Abs(rAr) <= 1.0e-10) break;
                double a = (double) rr / (double) rAr;
                x = x + (a * r);
            }
            return x;
        }

        // https://www.coursera.org/learn/predictive-analytics/lecture/RhkFB/parallelizing-gradient-descent
        // "Hogwild! A lock-free approach to parallelizing stochastic gradient descent"
        // https://arxiv.org/abs/1106.5730

        // Parallelize vector and matrix operations
        // http://www.dcs.warwick.ac.uk/pmbs/pmbs14/PMBS14/Workshop_Schedule_files/8-CUDAHPCG.pdf

        // An introduction to the conjugate gradient method without the agonizing pain
        // https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf

        // https://github.com/gmarkall/cuda_cg/blob/master/gpu_solve.cu
    }
}