Skip to content

Instantly share code, notes, and snippets.

@taumuon
Created October 1, 2014 20:15
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save taumuon/5d4db567c57f9846971d to your computer and use it in GitHub Desktop.
SIMD Vectorisation in .NET
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace Burgers
{
class Program
{
static Vector<double> dummy;
static Program() { dummy = Vector<double>.One; }
private static double BurgersAnalytical(double t, double x, double nu)
{
return -2 * nu * (-(-8 * t + 2 * x) * Math.Exp(-Math.Pow((-4 * t + x), 2) / (4 * nu * (t + 1))) / (4 * nu * (t + 1)) - (-8 * t + 2 * x - 12.5663706143592) * Math.Exp(-Math.Pow(-4 * t + x - 6.28318530717959, 2) / (4 * nu * (t + 1))) / (4 * nu * (t + 1))) / (Math.Exp(-Math.Pow(-4 * t + x - 6.28318530717959, 2) / (4 * nu * (t + 1))) + Math.Exp(-Math.Pow(-4 * t + x, 2) / (4 * nu * (t + 1)))) + 4;
}
private const string path = @"C:\FILES\Learning\mooc\numerical methods\BurgersCode\output";
private static double[] linspace(double first, double last, int num)
{
var step = (last - first) / (double)num;
return Enumerable.Range(0, num).Select(v => (v * step) + first).ToArray();
}
private static double[] GetAnalytical(double[] x, double t, double nu)
{
double[] u = new double[x.Length];
for (int i = 0; i < x.Length; ++i)
{
u[i] = BurgersAnalytical(t, x[i], nu);
}
return u;
}
private static double[] GetCalculated(int nt, int nx, double dx, double dt, double nu, double[] u)
{
for (int tStep = 0; tStep < nt; tStep++)
{
double[] un = new double[nx];
Array.Copy(u, un, u.Length);
for (int i = 1; i < nx - 1; i++)
{
u[i] = un[i] - un[i] * dt / dx * (un[i] - un[i - 1]) + Math.Pow(nu * dt / dx, 2.0) *
(un[i + 1] - 2 * un[i] + un[i - 1]);
}
u[0] = un[0] - un[0] * dt / dx * (un[0] - un[nx - 1]) + Math.Pow(nu * dt / dx, 2.0) *
(un[1] - 2 * un[0] + un[nx - 1]);
u[nx - 1] = un[nx - 1] - un[nx - 1] * dt / dx * (un[nx - 1] - un[nx - 2]) + Math.Pow(nu * dt / dx, 2.0) *
(un[0] - 2 * un[nx - 1] + un[nx - 2]);
}
return u;
}
// Reduce new array allocation and copying, ping-pong between them
private static double[] GetCalculated2(int nt, int nx, double dx, double dt, double nu, double[] initial)
{
double[] u = new double[nx];
double[] un = new double[nx];
Array.Copy(initial, un, un.Length);
for (int tStep = 0; tStep < nt; tStep++)
{
for (int i = 1; i < nx - 1; i++)
{
u[i] = un[i] - un[i] * dt / dx * (un[i] - un[i - 1]) + Math.Pow(nu * dt / dx, 2.0) *
(un[i + 1] - 2 * un[i] + un[i - 1]);
}
u[0] = un[0] - un[0] * dt / dx * (un[0] - un[nx - 1]) + Math.Pow(nu * dt / dx, 2.0) *
(un[1] - 2 * un[0] + un[nx - 1]);
u[nx - 1] = un[nx - 1] - un[nx - 1] * dt / dx * (un[nx - 1] - un[nx - 2]) + Math.Pow(nu * dt / dx, 2.0) *
(un[0] - 2 * un[nx - 1] + un[nx - 2]);
double[] swap = u;
u = un;
un = swap;
}
return un;
}
// Pull calculation of (nu * dt / dx)^2 out into a variable
private static double[] GetCalculated3(int nt, int nx, double dx, double dt, double nu, double[] initial)
{
double[] u = new double[nx];
double[] un = new double[nx];
Array.Copy(initial, un, un.Length);
double factor = Math.Pow(nu * dt / dx, 2.0);
for (int tStep = 0; tStep < nt; tStep++)
{
for (int i = 1; i < nx - 1; i++)
{
u[i] = un[i] - un[i] * dt / dx * (un[i] - un[i - 1]) + factor *
(un[i + 1] - 2 * un[i] + un[i - 1]);
}
u[0] = un[0] - un[0] * dt / dx * (un[0] - un[nx - 1]) + factor *
(un[1] - 2 * un[0] + un[nx - 1]);
u[nx - 1] = un[nx - 1] - un[nx - 1] * dt / dx * (un[nx - 1] - un[nx - 2]) + factor *
(un[0] - 2 * un[nx - 1] + un[nx - 2]);
double[] swap = u;
u = un;
un = swap;
}
return un;
}
private static double[] GetCalculated4(int nt, int nx, double dx, double dt, double nu, double[] initial)
{
var nx2 = nx + (Vector<double>.Length - (nx % Vector<double>.Length));
double[] u = new double[nx2];
double[] un = new double[nx2];
Array.Copy(initial, un, initial.Length);
double factor = Math.Pow(nu * dt / dx, 2.0);
for (int tStep = 0; tStep < nt; tStep++)
{
for (int i = 1; i < nx2 - 1; i += Vector<double>.Length)
{
var vectorIn0 = new Vector<double>(un, i);
var vectorInPrev = new Vector<double>(un, i - 1);
var vectorInNext = new Vector<double>(un, i + 1);
var vectorOut = vectorIn0 - vectorIn0 * (dt / dx) * (vectorIn0 - vectorInPrev) + factor *
(vectorInNext - 2.0 * vectorIn0 + vectorInPrev);
vectorOut.CopyTo(u, i);
}
//for (int i = 1; i < nx - 1; i++)
//{
// u[i] = un[i] - un[i] * dt / dx * (un[i] - un[i - 1]) + factor *
// (un[i + 1] - 2 * un[i] + un[i - 1]);
//}
u[0] = un[0] - un[0] * dt / dx * (un[0] - un[nx - 1]) + factor *
(un[1] - 2 * un[0] + un[nx - 1]);
u[nx - 1] = un[nx - 1] - un[nx - 1] * dt / dx * (un[nx - 1] - un[nx - 2]) + factor *
(un[0] - 2 * un[nx - 1] + un[nx - 2]);
double[] swap = u;
u = un;
un = swap;
}
return un;
}
static void Main(string[] args)
{
if (!VectorMath.IsHardwareAccelerated)
{
Console.WriteLine("Not hardware accelerated!");
}
else
{
Console.WriteLine("Vector<double>.Length: " + Vector<double>.Length);
}
int nx = 10001;
int nt = 1000;
double dx = 2.0 * Math.PI / (nx - 1.0);
double nu = 0.07;
double dt = dx * nu;
double[] x = linspace(0.0, 2.0 * Math.PI, nx);
double[] initial = GetAnalytical(x, 0.0, nu);
// double[] justToWarmupJit = GetCalculated(1, nx, dx, dt, nu, initial);
double[] justToWarmupJit = GetCalculated4(1, nx, dx, dt, nu, initial);
var stopwatch = new System.Diagnostics.Stopwatch();
stopwatch.Start();
// double[] calculated = GetCalculated(nt, nx, dx, dt, nu, initial);
double[] calculated = GetCalculated4(nt, nx, dx, dt, nu, initial);
stopwatch.Stop();
double[] analytical = GetAnalytical(x, nt * dt, nu);
var fileName = System.IO.Path.Combine(path, "out.dat");
using (var writer = new StreamWriter(fileName))
{
writer.WriteLine(string.Format("#x\tanalytical\tcalculated"));
for (int i = 0; i < nx; ++i)
{
writer.WriteLine(string.Format("{0}\t{1}\t{2}", x[i], analytical[i], calculated[i]));
}
}
Console.WriteLine("Elapsed: " + stopwatch.ElapsedMilliseconds);
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment