Friday, April 11, 2014

NativeInterop is live on NuGet

Today I released a first version of my NativeInterop package on NuGet. You can find a description of the purpose and usage of the package as well as the source code on BitBucket.

The motivation to build and release this package really evolved from two practical issues I encountered when building performance critical C# code:

  1. Creating a native, high-performance generic array data structure in C# seems to be impossible (see A minimalistic native 64 bit array ...).
  2. Reading structured binary data from a byte[] requires some ugly hacks to get both decent performance and genericity (see e.g. Reading Unmanaged Data Into Structures).
The reason for both of these issues is the fact, that C# lacks a "unmanaged" type constraint and thus you cannot express something like
static unsafe T Read<T>(T* p) {
    return *p;
}
But in F#, you can; you'd simply encode this as
open NativeInterop
let pVal = NativePtr.read ptr
where ptr is of type nativeptr<'T> and 'T is constrained to unmanaged types.

The performance offered by the NativeInterop package should be on par with non-generic unsafe C# code. The NativeInterop package also contains an implementation of NativeArray64, but this time without using inline IL. It turned out that in newer versions of the .NET framework, the AGUs are utilized correctly for the address (offset) computation (instead of emitting IMUL instructions): Calling NativePtr.set<'T>/get<'T>/set64<'T>/get64<'T> (or NativePtr.Get/Set or IntPtr.Get/Set or NativePtr.get_Item/set_Item, respectively) should all result in the generation of a single mov instruction.

Monday, April 7, 2014

A first look at RyuJIT CTP3 and SIMD (SSE2) support for .NET

Not being able to exploit today's processors SIMD processing capabilities is a major culprit when implementing high-performance (e. g. numerical) applications in C# (or any other CLI language). While there is Mono.Simd, there is no solutions for applications running on top of Microsoft's own runtime (CLR), despite popular demand ... until now!

With BUILD 2014, Microsoft released a new preview version of the next generation JIT compiler "RyuJIT" that, combined with a special SIMD library that can be installed via NuGet, supports SIMD intrinsics (only SSE2 for now, but AVX is in the works).

Finally! I couldn't wait to try out the new bits; thus I modified the C# version of my existing XRaySimulator* to make use of SSE2 by implementing a simple packet ray tracing technique, i. e. instead of tracing individual rays, this version traces bundles of 2x2 (SSE2) or 4x2 (AVX) rays. Because the rays are largely "coherent" they typically hit the same objects (cache hit rate!).

The contenders

Currently there are a total of six different variants of the XRaySimulator:
  •  "C#": This is the baseline, scalar managed implementation.
  •  "C# adj. trav.": A further optimized version that exploits the fact that once a ray is inside a volume (finite element) mesh, it must hit a face of an adjacent element (hexahedron).
  •  "C#/SSE2": Like "C#", but using 2x2 (X-)ray packets; doesn't use "adjacency traversal" due to the high branching factor
  •  "C++": A C++11 reimplementation of "C#"; I tried to stay as close as possible to "C#" while still using at least half-way decent, idiomatic C++.
  •  "C++ adj. trav.": Corresponds to "C# adj. trav."
  •  "C++/AVX": Vectorized version of "C++" using 4x2 ray bundles thanks to AVX
Note that most of these implementations are to be considered "quick-and-dirty, yet somehow working hacks..." If you don't mind the ugliness, though, you may follow the links to BitBucket and have a look at the code (Visual Studio 2013 projects; you also need the latest Roslyn CTP in order to compile the C#/SSE2 branch).

Performance analysis

So, who wins? The following figure shows the performance of the different versions in million rays per second (MRay/s) rendering an FE model consisting of 28672 hexahedral elements (344064 triangles) at a resolution of 6400 x 4800 pixels on an Intel Core i7-2600K (3.4 - 4.2 GHz) with 32 GB DDR4 RAM running under Windows 8.1 Pro:


As expected "C++/AVX" blasts away the rest of the pack with a stunning 13 MRay/s. And while "C++/AVX" delivers a speed-up of 5.2 over "C++", "C#/SSE2" only improves by a factor of 2.5 compared to "C#" and displays only insignificant performance gains over the optimized scalar version "C# adj. trav."

Now, given that SSE2 uses only 128-bit-wide vector lanes compared to AVX's generous 256 bit and the generally much more aggressive optimizer of the Visual C++ compiler, it's not exactly surprising to see an obvious performance difference between the "C++/AVX" and "C#/SSE2" case. Yet, I still would have expected the speed-up of  "C#/SSE2" to reach a value a little closer to 4x instead of 2.5. What's going on there?

According to Visual Studio's built-in profiler all of the implementations spend the majority of their time in the intersection routine of the AABB (axis-aligned bounding box) - which is a good thing, because this intersection test is very fast compared to a triangle intersection test. Thus the quality of the generated machine code for this method/function is critical for the overall performance of the renderer.

The source code of the C#/SSE2 version looks like this:
Loading ....

And here's the source for the C++/AVX version:
Loading ....

(Note: The C++ code uses a hard-coded vector lane width of 8 floats.)

Almost identical; yet, if you compare what both RyuJIT and Visual C++ make of these sources, you'll first notice that the machine code emitted by RyuJIT is much more convoluted and thus longer:
Most of that "additional stuff" that's going on in the C#/RyuJIT version seems to be related to null-pointer checks (both AABB and RayPack are classes and thus reference types). Still, I wonder if all those load/store operations are truly neccessary.

Preliminary conclusions

It seems like Microsoft has finally awakend and makes the long overdue investments in .NET performance. Thanks Google and Apple! Although RyuJIT will still require a lot of optimizations, in particular with respect to the generated SIMD code, Redmond's latest moves are promising. A next generation JIT, SIMD support, AOT compilation using the Visual C++ optimizer backend... What will come next? GPGPU support? Large arrays? A decent, modern, performant desktop UI framework? True first-class support for F#?

The future is bright!

*XRaySimulator is a visualization tool that renders X-ray-like images of finite element models. It uses a modified ray tracing algorithm to compute the energy absorption within each intersected element based on the element's material properties. A BVH (bounding volume hierarchy) is used to speed-up the intersection computation.
Details (German): http://www.uni-ulm.de/fileadmin/website_uni_ulm/uzwr/projekte/p10-2.pdf

**The C# versions of XRaySimulator on BitBucket currently don't support saving the rendered image to a file. In older versions, I used to use Tao.DevIL, but that only works on x86 and the preview releases of RyuJIT only emit x64 machine code. The C++ versions use a custom TGA output filter.

Sunday, April 25, 2010

A minimalistic native 64 bit array implementation for .NET (updated code)

If you ever felt the need to process huge amounts of data via a algorithm implemented using .NET/the CLR, you’ve surely ran into the 2^31-items-limit of the CLR’s current array implementation that only supports Int32 array indices (this also affects other collections like List<T> as those use standard arrays for storage internally).
You can try to circumvent this limitation by implementing your own array-like data-structure, either by emulating continuous storage via a collection of standard .NET-arrays (partition your data in chunks with 2^31 items each), or you can use native APIs and some evil pointer arithmetic to get maximum performance.
A while ago I tried to implement the latter approach in C#, which isn’t a big deal, only a matter of some unsafe-blocks for pointer arithmetic and a call to Marshal.AllocHGlobal for allocating memory on the unmanged heap. However, when I tried to make that custom collection into a generic one, I ran into an unsolvable problem:
public unsafe T this[long index]
{
    get
    {                
        return *((T*)pBase + index);
    }
    set
    {
        *((T*)pBase + index) = value;
    }
}

This code does not compile. The reason for that is, that there is no way to tell the C# compiler that T shall be constrained to unmanaged types.
Interestingly, F# 2.0 does feature such a constraint! This is how a minimalistic F# implementation of such an native 64 bit array could look like:
namespace NativeTools
 
#nowarn "9"
#nowarn "42"
 
open System
open System.Runtime
open Microsoft.FSharp.NativeInterop
 
module internal PointerArithmetic =
    [<CompiledName("AddIntPtrToIntPtr")>]
    [<Unverifiable>]
    let inline addNativeInt (x: nativeptr<'T>) (n: nativeint) : nativeptr<'T> = 
        (NativePtr.toNativeInt x) + n * (# "sizeof !0" type('T) : nativeint #) |> NativePtr.ofNativeInt
    
    // "reinterpret_cast<IntPtr>(x)"... EVIL!
    [<CompiledName("Int64ToIntPtr")>]
    [<Unverifiable>]
    let inline int64ToNativeint (x: int64) = (# "" x : nativeint #)
 
    [<CompiledName("AddInt64ToIntPtr")>]
    [<Unverifiable>]
    let inline addInt64 (x: nativeptr<'a>) (o: int64) : nativeptr<'a> = addNativeInt x (int64ToNativeint o)
    
[<Sealed>]
type NativeArray64<'T when 'T: unmanaged>(length: int64) =
    let itemSize: int64 = (int64)(InteropServices.Marshal.SizeOf(typeof<'T>))
    let mutable isDisposed = false
    let allocatedBytes = length * itemSize
    let blob = InteropServices.Marshal.AllocHGlobal(nativeint allocatedBytes)
    let pBlobBase: nativeptr<'T> = NativePtr.ofNativeInt blob
    let disposeLock = new Object()
 
    member this.Length = length
    member this.BaseAddress = pBlobBase
    member this.ItemSize = itemSize
    member this.IsDisposed = isDisposed
    member this.AllocatedBytes = allocatedBytes
    
    member private this.Free () =
        lock disposeLock (fun () ->
            if isDisposed
                then ()
                else InteropServices.Marshal.FreeHGlobal blob
                     isDisposed <- true
        )
           
    member this.Item
        with get (idx: int64) =
                        NativePtr.read (PointerArithmetic.addInt64 pBlobBase idx)                    
        and  set (idx: int64) (value: 'T) =
                        NativePtr.write (PointerArithmetic.addInt64 pBlobBase idx) value
        
    member private this.Items = seq {
            for i in 0L .. length - 1L do
                yield this.[i]
        }
 
    override this.Finalize () = this.Free()
    
    interface IDisposable with
        member this.Dispose () =
            GC.SuppressFinalize this
            this.Free()
 
    interface Collections.Generic.IEnumerable<'T> with
        member this.GetEnumerator () : Collections.Generic.IEnumerator<'T> =
            this.Items.GetEnumerator()
        member this.GetEnumerator () : Collections.IEnumerator =
            this.Items.GetEnumerator() :> Collections.IEnumerator

UPDATE 2010-04-25: Removed a few bugs.

You can use this data structure in your C# code like a normal array:
var length = 8L * 1024L * 1024L * 1024L;

// allocate a byte-array of 8 GiB

using(arr = new NativeTools.NativeArray64<byte>(length))

{

    arr[0] = 123;

    arr[length-1] = 222;

    Console.WriteLine("Allocated " + arr.AllocatedBytes);

}

// auto-disposed ...

Saturday, April 11, 2009

Generic .NET Math — First Benchmarks

As promised, I did some first benchmarking of different methods to implement generic arithmetic for .NET. I implemented simple vector classes using the Operator<T> class found in the MiscUtil library, another generic vector class using the approach described in my last posting (using INumeric/Calculators) and compared that to a completely non-generic vector implementation. The benchmarking method solved the following recursive problem for i = 100,000,000 on a Core 2 Quad Q9550:
sum_i := sum_(i-1) + <u_i, v> + <u_i, t>
u_i := u_(i-1) * sum_i
As the figure above shows, the non-generic version is far superior to the other implementations, especially in Long Mode (x64). The x86 version of the interface constraint approach (INumeric/Calculator) is only a tiny bit faster than MiscUtil’s Operator<T> class using runtime code generation. For x64 however, it’s about twice as fast.

References:
  1. Rüdiger Klaehn: Using generics for calculations
  2. Keith Farmer: Operator Overloading with Generics
  3. Bill Fugina: Arithmetic in Generic Classes
  4. Roger Alsing: Linq Expressions - Calculating with generics

Friday, April 10, 2009

Generic Arithmetic for .NET

To me, one of the major design flaws and sources of frustration in .NET’s generics system is the absence of suppport for some sort of common arithmetic interface for numeric types (float, double, int, etc.).

But why should anyone care about such things? Well, the problem becomes obvious as soon as you try to implement—say—a generic vector class:

(Note: I'll only show the relavant parts of the code here to avoid clutter and confusion)

   1:  // a vector in R^n
   2:  class Vector<T>
   3:  {
   4:      T[] components;
   5:    
   6:      // ctors, properties etc. ...
   7:  
   8:      // operator defintion, here: vector addition
   9:      public staticVector<T> operator +(Vector<T> lhs,
  10:                                        Vector<T> rhs)
  11:      {
  12:          Vector<T> res = new Vector<T>(lhs.Length);
  13:  
  14:          for (int i = 0; i < res.Length; ++i)
  15:          {
  16:                       // compiler error!
  17:              res[i] = lhs[i] + rhs[i];
  18:          }
  19:  
  20:          return res;
  21:      }
  22:  }

The reason for the compiler error is, that there’s simply not enough information available to the compiler to check, whether T supports op_Addition. Now, if the CLI designers had included some interface definition like

   1:  interface INumeric<T>
   2:  {
   3:      T Add(T rhs);
   4:      T Sub(T rhs);
   5:      T Mul(T rhs);
   6:      T Div(T rhs);
   7:      // etc.

... (or similar) and if they had implemented this interface for every numeric primitive type, we could simply use an interface constraint to solve our problem:

   1:  // a vector in R^n
   2:  class Vector<T> where T: INumeric<T>
   3:  {
   4:      T[] components;
   5:      
   6:      // ctors, properties ...
   7:      
   8:      public static Vector<T> operator +(Vector<T> lhs, Vector<T> rhs)
   9:      {
  10:          Vector<T> res = new Vector<T>(lhs.Length);
  11:      
  12:          for (int i = 0; i < res.Length; ++i)
  13:          {
  14:              res[i] = lhs[i].Add(rhs[i]);       // OK, T implements INumeric<T>
  15:          }
  16:          
  17:          return res;
  18:      }
  19:  }

(Note: Operators are always statically defined in .NET, hence we can only use normal prefix method calls as static methods can never be virtual, thus neither abstract nor part of an interface defintion.)

Unfortunately, such an interface does not exist. Several people already tried to circumvent this limitation. One of the first approaches uses special numerical interfaces that implement arithmetic operations for concrete types. This technique works reasonably well, besides having the disadvantage of forcing the user to explicitly provide an implementation of such an interfaces (or at least a type name of a type implementing the according interface).

More recently, a new solution came up, which uses runtime code generation by building and compiling expression trees on the fly (or using LCG for older .NET version). This indeed is a very elegant solution to the problem. Using an implementation of this concept like the one found in the MiscUtil library you can simply write something like

   1:  T sum = Operator<T>.Add(lhs, rhs);

... inside your generic math algorithms.

The appropriate "Add" method for type T is generated on-the-fly upon the first call and cached for later reuse. The drawback of this method is the need for delegate invocations, which tend to be a lot more expensive than interface calls (at least my first benchmarks indicate that, especially on x64; curiously, the effect is much less severe on x86).

Interestingly, the latest F# CTP also provides generic vector and matrix classes (namespace Microsoft.FSharp.Math in assembly FSharp.PowerPack.dll). You can instantiate e.g. matrices like that:

   1:  Vector<float> v = VectorModule.Generic.create<float>(3, 3, 0f);

This returns a 3×3 matrix of floats initialized to 0.0f. So how did they do that? As far as I understand, their approach is very similar to the first method I described (using interfaces to define "calculator" types). To avoid having the user explicitly specify the operators’ implementation they use some kind of “type registry”, mapping each known numeric type to an appropriate implementation of an interface that defines basic arithmetic operations.

I’ve tried to do a quick n’ dirty implementation of this idea in C# which resulted in something like this (I haven’t yet checked if this actually runs...):

   1:  namespace GenericMath
   2:  {
   3:      public interface INumeric<T>
   4:      {
   5:          T Add(T lhs, T rhs);
   6:          T Sub(T lhs, T rhs);
   7:          T Mul(T lhs, T rhs);
   8:          T Div(T lhs, T rhs);
   9:          T Abs(T x);
  10:          T Zero { get; }
  11:          T One { get; }
  12:      }
  13:   
  14:      class FloatNumerics: INumeric<float>
  15:      {
  16:          public float Add(float lhs, float rhs)
  17:          {
  18:              return lhs + rhs;
  19:          }
  20:   
  21:          public float Sub(float lhs, float rhs)
  22:          {
  23:              return lhs - rhs;
  24:          }
  25:   
  26:          public float Mul(float lhs, float rhs)
  27:          {
  28:              return lhs * rhs;
  29:          }
  30:   
  31:          public float Div(float lhs, float rhs)
  32:          {
  33:              return lhs / rhs;
  34:          }
  35:   
  36:          public float Abs(float x)
  37:          {
  38:              return Math.Abs(x);
  39:          }
  40:   
  41:          public float Zero
  42:          {
  43:              get { return 0f; }
  44:          }
  45:   
  46:          public float One
  47:          {
  48:              get { return 1f; }
  49:          }
  50:      }
  51:   
  52:      // implementations for other types ...
  53:   
  54:      public static class ArithmeticAssociations
  55:      {
  56:          static readonly IDictionary<Type, object> operationsDict = new Dictionary<Type, object>
  57:          {
  58:              { typeof(float), new FloatNumerics() },
  59:              { typeof(double), new DoubleNumerics() },
  60:              { typeof(int), new Int32Numerics() }
  61:          };
  62:   
  63:          static string numericIfaceName = typeof(INumeric<>).FullName;
  64:          static string numericIfaceNameShort = typeof(INumeric<>).Name;
  65:   
  66:          public static bool UnregisterType<T>()
  67:          {
  68:              return operationsDict.Remove(typeof(T));
  69:          }
  70:   
  71:          public static void RegisterType<T>(INumeric<T> arithmetics)
  72:          {
  73:              operationsDict[typeof(T)] = arithmetics;
  74:          }
  75:   
  76:          public static INumeric<T> TryGetNumericOps<T>()
  77:          {
  78:              Type t = typeof(T);
  79:              
  80:              if (!operationsDict.ContainsKey(t))
  81:              {
  82:                  throw new ArgumentException("No implementation of " 
  83:                      + numericIfaceNameShort + "<" + t.Name + "> registered.");
  84:              }
  85:   
  86:              Type opsType = operationsDict[t].GetType();
  87:              
  88:              if (opsType.GetInterface(numericIfaceName) == null)
  89:              {
  90:                  throw new ArgumentException("Arithmetic object associated with " 
  91:                      + t.Name + " does not implement " + numericIfaceNameShort);
  92:              }
  93:   
  94:              return (INumeric<T>)operationsDict[t];
  95:          }
  96:      }
  97:   
  98:      // a sample class that uses INumeric + ArithmeticAssociations
  99:      public sealed class Vector<T>
 100:      {
 101:          static INumeric<T> ops = ArithmeticAssociations.TryGetNumericOps<T>();
 102:   
 103:          T[] components;
 104:   
 105:          public T this[int i]
 106:          {
 107:              get { return components[i]; }
 108:              set { components[i] = value; }
 109:          }
 110:   
 111:          public int Length { get { return components.Length; } }
 112:   
 113:          public Vector(int size)
 114:          {
 115:              components = new T[size];
 116:          }
 117:   
 118:          public static T operator *(Vector<T> lhs, Vector<T> rhs)
 119:          {
 120:              T res = ops.Zero;
 121:   
 122:              for (int i = 0; i < lhs.Length; ++i)
 123:              {
 124:                  res = ops.Add(res, ops.Mul(lhs[i], rhs[i]));
 125:              }
 126:   
 127:              return res;
 128:          }
 129:   
 130:          public static Vector<T> operator *(T lambda, Vector<T> vec)
 131:          {
 132:              Vector<T> result = new Vector<T>(vec.Length);
 133:   
 134:              for (int i = 0; i < result.Length; ++i)
 135:              {
 136:                  result[i] = ops.Mul(vec[i], lambda);
 137:              }
 138:   
 139:              return result;
 140:          }
 141:   
 142:          // other properties, methods and operators ...
 143:      }
 144:   
 145:      // ...
 146:  }

Objects of the type Vector can now be created by this simple line:

   1:  Vector<float> v = new Vector<float>(3);

As you can see, this is totally transparent to the user (no need to specifiy a "Calculator" type parameter) and doesn’t have the potential performance problems of the delegate/Func method.

Maybe someone finds this useful and can even improve the idea. I also plan to do a performance comparision of the different methods soon.


Update 2009-04-11:
I found another nice blog entry that covers the same topic and also presents a solution similar to the one showed above.