diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs
new file mode 100644
index 00000000..fec0e4d8
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs
@@ -0,0 +1,566 @@
+// <copyright file="ManagedLinearAlgebraProvider.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+// Copyright (c) 2009 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// The managed linear algebra provider.
+    /// </summary>
+    public class AtlasLinearAlgebraProvider : ILinearAlgebraProvider
+    {
+        #region ILinearAlgebraProvider<double> Members
+
+        /// <summary>
+        /// Adds a scaled vector to another: <c>y += alpha*x</c>.
+        /// </summary>
+        /// <param name="y">The vector to update.</param>
+        /// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
+        /// <param name="x">The vector to add to <paramref name="y"/>.</param>
+        /// <remarks>This equivalent to the AXPY BLAS routine.</remarks>
+        public void AddVectorToScaledVector(double[] y, double alpha, double[] x)
+        {
+            if (y == null)
+            {
+                throw new ArgumentNullException("y");
+            }
+
+            if (x == null)
+            {
+                throw new ArgumentNullException("x");
+            }
+
+            if (y.Length != x.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (alpha == 0.0)
+            {
+                return;
+            }
+
+            SafeNativeMethods.d_axpy(y.Length, alpha, x, y);
+        }
+
+        /// <summary>
+        /// Scales an array. Can be used to scale a vector and a matrix.
+        /// </summary>
+        /// <param name="alpha">The scalar.</param>
+        /// <param name="x">The values to scale.</param>
+        /// <remarks>This is equivalent to the SCAL BLAS routine.</remarks>
+        public void ScaleArray(double alpha, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public int QueryWorkspaceBlockSize(string methodName)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double DotProduct(double[] x, double[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double MatrixNorm(Norm norm, double[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double MatrixNorm(Norm norm, double[] matrix, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, double[] b, double beta, double[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(double[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(double[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(double[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(double[] a, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(double[] a, int[] ipiv, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(double[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(double[] r, double[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(double[] r, double[] q, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, double[] s, double[] u, double[] vt, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+
+        #region ILinearAlgebraProvider<float> Members
+
+
+        public void AddVectorToScaledVector(float[] y, float alpha, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void ScaleArray(float alpha, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float DotProduct(float[] x, float[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float MatrixNorm(Norm norm, float[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float MatrixNorm(Norm norm, float[] matrix, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, float[] b, float beta, float[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(float[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(float[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(float[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(float[] a, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(float[] a, int[] ipiv, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, float[] a, int ipiv, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int ipiv, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(float[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(float[] r, float[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(float[] r, float[] q, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, float[] q, float[] r, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, float[] s, float[] u, float[] vt, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+
+        #region ILinearAlgebraProvider<Complex> Members
+
+
+        public void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void ScaleArray(Complex alpha, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex DotProduct(Complex[] x, Complex[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex MatrixNorm(Norm norm, Complex[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex MatrixNorm(Norm norm, Complex[] matrix, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, Complex[] b, Complex beta, Complex[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(Complex[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(Complex[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(Complex[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(Complex[] a, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(Complex[] a, int[] ipiv, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, Complex[] a, int ipiv, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int ipiv, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(Complex[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(Complex[] r, Complex[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(Complex[] r, Complex[] q, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, Complex[] q, Complex[] r, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+    }
+}
\ No newline at end of file
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.tt b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.tt
new file mode 100644
index 00000000..00754c05
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.tt
@@ -0,0 +1,4 @@
+<#@ template language="C#v3.5" debug="true" #>
+<#@ output extenstion="cs" #>
+<# string library = "Atlas";#>
+<#@ include file="..\NativeAlgebraProvider.include" #>
\ No newline at end of file
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.cs
new file mode 100644
index 00000000..c7915608
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.cs
@@ -0,0 +1,55 @@
+// <copyright file="Constants.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System.Runtime.InteropServices;
+using System.Security;
+
+namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas
+{
+    [SuppressUnmanagedCodeSecurity]
+    internal static class SafeNativeMethods
+    {
+        private const string DllName = "MathNET.Numerics.ATLAS.dll";
+
+        #region BLAS
+        
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void s_axpy(int n, float alpha, float[] x, [In, Out] float[] y);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void d_axpy(int n, double alpha, double[] x, [In, Out] double[] y);
+
+        //[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        //internal static extern void c_axpy(int n, ref Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void z_axpy(int n, ref Complex alpha, Complex[] x, [In, Out] Complex[] y);
+        
+        #endregion BLAS
+	}
+}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.tt b/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.tt
new file mode 100644
index 00000000..5f17c830
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/SafeNativeMethods.tt
@@ -0,0 +1,6 @@
+<#@ template language="C#v3.5" debug="true" #>
+<#@ output extenstion="cs" #>
+<# string namespaceSuffix = "Atlas";
+   string library = "ATLAS"; 
+#>
+<#@ include file="..\SafeNativeMethods.include" #>
\ No newline at end of file
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs
new file mode 100644
index 00000000..fad34bec
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs
@@ -0,0 +1,566 @@
+// <copyright file="ManagedLinearAlgebraProvider.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+// Copyright (c) 2009 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// The managed linear algebra provider.
+    /// </summary>
+    public class MklLinearAlgebraProvider : ILinearAlgebraProvider
+    {
+        #region ILinearAlgebraProvider<double> Members
+
+        /// <summary>
+        /// Adds a scaled vector to another: <c>y += alpha*x</c>.
+        /// </summary>
+        /// <param name="y">The vector to update.</param>
+        /// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
+        /// <param name="x">The vector to add to <paramref name="y"/>.</param>
+        /// <remarks>This equivalent to the AXPY BLAS routine.</remarks>
+        public void AddVectorToScaledVector(double[] y, double alpha, double[] x)
+        {
+            if (y == null)
+            {
+                throw new ArgumentNullException("y");
+            }
+
+            if (x == null)
+            {
+                throw new ArgumentNullException("x");
+            }
+
+            if (y.Length != x.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (alpha == 0.0)
+            {
+                return;
+            }
+
+            SafeNativeMethods.d_axpy(y.Length, alpha, x, y);
+        }
+
+        /// <summary>
+        /// Scales an array. Can be used to scale a vector and a matrix.
+        /// </summary>
+        /// <param name="alpha">The scalar.</param>
+        /// <param name="x">The values to scale.</param>
+        /// <remarks>This is equivalent to the SCAL BLAS routine.</remarks>
+        public void ScaleArray(double alpha, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public int QueryWorkspaceBlockSize(string methodName)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double DotProduct(double[] x, double[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double MatrixNorm(Norm norm, double[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double MatrixNorm(Norm norm, double[] matrix, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, double[] b, double beta, double[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(double[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(double[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(double[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(double[] a, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(double[] a, int[] ipiv, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(double[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(double[] r, double[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(double[] r, double[] q, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, double[] s, double[] u, double[] vt, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+
+        #region ILinearAlgebraProvider<float> Members
+
+
+        public void AddVectorToScaledVector(float[] y, float alpha, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void ScaleArray(float alpha, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float DotProduct(float[] x, float[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float MatrixNorm(Norm norm, float[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float MatrixNorm(Norm norm, float[] matrix, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, float[] b, float beta, float[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(float[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(float[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(float[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(float[] a, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(float[] a, int[] ipiv, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, float[] a, int ipiv, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int ipiv, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(float[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(float[] r, float[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(float[] r, float[] q, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, float[] q, float[] r, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, float[] s, float[] u, float[] vt, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+
+        #region ILinearAlgebraProvider<Complex> Members
+
+
+        public void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void ScaleArray(Complex alpha, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex DotProduct(Complex[] x, Complex[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex MatrixNorm(Norm norm, Complex[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex MatrixNorm(Norm norm, Complex[] matrix, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, Complex[] b, Complex beta, Complex[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(Complex[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(Complex[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(Complex[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(Complex[] a, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(Complex[] a, int[] ipiv, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, Complex[] a, int ipiv, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int ipiv, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(Complex[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(Complex[] r, Complex[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(Complex[] r, Complex[] q, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, Complex[] q, Complex[] r, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+    }
+}
\ No newline at end of file
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.tt b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.tt
new file mode 100644
index 00000000..2eb924ad
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.tt
@@ -0,0 +1,4 @@
+<#@ template language="C#v3.5" debug="true" #>
+<#@ output extenstion="cs" #>
+<# string library = "Mkl";#>
+<#@ include file="..\NativeAlgebraProvider.include" #>
\ No newline at end of file
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs
new file mode 100644
index 00000000..b605cef1
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.cs
@@ -0,0 +1,55 @@
+// <copyright file="Constants.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System.Runtime.InteropServices;
+using System.Security;
+
+namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl
+{
+    [SuppressUnmanagedCodeSecurity]
+    internal static class SafeNativeMethods
+    {
+        private const string DllName = "MathNET.Numerics.MKL.dll";
+
+        #region BLAS
+        
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void s_axpy(int n, float alpha, float[] x, [In, Out] float[] y);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void d_axpy(int n, double alpha, double[] x, [In, Out] double[] y);
+
+        //[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        //internal static extern void c_axpy(int n, ref Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void z_axpy(int n, ref Complex alpha, Complex[] x, [In, Out] Complex[] y);
+        
+        #endregion BLAS
+	}
+}
diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.tt b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.tt
new file mode 100644
index 00000000..993a9ac7
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/SafeNativeMethods.tt
@@ -0,0 +1,6 @@
+<#@ template language="C#v3.5" debug="true" #>
+<#@ output extenstion="cs" #>
+<# string namespaceSuffix = "Mkl";
+   string library = "MKL"; 
+#>
+<#@ include file="..\SafeNativeMethods.include" #>
\ No newline at end of file
diff --git a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include
new file mode 100644
index 00000000..b35c43e6
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include
@@ -0,0 +1,566 @@
+// <copyright file="ManagedLinearAlgebraProvider.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+// Copyright (c) 2009 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#>
+{
+    using System;
+    using Properties;
+
+    /// <summary>
+    /// The managed linear algebra provider.
+    /// </summary>
+    public class <#=library#>LinearAlgebraProvider : ILinearAlgebraProvider
+    {
+        #region ILinearAlgebraProvider<double> Members
+
+        /// <summary>
+        /// Adds a scaled vector to another: <c>y += alpha*x</c>.
+        /// </summary>
+        /// <param name="y">The vector to update.</param>
+        /// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
+        /// <param name="x">The vector to add to <paramref name="y"/>.</param>
+        /// <remarks>This equivalent to the AXPY BLAS routine.</remarks>
+        public void AddVectorToScaledVector(double[] y, double alpha, double[] x)
+        {
+            if (y == null)
+            {
+                throw new ArgumentNullException("y");
+            }
+
+            if (x == null)
+            {
+                throw new ArgumentNullException("x");
+            }
+
+            if (y.Length != x.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (alpha == 0.0)
+            {
+                return;
+            }
+
+            SafeNativeMethods.d_axpy(y.Length, alpha, x, y);
+        }
+
+        /// <summary>
+        /// Scales an array. Can be used to scale a vector and a matrix.
+        /// </summary>
+        /// <param name="alpha">The scalar.</param>
+        /// <param name="x">The values to scale.</param>
+        /// <remarks>This is equivalent to the SCAL BLAS routine.</remarks>
+        public void ScaleArray(double alpha, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public int QueryWorkspaceBlockSize(string methodName)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double DotProduct(double[] x, double[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double MatrixNorm(Norm norm, double[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public double MatrixNorm(Norm norm, double[] matrix, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(double[] x, double[] y, double[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, double[] b, double beta, double[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(double[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(double[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(double[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(double[] a, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(double[] a, int[] ipiv, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(double[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, double[] a, double[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(double[] r, double[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(double[] r, double[] q, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x, double[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, double[] s, double[] u, double[] vt, double[] b, double[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+
+        #region ILinearAlgebraProvider<float> Members
+
+
+        public void AddVectorToScaledVector(float[] y, float alpha, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void ScaleArray(float alpha, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float DotProduct(float[] x, float[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float MatrixNorm(Norm norm, float[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public float MatrixNorm(Norm norm, float[] matrix, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(float[] x, float[] y, float[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, float[] b, float beta, float[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(float[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(float[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(float[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(float[] a, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(float[] a, int[] ipiv, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, float[] a, int ipiv, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int ipiv, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(float[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, float[] a, float[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(float[] r, float[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(float[] r, float[] q, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, float[] q, float[] r, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x, float[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, float[] s, float[] u, float[] vt, float[] b, float[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+
+        #region ILinearAlgebraProvider<Complex> Members
+
+
+        public void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void ScaleArray(Complex alpha, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex DotProduct(Complex[] x, Complex[] y)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void AddArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex MatrixNorm(Norm norm, Complex[] matrix)
+        {
+            throw new NotImplementedException();
+        }
+
+        public Complex MatrixNorm(Norm norm, Complex[] matrix, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiply(Complex[] x, Complex[] y, Complex[] result)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, Complex[] b, Complex beta, Complex[] c)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUFactor(Complex[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(Complex[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(Complex[] a, int[] ipiv)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverse(Complex[] a, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUInverseFactored(Complex[] a, int[] ipiv, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(int columnsOfB, Complex[] a, int ipiv, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int ipiv, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskyFactor(Complex[] a)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolve(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void CholeskySolveFactored(int columnsOfB, Complex[] a, Complex[] b)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(Complex[] r, Complex[] q)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRFactor(Complex[] r, Complex[] q, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void QRSolveFactored(int columnsOfB, Complex[] q, Complex[] r, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SinguarValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x, Complex[] work)
+        {
+            throw new NotImplementedException();
+        }
+
+        public void SvdSolveFactored(int columnsOfB, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x)
+        {
+            throw new NotImplementedException();
+        }
+
+        #endregion
+    }
+}
\ No newline at end of file
diff --git a/src/Numerics/Algorithms/LinearAlgebra/SafeNativeMethods.include b/src/Numerics/Algorithms/LinearAlgebra/SafeNativeMethods.include
new file mode 100644
index 00000000..ac6f11fc
--- /dev/null
+++ b/src/Numerics/Algorithms/LinearAlgebra/SafeNativeMethods.include
@@ -0,0 +1,55 @@
+// <copyright file="Constants.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System.Runtime.InteropServices;
+using System.Security;
+
+namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
+{
+    [SuppressUnmanagedCodeSecurity]
+    internal static class SafeNativeMethods
+    {
+        private const string DllName = "MathNET.Numerics.<#=library#>.dll";
+
+        #region BLAS
+        
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void s_axpy(int n, float alpha, float[] x, [In, Out] float[] y);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void d_axpy(int n, double alpha, double[] x, [In, Out] double[] y);
+
+        //[DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        //internal static extern void c_axpy(int n, ref Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern void z_axpy(int n, ref Complex alpha, Complex[] x, [In, Out] Complex[] y);
+        
+        #endregion BLAS
+	}
+}
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 42efec89..70730d90 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -44,9 +44,30 @@
     <Reference Include="System.Xml" />
   </ItemGroup>
   <ItemGroup>
+    <Compile Include="Algorithms\LinearAlgebra\Atlas\AtlasLinearAlgebraProvider.cs">
+      <DependentUpon>AtlasLinearAlgebraProvider.tt</DependentUpon>
+      <AutoGen>True</AutoGen>
+      <DesignTime>True</DesignTime>
+    </Compile>
+    <Compile Include="Algorithms\LinearAlgebra\Atlas\SafeNativeMethods.cs">
+      <DependentUpon>SafeNativeMethods.tt</DependentUpon>
+      <AutoGen>True</AutoGen>
+      <DesignTime>True</DesignTime>
+    </Compile>
+    <None Include="Algorithms\LinearAlgebra\NativeAlgebraProvider.include" />
     <Compile Include="Algorithms\LinearAlgebra\ILinearAlgebraProvider.cs" />
     <Compile Include="Algorithms\LinearAlgebra\ILinearAlgebraProviderOfT.cs" />
     <Compile Include="Algorithms\LinearAlgebra\ManagedLinearAlgebraProvider.cs" />
+    <Compile Include="Algorithms\LinearAlgebra\Mkl\MklLinearAlgebraProvider.cs">
+      <DependentUpon>MklLinearAlgebraProvider.tt</DependentUpon>
+      <AutoGen>True</AutoGen>
+      <DesignTime>True</DesignTime>
+    </Compile>
+    <Compile Include="Algorithms\LinearAlgebra\Mkl\SafeNativeMethods.cs">
+      <AutoGen>True</AutoGen>
+      <DesignTime>True</DesignTime>
+      <DependentUpon>SafeNativeMethods.tt</DependentUpon>
+    </Compile>
     <Compile Include="Combinatorics.cs" />
     <Compile Include="Complex.cs" />
     <Compile Include="Constants.cs" />
@@ -144,6 +165,28 @@
     <None Include="..\MathNet.Numerics.snk">
       <Link>MathNet.Numerics.snk</Link>
     </None>
+    <None Include="Algorithms\LinearAlgebra\Atlas\AtlasLinearAlgebraProvider.tt">
+      <Generator>TextTemplatingFileGenerator</Generator>
+      <LastGenOutput>AtlasLinearAlgebraProvider.cs</LastGenOutput>
+    </None>
+    <None Include="Algorithms\LinearAlgebra\Atlas\SafeNativeMethods.tt">
+      <Generator>TextTemplatingFileGenerator</Generator>
+      <LastGenOutput>SafeNativeMethods.cs</LastGenOutput>
+    </None>
+    <None Include="Algorithms\LinearAlgebra\Mkl\MklLinearAlgebraProvider.tt">
+      <Generator>TextTemplatingFileGenerator</Generator>
+      <LastGenOutput>MklLinearAlgebraProvider.cs</LastGenOutput>
+    </None>
+    <None Include="Algorithms\LinearAlgebra\Mkl\SafeNativeMethods.tt">
+      <Generator>TextTemplatingFileGenerator</Generator>
+      <LastGenOutput>SafeNativeMethods.cs</LastGenOutput>
+    </None>
+    <None Include="Algorithms\LinearAlgebra\SafeNativeMethods.include">
+      <LastGenOutput>SafeNativeMethods.cs</LastGenOutput>
+    </None>
+  </ItemGroup>
+  <ItemGroup>
+    <Service Include="{B4F97281-0DBD-4835-9ED8-7DFB966E87FF}" />
   </ItemGroup>
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it.