Add direct sparse solver

src/Benchmark/Program.cs

@@ -16,6 +16,7 @@ namespace Benchmark
                         typeof(Transforms.FFT),
                         typeof(LinearAlgebra.DenseMatrixProduct),
                         typeof(LinearAlgebra.DenseVector),
+                        typeof(SparseSolver.DirectSparseSolver),
                     });
 
             switcher.Run(args);

src/Benchmark/SparseSolver/DirectSparseSolver.cs

@@ -0,0 +1,284 @@
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Configs;
+using BenchmarkDotNet.Environments;
+using BenchmarkDotNet.Jobs;
+using MathNet.Numerics;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Providers.Common.Mkl;
+using MathNet.Numerics.Providers.SparseSolver;
+using System;
+using System.Collections.Generic;
+
+namespace Benchmark.SparseSolver
+{
+    [Config(typeof(Config))]
+    public class DirectSparseSolver
+    {
+        class Config : ManualConfig
+        {
+            public Config()
+            {
+                Add(Job.Clr.With(Platform.X64).With(Jit.RyuJit));
+                Add(Job.Clr.With(Platform.X86).With(Jit.LegacyJit));
+#if !NET461
+                Add(Job.Core.With(Platform.X64).With(Jit.RyuJit));
+#endif
+            }
+        }
+
+        public enum ProviderId
+        {
+            NativeMKL,
+        }
+
+        [Params(32, 128, 1024)]
+        public int N { get; set; }
+
+        [Params(ProviderId.NativeMKL)]
+        public ProviderId Provider { get; set; }
+
+        int rowCount;
+        int columnCount;
+        int valueCount;
+        double[] values;
+        int[] rowPointers;
+        int[] columnIndices;
+        double[] rhs;
+
+        [GlobalSetup]
+        public void GlobalSetup()
+        {
+            switch (Provider)
+            {
+                case ProviderId.NativeMKL:
+                    Control.UseNativeMKL(MklConsistency.Auto, MklPrecision.Double, MklAccuracy.High);
+                    break;
+            }
+
+            var domain = new Domain(N);
+            domain.DefineProblem();
+
+            // Kmatrix is symmetric so, we need only upper triangle entries.         
+            var storage = domain.Kmatrix.UpperTriangle().Storage as SparseCompressedRowMatrixStorage<double>;
+            rowCount = storage.RowCount;
+            columnCount = storage.ColumnCount;
+            valueCount = storage.ValueCount;
+            values = storage.Values;
+            rowPointers = storage.RowPointers;
+            columnIndices = storage.ColumnIndices;
+
+            rhs = domain.Rhs.ToArray();
+        }
+
+        [Benchmark(OperationsPerInvoke = 1)]
+        public double[] SolveProblem()
+        {
+            double[] solution = new double[rowCount];
+            SparseSolverControl.Provider.Solve(DssMatrixStructure.Symmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, rhs, solution);
+            return solution;
+        }
+
+        #region Finite element method to solve Poisson's equation
+
+        class Node
+        {
+            public int ID;
+            public double X;
+
+            public Node(int id, double x)
+            {
+                ID = id;
+                X = x;
+            }
+        }
+
+        class Element
+        {
+            public int ID;
+            public Node[] Nodes;
+            public double Alpha;
+            public double Gamma;
+            public Matrix<double> Kmatrix;
+            public Vector<double> Rhs;
+
+            public Element(int id, Node node1, Node node2, double alpha, double gamma)
+            {
+                ID = id;
+                Nodes = new[] { node1, node2 };
+                Alpha = alpha;
+                Gamma = gamma;
+            }
+
+            public void ComputeMatrices()
+            {
+                var length = Nodes[1].X - Nodes[0].X; // length
+
+                Kmatrix = Matrix<double>.Build.Dense(2, 2);
+                Kmatrix[0, 0] = Alpha / length;
+                Kmatrix[0, 1] = -Alpha / length;
+                Kmatrix[1, 0] = -Alpha / length;
+                Kmatrix[1, 1] = Alpha / length;
+
+                Rhs = Vector<double>.Build.Dense(2);
+                Rhs[0] = -Gamma * length / 2.0;
+                Rhs[1] = -Gamma * length / 2.0;
+            }
+        }
+
+        class Domain
+        {
+            // Length of the domain in [m]
+            public double Length;
+
+            // Dielectric constant of the domain
+            public double Permittivity;
+
+            // Charge density in [C/m^3]
+            public double ChargeDensity;
+
+            public double VoltageAtLeft;
+            public double VoltageAtRight;
+
+            Node[] nodes;
+            Element[] elements;
+            List<Tuple<Node, double>> boundaries;
+
+            public Matrix<double> Kmatrix;
+            public Vector<double> Rhs;
+
+            public Domain(int elementCount = 4, double length = 0.08, double relativePermittivity = 1, double chargeDensity = 1E-8)
+            {
+                // Boundary value problem: 
+                //
+                // [0] ------ [1] ------ ... ------ [N]
+                //
+                // each element is characterized by an electron charge density and a dielectric constant
+                //
+                // V at node1 = 1V 
+                // V at node5 = 0V (ground)
+
+                Length = length;
+
+                Permittivity = Constants.ElectricPermittivity * relativePermittivity;
+                ChargeDensity = chargeDensity;
+
+                // Create nodes and elements 
+                nodes = new Node[elementCount + 1];
+                elements = new Element[elementCount];
+
+                var dx = Length / elements.Length;
+                for (int i = 0; i < nodes.Length; i++)
+                {
+                    nodes[i] = new Node(i, dx * i);
+                }
+                for (int i = 0; i < elements.Length; i++)
+                {
+                    elements[i] = new Element(i, nodes[i], nodes[i + 1], Permittivity, ChargeDensity);
+                }
+
+                // Initialization of the global K matrix and right-hand side vector
+                Kmatrix = Matrix<double>.Build.Sparse(nodes.Length, nodes.Length);
+                Rhs = Vector<double>.Build.Dense(nodes.Length);
+            }
+
+            public void DefineProblem(double Va = 1.0, double Vb = 0.0)
+            {
+                VoltageAtLeft = Va; // Boundary condition at the leftmost node
+                VoltageAtRight = Vb; // Boundary condition at the rightmost node
+
+                // Apply Dirichlet boundary conditions
+                boundaries = new List<Tuple<Node, double>>();
+                foreach (var node in nodes)
+                {
+                    if (node.X == 0)
+                    {
+                        boundaries.Add(new Tuple<Node, double>(node, VoltageAtLeft));
+                    }
+                    else if (node.X == Length)
+                    {
+                        boundaries.Add(new Tuple<Node, double>(node, VoltageAtRight));
+                    }
+                }
+
+                // Form the element matrices and assemble to the global matrix
+                foreach (var element in elements)
+                {
+                    element.ComputeMatrices();
+
+                    // Assemble element matrix into the global K matrix
+                    for (int i = 0; i < element.Nodes.Length; i++)
+                    {
+                        var row = element.Nodes[i].ID;
+                        for (int j = 0; j < element.Nodes.Length; j++)
+                        {
+                            var col = element.Nodes[j].ID;
+                            Kmatrix[row, col] += element.Kmatrix[i, j];
+                        }
+                        Rhs[row] += element.Rhs[i];
+                    }
+                }
+
+                // Imposition of Dirichlet boundary conditions
+                foreach (var boundary in boundaries)
+                {
+                    var node = boundary.Item1;
+                    var i = node.ID;
+                    var val = boundary.Item2;
+
+                    for (int j = 0; j < nodes.Length; j++)
+                    {
+                        if (nodes[j].ID != i)
+                            Rhs[j] = Rhs[j] - Kmatrix[j, i] * val;
+                    }
+
+                    Kmatrix.SetColumn(i, new double[Kmatrix.RowCount]);
+                    Kmatrix.SetRow(i, new double[Kmatrix.ColumnCount]);
+                    Kmatrix[i, i] = 1.0;
+                    Rhs[i] = val;
+                }
+            }
+
+            public double[] SolveProblem()
+            {
+                // Kmatrix is symmetric so, we need only upper triangle entries.
+                var storage = Kmatrix.UpperTriangle().Storage as SparseCompressedRowMatrixStorage<double>;
+                var rowCount = storage.RowCount;
+                var columnCount = storage.ColumnCount;
+                var valueCount = storage.ValueCount;
+                var values = storage.Values;
+                var rowPointers = storage.RowPointers;
+                var columnIndices = storage.ColumnIndices;
+
+                var rhs = Rhs.ToArray();
+                var solution = new double[rowCount];
+
+                SparseSolverControl.Provider.Solve(DssMatrixStructure.Symmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                    rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                    1, rhs, solution);
+
+                return solution;
+            }
+
+            public double[] GetExactSolution()
+            {
+                // Poisson's equation : ∇(ε∇V) = ρ
+                // solution : V(x) = ρ/ε/2*x^2 - (ρ/ε/2*d + (Va - Vb)/d)*x + Va, where d = length
+
+                double[] Vexact = new double[nodes.Length];
+                var factor = ChargeDensity / Permittivity * 0.5;
+                for (int i = 0; i < nodes.Length; i++)
+                {
+                    var x = nodes[i].X;
+                    Vexact[i] = factor * x * x - factor * Length * x - (VoltageAtLeft - VoltageAtRight) / Length * x + VoltageAtLeft;
+                }
+
+                return Vexact;
+            }
+        }
+
+        #endregion
+    }
+}

src/NativeProviders/MKL/capabilities.cpp

@@ -84,6 +84,10 @@ extern "C" {
 		case 384: return 1; // basic FFT (major - breaking)
 		case 385: return 0; // basic FFT (minor - non-breaking)
 
+        // SPARSE SOLVER
+        case 512: return 1;
+        case 513: return 0;
+
 		default: return 0; // unknown or not supported
 
 		}

src/NativeProviders/MKL/dss.c

@@ -0,0 +1,194 @@
+#include "wrapper_common.h"
+#include "dss.h"
+
+#if __cplusplus
+extern "C" {
+#endif
+
+    // Notes: zero-based indexing is used for rowIdx[] and colPtr[]. 
+
+    DLLEXPORT dss_int s_dss_solve(const dss_int matrixStructure, const dss_int matrixType, const dss_int systemType,
+        const dss_int nRows, const dss_int nCols, const dss_int nnz, const dss_int rowIdx[], const dss_int colPtr[], const float values[],
+        const dss_int nRhs, const float rhsValues[], float solValues[])
+    {
+        _MKL_DSS_HANDLE_t handle;
+        dss_int error;
+
+        dss_int opt = MKL_DSS_MSG_LVL_WARNING + MKL_DSS_TERM_LVL_ERROR + MKL_DSS_ZERO_BASED_INDEXING + MKL_DSS_AUTO_ORDER + MKL_DSS_SINGLE_PRECISION;
+        if (systemType) opt += MKL_DSS_TRANSPOSE_SOLVE; // solve a transposed system, A'x = b
+        
+        // Initialize the solver
+        error = dss_create(handle, opt);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Define the non-zero structure of the matrix
+        dss_int sym = (matrixStructure == 0)
+            ? MKL_DSS_SYMMETRIC_STRUCTURE
+            : (matrixStructure == 1)
+                ? MKL_DSS_SYMMETRIC
+                : MKL_DSS_NON_SYMMETRIC;
+        error = dss_define_structure(handle, sym, rowIdx, nRows, nCols, colPtr, nnz);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Reorder the matrix
+        error = dss_reorder(handle, opt, 0);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Factor the matrix
+        dss_int type = (matrixType == 0)
+            ? MKL_DSS_POSITIVE_DEFINITE
+            : MKL_DSS_INDEFINITE;
+        error = dss_factor_real(handle, type, values);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Get the solution vector
+        error = dss_solve_real(handle, opt, rhsValues, nRhs, solValues);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Deallocate solver storage
+        error = dss_delete(handle, opt);
+        return error;
+    }    
+
+    DLLEXPORT dss_int d_dss_solve(const dss_int matrixStructure, const dss_int matrixType, const dss_int systemType,
+        const dss_int nRows, const dss_int nCols, const dss_int nnz, const dss_int rowIdx[], const dss_int colPtr[], const double values[],
+        const dss_int nRhs, const double rhsValues[], double solValues[])
+    {
+        _MKL_DSS_HANDLE_t handle;
+        dss_int error;
+
+        dss_int opt = MKL_DSS_MSG_LVL_WARNING + MKL_DSS_TERM_LVL_ERROR + MKL_DSS_ZERO_BASED_INDEXING + MKL_DSS_AUTO_ORDER;
+        if (systemType) opt += MKL_DSS_TRANSPOSE_SOLVE; // solve a transposed system, A'x = b
+
+        // Initialize the solver
+        error = dss_create(handle, opt);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Define the non-zero structure of the matrix
+        dss_int sym = (matrixStructure == 0)
+            ? MKL_DSS_SYMMETRIC_STRUCTURE
+            : (matrixStructure == 1)
+                ? MKL_DSS_SYMMETRIC
+                : MKL_DSS_NON_SYMMETRIC;
+        error = dss_define_structure(handle, sym, rowIdx, nRows, nCols, colPtr, nnz);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Reorder the matrix
+        error = dss_reorder(handle, opt, 0);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Factor the matrix
+        dss_int type = (matrixType == 0)
+            ? MKL_DSS_POSITIVE_DEFINITE
+            : MKL_DSS_INDEFINITE;
+        error = dss_factor_real(handle, type, values);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Get the solution vector
+        error = dss_solve_real(handle, opt, rhsValues, nRhs, solValues);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Deallocate solver storage
+        error = dss_delete(handle, opt);
+        return error;
+    }
+
+    DLLEXPORT dss_int c_dss_solve(const dss_int matrixStructure, const dss_int matrixType, const dss_int systemType,
+        const dss_int nRows, const int nCols, const int nnz, const dss_int const rowIdx[], const dss_int colPtr[], const dss_complex_float values[],
+        const dss_int nRhs, const dss_complex_float rhsValues[], dss_complex_float solValues[])
+    {
+        _MKL_DSS_HANDLE_t handle;
+        dss_int error;
+
+        dss_int opt = MKL_DSS_MSG_LVL_WARNING + MKL_DSS_TERM_LVL_ERROR + MKL_DSS_ZERO_BASED_INDEXING + MKL_DSS_AUTO_ORDER + MKL_DSS_SINGLE_PRECISION;
+        if (systemType == 1) opt += MKL_DSS_CONJUGATE_SOLVE; // solve a conjugate transposed system, A¢Óx = b
+        else if(systemType == 2) opt += MKL_DSS_TRANSPOSE_SOLVE; // solve a transposed system, A'x = b
+
+        // Initialize the solver
+        error = dss_create(handle, opt);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Define the non-zero structure of the matrix
+        dss_int sym = (matrixStructure == 0)
+            ? MKL_DSS_SYMMETRIC_STRUCTURE_COMPLEX
+            : (matrixStructure == 1)
+                ? MKL_DSS_SYMMETRIC_COMPLEX
+                : MKL_DSS_NON_SYMMETRIC_COMPLEX;
+        error = dss_define_structure(handle, sym, rowIdx, nRows, nCols, colPtr, nnz);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Reorder the matrix
+        error = dss_reorder(handle, opt, 0);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Factor the matrix
+        dss_int type = (matrixType == 0)
+            ? MKL_DSS_POSITIVE_DEFINITE
+            : (matrixType == 1)
+                ? MKL_DSS_INDEFINITE
+                : (matrixType == 2)
+                    ? MKL_DSS_HERMITIAN_POSITIVE_DEFINITE
+                    : MKL_DSS_HERMITIAN_INDEFINITE;
+        error = dss_factor_complex(handle, type, values);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Get the solution vector
+        error = dss_solve_real(handle, opt, rhsValues, nRhs, solValues);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Deallocate solver storage
+        error = dss_delete(handle, opt);
+        return error;
+    }
+
+    DLLEXPORT dss_int z_dss_solve(const dss_int matrixStructure, const dss_int matrixType, const dss_int systemType,
+        const dss_int nRows, const dss_int nCols, const dss_int nnz, const dss_int rowIdx[], const dss_int colPtr[], const dss_complex_double values[],
+        const dss_int nRhs, const dss_complex_double rhsValues[], dss_complex_double solValues[])
+    {
+        _MKL_DSS_HANDLE_t handle;
+        dss_int error;
+
+        dss_int opt = MKL_DSS_MSG_LVL_WARNING + MKL_DSS_TERM_LVL_ERROR + MKL_DSS_ZERO_BASED_INDEXING + MKL_DSS_AUTO_ORDER;
+        if (systemType == 1) opt += MKL_DSS_CONJUGATE_SOLVE; // solve a conjugate transposed system, A¢Óx = b
+        else if (systemType == 2) opt += MKL_DSS_TRANSPOSE_SOLVE; // solve a transposed system, A'x = b
+
+        // Initialize the solver
+        error = dss_create(handle, opt);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Define the non-zero structure of the matrix
+        dss_int sym = (matrixStructure == 0)
+            ? MKL_DSS_SYMMETRIC_STRUCTURE_COMPLEX
+            : (matrixStructure == 1)
+                ? MKL_DSS_SYMMETRIC_COMPLEX
+                : MKL_DSS_NON_SYMMETRIC_COMPLEX;
+        error = dss_define_structure(handle, sym, rowIdx, nRows, nCols, colPtr, nnz);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Reorder the matrix
+        error = dss_reorder(handle, opt, 0);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Factor the matrix
+        dss_int type = (matrixType == 0)
+            ? MKL_DSS_POSITIVE_DEFINITE
+            : (matrixType == 1)
+                ? MKL_DSS_INDEFINITE
+                : (matrixType == 2)
+                    ? MKL_DSS_HERMITIAN_POSITIVE_DEFINITE
+                    : MKL_DSS_HERMITIAN_INDEFINITE;
+        error = dss_factor_complex(handle, type, values);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Get the solution vector
+        error = dss_solve_real(handle, opt, rhsValues, nRhs, solValues);
+        if (error != MKL_DSS_SUCCESS) return error;
+
+        // Deallocate solver storage
+        error = dss_delete(handle, opt);
+        return error;
+    }  
+
+#if __cplusplus
+}
+#endif

src/NativeProviders/MKL/dss.h

@@ -0,0 +1,9 @@
+#pragma once
+
+#include "mkl_dss.h"
+#include "mkl_types.h"
+
+#define dss_int             MKL_INT
+#define dss_complex_float   MKL_Complex8
+#define dss_complex_double  MKL_Complex16
+

src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj

@@ -194,6 +194,7 @@
     <ClCompile Include="..\..\Common\lapack.cpp" />
     <ClCompile Include="..\..\Common\WindowsDLL.cpp" />
     <ClCompile Include="..\..\MKL\capabilities.cpp" />
+    <ClCompile Include="..\..\MKL\dss.c" />
     <ClCompile Include="..\..\MKL\memory.c" />
     <ClCompile Include="..\..\MKL\vector_functions.c" />
     <ClCompile Include="..\..\MKL\fft.cpp" />
@@ -204,10 +205,11 @@
   <ItemGroup>
     <ClInclude Include="..\..\Common\lapack_common.h" />
     <ClInclude Include="..\..\MKL\blas.h" />
+    <ClInclude Include="..\..\MKL\dss.h" />
     <ClInclude Include="..\..\MKL\lapack.h" />
     <ClInclude Include="..\..\MKL\resource.h" />
   </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
   </ImportGroup>
-</Project>
+</Project>

src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters

@@ -36,6 +36,9 @@
     <ClCompile Include="..\..\MKL\fft.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="..\..\MKL\dss.c">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="..\..\MKL\resource.rc">
@@ -55,6 +58,9 @@
     <ClInclude Include="..\..\MKL\blas.h">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="..\..\MKL\dss.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <None Include="$(MSBuildThisFileDirectory)icon.png">

src/Numerics.Tests/Providers/SparseSolver/double/SparseSolverProviderTests.cs

@@ -0,0 +1,585 @@
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Providers.SparseSolver;
+using NUnit.Framework;
+using System;
+using System.Collections.Generic;
+
+namespace MathNet.Numerics.UnitTests.Providers.SparseSolver.Double
+{
+
+#if NATIVE
+#if MKL
+
+    /// <summary>
+    /// Base class for sparse solver provider tests.
+    /// </summary>
+    [TestFixture, Category("LAProvider")]
+    public class SparseSolverProviderTests
+    {
+        readonly double[] _b4 = { 1.0, 2.0, 3.0, 4.0};
+        readonly double[] _b5 = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+        
+        /// <summary>
+        /// Test matrix to use.
+        /// </summary>
+        readonly IDictionary<string, SparseMatrix> _matrices = new Dictionary<string, SparseMatrix>
+        {
+            {"SymmetricPositiveDefinite5x5", (SparseMatrix)Matrix<double>.Build.SparseOfColumnArrays(new [] {9.0, 1.5, 6.0, 0.75, 3.0}, new [] {1.5, 0.5, 0.0, 0.0, 0.0}, new [] {6.0, 0.0, 12.0, 0.0, 0.0 }, new [] { 0.75, 0.0, 0.0, 0.625, 0.0}, new [] {3.0, 0.0, 0.0, 0.0, 16.0})},
+            {"Triangle5x5", (SparseMatrix)Matrix<double>.Build.SparseOfColumnArrays(new [] {1.0, 0.0, 0.0, 0.0, 0.0}, new [] {5.0, 2.0, 0.0, 0.0, 0.0}, new [] {0.0, 8.0, 3.0, 0.0, 0.0 }, new [] { 0.0, 0.0, 9.0, 4.0, 0.0}, new [] {0.0, 0.0, 0.0, 10.0, 5.0})},
+            {"Square4x4", (SparseMatrix)Matrix<double>.Build.SparseOfColumnArrays(new [] {1.0, 1.0, 1.0, 2.0 },new [] {2.0, 0.0, 0.0, 2.0 },new [] {0.0, 0.0, 2.0, 1.0 },new [] {4.0, 1.0, 1.0, 0.0 })},
+        };
+
+        /// <summary>
+        /// Can solve Ax=b using direct sparse solver.
+        /// </summary>
+        [Test]
+        public void CanSolveSymmetricPositiveDefiniteMatrix()
+        {
+            var A = _matrices["SymmetricPositiveDefinite5x5"].UpperTriangle();
+
+            var csr = A.Storage as SparseCompressedRowMatrixStorage<double>;
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var xactual = new double[rowCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Symmetric, DssMatrixType.PositiveDefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, _b5, xactual);
+
+            Assert.That(error, Is.EqualTo(DssStatus.MKL_DSS_SUCCESS));
+
+            var xtrue = new double[] { -979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0 };
+
+            for (int i = 0; i < xtrue.Length; i++)
+                AssertHelpers.AlmostEqualRelative(xtrue[i], xactual[i], 13);
+        }
+
+        /// <summary>
+        /// Can solve Ax=b using direct sparse solver.
+        /// </summary>
+        [Test]
+        public void CanSolveUpperTriangularMatrix()
+        {
+            var A = _matrices["Triangle5x5"].UpperTriangle();
+
+            var csr = A.Storage as SparseCompressedRowMatrixStorage<double>;
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var xactual = new double[rowCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, _b5, xactual);
+
+            Assert.That(error, Is.EqualTo(DssStatus.MKL_DSS_SUCCESS));
+
+            var xtrue = new double[] { 106.0, -21.0, 5.5, -1.5, 1.0 };
+
+            for (int i = 0; i < xtrue.Length; i++)
+                AssertHelpers.AlmostEqualRelative(xtrue[i], xactual[i], 13);
+        }
+
+        /// <summary>
+        /// Can solve Ax=b using direct sparse solver.
+        /// </summary>
+        [Test]
+        public void CanSolveSquareMatrix()
+        {
+            var A = _matrices["Square4x4"];
+
+            var csr = A.Storage as SparseCompressedRowMatrixStorage<double>;
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var xactual = new double[rowCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, _b4, xactual);
+
+            Assert.That(error, Is.EqualTo(DssStatus.MKL_DSS_SUCCESS));
+
+            var xtrue = new double[] { 2.1, -0.35, 0.5, -0.1 };
+
+            for (int i = 0; i < xtrue.Length; i++)
+                AssertHelpers.AlmostEqualRelative(xtrue[i], xactual[i], 10);
+        }
+
+        /// <summary>
+        /// Can inverse A by using AX = I.
+        /// </summary>
+        [Test]
+        public void CanInverseSquareMatrix()
+        {
+            var A = _matrices["SymmetricPositiveDefinite5x5"];
+            var Atr = A.UpperTriangle();
+
+            var csr = Atr.Storage as SparseCompressedRowMatrixStorage<double>;
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var Identity = Matrix<double>.Build.DenseIdentity(columnCount);
+            var b = Identity.ToColumnMajorArray();
+            var Xactual = new double[rowCount * columnCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Symmetric, DssMatrixType.PositiveDefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                Identity.ColumnCount, b, Xactual);
+
+            Assert.That(error, Is.EqualTo(DssStatus.MKL_DSS_SUCCESS));
+
+            var Ainverse_actual = Matrix<double>.Build.SparseOfColumnMajor(rowCount, columnCount, Xactual);
+
+            var Ainverse_expected = Matrix<double>.Build.DenseOfColumnArrays(
+                new[] { 80.0/3.0, -80.0, -40.0/3.0, -32.0, -5.0 },
+                new[] { -80.0,    242.0,      40.0,  96.0, 15.0 },
+                new[] { -40.0/3.0, 40.0,      6.75,  16.0,  2.5 },
+                new[] { -32.0,     96.0,      16.0,  40.0,  6.0 },
+                new[] { -5.0,      15.0,       2.5,   6.0,  1.0});
+
+            for (int i = 0; i < Ainverse_actual.RowCount; i++)
+                for (int j = 0; j < Ainverse_actual.ColumnCount; j++)
+                    AssertHelpers.AlmostEqualRelative(Ainverse_actual[i, j], Ainverse_expected[i, j], 10);
+        }
+
+        /// <summary>
+        /// Can solve 1D boundary problem.
+        /// </summary>
+        [TestCase(4, 8001, 10, 1000)]
+        [TestCase(40, 8001, 0.1, 100)]
+        [TestCase(400, 8001, 0.001, 10)]
+        [TestCase(4000, 8001, 0.00001, 1)]
+        public void CanSolvePoissonEquation(int elementCount, int interpolationCount, double errV, double errE)
+        {
+            var domain = new Domain(elementCount, 0.08);
+            domain.DefineBoundaries(1.0, 0.0);
+            domain.DefineProblem();
+            domain.SolveProblem();
+
+            var xGrid = Generate.LinearSpaced(interpolationCount, 0.0, domain.Length);
+
+            var interpolated = domain.InterpolateAt(xGrid);
+            var Vactual = interpolated.Item1; // interpolated electric potential
+            var Eactual = interpolated.Item2; // interpolated electric field
+
+            var exact = domain.GetExactSolution(xGrid);
+            var Vexpected = exact.Item1; // expected electric potential
+            var Eexpected = exact.Item2; // expected electric field
+
+            var Vdiff = Vector<double>.Build.Dense(Vactual.Length, (i) => Vexpected[i] - Vactual[i]);
+            var Vnorm2 = Vdiff.L2Norm();
+            Assert.LessOrEqual(Vnorm2, errV);
+
+            var Ediff = Vector<double>.Build.Dense(Eactual.Length, (i) => Eexpected[i] - Eactual[i]);
+            var Enorm2 = Ediff.L2Norm();
+            Assert.LessOrEqual(Enorm2, errE);
+        }
+
+        #region Finite element method to solve Poisson's equation
+
+        class Node
+        {
+            public int ID;
+            public double X;
+            public double PrimaryValue; // electric potential at the node
+            public double SecondaryValue; // electric field at the node
+
+            public Node(double x)
+            {
+                X = x;
+            }
+        }
+
+        class Element // Linear element
+        {
+            public int ID;
+            public Node[] Nodes;
+            public double Alpha;
+            public double Gamma;
+            public Matrix<double> Kmatrix;
+            public Vector<double> Rhs;
+
+            public Element(Node node1, Node node2, double alpha, double gamma)
+            {
+                Nodes = new[] { node1, node2 };
+                Alpha = alpha;
+                Gamma = gamma;
+            }
+
+            public void ComputeMatrices()
+            {
+                // The master equation:
+                //   ∇(α∇V) + γ = 0
+                //
+                // Let's define a transformation from u to x
+                //    x = a * u + b
+                //
+                // for Node1, u = -1 gives -a + b = x1
+                // for Node2, u = +1 gives a + b = x2
+                //
+                // So, 
+                //    x = (x2 - x1)/2*u + (x2 + x1)/2
+                //    u = 2*(x - x1)/(x2 - x1) - 1
+                //
+                // This gives 
+                //    dx = l/2*du where l = x2 - x1 = x21 is length of the line
+                //
+                // Interpolation function, Ni(u) = ai + bi*u for i = 1, 2
+                //    N1(-1) = a1 - b1 = 1
+                //    N1(+1) = a1 + b1 = 0 -> a1 = 1/2, b1 = -1/2
+                //    N2(-1) = a2 - b2 = 0
+                //    N2(+1) = a2 + b2 = 1 -> a2 = 1/2, b2 = 1/2
+                // so, 
+                //    N1 = (1 - u) / 2
+                //    N2 = (1 + u) / 2
+                //
+                // Using the chain rule of differentiation, 
+                //    ∂Ni(x)/∂u = ∂Ni/∂x ∂x/∂u
+                //
+                // Here,
+                //    J = ∂x/∂u = x21/2 = l/2
+                //
+                //    |J| = l/2  where l is the length of the line
+                //
+                // Therefore, we can get
+                //    ∂N1/∂x = 2/l ∂N1/∂u = - 1/l
+                //    ∂N2/∂x = 2/l ∂N2/∂u = 1/l
+                //-----------------------------------------------------------------------------
+                // By using V(u) = Σ Vi*Ni(u) and ω = Ni for i = 1, 2 
+                // the weak form of the master equation is given as
+                //
+                //    ∫ ω [d/dx(α dV/dx) + γ] dx = 0
+                //
+                //    K v = f + p
+                //
+                // where
+                //    Kij = ∫ (dNi/dx) α (dNj/dx) dx
+                //    fi  = ∫ Ni γ dx
+                //    pi  = -Ni(x2)D(x2) + Ni(x1)D(x1)  where D(x) = -α dV(x)/dx
+                //-----------------------------------------------------------------------------
+                //    Kij = (l/2) ∫ α (∂Ni/∂x) (∂Nj/∂x) du, where α = const.
+                //        = α*(∂Ni/∂x)*(∂Nj/∂x)*l
+                //    K11 = α/l
+                //    K12 = M21 = - α/l
+                //    K22 = α/l
+                //-----------------------------------------------------------------------------
+                //    fi = γ*(l/2) ∫ Ni du, where γ = const.
+                //       = γ*(l/2)
+                //    f1 = γ*l/2
+                //    f2 = γ*l/2
+                //-----------------------------------------------------------------------------
+                //    p1 = -N1(x2)D(x2) + N1(x1)D(x1) = D(x1) = D1
+                //    p2 = -N2(x2)D(x2) + N2(x1)D(x1) = -D(x2) = -D2
+                // For a sufficiently large number of elements in the domain, we can ignore p.
+
+                var x21 = Nodes[1].X - Nodes[0].X; // element length
+
+                Kmatrix = Matrix<double>.Build.Dense(2, 2);
+                Kmatrix[0, 0] = Alpha / x21;
+                Kmatrix[0, 1] = -Alpha / x21;
+                Kmatrix[1, 0] = -Alpha / x21;
+                Kmatrix[1, 1] = Alpha / x21;
+
+                Rhs = Vector<double>.Build.Dense(2);
+                Rhs[0] = Gamma * x21 / 2.0;
+                Rhs[1] = Gamma * x21 / 2.0;
+            }
+
+            public bool Contains(Node point)
+            {
+                // A point P can be described with a line A-B
+                //    P = A + s*(B - A)
+                //    s = PA/BA
+                // if (s >= 0 && s <= 1), then P is inside the line A-B.
+
+                var s = (point.X - Nodes[0].X) / (Nodes[1].X - Nodes[0].X);
+                return s >= 0.0 && (1.0 - s) >= 0.0;
+            }
+
+            public void InterpolateAt(Node point)
+            {
+                // V(u) can be described with interpolation functions
+                //    V(u) = V1*N1(u) + V2*N2(u)
+                // where 
+                //    N1 = (1 - u) / 2
+                //    N2 = (1 + u) / 2
+                // The transformation from x to u,
+                //    u = 2*(x - x1)/(x2 - x1) - 1
+                // gives
+                //    V(x) = V1*(x2 - x)/(x2 - x1) + V2*(x - x1)/(x2 - x1)
+
+                var V1 = Nodes[0].PrimaryValue;
+                var V2 = Nodes[1].PrimaryValue;
+
+                var x21 = Nodes[1].X - Nodes[0].X;
+                var xp1 = point.X - Nodes[0].X;
+                var u = 2.0 * xp1 / x21 - 1.0;
+
+                var Vx = V1 * (1 - u) * 0.5 + V2 * (1 + u) * 0.5;
+
+                // Electric field, 
+                //    E = -∇V
+                // where
+                //    ∇V = [ (∂/∂x) ∑ViNi ]
+                //       = [ (∂/∂x)(V1*N1 + V2*N2) ]
+                //       = [ V1*(∂N1/∂x) + V2*(∂N2/∂x) ]
+                //
+                // The derivatives of Ni are
+                //    ∂N1/∂x = 2/l ∂N1/∂u = - 1/l where l = x21
+                //    ∂N2/∂x = 2/l ∂N2/∂u = 1/l
+                //
+                // Therefore,
+                //    E = (V1 - V2) / x21
+
+                var Ex = (V1 - V2) / x21;
+
+                point.PrimaryValue = Vx;
+                point.SecondaryValue = Ex;
+            }
+        }
+
+        class Domain
+        {
+            // Length of the domain in [m]
+            public double Length;
+
+            // Dielectric constant of the domain
+            public double Permittivity;
+
+            // Charge density in [C/m^3]
+            public double ChargeDensity;
+
+            // Boundary conditions
+            public double VoltageAtLeft;
+            public double VoltageAtRight;
+            public List<Tuple<Node, double>> Boundaries;
+
+            public Node[] Nodes;
+            public Element[] Elements;
+
+            public Matrix<double> Kmatrix; 
+            public Vector<double> Rhs;
+
+            public Domain(int elementCount = 4, double length = 0.08, double relativePermittivity = 1, double chargeDensity = 1E-8)
+            {
+                Length = length;
+
+                Permittivity = Constants.ElectricPermittivity * relativePermittivity;
+                ChargeDensity = chargeDensity;
+
+                // Create nodes and elements 
+                Nodes = new Node[elementCount + 1];
+                Elements = new Element[elementCount];
+
+                var dx = length / (double)elementCount;
+
+                // nodes:    0   1   2   3 ...    n+1
+                // elements:   0   1   2 ...    n
+                //           +---+---+---+ ... ---+
+                // x axis:   0                    length
+
+                for (int i = 0; i < Nodes.Length; i++)
+                {
+                    Nodes[i] = new Node(dx * i) { ID = i };
+                }
+                for (int i = 0; i < Elements.Length; i++)
+                {
+                    Elements[i] = new Element(Nodes[i], Nodes[i + 1], Permittivity, -ChargeDensity) { ID = i };
+                }
+
+                // Initialization of the global K matrix and right-hand side vector
+                // We know the K matrix is a symmetric matrix, so we will only handle the upper triangular parts.
+                Kmatrix = Matrix<double>.Build.Sparse(Nodes.Length, Nodes.Length);
+                Rhs = Vector<double>.Build.Dense(Nodes.Length);
+            }
+
+            public void DefineBoundaries(double Vleft, double Vright)
+            {
+                VoltageAtLeft = Vleft;
+                VoltageAtRight = Vright;
+
+                // Dirichlet boundary conditions
+                Boundaries = new List<Tuple<Node, double>>();
+                foreach (var node in Nodes)
+                {
+                    if (node.X == 0d)
+                    {
+                        Boundaries.Add(new Tuple<Node, double>(node, Vleft));
+                    }
+                    else if (node.X == Length)
+                    {
+                        Boundaries.Add(new Tuple<Node, double>(node, Vright));
+                    }
+                }
+            }
+
+            public void DefineProblem()
+            {
+                // Form the element matrices and assemble to the global matrix
+                foreach (var element in Elements)
+                {
+                    element.ComputeMatrices();
+
+                    // Assemble element matrix into the global K matrix
+                    for (int i = 0; i < element.Nodes.Length; i++)
+                    {
+                        var row = element.Nodes[i].ID;
+                        for (int j = 0; j < element.Nodes.Length; j++)
+                        {
+                            var col = element.Nodes[j].ID;
+                            if (row <= col) // only upper triangular parts are handled
+                            {
+                                Kmatrix[row, col] += element.Kmatrix[i, j];
+                            }
+                        }
+                        Rhs[row] += element.Rhs[i];
+                    }
+                }
+
+                // Imposition of Dirichlet boundary conditions
+                //
+                // If xn is given, i.e. x2 = b0, then the linear equations, 
+                //    [ K11 K12 K13 K14 ][ x1 ] = [ b1 ]
+                //    [ K21 K22 K23 K24 ][ x2 ]   [ b2 ]
+                //    [ K31 K32 K33 K34 ][ x3 ]   [ b3 ]
+                //    [ K41 K42 K43 K44 ][ x4 ]   [ b4 ]
+                // can be changed to
+                //    [ K11  0  K13 K14 ][ x1 ] = [ b1 - K12*b0 ]
+                //    [  0   1   0   0  ][ x2 ]   [ b0          ]
+                //    [ K31  0  K33 K34 ][ x3 ]   [ b3 - K32*b0 ]
+                //    [ K41  0  K43 K44 ][ x4 ]   [ b4 - K42*b0 ]
+
+                foreach (var boundary in Boundaries)
+                {
+                    var node = boundary.Item1;
+                    var i = node.ID;
+                    var val = boundary.Item2;
+
+                    for (int j = 0; j < Nodes.Length; j++)
+                    {
+                        if (Nodes[j].ID != i)
+                        {
+                            Rhs[j] -= (j <= i)
+                                ? Kmatrix[j, i] * val
+                                : Kmatrix[i, j] * val; // Kmatrix has only upper triangular parts
+                        }
+                    }
+
+                    var storage = Kmatrix.Storage as SparseCompressedRowMatrixStorage<double>;
+                    storage.MapIndexedInplace(
+                        (row, col, x) => (row == i && col == i)
+                                            ? 1d
+                                            : (row == i) || (col == i)
+                                                ? 0d
+                                                : x,
+                        Zeros.AllowSkip);
+
+                    Rhs[i] = val;
+                }
+            }
+
+            public void SolveProblem()
+            {
+                // Note that Kmatrix is actually an upper triangular, but considered as a symmetric.
+                var storage = Kmatrix.Storage as SparseCompressedRowMatrixStorage<double>;
+                var rowCount = storage.RowCount;
+                var columnCount = storage.ColumnCount;
+                var valueCount = storage.ValueCount;
+                var values = storage.Values;
+                var rowPointers = storage.RowPointers;
+                var columnIndices = storage.ColumnIndices;
+
+                var rhs = Rhs.ToArray();
+                var solution = new double[rowCount];
+
+                SparseSolverControl.Provider.Solve(DssMatrixStructure.Symmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                    rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                    1, rhs, solution);
+
+                for (int i = 0; i < solution.Length; i++)
+                {
+                    Nodes[i].PrimaryValue = solution[i];
+                }
+            }
+
+            public Tuple<double[], double[]> InterpolateAt(double[] xGrid)
+            {
+                var Vactual = new double[xGrid.Length];
+                var Eactual = new double[xGrid.Length];
+
+                for (int i = 0; i < xGrid.Length; i++)
+                {
+                    var point = new Node(xGrid[i]);
+                    for (int j = 0; j < Elements.Length; j++)
+                    {
+                        if (Elements[j].Contains(point))
+                        {
+                            Elements[j].InterpolateAt(point);
+                            Vactual[i] = point.PrimaryValue;
+                            Eactual[i] = point.SecondaryValue;
+                            break;
+                        }
+                    }
+                }
+
+                return new Tuple<double[], double[]>(Vactual, Eactual);
+            }
+
+            public Tuple<double[], double[]> GetExactSolution(double[] xGrid)
+            {
+                // Poisson's equation: ∇(ε∇V) = ρ
+                // Solution:
+                //    V(x) = ρ/ε/2*x^2 - (ρ/ε/2*d + (Va - Vb)/d)*x + Va, where d = length
+                //    E(x) = -∇V = ρ/ε*x - (ρ/ε/2*d + (Va - Vb)/d)
+
+                double[] Vexact = new double[xGrid.Length]; // electric potentials
+                double[] Eexact = new double[xGrid.Length]; // electric fields
+                var factor = ChargeDensity / Permittivity * 0.5; // ρ/ε/2
+
+                for (int i = 0; i < xGrid.Length; i++)
+                {
+                    var x = xGrid[i];
+                    Vexact[i] = factor * x * x - factor * Length * x - (VoltageAtLeft - VoltageAtRight) / Length * x + VoltageAtLeft;
+                    Eexact[i] = -2d * factor * x + factor * Length + (VoltageAtLeft - VoltageAtRight) / Length;
+                }
+
+                return new Tuple<double[], double[]>(Vexact, Eexact);
+            }
+        }
+
+        #endregion
+    }
+
+#endif
+#endif
+
+}
+

src/Numerics/Control.cs

@@ -33,6 +33,7 @@ using System.Reflection;
 using System.Runtime.InteropServices;
 using System.Text;
 using System.Threading.Tasks;
+using MathNet.Numerics.Providers.SparseSolver;
 using MathNet.Numerics.Providers.FourierTransform;
 using MathNet.Numerics.Providers.LinearAlgebra;
 
@@ -70,12 +71,14 @@ namespace MathNet.Numerics
         {
             LinearAlgebraControl.UseManaged();
             FourierTransformControl.UseManaged();
+            SparseSolverControl.UseManaged();
         }
 
         public static void UseManagedReference()
         {
             LinearAlgebraControl.UseManagedReference();
             FourierTransformControl.UseManaged();
+            SparseSolverControl.UseManaged();
         }
 
         /// <summary>
@@ -86,6 +89,7 @@ namespace MathNet.Numerics
         {
             LinearAlgebraControl.UseDefault();
             FourierTransformControl.UseDefault();
+            SparseSolverControl.UseDefault();
         }
 
         /// <summary>
@@ -95,6 +99,7 @@ namespace MathNet.Numerics
         {
             LinearAlgebraControl.UseBest();
             FourierTransformControl.UseBest();
+            SparseSolverControl.UseBest();
         }
 
 #if NATIVE
@@ -107,6 +112,7 @@ namespace MathNet.Numerics
         {
             LinearAlgebraControl.UseNativeMKL();
             FourierTransformControl.UseNativeMKL();
+            SparseSolverControl.UseNativeMKL();
         }
 
         /// <summary>
@@ -121,6 +127,7 @@ namespace MathNet.Numerics
         {
             LinearAlgebraControl.UseNativeMKL(consistency, precision, accuracy);
             FourierTransformControl.UseNativeMKL();
+            SparseSolverControl.UseNativeMKL();
         }
 
         /// <summary>
@@ -134,7 +141,8 @@ namespace MathNet.Numerics
         {
             bool linearAlgebra = LinearAlgebraControl.TryUseNativeMKL();
             bool fourierTransform = FourierTransformControl.TryUseNativeMKL();
-            return linearAlgebra || fourierTransform;
+            bool directSparseSolver = SparseSolverControl.TryUseNativeMKL();
+            return linearAlgebra || fourierTransform || directSparseSolver;
         }
 
         /// <summary>
@@ -200,6 +208,7 @@ namespace MathNet.Numerics
         {
             LinearAlgebraControl.FreeResources();
             FourierTransformControl.FreeResources();
+            SparseSolverControl.FreeResources();
         }
 
         public static void UseSingleThread()
@@ -209,6 +218,7 @@ namespace MathNet.Numerics
 
             LinearAlgebraControl.Provider.InitializeVerify();
             FourierTransformControl.Provider.InitializeVerify();
+            SparseSolverControl.Provider.InitializeVerify();
         }
 
         public static void UseMultiThreading()
@@ -218,6 +228,7 @@ namespace MathNet.Numerics
 
             LinearAlgebraControl.Provider.InitializeVerify();
             FourierTransformControl.Provider.InitializeVerify();
+            SparseSolverControl.Provider.InitializeVerify();
         }
 
         /// <summary>
@@ -247,6 +258,7 @@ namespace MathNet.Numerics
                 _nativeProviderHintPath = value;
                 LinearAlgebraControl.HintPath = value;
                 FourierTransformControl.HintPath = value;
+                SparseSolverControl.HintPath = value;
             }
         }
 
@@ -265,6 +277,7 @@ namespace MathNet.Numerics
                 // Reinitialize providers:
                 LinearAlgebraControl.Provider.InitializeVerify();
                 FourierTransformControl.Provider.InitializeVerify();
+                SparseSolverControl.Provider.InitializeVerify();
             }
         }
 
@@ -323,6 +336,7 @@ namespace MathNet.Numerics
 #endif
             sb.AppendLine($"Linear Algebra Provider: {LinearAlgebraControl.Provider}");
             sb.AppendLine($"Fourier Transform Provider: {FourierTransformControl.Provider}");
+            sb.AppendLine($"Sparse Solver Provider: {SparseSolverControl.Provider}");
             sb.AppendLine($"Max Degree of Parallelism: {MaxDegreeOfParallelism}");
             sb.AppendLine($"Parallelize Elements: {ParallelizeElements}");
             sb.AppendLine($"Parallelize Order: {ParallelizeOrder}");

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -291,6 +291,75 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             MapInplace(x => x, Zeros.AllowSkip);
         }
 
+        /// <summary>
+        /// Fill zeros explicitly on the diagonal entries as required by the Intel MKL direct sparse solver. 
+        /// </summary>
+        public void PopulateExplicitZerosOnDiagonal()
+        {
+            var delta = 0; // number of missing diagonal entries
+
+            for (int row = 0; row < RowCount; row++)
+            {
+                var found = false;
+                for (int j = RowPointers[row]; j < RowPointers[row + 1]; j++)
+                {
+                    if (ColumnIndices[j] == row)
+                    {
+                        found = true;
+                        break;
+                    }
+                }
+                if (!found) delta++;
+            }
+
+            if (delta > 0)
+            {
+                var size = Values.Length + delta;
+                if (size > int.MaxValue)
+                {
+                    throw new NotSupportedException(Resources.TooManyElements);
+                }
+
+                var newRowPointers = new int[RowCount + 1];
+                var newColumnIndices = new int[size];
+                var newValues = new T[size];
+
+                delta = 0;
+                for (int row = 0; row < RowCount; row++)
+                {
+                    var found = false;
+                    for (int j = RowPointers[row]; j < RowPointers[row + 1]; j++)
+                    {
+                        newColumnIndices[j + delta] = ColumnIndices[j];
+                        newValues[j + delta] = Values[j];
+                        if (ColumnIndices[j] == row)
+                        {
+                            found = true;
+                        }
+                    }
+                    if (!found)
+                    {
+                        var start = RowPointers[row] + delta;
+                        var end = RowPointers[row + 1] + delta;
+                        var count = end - start + 1;
+
+                        newColumnIndices[end] = row;
+                        newValues[end] = Zero;
+
+                        // Ordering may be not necessary
+                        Sorting.Sort(newColumnIndices, newValues, start, count);
+
+                        delta++;
+                    }
+                    newRowPointers[row + 1] = RowPointers[row + 1] + delta;
+                }
+
+                Array.Copy(newRowPointers, RowPointers, RowCount + 1);
+                ColumnIndices = newColumnIndices;
+                Values = newValues;
+            }
+        }
+
         /// <summary>
         /// Returns a hash code for this instance.
         /// </summary>

src/Numerics/Providers/Common/Mkl/MklProviderCapabilities.cs

@@ -56,7 +56,9 @@ namespace MathNet.Numerics.Providers.Common.Mkl
         VectorFunctionsMajor = 130,
         VectorFunctionsMinor = 131,
         FourierTransformMajor = 384,
-        FourierTransformMinor = 385
+        FourierTransformMinor = 385,
+        SparseSolverMajor = 512,
+        SparseSolverMinor = 513,
     }
 }
 

src/Numerics/Providers/Common/Mkl/SafeNativeMethods.cs

@@ -430,6 +430,30 @@ namespace MathNet.Numerics.Providers.Common.Mkl
 
         #endregion FFT
 
+        #region Direct Sparse Solver
+
+        [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int s_dss_solve(int matrixStructure, int matrixType, int systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, float[] values,
+            int nRhs, [In, Out] float[] rhs, [In, Out] float[] solution);
+
+        [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int d_dss_solve(int matrixStructure, int matrixType, int systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, double[] values,
+            int nRhs, [In, Out] double[] rhs, [In, Out] double[] solution);
+
+        [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int c_dss_solve(int matrixStructure, int matrixType, int systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, Complex32[] values,
+            int nRhs, [In, Out] Complex32[] rhs, [In, Out] Complex32[] solution);
+
+        [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int z_dss_solve(int matrixStructure, int matrixType, int systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, Complex[] values,
+            int nRhs, [In, Out] Complex[] rhs, [In, Out] Complex[] solution);
+
+        #endregion Direct Sparse Solver
+
         // ReSharper restore InconsistentNaming
     }
 }

src/Numerics/Providers/Experimental.cs

@@ -0,0 +1,439 @@
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Storage;
+using MathNet.Numerics.Properties;
+using MathNet.Numerics.Providers.SparseSolver;
+using System;
+
+namespace MathNet.Numerics
+{
+    using Complex = System.Numerics.Complex;
+
+    public static class Experimental
+    {
+        public static DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, float[] values,
+            int nRhs, float[] rhs, float[] solution)
+        {
+            return SparseSolverControl.Provider.Solve(matrixStructure, matrixType, systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+        }
+
+        public static DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, double[] values,
+            int nRhs, double[] rhs, double[] solution)
+        {
+            return SparseSolverControl.Provider.Solve(matrixStructure, matrixType, systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+        }
+
+        public static DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, Complex32[] values,
+            int nRhs, Complex32[] rhs, Complex32[] solution)
+        {
+            return SparseSolverControl.Provider.Solve(matrixStructure, matrixType, systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+        }
+
+        public static DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, Complex[] values,
+            int nRhs, Complex[] rhs, Complex[] solution)
+        {
+            return SparseSolverControl.Provider.Solve(matrixStructure, matrixType, systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+        }
+
+        // solve A x = b
+        // The symmetricity or definiteness of A is not checked.
+
+        public static DssStatus Solve(this Matrix<float> matrix, Vector<float> input, Vector<float> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+            if (result.Count != input.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+            if (input.Count != matrix.RowCount)
+            {
+                throw LinearAlgebra.Single.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<float>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, for example, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var rhs = input.ToArray();
+            var solution = new float[rowCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result.SetValues(solution);
+
+            return error;
+        }
+
+        public static DssStatus Solve(this Matrix<double> matrix, Vector<double> input, Vector<double> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+
+            if (result.Count != input.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (input.Count != matrix.RowCount)
+            {
+                throw LinearAlgebra.Double.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<double>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values; 
+            var rowPointers = csr.RowPointers; 
+            var columnIndices = csr.ColumnIndices;
+
+            var rhs = input.ToArray();
+            var solution = new double[rowCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result.SetValues(solution);
+
+            return error;
+        }
+
+        public static DssStatus Solve(this Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+
+            if (result.Count != input.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (input.Count != matrix.RowCount)
+            {
+                throw MathNet.Numerics.LinearAlgebra.Complex32.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<Complex32>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var rhs = input.ToArray();
+            var solution = new Complex32[rowCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result.SetValues(solution);
+
+            return error;
+        }
+        
+        public static DssStatus Solve(this Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+
+            if (result.Count != input.Count)
+            {
+                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
+            }
+
+            if (input.Count != matrix.RowCount)
+            {
+                throw LinearAlgebra.Complex.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<Complex>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var rhs = input.ToArray();
+            var solution = new Complex[rowCount];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                1, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result.SetValues(solution);
+
+            return error;
+        }
+
+        // Solve A X = B
+        // The symmetricity or definiteness of A is not checked.
+
+        public static DssStatus Solve(this Matrix<float> matrix, Matrix<float> input, Matrix<float> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+            if (result.ColumnCount != input.ColumnCount || result.RowCount != input.RowCount)
+            {
+                throw LinearAlgebra.Single.Matrix.DimensionsDontMatch<ArgumentException>(input, result);
+            }
+            if (input.RowCount != matrix.RowCount)
+            {
+                throw LinearAlgebra.Single.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<float>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, for example, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var nRhs = input.ColumnCount;
+            var rhs = new float[rowCount * nRhs];
+            Array.Copy(input.ToColumnMajorArray(), rhs, rhs.Length);
+
+            var solution = new float[rowCount * nRhs];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result = Matrix<float>.Build.DenseOfColumnMajor(rowCount, nRhs, solution);
+
+            return error;
+        }
+
+        public static DssStatus Solve(this Matrix<double> matrix, Matrix<double> input, Matrix<double> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+            if (result.ColumnCount != input.ColumnCount || result.RowCount != input.RowCount)
+            {
+                throw LinearAlgebra.Double.Matrix.DimensionsDontMatch<ArgumentException>(input, result);
+            }
+            if (input.RowCount != matrix.RowCount)
+            {
+                throw LinearAlgebra.Double.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<double>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var nRhs = input.ColumnCount;
+            var rhs = new double[rowCount * nRhs];
+            Array.Copy(input.ToColumnMajorArray(), rhs, rhs.Length);
+
+            var solution = new double[rowCount * nRhs];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result = Matrix<double>.Build.DenseOfColumnMajor(rowCount, nRhs, solution);
+
+            return error;
+        }
+
+        public static DssStatus Solve(this Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+            if (result.ColumnCount != input.ColumnCount || result.RowCount != input.RowCount)
+            {
+                throw LinearAlgebra.Complex32.Matrix.DimensionsDontMatch<ArgumentException>(input, result);
+            }
+            if (input.RowCount != matrix.RowCount)
+            {
+                throw LinearAlgebra.Complex32.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<Complex32>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var nRhs = input.ColumnCount;
+            var rhs = new Complex32[rowCount * nRhs];
+            Array.Copy(input.ToColumnMajorArray(), rhs, rhs.Length);
+
+            var solution = new Complex32[rowCount * nRhs];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result = Matrix<Complex32>.Build.DenseOfColumnMajor(rowCount, nRhs, solution);
+
+            return error;
+        }
+
+        public static DssStatus Solve(this Matrix<Complex> matrix, Matrix<Complex> input, Matrix<Complex> result)
+        {
+            if (matrix.RowCount != matrix.ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare, nameof(matrix));
+            }
+            if (result.ColumnCount != input.ColumnCount || result.RowCount != input.RowCount)
+            {
+                throw LinearAlgebra.Complex.Matrix.DimensionsDontMatch<ArgumentException>(input, result);
+            }
+            if (input.RowCount != matrix.RowCount)
+            {
+                throw LinearAlgebra.Complex.Matrix.DimensionsDontMatch<ArgumentException>(input, matrix);
+            }
+
+            var csr = matrix.Storage as SparseCompressedRowMatrixStorage<Complex>;
+            if (csr == null)
+            {
+                throw new ArgumentException(Resources.MatrixMustBeSparse, nameof(matrix));
+            }
+
+            // No diagonal element can be omitted from the values array.
+            // If there is a zero value on the diagonal, that element nonetheless must be explicitly represented.
+            csr.PopulateExplicitZerosOnDiagonal();
+
+            var rowCount = csr.RowCount;
+            var columnCount = csr.ColumnCount;
+            var valueCount = csr.ValueCount;
+
+            var values = csr.Values;
+            var rowPointers = csr.RowPointers;
+            var columnIndices = csr.ColumnIndices;
+
+            var nRhs = input.ColumnCount;
+            var rhs = new Complex[rowCount * nRhs];
+            Array.Copy(input.ToColumnMajorArray(), rhs, rhs.Length);
+
+            var solution = new Complex[rowCount * nRhs];
+
+            var error = SparseSolverControl.Provider.Solve(DssMatrixStructure.Nonsymmetric, DssMatrixType.Indefinite, DssSystemType.NonTransposed,
+                rowCount, columnCount, valueCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+
+            if (error == DssStatus.MKL_DSS_SUCCESS)
+                result = Matrix<Complex>.Build.DenseOfColumnMajor(rowCount, nRhs, solution);
+
+            return error;
+        }
+    }
+}

src/Numerics/Providers/SparseSolver/ISparseSolverProvider.cs

@@ -0,0 +1,100 @@
+using Complex = System.Numerics.Complex;
+
+namespace MathNet.Numerics.Providers.SparseSolver
+{
+    public enum DssMatrixStructure : int
+    {
+        SymmetricStructure = 0,
+        Symmetric = 1,
+        Nonsymmetric = 2
+    }
+
+    public enum DssMatrixType : int
+    {
+        PositiveDefinite = 0,
+        Indefinite = 1,
+        HermitianPositiveDefinite = 2,
+        HermitianIndefinite = 3
+    }
+
+    public enum DssSystemType : int
+    {
+        /// <summary>
+        /// Solve a system, Ax = b.
+        /// </summary>
+        NonTransposed = 0,
+        /// <summary>
+        /// Solve a conjugate transposed system, A†x = b
+        /// </summary>
+        ConjugateTransposed = 1,
+        /// <summary>
+        /// Solve a transposed system, A'x = b
+        /// </summary>
+        Transposed = 2
+    }
+
+    public enum DssStatus : int
+    {
+        /// <summary>
+        /// The operation was successful.
+        /// </summary>
+        MKL_DSS_SUCCESS = 0, 
+        MKL_DSS_ZERO_PIVOT = -1,
+        MKL_DSS_OUT_OF_MEMORY = -2,
+        MKL_DSS_FAILURE = -3,
+        MKL_DSS_ROW_ERR = -4,
+        MKL_DSS_COL_ERR = -5,
+        MKL_DSS_TOO_FEW_VALUES = -6,
+        MKL_DSS_TOO_MANY_VALUES = -7,
+        MKL_DSS_NOT_SQUARE = -8,
+        MKL_DSS_STATE_ERR = -9,
+        MKL_DSS_INVALID_OPTION = -10,
+        MKL_DSS_OPTION_CONFLICT = -11,
+        MKL_DSS_MSG_LVL_ERR = -12,
+        MKL_DSS_TERM_LVL_ERR = -13,
+        MKL_DSS_STRUCTURE_ERR = -14,
+        MKL_DSS_REORDER_ERR = -15,
+        MKL_DSS_VALUES_ERR = -16,
+        MKL_DSS_STATISTICS_INVALID_MATRIX = -17,
+        MKL_DSS_STATISTICS_INVALID_STATE = -18,
+        MKL_DSS_STATISTICS_INVALID_STRING = -19,
+        MKL_DSS_REORDER1_ERR = -20,
+        MKL_DSS_PREORDER_ERR = -21,
+        MKL_DSS_DIAG_ERR = -22,
+        MKL_DSS_I32BIT_ERR = -23,
+        MKL_DSS_OOC_MEM_ERR = -24,
+        MKL_DSS_OOC_OC_ERR = -25,
+        MKL_DSS_OOC_RW_ERR = -26,
+    }
+
+    public interface ISparseSolverProvider :
+        ISparseSolverProvider<double>,
+        ISparseSolverProvider<float>,
+        ISparseSolverProvider<Complex>,
+        ISparseSolverProvider<Complex32>
+    {
+        /// <summary>
+        /// Try to find out whether the provider is available, at least in principle.
+        /// Verification may still fail if available, but it will certainly fail if unavailable.
+        /// </summary>
+        bool IsAvailable();
+
+        /// <summary>
+        /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider
+        /// </summary>
+        void InitializeVerify();
+
+        /// <summary>
+        /// Frees memory buffers, caches and handles allocated in or to the provider.
+        /// Does not unload the provider itself, it is still usable afterwards.
+        /// </summary>
+        void FreeResources();
+    }
+
+    public interface ISparseSolverProvider<T>
+        where T : struct
+    {
+        DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType, int rows, int cols, int nnz, int[] rowIdx, int[] colPtr, T[] values, int nRhs, T[] rhs, T[] solution);
+    }
+}
+

src/Numerics/Providers/SparseSolver/Managed/ManagedSparseSolverProvider.cs

@@ -0,0 +1,68 @@
+using System;
+using System.Numerics;
+
+namespace MathNet.Numerics.Providers.SparseSolver.Managed
+{
+    /// <summary>
+    /// The managed sparse solver provider
+    /// </summary>
+    internal partial class ManagedSparseSolverProvider : ISparseSolverProvider
+    {
+        /// <summary>
+        /// Try to find out whether the provider is available, at least in principle.
+        /// Verification may still fail if available, but it will certainly fail if unavailable.
+        /// </summary>
+        public virtual bool IsAvailable()
+        {
+            return true;
+        }
+
+        /// <summary>
+        /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider
+        /// </summary>
+        public virtual void InitializeVerify()
+        {
+        }
+
+        /// <summary>
+        /// Frees memory buffers, caches and handles allocated in or to the provider.
+        /// Does not unload the provider itself, it is still usable afterwards.
+        /// </summary>
+        public virtual void FreeResources()
+        {
+        }
+
+        public override string ToString()
+        {
+            return "Managed";
+        }
+
+        public virtual DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] ColumnIndices, float[] values,
+            int nRhs, float[] rhs, float[] solution)
+        {
+            throw new NotImplementedException();
+        }
+
+        public virtual DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] ColumnIndices, double[] values,
+            int nRhs, double[] rhs, double[] solution)
+        {
+            throw new NotImplementedException();
+        }
+
+        public virtual DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] ColumnIndices, Complex32[] values,
+            int nRhs, Complex32[] rhs, Complex32[] solution)
+        {
+            throw new NotImplementedException();
+        }
+
+        public virtual DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] ColumnIndices, Complex[] values,
+            int nRhs, Complex[] rhs, Complex[] solution)
+        {
+            throw new NotImplementedException();
+        }
+    }
+}

src/Numerics/Providers/SparseSolver/Mkl/MklSparseSolverProvider.Complex.cs

@@ -0,0 +1,69 @@
+#if NATIVE
+
+using MathNet.Numerics.Properties;
+using MathNet.Numerics.Providers.Common.Mkl;
+using System;
+using System.Security;
+using Complex = System.Numerics.Complex;
+
+namespace MathNet.Numerics.Providers.SparseSolver.Mkl
+{
+    /// <summary>
+    /// Intel's Math Kernel Library (MKL) direct sparse solver provider.
+    /// </summary>
+    internal partial class MklSparseSolverProvider
+    {
+        [SecuritySafeCritical]
+        public override DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, Complex[] values,
+            int nRhs, Complex[] rhs, Complex[] solution)
+        {
+            if (rowCount != columnCount)
+            {
+                throw new ArgumentNullException(Resources.ArgumentMatrixSymmetric);
+            }
+
+            if (rowPointers == null)
+            {
+                throw new ArgumentNullException(nameof(rowPointers));
+            }
+
+            if (columnIndices == null)
+            {
+                throw new ArgumentNullException(nameof(columnIndices));
+            }
+
+            if (values == null)
+            {
+                throw new ArgumentNullException(nameof(values));
+            }
+
+            if (rhs == null)
+            {
+                throw new ArgumentNullException(nameof(rhs));
+            }
+
+            if (solution == null)
+            {
+                throw new ArgumentNullException(nameof(solution));
+            }
+
+            if (rowCount * nRhs != rhs.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(rhs));
+            }
+
+            if (columnCount * nRhs != solution.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(solution));
+            }
+
+            var error = SafeNativeMethods.z_dss_solve((int)matrixStructure, (int)matrixType, (int)systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+            return (DssStatus)error;
+        }
+    }
+}
+
+#endif

src/Numerics/Providers/SparseSolver/Mkl/MklSparseSolverProvider.Complex32.cs

@@ -0,0 +1,68 @@
+#if NATIVE
+
+using MathNet.Numerics.Properties;
+using MathNet.Numerics.Providers.Common.Mkl;
+using System;
+using System.Security;
+
+namespace MathNet.Numerics.Providers.SparseSolver.Mkl
+{
+    /// <summary>
+    /// Intel's Math Kernel Library (MKL) direct sparse solver provider.
+    /// </summary>
+    internal partial class MklSparseSolverProvider
+    {
+        [SecuritySafeCritical]
+        public override DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, Complex32[] values,
+            int nRhs, Complex32[] rhs, Complex32[] solution)
+        {
+            if (rowCount != columnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSymmetric);
+            }
+
+            if (rowPointers == null)
+            {
+                throw new ArgumentNullException(nameof(rowPointers));
+            }
+
+            if (columnIndices == null)
+            {
+                throw new ArgumentNullException(nameof(columnIndices));
+            }
+
+            if (values == null)
+            {
+                throw new ArgumentNullException(nameof(values));
+            }
+
+            if (rhs == null)
+            {
+                throw new ArgumentNullException(nameof(rhs));
+            }
+
+            if (solution == null)
+            {
+                throw new ArgumentNullException(nameof(solution));
+            }
+
+            if (rowCount * nRhs != rhs.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(rhs));
+            }
+
+            if (columnCount * nRhs != solution.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(solution));
+            }
+
+            var error = SafeNativeMethods.c_dss_solve((int)matrixStructure, (int)matrixType, (int)systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+            return (DssStatus)error;
+        }
+    }
+}
+
+#endif

src/Numerics/Providers/SparseSolver/Mkl/MklSparseSolverProvider.Double.cs

@@ -0,0 +1,68 @@
+#if NATIVE
+
+using MathNet.Numerics.Properties;
+using MathNet.Numerics.Providers.Common.Mkl;
+using System;
+using System.Security;
+
+namespace MathNet.Numerics.Providers.SparseSolver.Mkl
+{
+    /// <summary>
+    /// Intel's Math Kernel Library (MKL) direct sparse solver provider.
+    /// </summary>
+    internal partial class MklSparseSolverProvider
+    {
+        [SecuritySafeCritical]
+        public override DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, double[] values,
+            int nRhs, double[] rhs, double[] solution)
+        {
+            if (rowCount != columnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSymmetric);
+            }
+
+            if (rowPointers == null)
+            {
+                throw new ArgumentNullException(nameof(rowPointers));
+            }
+
+            if (columnIndices == null)
+            {
+                throw new ArgumentNullException(nameof(columnIndices));
+            }
+
+            if (values == null)
+            {
+                throw new ArgumentNullException(nameof(values));
+            }
+
+            if (rhs == null)
+            {
+                throw new ArgumentNullException(nameof(rhs));
+            }
+
+            if (solution == null)
+            {
+                throw new ArgumentNullException(nameof(solution));
+            }
+
+            if (rowCount * nRhs != rhs.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(rhs));
+            }
+
+            if (columnCount * nRhs != solution.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(solution));
+            }
+
+            var error = SafeNativeMethods.d_dss_solve((int)matrixStructure, (int)matrixType, (int)systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+            return (DssStatus)error;
+        }
+    }
+}
+
+#endif

src/Numerics/Providers/SparseSolver/Mkl/MklSparseSolverProvider.Single.cs

@@ -0,0 +1,67 @@
+#if NATIVE
+
+using MathNet.Numerics.Properties;
+using MathNet.Numerics.Providers.Common.Mkl;
+using System;
+using System.Security;
+
+namespace MathNet.Numerics.Providers.SparseSolver.Mkl
+{
+    /// <summary>
+    /// Intel's Math Kernel Library (MKL) direct sparse solver provider.
+    /// </summary>
+    internal partial class MklSparseSolverProvider
+    {
+        [SecuritySafeCritical]
+        public override DssStatus Solve(DssMatrixStructure matrixStructure, DssMatrixType matrixType, DssSystemType systemType,
+            int rowCount, int columnCount, int nonZerosCount, int[] rowPointers, int[] columnIndices, float[] values,
+            int nRhs, float[] rhs, float[] solution)
+        {
+            if (rowCount != columnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSymmetric);
+            }
+
+            if (rowPointers == null)
+            {
+                throw new ArgumentNullException(nameof(rowPointers));
+            }
+
+            if (columnIndices == null)
+            {
+                throw new ArgumentNullException(nameof(columnIndices));
+            }
+
+            if (values == null)
+            {
+                throw new ArgumentNullException(nameof(values));
+            }
+
+            if (rhs == null)
+            {
+                throw new ArgumentNullException(nameof(rhs));
+            }
+
+            if (solution == null)
+            {
+                throw new ArgumentNullException(nameof(solution));
+            }
+
+            if (rowCount * nRhs != rhs.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(rhs));
+            }
+
+            if (columnCount * nRhs != solution.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentArraysSameLength, nameof(solution));
+            }
+
+            var error = SafeNativeMethods.s_dss_solve((int)matrixStructure, (int)matrixType, (int)systemType,
+                rowCount, columnCount, nonZerosCount, rowPointers, columnIndices, values,
+                nRhs, rhs, solution);
+            return (DssStatus)error;
+        }
+    }
+}
+#endif

src/Numerics/Providers/SparseSolver/Mkl/MklSparseSolverProvider.cs

@@ -0,0 +1,77 @@
+#if NATIVE
+
+using MathNet.Numerics.Providers.Common.Mkl;
+using System;
+
+namespace MathNet.Numerics.Providers.SparseSolver.Mkl
+{
+    /// <summary>
+    /// Intel's Math Kernel Library (MKL) sparse solver provider.
+    /// </summary>
+    internal partial class MklSparseSolverProvider : Managed.ManagedSparseSolverProvider, IDisposable
+    {
+        const int MinimumCompatibleRevision = 12;
+
+        readonly string _hintPath;
+
+        int sparseSolverMajor;
+        int sparseSolverMinor;
+
+        /// <param name="hintPath">Hint path where to look for the native binaries</param>
+        internal MklSparseSolverProvider(string hintPath)
+        {
+            _hintPath = hintPath;
+        }
+
+        /// <summary>
+        /// Try to find out whether the provider is available, at least in principle.
+        /// Verification may still fail if available, but it will certainly fail if unavailable.
+        /// </summary>
+        public override bool IsAvailable()
+        {
+            return MklProvider.IsAvailable(hintPath: _hintPath);
+        }
+
+        /// <summary>
+        /// Initialize and verify that the provided is indeed available.
+        /// If calling this method fails, consider to fall back to alternatives like the managed provider.
+        /// </summary>
+        public override void InitializeVerify()
+        {
+            int revision = MklProvider.Load(_hintPath);
+            if (revision < MinimumCompatibleRevision)
+            {
+                throw new NotSupportedException($"MKL Native Provider revision r{revision} is too old. Consider upgrading to a newer version. Revision r{MinimumCompatibleRevision} and newer are supported.");
+            }
+                        
+            sparseSolverMajor = SafeNativeMethods.query_capability((int)ProviderCapability.SparseSolverMajor);
+            sparseSolverMinor = SafeNativeMethods.query_capability((int)ProviderCapability.SparseSolverMinor);
+            if (!(sparseSolverMajor == 1 && sparseSolverMinor >= 0))
+            {
+                throw new NotSupportedException(string.Format("MKL Native Provider not compatible. Expecting sparse solver v1 but provider implements v{0}.", sparseSolverMajor));
+            }
+        }
+
+        /// <summary>
+        /// Frees memory buffers, caches and handles allocated in or to the provider.
+        /// Does not unload the provider itself, it is still usable afterwards.
+        /// </summary>
+        public override void FreeResources()
+        {
+            MklProvider.FreeResources();
+        }
+
+        public override string ToString()
+        {
+            return MklProvider.Describe();
+        }
+
+        public void Dispose()
+        {
+            FreeResources();
+        }
+    }
+}
+
+#endif
+

src/Numerics/Providers/SparseSolver/SparseSolverControl.cs

@@ -0,0 +1,165 @@
+using System;
+
+namespace MathNet.Numerics.Providers.SparseSolver
+{
+    public static class SparseSolverControl
+    {
+        const string EnvVarSSProvider = "MathNetNumericsSSProvider";
+        const string EnvVarSSProviderPath = "MathNetNumericsSSProviderPath";
+
+        static ISparseSolverProvider _sparseSolverProvider;
+        static readonly object StaticLock = new object();
+
+        /// <summary>
+        /// Gets or sets the sparse solver provider. Consider to use UseNativeMKL or UseManaged instead.
+        /// </summary>
+        /// <value>The linear algebra provider.</value>
+        public static ISparseSolverProvider Provider
+        {
+            get
+            {
+                if (_sparseSolverProvider == null)
+                {
+                    lock (StaticLock)
+                    {
+                        if (_sparseSolverProvider == null)
+                        {
+                            UseDefault();
+                        }
+                    }
+                }
+
+                return _sparseSolverProvider;
+            }
+            set
+            {
+                value.InitializeVerify();
+
+                // only actually set if verification did not throw
+                _sparseSolverProvider = value;
+            }
+        }
+
+        /// <summary>
+        /// Optional path to try to load native provider binaries from.
+        /// If not set, Numerics will fall back to the environment variable
+        /// `MathNetNumericsSSProviderPath` or the default probing paths.
+        /// </summary>
+        public static string HintPath { get; set; }
+
+        public static ISparseSolverProvider CreateManaged()
+        {
+            return new Managed.ManagedSparseSolverProvider();
+        }
+
+        public static void UseManaged()
+        {
+            Provider = CreateManaged();
+        }
+
+#if NATIVE
+        public static ISparseSolverProvider CreateNativeMKL()
+        {
+            return new Mkl.MklSparseSolverProvider(GetCombinedHintPath());
+        }
+
+        public static void UseNativeMKL()
+        {
+            Provider = CreateNativeMKL();
+        }
+
+        public static bool TryUseNativeMKL()
+        {
+            return TryUse(CreateNativeMKL());
+        }
+
+        /// <summary>
+        /// Try to use a native provider, if available.
+        /// </summary>
+        public static bool TryUseNative()
+        {
+            return TryUseNativeMKL();
+        }
+#endif
+
+        static bool TryUse(ISparseSolverProvider provider)
+        {
+            try
+            {
+                if (!provider.IsAvailable())
+                {
+                    return false;
+                }
+
+                Provider = provider;
+                return true;
+            }
+            catch
+            {
+                // intentionally swallow exceptions here - use the explicit variants if you're interested in why
+                return false;
+            }
+        }
+
+        /// <summary>
+        /// Use the best provider available.
+        /// </summary>
+        public static void UseBest()
+        {
+#if NATIVE
+            if (!TryUseNative())
+            {
+                UseManaged();
+            }
+#else
+            UseManaged();
+#endif
+        }
+
+        /// <summary>
+        /// Use a specific provider if configured, e.g. using the
+        /// "MathNetNumericsDSSProvider" environment variable,
+        /// or fall back to the best provider.
+        /// </summary>
+        public static void UseDefault()
+        {
+#if NATIVE
+            var value = Environment.GetEnvironmentVariable(EnvVarSSProvider);
+            switch (value != null ? value.ToUpperInvariant() : string.Empty)
+            {
+
+                case "MKL":
+                    UseNativeMKL();
+                    break;
+
+                default:
+                    UseBest();
+                    break;
+            }
+#else
+            UseBest();
+#endif
+        }
+
+        public static void FreeResources()
+        {
+            Provider.FreeResources();
+        }
+
+        static string GetCombinedHintPath()
+        {
+            if (!String.IsNullOrEmpty(HintPath))
+            {
+                return HintPath;
+            }
+
+            var value = Environment.GetEnvironmentVariable(EnvVarSSProviderPath);
+            if (!String.IsNullOrEmpty(value))
+            {
+                return value;
+            }
+
+            return null;
+        }
+    }
+}