Merge remote-tracking branch 'PohlodSlavik/master' into la-knuth

interesting Knuth algorithm for work with sparse matrixes.
Ill fix it.
PS: Sorry for my bad english, I just learn it)

Backup/MathNet.Numerics.sln

@@ -0,0 +1,88 @@
+
+Microsoft Visual Studio Solution File, Format Version 11.00
+# Visual Studio 2010
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Examples", "src\Examples\Examples.csproj", "{8239A6FF-1EF3-4DA4-A860-95C392DD6899}"
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Numerics", "src\Numerics\Numerics.csproj", "{B7CAE5F4-A23F-4438-B5BE-41226618B695}"
+EndProject
+Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "FSharp", "src\FSharp\FSharp.fsproj", "{37E8E802-A354-4114-BFC1-6E1357DA605B}"
+EndProject
+Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "FSharpExamples", "src\FSharpExamples\FSharpExamples.fsproj", "{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1}"
+EndProject
+Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "FSharpUnitTests", "src\FSharpUnitTests\FSharpUnitTests.fsproj", "{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D}"
+EndProject
+Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Tests", "Tests", "{4D50FB34-10BC-495A-8B2F-482E34B4D771}"
+EndProject
+Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Samples", "Samples", "{49EE74BD-301F-4C3B-B76A-07F90CC88CE7}"
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests", "src\UnitTests\UnitTests.csproj", "{DAF07AA8-C5C9-4963-98F7-2C3285064DAD}"
+EndProject
+Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Readme", "Readme", "{C2F37492-38AE-4186-8A7F-17B0B080942C}"
+	ProjectSection(SolutionItems) = preProject
+		AUTHORS.markdown = AUTHORS.markdown
+		COPYRIGHT.markdown = COPYRIGHT.markdown
+		README.markdown = README.markdown
+	EndProjectSection
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Numerics.IO", "src\Numerics.IO\Numerics.IO.csproj", "{EB1A5D32-F264-4BCE-BEB7-0B97085075BE}"
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|Any CPU = Debug|Any CPU
+		Release|Any CPU = Release|Any CPU
+		Release-Signed|Any CPU = Release-Signed|Any CPU
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{8239A6FF-1EF3-4DA4-A860-95C392DD6899}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{8239A6FF-1EF3-4DA4-A860-95C392DD6899}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{8239A6FF-1EF3-4DA4-A860-95C392DD6899}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{8239A6FF-1EF3-4DA4-A860-95C392DD6899}.Release|Any CPU.Build.0 = Release|Any CPU
+		{8239A6FF-1EF3-4DA4-A860-95C392DD6899}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
+		{8239A6FF-1EF3-4DA4-A860-95C392DD6899}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
+		{B7CAE5F4-A23F-4438-B5BE-41226618B695}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{B7CAE5F4-A23F-4438-B5BE-41226618B695}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release|Any CPU.Build.0 = Release|Any CPU
+		{B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
+		{B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
+		{37E8E802-A354-4114-BFC1-6E1357DA605B}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{37E8E802-A354-4114-BFC1-6E1357DA605B}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{37E8E802-A354-4114-BFC1-6E1357DA605B}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{37E8E802-A354-4114-BFC1-6E1357DA605B}.Release|Any CPU.Build.0 = Release|Any CPU
+		{37E8E802-A354-4114-BFC1-6E1357DA605B}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
+		{37E8E802-A354-4114-BFC1-6E1357DA605B}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
+		{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1}.Release|Any CPU.Build.0 = Release|Any CPU
+		{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
+		{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
+		{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D}.Release|Any CPU.Build.0 = Release|Any CPU
+		{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
+		{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
+		{DAF07AA8-C5C9-4963-98F7-2C3285064DAD}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{DAF07AA8-C5C9-4963-98F7-2C3285064DAD}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{DAF07AA8-C5C9-4963-98F7-2C3285064DAD}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{DAF07AA8-C5C9-4963-98F7-2C3285064DAD}.Release|Any CPU.Build.0 = Release|Any CPU
+		{DAF07AA8-C5C9-4963-98F7-2C3285064DAD}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
+		{DAF07AA8-C5C9-4963-98F7-2C3285064DAD}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
+		{EB1A5D32-F264-4BCE-BEB7-0B97085075BE}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{EB1A5D32-F264-4BCE-BEB7-0B97085075BE}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{EB1A5D32-F264-4BCE-BEB7-0B97085075BE}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{EB1A5D32-F264-4BCE-BEB7-0B97085075BE}.Release|Any CPU.Build.0 = Release|Any CPU
+		{EB1A5D32-F264-4BCE-BEB7-0B97085075BE}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
+		{EB1A5D32-F264-4BCE-BEB7-0B97085075BE}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+	GlobalSection(NestedProjects) = preSolution
+		{BC81EA37-8EE6-4BF9-B8A9-B30497AEF8B1} = {49EE74BD-301F-4C3B-B76A-07F90CC88CE7}
+		{8239A6FF-1EF3-4DA4-A860-95C392DD6899} = {49EE74BD-301F-4C3B-B76A-07F90CC88CE7}
+		{F2F8032B-A31D-4E33-A05E-F2CDCBFAA75D} = {4D50FB34-10BC-495A-8B2F-482E34B4D771}
+		{DAF07AA8-C5C9-4963-98F7-2C3285064DAD} = {4D50FB34-10BC-495A-8B2F-482E34B4D771}
+	EndGlobalSection
+EndGlobal

Backup/src/FSharp/AssemblyInfo.fs

@@ -0,0 +1,51 @@
+// <copyright file="AssemblyInfo.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics
+
+open System.Reflection
+open System.Resources;
+open System.Runtime.CompilerServices
+open System.Runtime.InteropServices
+
+[<assembly: AssemblyTitle("Math.NET Numerics F# Modules")>]
+[<assembly: AssemblyDescription("F# Modules for Math.NET Numerics, providing methods and algorithms for numerical computations in science, engineering and every day use.")>]
+[<assembly: AssemblyConfiguration("")>]
+[<assembly: AssemblyCompany("Math.NET Project")>]
+[<assembly: AssemblyProduct("Math.NET Numerics")>]
+[<assembly: AssemblyCopyright("Copyright © Math.NET Project")>]
+[<assembly: AssemblyTrademark("")>]
+[<assembly: AssemblyCulture("")>]
+[<assembly: ComVisible(false)>]
+[<assembly: Guid("048BC4EB-CE2B-4040-9967-4784F5405B0F")>]
+[<assembly: NeutralResourcesLanguage("en")>]
+[<assembly: AssemblyVersion("1.0.0.0")>]
+[<assembly: AssemblyFileVersion("1.0.0.0")>]
+()

Backup/src/FSharp/DenseMatrix.fs

@@ -0,0 +1,101 @@
+// <copyright file="DenseMatrix.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.LinearAlgebra
+open MathNet.Numerics.LinearAlgebra.Generic
+
+/// A module which implements functional dense vector operations.
+module DenseMatrix =
+
+    /// Initialize a matrix by calling a construction function for every element.
+    let inline init (n: int) (m: int) f =
+        let A = new DenseMatrix(n,m)
+        for i=0 to n-1 do
+            for j=0 to m-1 do
+                A.[i,j] <- f i j
+        A
+    
+    /// Create a matrix from a list of float lists. Every list in the master list specifies a row.
+    let inline ofList (fll: float list list) =
+        let n = List.length fll
+        let m = List.length (List.head fll)
+        let A = DenseMatrix(n,m)
+        fll |> List.iteri (fun i fl ->
+                            if (List.length fl) <> m then failwith "Each subrow must be of the same length." else
+                            List.iteri (fun j f -> A.[i,j] <- f) fl)
+        A
+    
+    /// Create a matrix from a list of sequences. Every sequence in the master sequence specifies a row.
+    let inline ofSeq (fss: #seq<#seq<float>>) =
+        let n = Seq.length fss
+        let m = Seq.length (Seq.head fss)
+        let A = DenseMatrix(n,m)
+        fss |> Seq.iteri (fun i fs ->
+                            if (Seq.length fs) <> m then failwith "Each subrow must be of the same length." else
+                            Seq.iteri (fun j f -> A.[i,j] <- f) fs)
+        A
+
+    /// Create a matrix from a 2D array of floating point numbers.
+    let inline ofArray2 (arr: float[,]) = new DenseMatrix(arr)
+    
+    /// Create a matrix with the given entries.
+    let inline initDense (n: int) (m: int) (es: #seq<int * int * float>) =
+        let A = new DenseMatrix(n,m)
+        Seq.iter (fun (i,j,f) -> A.[i,j] <- f) es
+        A
+    
+    /// Create a square matrix with constant diagonal entries.
+    let inline constDiag (n: int) (f: float) =
+        let A = new DenseMatrix(n,n)
+        for i=0 to n-1 do
+            A.[i,i] <- f
+        A
+    
+    /// Create a square matrix with the vector elements on the diagonal.
+    let inline diag (v: #Vector<float>) =
+        let n = v.Count
+        let A = new DenseMatrix(n,n)
+        for i=0 to n-1 do
+            A.[i,i] <- v.Item(i)
+        A
+
+    /// Initialize a matrix by calling a construction function for every row.
+    let inline initRow (n: int) (m: int) (f: int -> #Vector<float>) =
+        let A = new DenseMatrix(n,m)
+        for i=0 to n-1 do A.SetRow(i, f i)
+        A
+
+    /// Initialize a matrix by calling a construction function for every column.
+    let inline initCol (n: int) (m: int) (f: int -> #Vector<float>) =
+        let A = new DenseMatrix(n,m)
+        for i=0 to m-1 do A.SetColumn(i, f i)
+        A

Backup/src/FSharp/DenseVector.fs

@@ -0,0 +1,69 @@
+// <copyright file="DenseVector.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+
+open MathNet.Numerics.LinearAlgebra
+
+/// A module which implements functional dense vector operations.
+module DenseVector =
+
+    /// Initialize a vector by calling a construction function for every element.
+    let inline init (n: int) f =
+        let v = new Double.DenseVector(n)
+        for i=0 to n-1 do
+            v.[i] <- f i
+        v
+
+    /// Create a vector from a float list.
+    let inline ofList (fl: float list) =
+        let n = List.length fl
+        let v = Double.DenseVector(n)
+        fl |> List.iteri (fun i f -> v.[i] <- f)
+        v
+    
+    /// Create a vector from a sequences.
+    let inline ofSeq (fs: #seq<float>) =
+        let n = Seq.length fs
+        let v = DenseVector(n)
+        fs |> Seq.iteri (fun i f -> v.[i] <- f)
+        v
+    
+    /// Create a vector with evenly spaced entries: e.g. rangef -1.0 0.5 1.0 = [-1.0 -0.5 0.0 0.5 1.0]
+    let inline rangef (start: float) (step: float) (stop: float) =
+        let n = (int ((stop - start) / step)) + 1
+        let v = new DenseVector(n)
+        for i=0 to n-1 do
+            v.[i] <- (float i) * step + start
+        v
+    
+    /// Create a vector with integer entries in the given range.
+    let inline range (start: int) (stop: int) =
+        new DenseVector([| for i in [start .. stop] -> float i |])

Backup/src/FSharp/Extensions.fs

@@ -0,0 +1,68 @@
+// <copyright file="Extensions.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Generic
+
+// Module that contains implementation of useful F#-specific
+// extension members for generic Matrix and Vector types
+[<AutoOpen>]
+module FSharpExtensions = 
+
+  // A type extension for the generic vector type that 
+  // adds the 'GetSlice' method to allow vec.[a .. b] syntax
+  type MathNet.Numerics.LinearAlgebra.Generic.
+      Vector<'T when 'T : struct and 'T : (new : unit -> 'T) 
+                 and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable 
+                 and 'T :> System.ValueType> with
+
+    /// Gets a slice of a vector starting at a specified index
+    /// and ending at a specified index (both indices are optional)
+    /// This method can be used via the x.[start .. finish] syntax
+    member x.GetSlice(start, finish) = 
+      let start = defaultArg start 0
+      let finish = defaultArg finish (x.Count - 1)
+      x.SubVector(start, finish - start + 1)
+
+  // A type extension for the generic matrix type that
+  // adds the 'GetSlice' method to allow m.[r1 .. r2, c1 .. c2] syntax
+  type MathNet.Numerics.LinearAlgebra.Generic.
+      Matrix<'T when 'T : struct and 'T : (new : unit -> 'T) 
+                 and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable 
+                 and 'T :> System.ValueType> with
+
+    /// Gets a submatrix using a specified column range and 
+    /// row range (all indices are optional)
+    /// This method can be used via the x.[r1 .. r2, c1 .. c2 ] syntax
+    member x.GetSlice(rstart, rfinish, cstart, cfinish) = 
+      let cstart = defaultArg cstart 0
+      let rstart = defaultArg rstart 0
+      let cfinish = defaultArg cfinish (x.ColumnCount - 1)
+      let rfinish = defaultArg rfinish (x.RowCount - 1)
+      x.SubMatrix(rstart, rfinish - rstart + 1, cstart, cfinish - cstart + 1)

Backup/src/FSharp/FSharp.fsproj

@@ -0,0 +1,80 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="4.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <PropertyGroup>
+    <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
+    <Platform Condition=" '$(Platform)' == '' ">AnyCPU</Platform>
+    <ProductVersion>8.0.30703</ProductVersion>
+    <SchemaVersion>2.0</SchemaVersion>
+    <ProjectGuid>{37e8e802-a354-4114-bfc1-6e1357da605b}</ProjectGuid>
+    <OutputType>Library</OutputType>
+    <RootNamespace>FSharp</RootNamespace>
+    <AssemblyName>MathNet.Numerics.FSharp</AssemblyName>
+    <TargetFrameworkVersion>v4.0</TargetFrameworkVersion>
+    <Name>FSharp</Name>
+    <TargetFrameworkProfile />
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
+    <DebugSymbols>true</DebugSymbols>
+    <DebugType>full</DebugType>
+    <Optimize>false</Optimize>
+    <Tailcalls>false</Tailcalls>
+    <OutputPath>..\..\out\debug\Net40\</OutputPath>
+    <DefineConstants>DEBUG;TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <Tailcalls>true</Tailcalls>
+    <OutputPath>..\..\out\lib\Net40\</OutputPath>
+    <DefineConstants>TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+    <DocumentationFile>..\..\out\lib\Net40\MathNet.Numerics.FSharp.xml</DocumentationFile>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)' == 'Release' ">
+    <DocumentationFile>MathNet.Numerics.FSharp.XML</DocumentationFile>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release-Signed|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <Tailcalls>true</Tailcalls>
+    <DefineConstants>TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+    <DocumentationFile>..\..\out\lib\Net40\MathNet.Numerics.FSharp.xml</DocumentationFile>
+    <OutputPath>bin\Release-Signed\</OutputPath>
+  </PropertyGroup>
+  <Import Project="$(MSBuildExtensionsPath32)\FSharp\1.0\Microsoft.FSharp.Targets" Condition="!Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <Import Project="$(MSBuildExtensionsPath32)\..\Microsoft F#\v4.0\Microsoft.FSharp.Targets" Condition="Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <ItemGroup>
+    <Compile Include="Extensions.fs" />
+    <Compile Include="DenseVector.fs" />
+    <Compile Include="DenseMatrix.fs" />
+    <Compile Include="SparseMatrix.fs" />
+    <Compile Include="SparseVector.fs" />
+    <Compile Include="Vector.fs" />
+    <Compile Include="Matrix.fs" />
+    <Compile Include="Main.fs" />
+    <Compile Include="AssemblyInfo.fs" />
+  </ItemGroup>
+  <ItemGroup>
+    <Reference Include="FSharp.Core" />
+    <Reference Include="mscorlib" />
+    <Reference Include="System" />
+    <Reference Include="System.Core">
+      <RequiredTargetFramework>3.5</RequiredTargetFramework>
+    </Reference>
+    <Reference Include="System.Numerics" />
+    <ProjectReference Include="..\Numerics\Numerics.csproj">
+      <Name>Numerics</Name>
+      <Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
+      <Private>True</Private>
+    </ProjectReference>
+  </ItemGroup>
+  <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
+	     Other similar extension points exist, see Microsoft.Common.targets.
+	<Target Name="BeforeBuild">
+	</Target>
+	<Target Name="AfterBuild">
+	</Target>
+	-->
+</Project>

Backup/src/FSharp/Main.fs

@@ -0,0 +1,43 @@
+// <copyright file="Main.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics
+
+open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.LinearAlgebra.Generic
+
+/// A module which implements some F# utility functions.
+module FSharp =
+
+    /// Construct a dense matrix from a list of floating point numbers.
+    let inline matrix (lst: list<list<float>>) = DenseMatrix.ofList lst :> Matrix<float>
+
+    /// Construct a dense vector from a list of floating point numbers.
+    let inline vector (lst: list<float>) = DenseVector.ofList lst :> Vector<float>

Backup/src/FSharp/Matrix.fs

@@ -0,0 +1,243 @@
+// <copyright file="Matrix.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+
+open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.LinearAlgebra.Generic
+
+/// A module which implements functional matrix operations.
+module Matrix =
+    
+    /// Fold a function over all matrix elements.
+    let inline fold (f: 'a -> float -> 'a) (acc0: 'a) (A: #Matrix<float>) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        let mutable acc = acc0
+        for i=0 to n-1 do
+            for j=0 to m-1 do
+                acc <- f acc (A.Item(i,j))
+        acc
+
+    /// Fold a function over all matrix elements in reverse order.
+    let inline foldBack (f: float -> 'a -> 'a) (acc0: 'a) (A: #Matrix<float>) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        let mutable acc = acc0
+        for i in n-1 .. -1 .. 0 do
+            for j in m-1 .. -1 .. 0 do
+                acc <- f (A.Item(i,j)) acc
+        acc
+
+    /// Fold a matrix by applying a given function to all matrix elements.
+    let inline foldi (f: int -> int -> 'a -> float -> 'a) (acc0: 'a) (A: #Matrix<float>) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        let mutable acc = acc0
+        for i=0 to n-1 do
+            for j=0 to m-1 do
+                acc <- f i j acc (A.Item(i,j))
+        acc
+
+    /// Create a 2D array from a matrix.
+    let inline toArray2 (A: #Matrix<float>) =
+        let n = A.RowCount
+        let m = A.ColumnCount
+        Array2D.init n m (fun i j -> (A.Item(i,j)))
+    
+    /// Checks whether a predicate holds for all elements of a matrix.  
+    let inline forall (p: float -> bool) (A: #Matrix<float>) =
+        let mutable b = true
+        let mutable i = 0
+        let mutable j = 0
+        while b && i < A.RowCount do
+            b <- b && (p (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Chechks whether a predicate holds for at least one element of a matrix.
+    let inline exists (p: float -> bool) (A: #Matrix<float>) =
+        let mutable b = false
+        let mutable i = 0
+        let mutable j = 0
+        while not(b) && i < A.RowCount do
+            b <- b || (p (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Checks whether a position dependent predicate holds for all elements of a matrix.
+    let inline foralli (p: int -> int -> float -> bool) (A: #Matrix<float>) =
+        let mutable b = true
+        let mutable i = 0
+        let mutable j = 0
+        while b && i < A.RowCount do
+            b <- b && (p i j (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Checks whether a position dependent predicate holds for at least one element of a matrix.
+    let inline existsi (p: int -> int -> float -> bool) (A: #Matrix<float>) =
+        let mutable b = false
+        let mutable i = 0
+        let mutable j = 0
+        while not(b) && i < A.RowCount do
+            b <- b || (p i j (A.Item(i,j)))
+            j <- j+1
+            if j = A.ColumnCount then i <- i+1; j <- 0
+        b
+    
+    /// Map every matrix element using the given function.
+    let inline map (f: float -> float) (A: #Matrix<float>) =
+        let N = A.RowCount
+        let M = A.ColumnCount
+        let C = A.Clone()
+        for i=0 to N-1 do
+            for j=0 to M-1 do
+                C.[i,j] <- f (C.Item(i,j))
+        C
+    
+    /// Map every matrix element using the given position dependent function.
+    let inline mapi (f: int -> int -> float -> float) (A: #Matrix<float>) =
+        let N = A.RowCount
+        let M = A.ColumnCount
+        let C = A.Clone()
+        for i=0 to N-1 do
+            for j=0 to M-1 do
+                C.[i,j] <- f i j (C.Item(i,j))
+        C
+        
+    /// In-place map every matrix column using the given position dependent function.
+    let inline inplaceMapCols (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        for j = 0 to A.ColumnCount-1 do
+            A.SetColumn(j, f j (A.Column(j)))
+        ()
+
+    /// In-place map every matrix row using the given position dependent function.
+    let inline inplaceMapRows (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        for i = 0 to A.RowCount-1 do
+            A.SetRow(i, f i (A.Row(i)))
+        ()
+
+    /// Map every matrix column using the given position dependent function.
+    let inline mapCols (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        let A = A.Clone()
+        inplaceMapCols f A
+        A
+
+    /// Map every matrix row using the given position dependent function.
+    let inline mapRows (f: int -> Vector<float> -> Vector<float>) (A: #Matrix<float>) =
+        let A = A.Clone()
+        inplaceMapRows f A
+        A
+
+    /// In-place assignment.
+    let inline inplaceAssign (f: int -> int -> float) (A: #Matrix<float>) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                A.Item(i,j) <- f i j
+    
+    /// In-place map of every matrix element using a position dependent function.
+    let inline inplaceMapi (f: int -> int -> float -> float) (A: #Matrix<float>) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                A.Item(i,j) <- f i j (A.Item(i,j))
+        
+    /// Creates a sequence that iterates the non-zero entries in the matrix.
+    let inline nonZeroEntries (A: #Matrix<float>) =
+        seq { for i in 0 .. A.RowCount-1 do
+                for j in 0 .. A.ColumnCount-1 do
+                  if A.Item(i,j) <> 0.0 then yield (i,j, A.Item(i,j)) }
+    
+    /// Returns the sum of all elements of a matrix.
+    let inline sum (A: #Matrix<float>) =
+        let mutable f = 0.0
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f <- f + A.Item(i,j)
+        f
+
+    /// Returns the sum of the results generated by applying a position dependent function to each column of the matrix.
+    let inline sumColsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<float>) =
+        A.ColumnEnumerator() |> Seq.map (fun (j,col) -> f j col) |> Seq.reduce (+)
+    
+    /// Returns the sum of the results generated by applying a position dependent function to each row of the matrix.
+    let inline sumRowsBy (f: int -> Vector<float> -> 'a) (A: #Matrix<float>) =
+        A.RowEnumerator() |> Seq.map (fun (i,row) -> f i row)  |> Seq.reduce (+)
+
+    /// Iterates over all elements of a matrix.
+    let inline iter (f: float -> unit) (A: #Matrix<float>) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f (A.Item(i,j))
+        ()
+    
+    /// Iterates over all elements of a matrix using the element indices.
+    let inline iteri (f: int -> int -> float -> unit) (A: #Matrix<float>) =
+        for i=0 to A.RowCount-1 do
+            for j=0 to A.ColumnCount-1 do
+                f i j (A.Item(i,j))
+        ()
+        
+    /// Fold one column.
+    let inline foldCol (f: 'a -> float -> 'a) acc (A: #Matrix<float>) k =
+        let mutable macc = acc
+        for i=0 to A.RowCount-1 do
+            macc <- f macc (A.Item(i,k))
+        macc
+    
+    /// Fold one row.
+    let inline foldRow (f: 'a -> float -> 'a) acc (A: #Matrix<float>) k =
+        let mutable macc = acc
+        for i=0 to A.ColumnCount-1 do
+            macc <- f macc (A.Item(k,i))
+        macc
+
+    /// Fold all columns into one row vector.
+    let inline foldByCol (f: float -> float -> float) acc (A: #Matrix<float>) =
+        let v = new DenseVector(A.ColumnCount)
+        for k=0 to A.ColumnCount-1 do
+            let mutable macc = acc
+            for i=0 to A.RowCount-1 do
+                macc <- f macc (A.Item(i,k))
+            v.[k] <- macc
+        v :> Vector<float>
+    
+    /// Fold all rows into one column vector.
+    let inline foldByRow (f: float -> float -> float) acc (A: #Matrix<float>) =
+        let v = new DenseVector(A.RowCount)
+        for k=0 to A.RowCount-1 do
+            let mutable macc = acc
+            for i=0 to A.ColumnCount-1 do
+                macc <- f macc (A.Item(k,i))
+            v.[k] <- macc
+        v :> Vector<float>

Backup/src/FSharp/SparseMatrix.fs

@@ -0,0 +1,75 @@
+// <copyright file="SparseMatrix.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+
+open MathNet.Numerics.LinearAlgebra
+open MathNet.Numerics.LinearAlgebra.Generic
+/// A module which implements functional sparse vector operations.
+module SparseMatrix =
+    
+    /// Create a matrix from a list of float lists. Every list in the master list specifies a row.
+    let inline ofList (rows: int) (cols: int) (fll: list<int * int * float>) =
+        let A = new Double.SparseMatrix(rows, cols)
+        fll |> List.iter (fun (i, j, x) -> A.[i,j] <- x)
+        A
+    
+    /// Create a matrix from a list of sequences. Every sequence in the master sequence specifies a row.
+    let inline ofSeq (rows: int) (cols: int) (fss: #seq<int * int * float>) =
+        let A = new Double.SparseMatrix(rows, cols)
+        fss |> Seq.iter (fun (i, j, x) -> A.[i,j] <- x)
+        A
+    
+    /// Create a square matrix with constant diagonal entries.
+    let inline constDiag (n: int) (f: float) =
+        let A = new Double.SparseMatrix(n,n)
+        for i=0 to n-1 do
+            A.[i,i] <- f
+        A
+    
+    /// Create a square matrix with the vector elements on the diagonal.
+    let inline diag (v: #Vector<float>) =
+        let n = v.Count
+        let A = new Double.SparseMatrix(n,n)
+        for i=0 to n-1 do
+            A.[i,i] <- v.Item(i)
+        A
+
+    /// Initialize a matrix by calling a construction function for every row.
+    let inline initRow (n: int) (m: int) (f: int -> #Vector<float>) =
+        let A = new Double.SparseMatrix(n,m)
+        for i=0 to n-1 do A.SetRow(i, f i)
+        A
+
+    /// Initialize a matrix by calling a construction function for every column.
+    let inline initCol (n: int) (m: int) (f: int -> #Vector<float>) =
+        let A = new Double.SparseMatrix(n,m)
+        for i=0 to m-1 do A.SetColumn(i, f i)
+        A

Backup/src/FSharp/SparseVector.fs

@@ -0,0 +1,48 @@
+// <copyright file="SparseVector.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+
+open MathNet.Numerics.LinearAlgebra
+
+/// A module which implements functional sparse vector operations.
+module SparseVector =
+    
+    /// Create a sparse vector with a given dimension from a list of entry, value pairs.
+    let inline ofList (dim: int) (fl: list<int * float>) =
+        let v = new Double.SparseVector(dim)
+        fl |> List.iter (fun (i, f) -> v.[i] <- f)
+        v
+    
+    /// Create a sparse vector with a given dimension from a sequence of entry, value pairs.
+    let inline ofSeq (dim: int) (fs: #seq<int * float>) =
+        let v = new Double.SparseVector(dim)
+        fs |> Seq.iter (fun (i, f) -> v.[i] <- f)
+        v

Backup/src/FSharp/Vector.fs

@@ -0,0 +1,187 @@
+// <copyright file="Vector.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+
+open MathNet.Numerics.LinearAlgebra
+open MathNet.Numerics.LinearAlgebra.Generic
+
+/// A module which implements functional vector operations.
+module Vector =
+        
+    /// Transform a vector into an array.
+    let inline toArray (v: #Vector<float>) =
+        let n = v.Count
+        Array.init n (fun i -> v.Item(i))
+        
+    /// Transform a vector into an array.
+    let inline toList (v: #Vector<float>) =
+        let n = v.Count
+        List.init n (fun i -> v.Item(i))
+
+    /// In-place mutation by applying a function to every element of the vector.
+    let inline mapInPlace (f: float -> float) (v: #Vector<float>) =
+        for i=0 to v.Count-1 do
+            v.Item(i) <- f (v.Item(i))
+        ()
+
+    /// In-place mutation by applying a function to every element of the vector.
+    let inline mapiInPlace (f: int -> float -> float) (v: #Vector<float>) =
+        for i=0 to v.Count-1 do
+            v.Item(i) <- f i (v.Item(i))
+        ()
+        
+    /// In-place vector addition.
+    let inline addInPlace (v: #Vector<float>) (w: #Vector<float>) = v.Add(w, v)
+        
+    /// In place vector subtraction.
+    let inline subInPlace (v: #Vector<float>) (w: #Vector<float>) = v.Subtract(w, v)
+    
+    /// Functional map operator for vectors.
+    /// <include file='../../../../FSharpExamples/DenseVector.xml' path='example'/> 
+    let inline map f (v: #Vector<float>) =
+        let w = v.Clone()
+        mapInPlace (fun x -> f x) w
+        w
+        
+    /// Applies a function to all elements of the vector.
+    let inline iter (f: float -> unit) (v: #Vector<float>) =
+        for i=0 to v.Count-1 do
+            f (v.Item i)
+        
+    /// Applies a function to all elements of the vector.
+    let inline iteri (f: int -> float -> unit) (v: #Vector<float>) =
+        for i=0 to v.Count-1 do
+            f i (v.Item i)
+            
+    /// Maps a vector to a new vector by applying a function to every element.
+    let inline mapi (f: int -> float -> float) (v: #Vector<float>) =
+        let w = v.Clone()
+        mapiInPlace f w
+        w
+
+    /// Fold all entries of a vector.
+    let inline fold (f: 'a -> float -> 'a) (acc0: 'a) (v: #Vector<float>) =
+        let mutable acc = acc0
+        for i=0 to v.Count-1 do
+            acc <- f acc (v.Item(i))
+        acc
+
+    /// Fold all entries of a vector in reverse order.
+    let inline foldBack (f: float -> 'a -> 'a) (acc0: 'a) (v: #Vector<float>) =
+        let mutable acc = acc0
+        for i=2 to v.Count do
+            acc <- f (v.Item(v.Count - i)) acc
+        acc
+
+    /// Fold all entries of a vector using a position dependent folding function.
+    let inline foldi (f: int -> 'a -> float -> 'a) (acc0: 'a) (v: #Vector<float>) =
+        let mutable acc = acc0
+        for i=0 to v.Count-1 do
+            acc <- f i acc (v.Item(i))
+        acc
+        
+    /// Checks whether a predicate is satisfied for every element in the vector.
+    let inline forall (p: float -> bool) (v: #Vector<float>) =
+        let mutable b = true
+        let mutable i = 0
+        while b && i < v.Count do
+            b <- b && (p (v.Item(i)))
+            i <- i+1
+        b
+    
+    /// Checks whether there is an entry in the vector that satisfies a given predicate.
+    let inline exists (p: float -> bool) (v: #Vector<float>) =
+        let mutable b = false
+        let mutable i = 0
+        while not(b) && i < v.Count do
+            b <- b || (p (v.Item(i)))
+            i <- i+1
+        b
+    
+    /// Checks whether a predicate is true for all entries in a vector.
+    let inline foralli (p: int -> float -> bool) (v: #Vector<float>) =
+        let mutable b = true
+        let mutable i = 0
+        while b && i < v.Count do
+            b <- b && (p i (v.Item(i)))
+            i <- i+1
+        b
+    
+    /// Checks whether there is an entry in the vector that satisfies a given position dependent predicate.
+    let inline existsi (p: int -> float -> bool) (v: #Vector<float>) =
+        let mutable b = false
+        let mutable i = 0
+        while not(b) && i < v.Count do
+            b <- b || (p i (v.Item(i)))
+            i <- i+1
+        b
+
+    /// Scans a vector; like fold but returns the intermediate result.
+    let inline scan (f: float -> float -> float) (v: #Vector<float>) =
+        let w = v.Clone()
+        let mutable p = v.Item(0)
+        for i=1 to v.Count-1 do
+            p <- f p (v.Item(i))
+            w.[i] <- p
+        w
+
+    /// Scans a vector in reverse order; like foldBack but returns the intermediate result.
+    let inline scanBack (f: float -> float -> float) (v: #Vector<float>) =
+        let w = v.Clone()
+        let mutable p = v.Item(v.Count-1)
+        for i=2 to v.Count do
+            p <- f (v.Item(v.Count - i)) p
+            w.[v.Count - i] <- p
+        w
+
+    /// Reduces a vector: the result of this function will be f(...f(f(v[0],v[1]), v[2]),..., v[n]).
+    let inline reduce (f: float -> float -> float) (v: #Vector<float>) =
+        let mutable p = v.Item(0)
+        for i=1 to v.Count-1 do
+            p <- f p (v.Item(i))
+        p
+
+    /// Reduces a vector in reverse order: the result of this function will be f(v[1], ..., f(v[n-2], f(v[n-1],v[n]))...).
+    let inline reduceBack (f: float -> float -> float) (v: #Vector<float>) =
+        let mutable p = v.Item(v.Count-1)
+        for i=2 to v.Count do
+            p <- f (v.Item(v.Count - i)) p
+        p
+
+    /// Creates a new vector and inserts the given value at the given index.
+    let inline insert index value (v: #Vector<float>) =        
+        let newV = new DenseVector(v.Count + 1)                
+        for i = 0 to index - 1 do 
+            newV.Item(i) <- v.Item(i)        
+        newV.Item(index) <- value        
+        for i = index + 1 to v.Count do 
+            newV.Item(i) <- v.Item(i - 1)        
+        newV

Backup/src/FSharpExamples/Apply.fs

@@ -0,0 +1,94 @@
+// <copyright file="Apply.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+module MathNet.Numerics.FSharp.Examples.Apply
+
+open System.Numerics
+open MathNet.Numerics
+open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.LinearAlgebra.Generic
+/// Flag to specify wether we want pretty printing or tab separated output.
+let prettyPrint = false
+
+/// The size of the vector we want to map things for.
+let N = 1000000
+/// The number of times we repeat a call.
+let T = 10
+/// The list of all functions we want to test.
+let FunctionList : (string * (float -> float) * (float -> float)) [] =
+    [| ("Cosine", cos, System.Math.Cos);
+       ("Sine", sin, System.Math.Sin);
+       ("Tangent", tan, System.Math.Tan);
+       ("Inverse Cosine", acos, System.Math.Acos);
+       ("Inverse Sine", asin, System.Math.Asin);
+       ("Inverse Tangent", atan, System.Math.Atan);
+       ("Hyperbolic Cosine", cosh, System.Math.Cosh);
+       ("Hyperbolic Sine", sinh, System.Math.Sinh);
+       ("Hyperbolic Tangent", tanh, System.Math.Tanh);
+       ("Abs", abs, System.Math.Abs);
+       ("Exp", exp, System.Math.Exp);
+       ("Log", log, System.Math.Log);
+       ("Sqrt", sqrt, System.Math.Sqrt);
+       ("Error Function", SpecialFunctions.Erf, SpecialFunctions.Erf);
+       ("Error Function Complement", SpecialFunctions.Erfc, SpecialFunctions.Erfc);
+       ("Inverse Error Function", SpecialFunctions.ErfInv, SpecialFunctions.ErfInv);
+       ("Inverse Error Function Complement", SpecialFunctions.ErfcInv, SpecialFunctions.ErfcInv) |]
+
+/// A vector with random entries.
+let w =
+    let rnd = new Random.MersenneTwister()
+    (new DenseVector(Array.init N (fun _ -> rnd.NextDouble() * 10.0))) :> Vector<float>
+
+/// A stopwatch to time the execution.
+let sw = new System.Diagnostics.Stopwatch()
+
+
+for (name, fs, dotnet) in FunctionList do
+    if prettyPrint then printfn "Running %s on an %d dimensional vector for %d iterations:" name N T
+    else printf "%s" name
+
+    /// Perform the standard F# map function.
+    do
+        let v = w.Clone()
+        sw.Start()
+        for t in 1 .. T do Vector.mapInPlace fs v
+        sw.Stop()
+        if prettyPrint then printfn "\tVector.map (F#): %d milliseconds." sw.ElapsedMilliseconds
+        else printf "\t%d" sw.ElapsedMilliseconds
+        sw.Reset()
+    (*
+    /// Perform the Apply.Map function.
+    do
+        let v = w.Clone()
+        sw.Start()
+        for t in 1 .. T do v.Map(fun x -> dotnet x)
+        sw.Stop()
+        if prettyPrint then printfn "\tApply.Map (MKL): %d milliseconds." sw.ElapsedMilliseconds
+        else printf "\t%d" sw.ElapsedMilliseconds
+        sw.Reset()*)
+    
+    printfn ""

Backup/src/FSharpExamples/DenseVector.fs

@@ -0,0 +1,54 @@
+// <copyright file="DenseVector.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+module MathNet.Numerics.FSharp.Examples.DenseVector
+
+open MathNet.Numerics.FSharp
+open MathNet.Numerics.LinearAlgebra
+
+// Create a new 100 dimensional dense vector.
+let v = Double.DenseVector.init 100 (fun i -> float i / 100.0)
+
+// Another way to create a 100 dimensional dense vector is using the vector function.
+let w = vector (List.init 100 (fun i -> float i ** 2.0))
+
+// Vectors can also be constructed from sequences.
+let t = Double.DenseVector.ofSeq (seq { for i in 1 .. 100 do yield float i })
+
+// We can now add two vectors together ...
+let z = v + w
+
+// ... or scale them in the process.
+let x = v + 3.0 * t
+
+
+
+// We can create a vector from an integer range (in this case, 5 and 10 inclusive) ...
+let s = Double.DenseVector.range 5 10
+
+// ... or we can create a vector from a double range with a particular step size.
+let r = Double.DenseVector.rangef 0.0 0.1 10.0

Backup/src/FSharpExamples/FSharpExamples.fsproj

@@ -0,0 +1,77 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="4.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <PropertyGroup>
+    <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
+    <Platform Condition=" '$(Platform)' == '' ">AnyCPU</Platform>
+    <ProductVersion>8.0.30703</ProductVersion>
+    <SchemaVersion>2.0</SchemaVersion>
+    <ProjectGuid>{bc81ea37-8ee6-4bf9-b8a9-b30497aef8b1}</ProjectGuid>
+    <OutputType>Library</OutputType>
+    <RootNamespace>FSharpExamples</RootNamespace>
+    <AssemblyName>FSharpExamples</AssemblyName>
+    <TargetFrameworkVersion>v4.0</TargetFrameworkVersion>
+    <Name>FSharpExamples</Name>
+    <TargetFrameworkProfile />
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
+    <DebugSymbols>true</DebugSymbols>
+    <DebugType>full</DebugType>
+    <Optimize>false</Optimize>
+    <Tailcalls>false</Tailcalls>
+    <OutputPath>bin\Debug\</OutputPath>
+    <DefineConstants>DEBUG;TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <Tailcalls>true</Tailcalls>
+    <OutputPath>bin\Release\</OutputPath>
+    <DefineConstants>TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release-Signed|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <Tailcalls>true</Tailcalls>
+    <DefineConstants>TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+    <OutputPath>bin\Release-Signed\</OutputPath>
+  </PropertyGroup>
+  <ItemGroup>
+    <Reference Include="FSharp.Core" />
+    <Reference Include="mscorlib" />
+    <Reference Include="System" />
+    <Reference Include="System.Core">
+      <RequiredTargetFramework>3.5</RequiredTargetFramework>
+    </Reference>
+    <Reference Include="System.Numerics" />
+  </ItemGroup>
+  <ItemGroup>
+    <ProjectReference Include="..\FSharp\FSharp.fsproj">
+      <Name>FSharp</Name>
+      <Project>{37e8e802-a354-4114-bfc1-6e1357da605b}</Project>
+      <Private>True</Private>
+    </ProjectReference>
+    <ProjectReference Include="..\Numerics\Numerics.csproj">
+      <Name>Numerics</Name>
+      <Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
+      <Private>True</Private>
+    </ProjectReference>
+  </ItemGroup>
+  <ItemGroup>
+    <Compile Include="DenseVector.fs" />
+    <Compile Include="Apply.fs" />
+    <Compile Include="Histogram.fs" />
+    <Compile Include="MCMC.fs" />
+  </ItemGroup>
+  <Import Project="$(MSBuildExtensionsPath32)\FSharp\1.0\Microsoft.FSharp.Targets" Condition="!Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <Import Project="$(MSBuildExtensionsPath32)\..\Microsoft F#\v4.0\Microsoft.FSharp.Targets" Condition="Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
+	     Other similar extension points exist, see Microsoft.Common.targets.
+	<Target Name="BeforeBuild">
+	</Target>
+	<Target Name="AfterBuild">
+	</Target>
+	-->
+</Project>

Backup/src/FSharpExamples/Histogram.fs

@@ -0,0 +1,44 @@
+// <copyright file="Histogram.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+module MathNet.Numerics.FSharp.Examples.Histogram
+
+open MathNet.Numerics.Statistics
+
+/// The number of buckets to use in our histogram.
+let B = 9
+
+/// Create a small dataset.
+let data = [| 0.5; 1.5; 2.5; 3.5; 4.5; 5.5; 6.5; 7.5; 8.5; 9.5 |]
+
+/// A histogram with 9 buckets for this dataset.
+let hist = new Histogram(data, B)
+
+// Print some histogram information.
+printfn "Histogram.ToString(): %O" hist
+for i in 0 .. B-1 do
+    printfn "Bucket %d contains %f datapoints." i hist.[i].Count

Backup/src/FSharpExamples/MCMC.fs

@@ -0,0 +1,201 @@
+// <copyright file="MCMC.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+module MathNet.Numerics.FSharp.Examples.Mcmc
+
+open MathNet.Numerics.Random
+open MathNet.Numerics.Statistics
+open MathNet.Numerics.Distributions
+open MathNet.Numerics.Statistics.Mcmc
+
+/// The number of samples to gather for each sampler.
+let N = 10000
+/// The random number generator we use for the examples.
+let rnd = new MersenneTwister()
+
+//
+// Example 1: Sampling a Beta distributed variable through rejection sampling.
+//
+// Target Distribution: Beta(2.7, 6.3)
+//
+// -----------------------------------------------------------------------------
+do
+    printfn "Rejection Sampling Example"
+    
+    /// The target distribution.
+    let beta = new Beta(2.7, 6.3)
+
+    /// Samples uniform distributed variables.
+    let uniform = new ContinuousUniform(0.0, 1.0, RandomSource = rnd)
+
+    /// Implements the rejection sampling procedure.
+    let rs = new RejectionSampler<float>( ( fun x -> x**(beta.A-1.0) * (1.0 - x)**(beta.B-1.0) ),
+                                          ( fun x -> 0.021 ),
+                                          ( fun () -> uniform.Sample()) )
+
+    /// An array of samples from the rejection sampler.
+    let arr = rs.Sample(N)
+
+    /// The true distribution.
+    printfn "\tEmpirical Mean = %f (should be %f)" (Statistics.Mean(arr)) beta.Mean
+    printfn "\tEmpirical StdDev = %f (should be %f)" (Statistics.StandardDeviation(arr)) beta.StdDev
+    printfn "\tAcceptance rate = %f" rs.AcceptanceRate
+    printfn ""
+
+
+
+//
+// Example 2: Sampling a normal distributed variable through Metropolis sampling.
+//
+// Target Distribution: Normal(1.0, 3.5)
+//
+// -----------------------------------------------------------------------------
+do
+    printfn "Metropolis Sampling Example"
+    
+    let mean, stddev = 1.0, 3.5
+    let normal = new Normal(mean, stddev)
+
+    /// Implements the rejection sampling procedure.
+    let ms = new MetropolisSampler<float>( 0.1, (fun x -> log(normal.Density(x))),
+                                                (fun x -> Normal.Sample(rnd, x, 0.3)), 20,
+                                                RandomSource = rnd )
+
+    /// An array of samples from the rejection sampler.
+    let arr = ms.Sample(N)
+
+    /// The true distribution.
+    printfn "\tEmpirical Mean = %f (should be %f)" (Statistics.Mean(arr)) normal.Mean
+    printfn "\tEmpirical StdDev = %f (should be %f)" (Statistics.StandardDeviation(arr)) normal.StdDev
+    printfn "\tAcceptance rate = %f" ms.AcceptanceRate
+    printfn ""
+    
+
+
+//
+// Example 3: Sampling a normal distributed variable through Metropolis-Hastings sampling
+//              with a symmetric proposal distribution.
+//
+// Target Distribution: Normal(1.0, 3.5)
+//
+// -----------------------------------------------------------------------------------------
+do
+    printfn "Metropolis Hastings Sampling Example (Symmetric Proposal)"
+    let mean, stddev = 1.0, 3.5
+    let normal = new Normal(mean, stddev)
+    
+    /// Evaluates the log normal distribution.
+    let npdf x m s = -0.5*(x-m)*(x-m)/(s*s) - 0.5 * log(2.0 * System.Math.PI * s * s)
+
+    /// Implements the rejection sampling procedure.
+    let ms = new MetropolisHastingsSampler<float>( 0.1, (fun x -> log(normal.Density(x))), 
+                    (fun x y -> npdf x y 0.3), (fun x -> Normal.Sample(rnd, x, 0.3)), 10,
+                                                RandomSource = rnd )
+
+    /// An array of samples from the rejection sampler.
+    let arr = ms.Sample(N)
+
+    /// The true distribution.
+    printfn "\tEmpirical Mean = %f (should be %f)" (Statistics.Mean(arr)) normal.Mean
+    printfn "\tEmpirical StdDev = %f (should be %f)" (Statistics.StandardDeviation(arr)) normal.StdDev
+    printfn "\tAcceptance rate = %f" ms.AcceptanceRate
+    printfn ""
+
+
+
+//
+// Example 4: Sampling a normal distributed variable through Metropolis-Hastings sampling
+//              with a asymmetric proposal distribution.
+//
+// Target Distribution: Normal(1.0, 3.5)
+//
+// -----------------------------------------------------------------------------------------
+do
+    printfn "Metropolis Hastings Sampling Example (Assymetric Proposal)"
+    let mean, stddev = 1.0, 3.5
+    let normal = new Normal(mean, stddev)
+    
+    /// Evaluates the logarithm of the normal distribution function.
+    let npdf x m s = -0.5*(x-m)*(x-m)/(s*s) - 0.5 * log(2.0 * System.Math.PI * s * s)
+    
+    /// Samples from a mixture that is biased towards samples larger than x.
+    let mixSample x =
+        if Bernoulli.Sample(rnd, 0.5) = 1 then
+            Normal.Sample(rnd, x, 0.3)
+        else
+            Normal.Sample(rnd, x + 0.1, 0.3)
+    
+    /// The transition kernel for the proposal above.
+    let krnl xnew x = log (0.5 * exp(npdf xnew x 0.3) + 0.5 * exp(npdf xnew (x+0.1) 0.3))
+
+    /// Implements the rejection sampling procedure.
+    let ms = new MetropolisHastingsSampler<float>( 0.1, (fun x -> log(normal.Density(x))), 
+                    (fun xnew x -> krnl xnew x), (fun x -> mixSample x), 10,
+                                                RandomSource = rnd )
+
+    /// An array of samples from the rejection sampler.
+    let arr = ms.Sample(N)
+
+    /// The true distribution.
+    printfn "\tEmpirical Mean = %f (should be %f)" (Statistics.Mean(arr)) normal.Mean
+    printfn "\tEmpirical StdDev = %f (should be %f)" (Statistics.StandardDeviation(arr)) normal.StdDev
+    printfn "\tAcceptance rate = %f" ms.AcceptanceRate
+    printfn ""
+
+
+
+//
+// Example 5: Slice sampling a normal distributed random variable.
+//
+// Target Distribution: Normal(1.0, 3.5)
+//
+// -----------------------------------------------------------------------------------------
+do
+    printfn "Slice Sampling Example"
+    let mean, stddev = 1.0, 3.5
+    let normal = new Normal(mean, stddev)
+    
+    /// Evaluates the unnormalized logarithm of the normal distribution function.
+    let npdf x m s = -0.5*(x-m)*(x-m)/(s*s)
+
+    /// Implements the rejection sampling procedure.
+    let ms = new UnivariateSliceSampler( 0.1, (fun x -> npdf x mean stddev), 5, 1.0, RandomSource = rnd )
+
+    /// An array of samples from the rejection sampler.
+    let arr = ms.Sample(N)
+
+    /// The true distribution.
+    printfn "\tEmpirical Mean = %f (should be %f)" (Statistics.Mean(arr)) normal.Mean
+    printfn "\tEmpirical StdDev = %f (should be %f)" (Statistics.StandardDeviation(arr)) normal.StdDev
+    printfn ""
+
+
+
+
+
+
+System.Console.ReadLine() |> ignore

Backup/src/FSharpUnitTests/App.config

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="utf-8"?>
+<configuration>
+  <startup>
+    <supportedRuntime version="v4.0" sku=".NETFramework,Version=v4.0" />
+  </startup>
+</configuration>

Backup/src/FSharpUnitTests/FSharpUnitTests.fsproj

@@ -0,0 +1,76 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="4.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <PropertyGroup>
+    <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
+    <Platform Condition=" '$(Platform)' == '' ">AnyCPU</Platform>
+    <ProductVersion>8.0.30703</ProductVersion>
+    <SchemaVersion>2.0</SchemaVersion>
+    <ProjectGuid>{f2f8032b-a31d-4e33-a05e-f2cdcbfaa75d}</ProjectGuid>
+    <OutputType>Exe</OutputType>
+    <RootNamespace>FSharpUnitTests</RootNamespace>
+    <AssemblyName>MathNet.Numerics.FSharp.UnitTests</AssemblyName>
+    <TargetFrameworkVersion>v4.0</TargetFrameworkVersion>
+    <Name>FSharpUnitTests</Name>
+    <TargetFrameworkProfile />
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
+    <DebugSymbols>true</DebugSymbols>
+    <DebugType>full</DebugType>
+    <Optimize>false</Optimize>
+    <Tailcalls>false</Tailcalls>
+    <OutputPath>..\..\out\test\debug\Net40\</OutputPath>
+    <DefineConstants>DEBUG;TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <Tailcalls>true</Tailcalls>
+    <OutputPath>..\..\out\test\Net40\</OutputPath>
+    <DefineConstants>TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release-Signed|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <Tailcalls>true</Tailcalls>
+    <DefineConstants>TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+    <OutputPath>bin\Release-Signed\</OutputPath>
+  </PropertyGroup>
+  <Import Project="$(MSBuildExtensionsPath32)\FSharp\1.0\Microsoft.FSharp.Targets" Condition="!Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <Import Project="$(MSBuildExtensionsPath32)\..\Microsoft F#\v4.0\Microsoft.FSharp.Targets" Condition="Exists('$(MSBuildToolsPath)\Microsoft.Build.Tasks.v4.0.dll')" />
+  <ItemGroup>
+    <Compile Include="FsUnit.fs" />
+    <Compile Include="Program.fs" />
+    <None Include="App.config">
+      <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+    </None>
+  </ItemGroup>
+  <ItemGroup>
+    <Reference Include="FSharp.Core" />
+    <Reference Include="mscorlib" />
+    <Reference Include="System" />
+    <Reference Include="System.Core">
+      <RequiredTargetFramework>3.5</RequiredTargetFramework>
+    </Reference>
+    <ProjectReference Include="..\FSharp\FSharp.fsproj">
+      <Name>FSharp</Name>
+      <Project>{37e8e802-a354-4114-bfc1-6e1357da605b}</Project>
+      <Private>True</Private>
+    </ProjectReference>
+    <ProjectReference Include="..\Numerics\Numerics.csproj">
+      <Name>Numerics</Name>
+      <Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
+      <Private>True</Private>
+    </ProjectReference>
+    <Reference Include="System.Numerics" />
+  </ItemGroup>
+  <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
+	     Other similar extension points exist, see Microsoft.Common.targets.
+	<Target Name="BeforeBuild">
+	</Target>
+	<Target Name="AfterBuild">
+	</Target>
+	-->
+</Project>

Backup/src/FSharpUnitTests/FsUnit.fs

@@ -0,0 +1,244 @@
+(* 
+Copyright (c) 2008, Raymond W. Vernagus (R.Vernagus@gmail.com)
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+  * Redistributions of source code must retain the above copyright notice,
+    this list of conditions and the following disclaimer.
+  * Redistributions in binary form must reproduce the above copyright notice,
+    this list of conditions and the following disclaimer in the documentation
+    and/or other materials provided with the distribution.
+  * Neither the name of Raymond W. Vernagus nor the names of FsUnit's
+    contributors may be used to endorse or promote products derived from this
+    software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*)
+namespace FsUnit
+
+open MathNet.Numerics
+open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.LinearAlgebra.Generic
+
+type Result =
+    | Pass
+    | Fail of string
+    | Error of string
+
+type Expectation =
+    | True
+    | False
+    | Empty
+    | NullOrEmpty
+    | Null
+    | SameAs of obj
+
+type Spec =
+    abstract Check : unit -> Result
+    abstract NegatedMessage : string
+
+[<AutoOpen>]
+module SpecOps =
+    let internal safeCheck f =
+        try
+            f()
+        with ex -> Error (ex.ToString())
+    
+    let make f nmsg =
+        { new Spec with
+            member this.Check() = safeCheck f
+            member this.NegatedMessage = nmsg }
+    
+    let check (s: Spec) = s.Check()
+    
+    let not' f x =
+        let s = f x
+        match check s with
+        | Pass      -> make (fun () -> Fail s.NegatedMessage) ""
+        | Fail msg  -> make (fun () -> Pass) msg
+        | Error msg -> make (fun () -> Error msg) msg
+    
+    let equal expected actual =
+        make (fun () ->
+            if actual.Equals(expected)
+                then Pass
+                else Fail (sprintf "Expected: %A\nActual: %A" expected actual))
+            (sprintf "NOT Expected: %A\nActual: %A" expected actual)
+
+    let have n lbl s =
+        make (fun() ->
+            let len = Seq.length s
+            if len = n
+                then Pass
+                else Fail (sprintf "Expected: %d %s\nActual: %d %s" n lbl len lbl))
+            (sprintf "NOT Expected: %d %s\nActual: %d %s" n lbl n lbl)
+    
+    let contain x s =
+        make (fun () ->
+            let exists = s |> Seq.exists (fun x' -> x' = x)
+            if exists
+                then Pass
+                else Fail (sprintf "Expected %A to contain %A" s x))
+            (sprintf "Did NOT expect %A to contain %A" s x)
+
+    let raise'<'a when 'a :> System.Exception> f =
+        let expType = typeof<'a>
+        make (fun () -> 
+            try
+                f()
+                Fail (sprintf "Expected %s but no exception was raised" expType.FullName)
+            with ex ->
+                let actualExType = ex.GetType()
+                if expType = actualExType
+                    then Pass
+                    else Fail (sprintf "Expected: %s\nActual: %s" expType.FullName actualExType.FullName))
+            (sprintf "Did NOT expect %s to be raised" expType.FullName)
+
+    let be expectation x =
+        let x = box x
+        let msg = sprintf "Expected: %A\nActual: %A" expectation x
+        let negmsg = sprintf "NOT Expected: %A\nActual: %A" expectation x
+        match expectation with
+        | True ->
+            make (fun () ->
+                if x = box true
+                    then Pass
+                    else Fail msg)
+                negmsg
+        | False ->
+            make (fun () ->
+                if x = box false
+                    then Pass
+                    else Fail msg)
+                negmsg
+        | Empty ->
+            make (fun () ->
+                if x = box System.String.Empty
+                    then Pass
+                    else Fail msg)
+                negmsg
+        | NullOrEmpty ->
+            make (fun () ->
+                if System.String.IsNullOrEmpty(x :?> string)
+                    then Pass
+                    else Fail msg)
+                negmsg
+        | Null ->
+            make (fun () ->
+                if x = null
+                    then Pass
+                    else Fail msg)
+                negmsg
+        | SameAs other ->
+            make (fun () ->
+                if System.Object.ReferenceEquals(x, other)
+                    then Pass
+                    else Fail (sprintf "Expected actual to be same reference as expected %A" other))
+                (sprintf "Expected %A to have different reference than %A" x other)
+
+    let array_equal (expected: float []) (actual: float []) =
+        make (fun () ->
+            let mutable f = true
+            for i=0 to expected.Length-1 do
+                f <- f && (expected.[i] = actual.[i])
+            if f
+                then Pass
+                else Fail (sprintf "Expected: %A\nActual: %A" expected actual))
+            (sprintf "NOT Expected: %A\nActual: %A" expected actual)
+
+    let array2_equal (expected: float [,]) (actual: float [,]) =
+        make (fun () ->
+            let mutable f = true
+            for i=0 to expected.GetLength(0)-1 do
+                for j=0 to expected.GetLength(1)-1 do
+                    f <- f && (expected.[i,j] = actual.[i,j])
+            if f
+                then Pass
+                else Fail (sprintf "Expected: %A\nActual: %A" expected actual))
+            (sprintf "NOT Expected: %A\nActual: %A" expected actual)
+
+    let approximately_equal (places : int) (expected: float) (actual: float) =
+        make (fun () ->
+            if Precision.AlmostEqualInDecimalPlaces(actual, expected, places)
+                then Pass
+                else Fail (sprintf "Expected: %A\nActual: %A" expected actual))
+            (sprintf "NOT Expected: %A\nActual: %A" expected actual)
+            
+    let approximately_vector_equal (places : int) (expected: #Vector<float>) (actual: #Vector<float>) =
+        make (fun () ->
+            let mutable f = true
+            for i=0 to expected.Count-1 do
+                f <- f && Precision.AlmostEqualInDecimalPlaces(expected.[i], actual.[i], places)
+            if f
+                then Pass
+                else Fail (sprintf "Expected: %A\nActual: %A" expected actual))
+            (sprintf "NOT Expected: %A\nActual: %A" expected actual)
+            
+    let approximately_matrix_equal (places : int) (expected: #Matrix<float>) (actual: #Matrix<float>) =
+        make (fun () ->
+            let mutable f = true
+            for i=0 to expected.RowCount-1 do
+              for j=0 to expected.ColumnCount-1 do
+                f <- f && Precision.AlmostEqualInDecimalPlaces(expected.[i,j], actual.[i,j], places)
+            if f
+                then Pass
+                else Fail (sprintf "Expected: %A\nActual: %A" expected actual))
+            (sprintf "NOT Expected: %A\nActual: %A" expected actual)
+
+module Results =
+    open Microsoft.FSharp.Text
+    
+    let internal currentResults = new ResizeArray<string * Result>()
+    
+    let add = currentResults.Add
+    
+    let passedCount () =
+        currentResults
+        |> Seq.filter (function _,Pass -> true | _ -> false)
+        |> Seq.length
+    
+    let failed () =
+        currentResults
+        |> Seq.filter (function _,Fail _ -> true | _ -> false)
+    
+    let failedCount () =
+        failed()
+        |> Seq.length
+    
+    let erred () =
+        currentResults
+        |> Seq.filter (function _,Error _ -> true | _ -> false)
+    
+    let erredCount () =
+        erred()
+        |> Seq.length
+        
+    let summary () =
+        let buff = new System.Text.StringBuilder()
+        buff.AppendFormat("{0} passed.\n{1} failed.\n{2} erred.", passedCount(), failedCount(), erredCount()) |> ignore
+        failed()
+        |> Seq.iter (function (lbl,Fail msg) -> buff.AppendFormat("\n----\nFailed: {0}\n{1}", lbl, msg) |> ignore | _ -> ())
+        erred()
+        |> Seq.iter (function (lbl,Error msg) -> buff.AppendFormat("\n----\nErred: {0}\n{1}", lbl, msg) |> ignore | _ -> ())
+        buff.ToString()
+    
+
+[<AutoOpen>]
+module SpecHelpers =
+    let spec lbl s = (lbl, check s)
+    
+    let should f x = f x
+    
+    let specs lbl (results: seq<string * Result>) =
+        results |> Seq.iter (fun x -> Results.add x)

Backup/src/FSharpUnitTests/Program.fs

@@ -0,0 +1,224 @@
+open FsUnit
+open MathNet.Numerics.FSharp
+open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.LinearAlgebra.Generic
+
+/// Unit tests for the dense vector type.
+let DenseVectorTests =
+
+    /// A small uniform vector.
+    let smallv = new DenseVector(5, 0.3 )
+
+    /// A large vector with increasingly large entries
+    let largev = new DenseVector( Array.init 100 (fun i -> float i / 100.0) )
+    
+    specs "DenseVector" [
+        spec "DenseVector.init"
+            (DenseVector.init 100 (fun i -> float i / 100.0) |> should equal largev)
+        spec "DenseVector.ofList"
+            (DenseVector.ofList [ for i in 0 .. 99 -> float i / 100.0 ] |> should equal largev)
+        spec "DenseVector.ofSeq"
+            (DenseVector.ofSeq (seq { for i in 0 .. 99 -> float i / 100.0 }) |> should equal largev)
+        spec "DenseVector.rangef"
+            (DenseVector.rangef 0.0 0.01 0.99 |> should equal (new DenseVector( [| for i in 0 .. 99 -> 0.01 * float i |] ) ))
+        spec "DenseVector.range"
+            (DenseVector.range 0 99 |> should equal (new DenseVector( [| for i in 0 .. 99 -> float i |] ) ))
+    ]
+
+
+/// Unit tests for the sparse vector type.
+let SparseVectorTests =
+
+    /// A small uniform vector.
+    let smallv = new DenseVector( [|0.0;0.3;0.0;0.0;0.0|] ) :> Vector<float>
+    
+    specs "SparseVector" [
+        spec "SparseVector.ofList"
+            ((SparseVector.ofList 5 [ (1,0.3) ] :> Vector<float>) |> should equal smallv)
+        spec "SparseVector.ofSeq"
+            ((SparseVector.ofSeq 5 (List.toSeq [ (1,0.3) ]) :> Vector<float>) |> should equal smallv)
+    ]
+
+
+/// Unit tests for the vector type.
+let VectorTests =
+
+    /// A small uniform vector.
+    let smallv = new DenseVector( [|0.3;0.3;0.3;0.3;0.3|] ) :> Vector<float>
+
+    /// A large vector with increasingly large entries
+    let largev = new DenseVector( Array.init 100 (fun i -> float i / 100.0) ) :> Vector<float>
+    
+    specs "Vector" [
+        spec "Vector.toArray"
+            (Vector.toArray smallv |> should array_equal [|0.3;0.3;0.3;0.3;0.3|])
+        spec "Vector.toList"
+            (Vector.toList smallv |> should equal [0.3;0.3;0.3;0.3;0.3])
+        spec "Vector.mapInPlace"
+            ( let w = smallv.Clone()
+              Vector.mapInPlace (fun x -> 2.0 * x) w
+              w |> should equal (2.0 * smallv))
+        spec "Vector.mapiInPlace"
+            ( let w = largev.Clone()
+              Vector.mapiInPlace (fun i x -> float i / 100.0) w
+              w |> should equal (largev))
+        spec "Vector.addInPlace"
+            ( let w = largev.Clone()
+              Vector.addInPlace w largev
+              w |> should equal (2.0 * largev))
+        spec "Vector.subInPlace"
+            ( let w = largev.Clone()
+              Vector.subInPlace w largev
+              w |> should equal (0.0 * largev))
+        spec "Vector.map"
+            (Vector.map (fun x -> 2.0 * x) largev |> should equal (2.0 * largev))
+        spec "Vector.mapi"
+            (Vector.mapi (fun i x -> float i / 100.0) largev |> should equal largev)
+        spec "Vector.fold"
+            (Vector.fold (fun a b -> a - b) 0.0 smallv |> should equal -1.5)
+        spec "Vector.foldBack"
+            (Vector.foldBack (fun a b -> a - b) 0.0 smallv |> should equal 0.0)
+        spec "Vector.foldi"
+            (Vector.foldi (fun i a b -> a + b) 0.0 smallv |> should equal 1.5)
+        spec "Vector.forall"
+            (Vector.forall (fun x -> x = 0.3) smallv |> should equal true)
+        spec "Vector.exists"
+            (Vector.exists (fun x -> x = 0.3) smallv |> should equal true)
+        spec "Vector.foralli"
+            (Vector.foralli (fun i x -> x = 0.3 && i < 5) smallv |> should equal true)
+        spec "Vector.existsi"
+            (Vector.existsi (fun i x -> x = 0.3 && i = 2) smallv |> should equal true)
+        spec "Vector.scan"
+            (Vector.scan (fun acc x -> acc + x) smallv |> should approximately_vector_equal 14 (new DenseVector( [|0.3;0.6;0.9;1.2;1.5|] ) :> Vector<float>) )
+        spec "Vector.scanBack"
+            (Vector.scanBack (fun x acc -> acc + x) smallv |> should approximately_vector_equal 14 (new DenseVector( [|1.5;1.2;0.9;0.6;0.3|] ) :> Vector<float>) )
+        spec "Vector.reduce"
+            (Vector.reduce (fun acc x -> acc ** x) smallv |> should approximately_equal 14 0.990295218585507)
+        spec "Vector.reduceBack"
+            (Vector.reduceBack (fun x acc -> x ** acc) smallv |> should approximately_equal 14 0.488911287726319)
+        spec "Vector.insert"
+            (Vector.insert 2 0.5 smallv |> should approximately_vector_equal 14 (new DenseVector ( [|0.3;0.3;0.5;0.3;0.3;0.3|] ) :> Vector<float>) )
+    ]
+
+
+/// Unit tests for the matrix type.
+let MatrixTests =
+
+    /// A small uniform vector.
+    let smallM = new DenseMatrix( Array2D.create 2 2 0.3 )
+    let failingFoldBackM = DenseMatrix.init 2 3 (fun i j -> 1.0)
+
+
+    /// A large vector with increasingly large entries
+    let largeM = new DenseMatrix( Array2D.init 100 100 (fun i j -> float i * 100.0 + float j) )
+    
+    specs "Matrix" [
+        spec "Matrix.fold"
+            (Matrix.fold (fun a b -> a - b) 0.0 smallM |> should equal -1.2)
+        spec "Matrix.foldBack"
+            (Matrix.foldBack (fun a b -> a - b) 0.0 smallM |> should equal 0.0)
+        spec "Matrix.foldBackSummation"
+            (Matrix.foldBack( fun a b -> a + b) 0.0 failingFoldBackM |> should equal 6.0)
+        spec "Matrix.foldi"
+            (Matrix.foldi (fun i j acc x -> acc + x + float (i+j)) 0.0 smallM |> should equal 5.2)
+        spec "Matrix.toArray2"
+            (Matrix.toArray2 smallM |> should array2_equal (Array2D.create 2 2 0.3))
+        spec "Matrix.forall"
+            (Matrix.forall (fun x -> x = 0.3) smallM |> should equal true)
+        spec "Matrix.exists"
+            (Matrix.exists (fun x -> x = 0.5) smallM |> should equal false)
+        spec "Matrix.foralli"
+            (Matrix.foralli (fun i j x -> x = float i * 100.0 + float j) largeM |> should equal true)
+        spec "Matrix.existsi"
+            (Matrix.existsi (fun i j x -> x = float i * 100.0 + float j) largeM |> should equal true)
+        spec "Matrix.map"
+            (Matrix.map (fun x -> 2.0 * x) smallM |> should equal (2.0 * smallM))
+        spec "Matrix.mapi"
+            (Matrix.mapi (fun i j x -> float i * 100.0 + float j + x) largeM |> should equal (2.0 * largeM))
+        spec "Matrix.mapCols"
+            (Matrix.mapCols (fun j col -> col.Add(float j)) smallM |> should approximately_matrix_equal 14 (matrix [[0.3;1.3];[0.3;1.3]]))
+        spec "Matrix.mapRows"
+            (Matrix.mapRows (fun i row -> row.Add(float i)) smallM |> should approximately_matrix_equal 14 (matrix [[0.3;0.3];[1.3;1.3]]))
+        spec "Matrix.inplaceAssign"
+            ( let N = smallM.Clone()
+              Matrix.inplaceAssign (fun i j -> 0.0) N
+              N |> should equal (0.0 * smallM))
+        spec "Matrix.inplaceMapi"
+            ( let N = largeM.Clone()
+              Matrix.inplaceMapi (fun i j x -> 2.0 * (float i * 100.0 + float j) + x) N
+              N |> should equal (3.0 * largeM))
+        spec "Matrix.nonZeroEntries"
+            (Seq.length (Matrix.nonZeroEntries smallM) |> should equal 4)
+        spec "Matrix.sum"
+            (Matrix.sum smallM |> should equal 1.2)
+        spec "Matrix.sumColsBy"
+            (Matrix.sumColsBy (fun j col -> col.[0] * col.[1]) (matrix [[1.0; 2.0]; [3.0; 4.0]]) |> should equal 11.0)
+        spec "Matrix.sumRowsBy"
+            (Matrix.sumRowsBy (fun i row -> row.[0] * row.[1]) (matrix [[1.0; 2.0]; [3.0; 4.0]]) |> should equal 14.0)
+        spec "Matrix.foldCol"
+            (Matrix.foldCol (+) 0.0 largeM 0 |> should equal 495000.0)
+        spec "Matrix.foldRow"
+            (Matrix.foldRow (+) 0.0 largeM 0 |> should equal 4950.0)
+        spec "Matrix.foldByCol"
+            (Matrix.foldByCol (+) 0.0 smallM |> should equal (DenseVector.ofList [0.6;0.6] :> Vector<float>))
+        spec "Matrix.foldByRow"
+            (Matrix.foldByRow (+) 0.0 smallM |> should equal (DenseVector.ofList [0.6;0.6] :> Vector<float>))
+    ]
+
+
+/// Unit tests for the dense matrix type.
+let DenseMatrixTests =
+
+    /// A small uniform vector.
+    let smallM = new DenseMatrix( Array2D.create 2 2 0.3 )
+
+    /// A large vector with increasingly large entries
+    let largeM = new DenseMatrix( Array2D.init 100 100 (fun i j -> float i * 100.0 + float j) )
+    
+    specs "DenseMatrix" [
+        spec "DenseMatrix.init"
+            (DenseMatrix.init 100 100 (fun i j -> float i * 100.0 + float j) |> should equal largeM)
+        spec "DenseMatrix.ofList"
+            (DenseMatrix.ofList [[0.3;0.3];[0.3;0.3]] |> should equal smallM)
+        spec "DenseMatrix.ofSeq"
+            (DenseMatrix.ofSeq (Seq.ofList [[0.3;0.3];[0.3;0.3]]) |> should equal smallM)
+        spec "DenseMatrix.ofArray2"
+            (DenseMatrix.ofArray2 (Array2D.create 2 2 0.3) |> should equal smallM)
+        spec "DenseMatrix.initDense"
+            (DenseMatrix.initDense 100 100 (seq { for i in 0 .. 99 do
+                                                   for j in 0 .. 99 -> (i,j, float i * 100.0 + float j)}) |> should equal largeM)
+        spec "DenseMatrix.constDiag"
+            (DenseMatrix.constDiag 100 2.0 |> should equal (2.0 * (DenseMatrix.Identity 100)))
+        spec "DenseMatrix.diag"
+            (DenseMatrix.diag (new DenseVector(100, 2.0)) |> should equal (2.0 * (DenseMatrix.Identity 100)))
+        spec "DenseMatrix.init_row"
+            (DenseMatrix.initRow 100 100 (fun i -> (DenseVector.init 100 (fun j -> float i * 100.0 + float j))) |> should equal largeM)
+        spec "DenseMatrix.init_col"
+            (DenseMatrix.initCol 100 100 (fun j -> (DenseVector.init 100 (fun i -> float i * 100.0 + float j))) |> should equal largeM)
+    ]
+    
+
+/// Unit tests for the sparse matrix type.
+let SparseMatrixTests =
+
+    /// A small uniform vector.
+    let smallM = DenseMatrix.init 4 4 (fun i j -> if i = 1 && j = 2 then 1.0 else 0.0) :> Matrix<float>
+    
+    specs "SparseMatrix" [
+        spec "SparseMatrix.ofList"
+            ((SparseMatrix.ofList 4 4 [(1,2,1.0)] :> Matrix<float>) |> should equal smallM)
+        spec "SparseMatrix.ofSeq"
+            ((SparseMatrix.ofSeq 4 4 (Seq.ofList [(1,2,1.0)]) :> Matrix<float>) |> should equal smallM)
+        spec "SparseMatrix.constDiag"
+            (SparseMatrix.constDiag 100 2.0 |> should equal (2.0 * (SparseMatrix.Identity 100)))
+        spec "SparseMatrix.diag"
+            (SparseMatrix.diag (new DenseVector(100, 2.0)) |> should equal (2.0 * (SparseMatrix.Identity 100)))
+    ]
+
+    
+/// Report on errors and success and exit.
+printfn "F# Test Results:"
+printfn "%s" (Results.summary())
+    
+let code = if Results.erredCount() > 0 || Results.failedCount() > 0 then -1 else 0;;
+exit code;;

_UpgradeReport_Files/UpgradeReport.css

@@ -0,0 +1,118 @@
+/* Body style, for the entire document */
+body
+{
+    background: #F3F3F4;
+    color: #1E1E1F;
+    font-family: "Segoe UI", Tahoma, Geneva, Verdana, sans-serif;
+    padding: 0;
+    margin: 0;
+}
+
+/* Header1 style, used for the main title */
+h1
+{
+    padding: 10px 0px 10px 10px;
+    font-size: 21pt;
+    background-color: #E2E2E2;
+    border-bottom: 1px #C1C1C2 solid; 
+    color: #201F20;
+    margin: 0;
+    font-weight: normal;
+}
+
+/* Header2 style, used for "Overview" and other sections */
+h2
+{
+    font-size: 18pt;
+    font-weight: normal;
+    padding: 15px 0 5px 0;
+    margin: 0;
+}
+
+/* Header3 style, used for sub-sections, such as project name */
+h3
+{
+    font-weight: normal;
+    font-size: 15pt;
+    margin: 0;
+    padding: 15px 0 5px 0;
+    background-color: transparent;
+}
+
+/* Color all hyperlinks one color */
+a
+{
+    color: #1382CE;
+}
+
+/* Table styles */ 
+table
+{
+    border-spacing: 0 0;
+    border-collapse: collapse;
+    font-size: 10pt;
+}
+
+table th
+{
+    background: #E7E7E8;
+    text-align: left;
+    text-decoration: none;
+    font-weight: normal;
+    padding: 3px 6px 3px 6px;
+}
+
+table td
+{
+    vertical-align: top;
+    padding: 3px 6px 5px 5px;
+    margin: 0px;
+    border: 1px solid #E7E7E8;
+    background: #F7F7F8;
+}
+
+/* Local link is a style for hyperlinks that link to file:/// content, there are lots so color them as 'normal' text until the user mouse overs */
+.localLink
+{
+    color: #1E1E1F;
+    background: #EEEEED;
+    text-decoration: none;
+}
+
+.localLink:hover
+{
+    color: #1382CE;
+    background: #FFFF99;
+    text-decoration: none;
+}
+
+/* Center text, used in the over views cells that contain message level counts */ 
+.textCentered
+{
+    text-align: center;
+}
+
+/* The message cells in message tables should take up all avaliable space */
+.messageCell
+{
+    width: 100%;
+}
+
+/* Padding around the content after the h1 */ 
+#content 
+{
+	padding: 0px 12px 12px 12px; 
+}
+
+/* The overview table expands to width, with a max width of 97% */ 
+#overview table
+{
+    width: auto;
+    max-width: 75%; 
+}
+
+/* The messages tables are always 97% width */
+#messages table
+{
+    width: 97%;
+}

_UpgradeReport_Files/UpgradeReport.xslt

@@ -0,0 +1,558 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<xsl:stylesheet version="1.0" xmlns:xsl="http://www.w3.org/1999/XSL/Transform" xmlns:msxsl='urn:schemas-microsoft-com:xslt'>
+  <xsl:output omit-xml-declaration="yes" />
+
+  <!-- Keys -->
+  <xsl:key name="ProjectKey" match="Event" use="@Project" />
+
+  <!-- String split template -->
+  <xsl:template name="SplitString">
+    <xsl:param name="source"      select="''" />
+    <xsl:param name="separator" select="','" />
+    <xsl:if test="not($source = '' or $separator = '')">
+      <xsl:variable name="head" select="substring-before(concat($source, $separator), $separator)" />
+      <xsl:variable name="tail" select="substring-after($source, $separator)" />
+      <part>
+        <xsl:value-of select="$head"/>
+      </part>
+      <xsl:call-template name="SplitString">
+        <xsl:with-param name="source" select="$tail" />
+        <xsl:with-param name="separator" select="$separator" />
+      </xsl:call-template>
+    </xsl:if>
+  </xsl:template>
+
+  <!-- Intermediate Templates -->
+  <xsl:template match="UpgradeReport" mode="ProjectOverviewXML">
+    <Projects>
+      <xsl:for-each select="Event[generate-id(.) = generate-id(key('ProjectKey', @Project))]">
+        <Project>
+          <xsl:variable name="pNode" select="current()" />
+          <xsl:variable name="errorCount" select="count(../Event[@Project = current()/@Project and @ErrorLevel=2])" />
+          <xsl:variable name="warningCount" select="count(../Event[@Project = current()/@Project and @ErrorLevel=1])" />
+          <xsl:variable name="messageCount" select="count(../Event[@Project = current()/@Project and @ErrorLevel=0])" />
+          <xsl:variable name="pathSplitSeparator">
+            <xsl:text>\</xsl:text>
+          </xsl:variable>
+          <xsl:variable name="projectName">
+            <xsl:choose>
+              <xsl:when test="@Project = ''">Solution</xsl:when>
+              <xsl:otherwise>
+                <xsl:value-of select="@Project"/>
+              </xsl:otherwise>
+            </xsl:choose>
+          </xsl:variable>
+          <xsl:attribute name="IsSolution">
+            <xsl:value-of select="@Project = ''"/>
+          </xsl:attribute>
+          <xsl:attribute name="Project">
+            <xsl:value-of select="$projectName"/>
+          </xsl:attribute>
+          <xsl:attribute name="ProjectDisplayName">
+
+            <xsl:variable name="localProjectName" select="@Project" />
+
+            <!-- Sometimes it is possible to have project name set to a path over a real project name,
+                 we split the string on '\' and if we end up with >1 part in the resulting tokens set
+                 we format the ProjectDisplayName as ..\prior\last -->
+            <xsl:variable name="pathTokens">
+              <xsl:call-template name="SplitString">
+                <xsl:with-param name="source" select="$localProjectName" />
+                <xsl:with-param name="separator" select="$pathSplitSeparator" />
+              </xsl:call-template>
+            </xsl:variable>
+
+            <xsl:choose>
+              <xsl:when test="count(msxsl:node-set($pathTokens)/part) > 1">
+                <xsl:value-of select="concat('..', $pathSplitSeparator, msxsl:node-set($pathTokens)/part[last() - 1], $pathSplitSeparator, msxsl:node-set($pathTokens)/part[last()])"/>
+              </xsl:when>
+              <xsl:otherwise>
+                <xsl:value-of select="$localProjectName"/>
+              </xsl:otherwise>
+            </xsl:choose>
+
+          </xsl:attribute>
+          <xsl:attribute name="ProjectSafeName">
+            <xsl:value-of select="translate($projectName, '\', '-')"/>
+          </xsl:attribute>
+          <xsl:attribute name="Solution">
+            <xsl:value-of select="/UpgradeReport/Properties/Property[@Name='Solution']/@Value"/>
+          </xsl:attribute>
+          <xsl:attribute name="Source">
+            <xsl:value-of select="@Source"/>
+          </xsl:attribute>
+          <xsl:attribute name="Status">
+            <xsl:choose>
+              <xsl:when test="$errorCount > 0">Error</xsl:when>
+              <xsl:when test="$warningCount > 0">Warning</xsl:when>
+              <xsl:otherwise>Success</xsl:otherwise>
+            </xsl:choose>
+          </xsl:attribute>
+          <xsl:attribute name="ErrorCount">
+            <xsl:value-of select="$errorCount" />
+          </xsl:attribute>
+          <xsl:attribute name="WarningCount">
+            <xsl:value-of select="$warningCount" />
+          </xsl:attribute>
+          <xsl:attribute name="MessageCount">
+            <xsl:value-of select="$messageCount" />
+          </xsl:attribute>
+          <xsl:attribute name="TotalCount">
+            <xsl:value-of select="$errorCount + $warningCount + $messageCount"/>
+          </xsl:attribute>
+          <xsl:for-each select="../Event[@Project = $pNode/@Project and @ErrorLevel=3]">
+            <ConversionStatus>
+              <xsl:value-of select="@Description"/>
+            </ConversionStatus>
+          </xsl:for-each>
+          <Messages>
+            <xsl:for-each select="../Event[@Project = $pNode/@Project and @ErrorLevel&lt;3]">
+              <Message>
+                <xsl:attribute name="Level">
+                  <xsl:value-of select="@ErrorLevel" />
+                </xsl:attribute>
+                <xsl:attribute name="Status">
+                  <xsl:choose>
+                    <xsl:when test="@ErrorLevel = 0">Message</xsl:when>
+                    <xsl:when test="@ErrorLevel = 1">Warning</xsl:when>
+                    <xsl:when test="@ErrorLevel = 2">Error</xsl:when>
+                    <xsl:otherwise>Message</xsl:otherwise>
+                  </xsl:choose>
+                </xsl:attribute>
+                <xsl:attribute name="Source">
+                  <xsl:value-of select="@Source"/>
+                </xsl:attribute>
+                <xsl:attribute name="Message">
+                  <xsl:value-of select="@Description"/>
+                </xsl:attribute>
+              </Message>
+            </xsl:for-each>
+          </Messages>
+        </Project>
+      </xsl:for-each>
+    </Projects>
+  </xsl:template>
+
+
+
+  <!-- Project Overview template -->
+  <xsl:template match="Projects" mode="ProjectOverview">
+
+    <table>
+      <tr>
+        <th></th>
+        <th _locID="ProjectTableHeader">Project</th>
+        <th _locID="PathTableHeader">Path</th>
+        <th _locID="ErrorsTableHeader">Errors</th>
+        <th _locID="WarningsTableHeader">Warnings</th>
+        <th _locID="MessagesTableHeader">Messages</th>
+      </tr>
+
+      <xsl:for-each select="Project">
+
+        <xsl:sort select="@ErrorCount" order="descending" />
+        <xsl:sort select="@WarningCount" order="descending" />
+        <!-- Always make solution last within a group -->
+        <xsl:sort select="@IsSolution" order="ascending" />
+        <xsl:sort select="@ProjectSafeName" order="ascending" />
+
+        <tr>
+          <td>
+            <img width="16" height="16">
+              <xsl:attribute name="src">
+                <xsl:choose>
+                  <xsl:when test="@Status = 'Error'">_UpgradeReport_Files\UpgradeReport_Error.png</xsl:when>
+                  <xsl:when test="@Status = 'Warning'">_UpgradeReport_Files\UpgradeReport_Warning.png</xsl:when>
+                  <xsl:when test="@Status = 'Success'">_UpgradeReport_Files\UpgradeReport_Success.png</xsl:when>
+                </xsl:choose>
+              </xsl:attribute>
+              <xsl:attribute name="alt">
+                <xsl:value-of select="@Status" />
+              </xsl:attribute>
+            </img>
+          </td>
+          <td>
+            <strong>
+              <a>
+                <xsl:attribute name="href">
+                  <xsl:value-of select="concat('#', @ProjectSafeName)"/>
+                </xsl:attribute>
+                <xsl:choose>
+                  <xsl:when test="@ProjectDisplayName = ''">
+                    <span _locID="OverviewSolutionSpan">Solution</span>
+                  </xsl:when>
+                  <xsl:otherwise>
+                    <xsl:value-of select="@ProjectDisplayName" />
+                  </xsl:otherwise>
+                </xsl:choose>
+              </a>
+            </strong>
+          </td>
+          <td>
+            <xsl:value-of select="@Source" />
+          </td>
+          <td class="textCentered">
+            <a>
+              <xsl:if test="@ErrorCount > 0">
+                <xsl:attribute name="href">
+                  <xsl:value-of select="concat('#', @ProjectSafeName, 'Error')"/>
+                </xsl:attribute>
+              </xsl:if>
+              <xsl:value-of select="@ErrorCount" />
+            </a>
+          </td>
+          <td class="textCentered">
+            <a>
+              <xsl:if test="@WarningCount > 0">
+                <xsl:attribute name="href">
+                  <xsl:value-of select="concat('#', @ProjectSafeName, 'Warning')"/>
+                </xsl:attribute>
+              </xsl:if>
+              <xsl:value-of select="@WarningCount" />
+            </a>
+          </td>
+          <td class="textCentered">
+            <a href="#">
+              <xsl:if test="@MessageCount > 0">
+                <xsl:attribute name="onclick">
+                  <xsl:variable name="apos">
+                    <xsl:text>'</xsl:text>
+                  </xsl:variable>
+                  <xsl:variable name="JS" select="concat('ScrollToFirstVisibleMessage(', $apos, @ProjectSafeName, $apos, ')')" />
+                  <xsl:value-of select="concat($JS, '; return false;')"/>
+                </xsl:attribute>
+              </xsl:if>
+              <xsl:value-of select="@MessageCount" />
+            </a>
+          </td>
+        </tr>
+      </xsl:for-each>
+    </table>
+  </xsl:template>
+
+  <!-- Show messages row -->
+  <xsl:template match="Project" mode="ProjectShowMessages">
+    <tr>
+      <xsl:attribute name="name">
+        <xsl:value-of select="concat('MessageRowHeaderShow', @ProjectSafeName)"/>
+      </xsl:attribute>
+      <td>
+        <img width="16" height="16" src="_UpgradeReport_Files\UpgradeReport_Information.png" />
+      </td>
+      <td class="messageCell">
+        <xsl:variable name="apos">
+          <xsl:text>'</xsl:text>
+        </xsl:variable>
+        <xsl:variable name="toggleRowsJS" select="concat('ToggleMessageVisibility(', $apos, @ProjectSafeName, $apos, ')')" />
+
+        <a _locID="ShowAdditionalMessages" href="#">
+          <xsl:attribute name="name">
+            <xsl:value-of select="concat(@ProjectSafeName, 'Message')"/>
+          </xsl:attribute>
+          <xsl:attribute name="onclick">
+            <xsl:value-of select="concat($toggleRowsJS, '; return false;')"/>
+          </xsl:attribute>
+          Show <xsl:value-of select="@MessageCount" /> additional messages
+        </a>
+      </td>
+    </tr>
+  </xsl:template>
+
+  <!-- Hide messages row -->
+  <xsl:template match="Project" mode="ProjectHideMessages">
+    <tr style="display: none">
+      <xsl:attribute name="name">
+        <xsl:value-of select="concat('MessageRowHeaderHide', @ProjectSafeName)"/>
+      </xsl:attribute>
+      <td>
+        <img width="16" height="16" src="_UpgradeReport_Files\UpgradeReport_Information.png" />
+      </td>
+      <td class="messageCell">
+        <xsl:variable name="apos">
+          <xsl:text>'</xsl:text>
+        </xsl:variable>
+        <xsl:variable name="toggleRowsJS" select="concat('ToggleMessageVisibility(', $apos, @ProjectSafeName, $apos, ')')" />
+
+        <a _locID="HideAdditionalMessages" href="#">
+          <xsl:attribute name="name">
+            <xsl:value-of select="concat(@ProjectSafeName, 'Message')"/>
+          </xsl:attribute>
+          <xsl:attribute name="onclick">
+            <xsl:value-of select="concat($toggleRowsJS, '; return false;')"/>
+          </xsl:attribute>
+          Hide <xsl:value-of select="@MessageCount" /> additional messages
+        </a>
+      </td>
+    </tr>
+  </xsl:template>
+
+  <!-- Message row templates -->
+  <xsl:template match="Message">
+    <tr>
+      <xsl:attribute name="name">
+        <xsl:value-of select="concat(@Status, 'RowClass', ../../@ProjectSafeName)"/>
+      </xsl:attribute>
+
+      <xsl:if test="@Level = 0">
+        <xsl:attribute name="style">display: none</xsl:attribute>
+      </xsl:if>
+      <td>
+        <a>
+          <xsl:attribute name="name">
+            <xsl:value-of select="concat(../../@ProjectSafeName, @Status)"/>
+          </xsl:attribute>
+        </a>
+        <img width="16" height="16">
+          <xsl:attribute name="src">
+            <xsl:choose>
+              <xsl:when test="@Status = 'Error'">_UpgradeReport_Files\UpgradeReport_Error.png</xsl:when>
+              <xsl:when test="@Status = 'Warning'">_UpgradeReport_Files\UpgradeReport_Warning.png</xsl:when>
+              <xsl:when test="@Status = 'Message'">_UpgradeReport_Files\UpgradeReport_Information.png</xsl:when>
+            </xsl:choose>
+          </xsl:attribute>
+          <xsl:attribute name="alt">
+            <xsl:value-of select="@Status" />
+          </xsl:attribute>
+        </img>
+      </td>
+      <td class="messageCell">
+        <strong>
+          <xsl:value-of select="@Source"/>:
+        </strong>
+        <span>
+          <xsl:value-of select="@Message"/>
+        </span>
+      </td>
+    </tr>
+  </xsl:template>
+
+  <!-- Project Details Template -->
+  <xsl:template match="Projects" mode="ProjectDetails">
+
+    <xsl:for-each select="Project">
+      <xsl:sort select="@ErrorCount" order="descending" />
+      <xsl:sort select="@WarningCount" order="descending" />
+      <!-- Always make solution last within a group -->
+      <xsl:sort select="@IsSolution" order="ascending" />
+      <xsl:sort select="@ProjectSafeName" order="ascending" />
+
+      <a>
+        <xsl:attribute name="name">
+          <xsl:value-of select="@ProjectSafeName"/>
+        </xsl:attribute>
+      </a>
+      <xsl:choose>
+        <xsl:when test="@ProjectDisplayName = ''">
+          <h3 _locID="ProjectDisplayNameHeader">Solution</h3>
+        </xsl:when>
+        <xsl:otherwise>
+          <h3>
+            <xsl:value-of select="@ProjectDisplayName"/>
+          </h3>
+        </xsl:otherwise>
+      </xsl:choose>
+
+      <table>
+        <tr>
+          <xsl:attribute name="id">
+            <xsl:value-of select="concat(@ProjectSafeName, 'HeaderRow')"/>
+          </xsl:attribute>
+          <th></th>
+          <th class="messageCell" _locID="MessageTableHeader">Message</th>
+        </tr>
+
+        <!-- Errors and warnings -->
+        <xsl:for-each select="Messages/Message[@Level &gt; 0]">
+          <xsl:sort select="@Level" order="descending" />
+          <xsl:apply-templates select="." />
+        </xsl:for-each>
+
+        <xsl:if test="@MessageCount > 0">
+          <xsl:apply-templates select="." mode="ProjectShowMessages" />
+        </xsl:if>
+
+        <!-- Messages -->
+        <xsl:for-each select="Messages/Message[@Level = 0]">
+          <xsl:apply-templates select="." />
+        </xsl:for-each>
+
+        <xsl:choose>
+          <!-- Additional row as a 'place holder' for 'Show/Hide' additional messages -->
+          <xsl:when test="@MessageCount > 0">
+            <xsl:apply-templates select="." mode="ProjectHideMessages" />
+          </xsl:when>
+          <!-- No messages at all, show blank row -->
+          <xsl:when test="@TotalCount = 0">
+            <tr>
+              <td>
+                <img width="16" height="16" src="_UpgradeReport_Files\UpgradeReport_Information.png" />
+              </td>
+              <xsl:choose>
+                <xsl:when test="@ProjectDisplayName = ''">
+                  <td class="messageCell" _locID="NoMessagesRow2">
+                    Solution logged no messages.
+                  </td>
+                </xsl:when>
+                <xsl:otherwise>
+                  <td class="messageCell" _locID="NoMessagesRow">
+                    <xsl:value-of select="@ProjectDisplayName" /> logged no messages.
+                  </td>
+                </xsl:otherwise>
+              </xsl:choose>
+            </tr>
+          </xsl:when>
+        </xsl:choose>
+      </table>
+    </xsl:for-each>
+  </xsl:template>
+
+  <!-- Document, matches "UpgradeReport" -->
+  <xsl:template match="UpgradeReport">
+    <!-- Output doc type the 'Mark of the web' which disabled prompting to run JavaScript from local HTML Files in IE -->
+    <!-- NOTE: The whitespace around the 'Mark of the web' is important it must be exact -->
+    <xsl:text disable-output-escaping="yes"><![CDATA[<!DOCTYPE html>
+<!-- saved from url=(0014)about:internet -->
+]]>
+    </xsl:text>
+    <html>
+      <head>
+        <meta content="en-us" http-equiv="Content-Language" />
+        <meta content="text/html; charset=utf-16" http-equiv="Content-Type" />
+        <link type="text/css" rel="stylesheet" href="_UpgradeReport_Files\UpgradeReport.css" />
+        <title _locID="ConversionReport0">
+          Migration Report
+        </title>
+
+        <script type="text/javascript" language="javascript">
+          <xsl:text disable-output-escaping="yes">
+          <![CDATA[
+          
+            // Startup 
+            // Hook up the the loaded event for the document/window, to linkify the document content
+            var startupFunction = function() { linkifyElement("messages"); };
+            
+            if(window.attachEvent)
+            {
+              window.attachEvent('onload', startupFunction);
+            }
+            else if (window.addEventListener) 
+            {
+              window.addEventListener('load', startupFunction, false);
+            }
+            else 
+            {
+              document.addEventListener('load', startupFunction, false);
+            } 
+            
+            // Toggles the visibility of table rows with the specified name 
+            function toggleTableRowsByName(name)
+            {
+               var allRows = document.getElementsByTagName('tr');
+               for (i=0; i < allRows.length; i++)
+               {
+                  var currentName = allRows[i].getAttribute('name');
+                  if(!!currentName && currentName.indexOf(name) == 0)
+                  {
+                      var isVisible = allRows[i].style.display == ''; 
+                      isVisible ? allRows[i].style.display = 'none' : allRows[i].style.display = '';
+                  }
+               }
+            }
+            
+            function scrollToFirstVisibleRow(name) 
+            {
+               var allRows = document.getElementsByTagName('tr');
+               for (i=0; i < allRows.length; i++)
+               {
+                  var currentName = allRows[i].getAttribute('name');
+                  var isVisible = allRows[i].style.display == ''; 
+                  if(!!currentName && currentName.indexOf(name) == 0 && isVisible)
+                  {
+                     allRows[i].scrollIntoView(true); 
+                     return true; 
+                  }
+               }
+               
+               return false;
+            }
+            
+            // Linkifies the specified text content, replaces candidate links with html links 
+            function linkify(text)
+            {
+                 if(!text || 0 === text.length)
+                 {
+                     return text; 
+                 }
+
+                 // Find http, https and ftp links and replace them with hyper links 
+                 var urlLink = /(http|https|ftp)\:\/\/[a-zA-Z0-9\-\.]+(:[a-zA-Z0-9]*)?\/?([a-zA-Z0-9\-\._\?\,\/\\\+&%\$#\=~;\{\}])*/gi;
+                 
+                 return text.replace(urlLink, '<a href="$&">$&</a>') ;
+            }
+            
+            // Linkifies the specified element by ID
+            function linkifyElement(id)
+            {
+                var element = document.getElementById(id);
+                if(!!element)
+                {
+                  element.innerHTML = linkify(element.innerHTML); 
+                }
+            }
+            
+            function ToggleMessageVisibility(projectName)
+            {
+              if(!projectName || 0 === projectName.length)
+              {
+                return; 
+              }
+              
+              toggleTableRowsByName("MessageRowClass" + projectName);
+              toggleTableRowsByName('MessageRowHeaderShow' + projectName);
+              toggleTableRowsByName('MessageRowHeaderHide' + projectName); 
+            }
+            
+            function ScrollToFirstVisibleMessage(projectName)
+            {
+              if(!projectName || 0 === projectName.length)
+              {
+                return; 
+              }
+              
+              // First try the 'Show messages' row
+              if(!scrollToFirstVisibleRow('MessageRowHeaderShow' + projectName))
+              {
+                // Failed to find a visible row for 'Show messages', try an actual message row 
+                scrollToFirstVisibleRow('MessageRowClass' + projectName); 
+              }
+            }
+          ]]>
+        </xsl:text>
+        </script>
+      </head>
+      <body>
+        <h1 _locID="ConversionReport">
+          Migration Report - <xsl:value-of select="Properties/Property[@Name='Solution']/@Value"/>
+        </h1>
+
+        <div id="content">
+          <h2 _locID="OverviewTitle">Overview</h2>
+          <xsl:variable name="projectOverview">
+            <xsl:apply-templates select="self::node()" mode="ProjectOverviewXML" />
+          </xsl:variable>
+
+          <div id="overview">
+            <xsl:apply-templates select="msxsl:node-set($projectOverview)/*" mode="ProjectOverview" />
+          </div>
+
+          <h2 _locID="SolutionAndProjectsTitle">Solution and projects</h2>
+
+          <div id="messages">
+            <xsl:apply-templates select="msxsl:node-set($projectOverview)/*" mode="ProjectDetails" />
+          </div>
+        </div>
+      </body>
+    </html>
+  </xsl:template>
+
+</xsl:stylesheet>

src/Numerics/LinearAlgebra/Complex/KnuthSparseMatrix.cs

@@ -0,0 +1,355 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+    using System;
+    using System.Collections.Generic;
+    using System.Linq;
+    using System.Text;
+    using MathNet.Numerics.LinearAlgebra.Storage;
+    using Generic;
+    using Properties;
+    using Threading;
+    using System.Numerics;
+    using Storage;
+
+    /// <summary>
+    /// A Matrix class with sparse storage. The underlying storage scheme is pointers to lefter and upper elements Format.
+    /// Algorithm was posted in "The Art of Computer Programming" - Knuth - 1968
+    /// </summary>
+    public class KnuthSparseMatrix : Matrix
+    {
+        /// <summary>
+        /// Storage element
+        /// </summary>
+        private KnuthSparseMatrixStorage<Complex> storage;
+
+        /// <summary>
+        /// Creates a <c>KnuthSparseMatrix</c> for the given number of rows and columns.
+        /// </summary>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
+        /// <returns>
+        /// A <c>SparseMatrix</c> with the given dimensions.
+        /// </returns>
+        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
+        {
+            return new KnuthSparseMatrix(numberOfRows, numberOfColumns);
+        }
+
+        /// <summary>
+        /// Creates a <see cref="SparseVector"/> with a the given dimension.
+        /// </summary>
+        /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
+        /// <returns>
+        /// A <see cref="SparseVector"/> with the given dimension.
+        /// </returns>
+        public override Vector<Complex> CreateVector(int size, bool fullyMutable = false)
+        {
+            return new SparseVector(size);
+        }
+
+        /// <summary>
+        /// Gets the number of non zero elements in the matrix.
+        /// </summary>
+        /// <value>The number of non zero elements.</value>
+        public int NonZerosCount()
+        {
+            int result = 0;
+            for (int i = 0; i < storage.Rows.Count; ++i)
+            {
+                KnuthNode<Complex> iter = storage.Rows[i];
+                while (iter.Left != storage.Rows[i])
+                {
+                    ++result;
+                    iter = iter.Left;
+                }
+            }
+            return result;
+        }
+
+        /// <summary>
+        /// Insert new value to current possition
+        /// </summary>
+        /// <param name="row">number of row</param>
+        /// <param name="col">number of column</param>
+        /// <param name="val">value to insert</param>
+        public void At(int row, int col, Complex val)
+        {
+            storage.At(row, col, val);
+        }
+
+        /// <summary>
+        /// Gets a value at current possition
+        /// </summary>
+        /// <param name="row">number of row</param>
+        /// <param name="col">number of column</param>
+        public Complex At(int row, int col)
+        {
+            return storage.At(row, col);
+        }
+
+        /// <summary>
+        /// Check if elements are equals
+        /// </summary>
+        /// <param name="obj">opject to compare</param>
+        /// <returns>if they are equals</returns>
+        public override bool Equals(object obj)
+        {
+            return base.Equals(obj);
+        }
+
+        /// <summary>
+        /// Gets hash code of element
+        /// </summary>
+        /// <returns>hash code</returns>
+        public override int GetHashCode()
+        {
+            return base.GetHashCode();
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="_storage">
+        /// Saved at KnuthSparseMatrixStorage matrix.
+        /// </param>
+        internal KnuthSparseMatrix(KnuthSparseMatrixStorage<Complex> _storage)
+            : base(_storage)
+        {
+            storage = _storage;
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="_storage">
+        /// Saved at KnuthSparseMatrixStorage matrix.
+        /// </param>
+        internal KnuthSparseMatrix(KnuthSparseMatrix mtx)
+            : base(mtx.storage)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="rows">
+        /// The number of rows.
+        /// </param>
+        /// <param name="columns">
+        /// The number of columns.
+        /// </param>
+        public KnuthSparseMatrix(int rows, int columns)
+            : this(new KnuthSparseMatrixStorage<Complex>(rows, columns, 0))
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class. This matrix is square with a given size.
+        /// </summary>
+        /// <param name="order">the size of the square matrix.</param>
+        /// <exception cref="ArgumentException">
+        /// If <paramref name="order"/> is less than one.
+        /// </exception>
+        public KnuthSparseMatrix(int order)
+            : this(order, order)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="SparseMatrix"/> class from a 2D array. 
+        /// </summary>
+        /// <param name="array">The 2D array to create this matrix from.</param>
+        public KnuthSparseMatrix(Complex[,] array)
+            : this(array.GetLength(0), array.GetLength(1))
+        {
+            for (var i = 0; i < storage.RowCount; i++)
+            {
+                for (var j = 0; j < storage.ColumnCount; j++)
+                {
+                    storage.At(i, j, array[i, j]);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Create new instanse of sparse matrix 
+        /// </summary>
+        /// <param name="mtx">matrix in list<list> presentation to set as sparse</param>
+        public KnuthSparseMatrix(List<List<Complex>> array)
+            : this(array.Count, array[0].Count)
+        {
+            for (var i = 0; i < storage.RowCount; i++)
+            {
+                for (var j = 0; j < storage.ColumnCount; j++)
+                {
+                    storage.At(i, j, array[i][j]);
+                }
+            }
+        }
+        /// <summary>
+        /// Change signs of all values to opposite
+        /// </summary>
+        public void Negate()
+        {
+            KnuthNode<Complex> rw;
+            for (int i = 0; i < storage.Rows.Count; ++i)
+            {
+                rw = storage.Rows[i];
+                while (rw.Left != storage.Rows[i])
+                {
+                    rw.Val *= -1;
+                    rw = rw.Left;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Multiply two sparse matrix
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>result of multiplying this matrix</returns>
+        public static KnuthSparseMatrix operator *(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            if (mtx1.storage.Cols.Count != mtx2.storage.Rows.Count)
+                throw new Exception("Wrong size of matrix to multiple");
+            KnuthSparseMatrix Res = new KnuthSparseMatrix(mtx1.storage.Rows.Count, mtx2.storage.Cols.Count);
+            Complex tmp;
+            KnuthNode<Complex> Iter1, Iter2;
+            for (int i = 0; i < mtx1.storage.Rows.Count; ++i)
+            {
+                for (int j = 0; j < mtx2.storage.Cols.Count; ++j)
+                {
+                    tmp = 0; //initalize mult of i-th row and j-th column
+                    Iter1 = mtx1.storage.Rows[i].Left;
+                    Iter2 = mtx2.storage.Cols[j].Up;
+                    while (Iter1 != mtx1.storage.Rows[i] && Iter2 != mtx2.storage.Cols[j])
+                    {
+                        if (Iter1.Col == Iter2.Row)
+                        {
+                            tmp += Iter1.Val * Iter2.Val;
+                            Iter1 = Iter1.Left;
+                            Iter2 = Iter2.Up;
+                        }
+                        else if (Iter1.Col < Iter2.Row)
+                        {
+                            Iter2 = Iter2.Up;
+                        }
+                        else //(Iter2.Col < Iter1.Col)
+                        {
+                            Iter1 = Iter1.Left;
+                        }
+                    }
+                    if (tmp != Complex.Zero) // paste not 0 element
+                    {
+                        Res.At(i, j, tmp);
+                    }
+
+                }
+            }
+            return Res;
+        }
+
+        /// <summary>
+        /// Adds two matrixes
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>result of summing this matrixes</returns>
+        public static KnuthSparseMatrix operator +(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            if (mtx1.storage.Rows.Count != mtx2.storage.Rows.Count || mtx1.storage.Cols.Count != mtx2.storage.Cols.Count)
+                throw new Exception("Wrong size of matrix to summ");
+            KnuthSparseMatrix Res = new KnuthSparseMatrix(mtx1.storage.Rows.Count, mtx1.storage.Cols.Count);
+            KnuthNode<Complex> Iter1, Iter2;
+            for (int i = 0; i < mtx1.storage.Rows.Count; ++i)
+            {
+                Iter1 = mtx1.storage.Rows[i].Left;
+                Iter2 = mtx2.storage.Rows[i].Left;
+                while (Iter1 != mtx1.storage.Rows[i] || Iter2 != mtx2.storage.Rows[i])
+                {
+                    if (Iter1 == mtx1.storage.Rows[i])
+                    {
+                        Res.At(Iter2.Row, Iter2.Col, Iter2.Val);
+                        Iter2 = Iter2.Left;
+                    }
+                    else if (Iter2 == mtx2.storage.Rows[i])
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val);
+                        Iter1 = Iter1.Left;
+                    }
+                    else if (Iter1.Col == Iter2.Col)
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val + Iter2.Val);
+                        Iter1 = Iter1.Left;
+                        Iter2 = Iter2.Left;
+                    }
+                    else if (Iter1.Col < Iter2.Col)
+                    {
+                        Res.At(Iter2.Row, Iter2.Col, Iter2.Val);
+                        Iter2 = Iter2.Left;
+                    }
+                    else //(Iter2.Col < Iter1.Col)
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val);
+                        Iter1 = Iter1.Left;
+                    }
+                }
+            }
+            return Res;
+        }
+
+        /// <summary>
+        /// Change signs of all values to opposite
+        /// </summary>
+        /// <param name="rightSide"></param>
+        /// <returns></returns>
+        public static KnuthSparseMatrix operator -(KnuthSparseMatrix rightSide)
+        {
+            if (rightSide == null)
+            {
+                throw new ArgumentNullException("rightSide");
+            }
+            KnuthSparseMatrix res = new KnuthSparseMatrix(rightSide);
+            res.Negate();
+            return res;
+        }
+
+        /// <summary>
+        /// Check if two matrixes is equals
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>true if matrixes are equal, else false</returns>
+        public static bool operator ==(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            return mtx1.Equals(mtx2);
+        }
+        /// <summary>
+        /// Check if two matrixes aren`t equals
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>true if matrixes are not equals, else - false</returns>
+        public static bool operator !=(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            return !mtx1.Equals(mtx2);
+        }
+
+        /// <summary>
+        /// Get current matrix in full standart form
+        /// </summary>
+        /// <returns>matrix in full form</returns>
+        public Complex[][] getFullMatrix()
+        {
+            return storage.getFullMatrix();
+        }
+    }
+}

src/Numerics/LinearAlgebra/Double/KnuthSparseMatrix.cs

@@ -0,0 +1,349 @@
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+{
+    using System;
+    using System.Collections.Generic;
+    using System.Linq;
+    using System.Text;
+    using MathNet.Numerics.LinearAlgebra.Storage;
+    using Generic;
+    using Properties;
+    using Threading;
+
+    /// <summary>
+    /// A Matrix class with sparse storage. The underlying storage scheme is pointers to lefter and upper elements Format.
+    /// Algorithm was posted in "The Art of Computer Programming" - Knuth - 1968
+    /// </summary>
+    public class KnuthSparseMatrix : Matrix
+    {
+        /// <summary>
+        /// Storage element
+        /// </summary>
+        private KnuthSparseMatrixStorage<double> storage;
+
+        /// <summary>
+        /// Creates a <c>KnuthSparseMatrix</c> for the given number of rows and columns.
+        /// </summary>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
+        /// <returns>
+        /// A <c>SparseMatrix</c> with the given dimensions.
+        /// </returns>
+        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
+        {
+            return new KnuthSparseMatrix(numberOfRows, numberOfColumns);
+        }
+
+        /// <summary>
+        /// Creates a <see cref="SparseVector"/> with a the given dimension.
+        /// </summary>
+        /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
+        /// <returns>
+        /// A <see cref="SparseVector"/> with the given dimension.
+        /// </returns>
+        public override Vector<double> CreateVector(int size, bool fullyMutable = false)
+        {
+            return new SparseVector(size);
+        }
+
+        /// <summary>
+        /// Gets the number of non zero elements in the matrix.
+        /// </summary>
+        /// <value>The number of non zero elements.</value>
+        public int NonZerosCount()
+        {
+            int result = 0;
+            for (int i = 0; i < storage.Rows.Count; ++i)
+            {
+                KnuthNode<double> iter = storage.Rows[i];
+                while (iter.Left != storage.Rows[i])
+                {
+                    ++result;
+                    iter = iter.Left;
+                }
+            }
+            return result;
+        }
+
+        /// <summary>
+        /// Insert new value to current possition
+        /// </summary>
+        /// <param name="row">number of row</param>
+        /// <param name="col">number of column</param>
+        /// <param name="val">value to insert</param>
+        public void At(int row, int col, double val)
+        {
+            storage.At(row, col, val);
+        }
+
+        /// <summary>
+        /// Gets a value at current possition
+        /// </summary>
+        /// <param name="row">number of row</param>
+        /// <param name="col">number of column</param>
+        public double At(int row, int col)
+        {
+            return storage.At(row, col);
+        }
+
+        /// <summary>
+        /// Check if elements are equals
+        /// </summary>
+        /// <param name="obj">opject to compare</param>
+        /// <returns>if they are equals</returns>
+        public override bool Equals(object obj)
+        {
+            return base.Equals(obj);
+        }
+
+        /// <summary>
+        /// Gets hash code of element
+        /// </summary>
+        /// <returns>hash code</returns>
+        public override int GetHashCode()
+        {
+            return base.GetHashCode();
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="_storage">
+        /// Saved at KnuthSparseMatrixStorage matrix.
+        /// </param>
+        internal KnuthSparseMatrix(KnuthSparseMatrixStorage<double> _storage)
+            : base(_storage)
+        {
+            storage = _storage;
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="_storage">
+        /// Saved at KnuthSparseMatrixStorage matrix.
+        /// </param>
+        internal KnuthSparseMatrix(KnuthSparseMatrix mtx)
+            : base(mtx.storage)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="rows">
+        /// The number of rows.
+        /// </param>
+        /// <param name="columns">
+        /// The number of columns.
+        /// </param>
+        public KnuthSparseMatrix(int rows, int columns)
+            : this(new KnuthSparseMatrixStorage<double>(rows, columns, 0d))
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class. This matrix is square with a given size.
+        /// </summary>
+        /// <param name="order">the size of the square matrix.</param>
+        /// <exception cref="ArgumentException">
+        /// If <paramref name="order"/> is less than one.
+        /// </exception>
+        public KnuthSparseMatrix(int order)
+            : this(order, order)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="SparseMatrix"/> class from a 2D array. 
+        /// </summary>
+        /// <param name="array">The 2D array to create this matrix from.</param>
+        public KnuthSparseMatrix(double[,] array)
+            : this(array.GetLength(0), array.GetLength(1))
+        {
+            for (var i = 0; i < storage.RowCount; i++)
+            {
+                for (var j = 0; j < storage.ColumnCount; j++)
+                {
+                    storage.At(i, j, array[i, j]);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Create new instanse of sparse matrix 
+        /// </summary>
+        /// <param name="mtx">matrix in list<list> presentation to set as sparse</param>
+        public KnuthSparseMatrix(List<List<double>> array)
+            : this(array.Count, array[0].Count)
+        {
+            for (var i = 0; i < storage.RowCount; i++)
+            {
+                for (var j = 0; j < storage.ColumnCount; j++)
+                {
+                    storage.At(i, j, array[i][j]);
+                }
+            }
+        }
+        /// <summary>
+        /// Change signs of all values to opposite
+        /// </summary>
+        public void Negate()
+        {
+            KnuthNode<double> rw;
+            for (int i = 0; i < storage.Rows.Count; ++i)
+            {
+                rw = storage.Rows[i];
+                while (rw.Left != storage.Rows[i])
+                {
+                    rw.Val *= -1;
+                    rw = rw.Left;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Multiply two sparse matrix
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>result of multiplying this matrix</returns>
+        public static KnuthSparseMatrix operator *(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            if (mtx1.storage.Cols.Count != mtx2.storage.Rows.Count)
+                throw new Exception("Wrong size of matrix to multiple");
+            KnuthSparseMatrix Res = new KnuthSparseMatrix(mtx1.storage.Rows.Count, mtx2.storage.Cols.Count);
+            double tmp;
+            KnuthNode<double> Iter1, Iter2;
+            for (int i = 0; i < mtx1.storage.Rows.Count; ++i)
+            {
+                for (int j = 0; j < mtx2.storage.Cols.Count; ++j)
+                {
+                    tmp = 0; //initalize mult of i-th row and j-th column
+                    Iter1 = mtx1.storage.Rows[i].Left;
+                    Iter2 = mtx2.storage.Cols[j].Up;
+                    while (Iter1 != mtx1.storage.Rows[i] && Iter2 != mtx2.storage.Cols[j])
+                    {
+                        if (Iter1.Col == Iter2.Row)
+                        {
+                            tmp += Iter1.Val * Iter2.Val;
+                            Iter1 = Iter1.Left;
+                            Iter2 = Iter2.Up;
+                        }
+                        else if (Iter1.Col < Iter2.Row)
+                        {
+                            Iter2 = Iter2.Up;
+                        }
+                        else //(Iter2.Col < Iter1.Col)
+                        {
+                            Iter1 = Iter1.Left;
+                        }
+                    }
+                    if (tmp != 0) // paste not 0 element
+                    {
+                        Res.At(i, j, tmp);
+                    }
+
+                }
+            }
+            return Res;
+        }
+
+        /// <summary>
+        /// Adds two matrixes
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>result of summing this matrixes</returns>
+        public static KnuthSparseMatrix operator +(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            if (mtx1.storage.Rows.Count != mtx2.storage.Rows.Count || mtx1.storage.Cols.Count != mtx2.storage.Cols.Count)
+                throw new Exception("Wrong size of matrix to summ");
+            KnuthSparseMatrix Res = new KnuthSparseMatrix(mtx1.storage.Rows.Count, mtx1.storage.Cols.Count);
+            KnuthNode<double> Iter1, Iter2;
+            for (int i = 0; i < mtx1.storage.Rows.Count; ++i)
+            {
+                Iter1 = mtx1.storage.Rows[i].Left;
+                Iter2 = mtx2.storage.Rows[i].Left;
+                while (Iter1 != mtx1.storage.Rows[i] || Iter2 != mtx2.storage.Rows[i])
+                {
+                    if (Iter1 == mtx1.storage.Rows[i])
+                    {
+                        Res.At(Iter2.Row, Iter2.Col, Iter2.Val);
+                        Iter2 = Iter2.Left;
+                    }
+                    else if (Iter2 == mtx2.storage.Rows[i])
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val);
+                        Iter1 = Iter1.Left;
+                    }
+                    else if (Iter1.Col == Iter2.Col)
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val + Iter2.Val);
+                        Iter1 = Iter1.Left;
+                        Iter2 = Iter2.Left;
+                    }
+                    else if (Iter1.Col < Iter2.Col)
+                    {
+                        Res.At(Iter2.Row, Iter2.Col, Iter2.Val);
+                        Iter2 = Iter2.Left;
+                    }
+                    else //(Iter2.Col < Iter1.Col)
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val);
+                        Iter1 = Iter1.Left;
+                    }
+                }
+            }
+            return Res;
+        }
+
+        /// <summary>
+        /// Change signs of all values to opposite
+        /// </summary>
+        /// <param name="rightSide"></param>
+        /// <returns></returns>
+        public static KnuthSparseMatrix operator -(KnuthSparseMatrix rightSide)
+        {
+            if (rightSide == null)
+            {
+                throw new ArgumentNullException("rightSide");
+            }
+            KnuthSparseMatrix res = new KnuthSparseMatrix(rightSide);
+            res.Negate();
+            return res;
+        }
+
+        /// <summary>
+        /// Check if two matrixes is equals
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>true if matrixes are equal, else false</returns>
+        public static bool operator ==(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            return mtx1.Equals(mtx2);
+        }
+        /// <summary>
+        /// Check if two matrixes aren`t equals
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>true if matrixes are not equals, else - false</returns>
+        public static bool operator !=(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            return !mtx1.Equals(mtx2);
+        }
+
+        /// <summary>
+        /// Get current matrix in full standart form
+        /// </summary>
+        /// <returns>matrix in full form</returns>
+        public double[][] getFullMatrix()
+        {
+            return storage.getFullMatrix();
+        }
+    }
+}

src/Numerics/LinearAlgebra/Single/KnuthSparseMatrix.cs

@@ -0,0 +1,353 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace MathNet.Numerics.LinearAlgebra.Single
+{
+    using System;
+    using System.Collections.Generic;
+    using System.Linq;
+    using System.Text;
+    using MathNet.Numerics.LinearAlgebra.Storage;
+    using Generic;
+    using Properties;
+    using Threading;
+
+    /// <summary>
+    /// A Matrix class with sparse storage. The underlying storage scheme is pointers to lefter and upper elements Format.
+    /// Algorithm was posted in "The Art of Computer Programming" - Knuth - 1968
+    /// </summary>
+    public class KnuthSparseMatrix : Matrix
+    {
+        /// <summary>
+        /// Storage element
+        /// </summary>
+        private KnuthSparseMatrixStorage<float> storage;
+
+        /// <summary>
+        /// Creates a <c>KnuthSparseMatrix</c> for the given number of rows and columns.
+        /// </summary>
+        /// <param name="numberOfRows">The number of rows.</param>
+        /// <param name="numberOfColumns">The number of columns.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable (e.g. not a diagonal matrix).</param>
+        /// <returns>
+        /// A <c>SparseMatrix</c> with the given dimensions.
+        /// </returns>
+        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
+        {
+            return new KnuthSparseMatrix(numberOfRows, numberOfColumns);
+        }
+
+        /// <summary>
+        /// Creates a <see cref="SparseVector"/> with a the given dimension.
+        /// </summary>
+        /// <param name="size">The size of the vector.</param>
+        /// <param name="fullyMutable">True if all fields must be mutable.</param>
+        /// <returns>
+        /// A <see cref="SparseVector"/> with the given dimension.
+        /// </returns>
+        public override Vector<float> CreateVector(int size, bool fullyMutable = false)
+        {
+            return new SparseVector(size);
+        }
+
+        /// <summary>
+        /// Gets the number of non zero elements in the matrix.
+        /// </summary>
+        /// <value>The number of non zero elements.</value>
+        public int NonZerosCount()
+        {
+            int result = 0;
+            for (int i = 0; i < storage.Rows.Count; ++i)
+            {
+                KnuthNode<float> iter = storage.Rows[i];
+                while (iter.Left != storage.Rows[i])
+                {
+                    ++result;
+                    iter = iter.Left;
+                }
+            }
+            return result;
+        }
+
+        /// <summary>
+        /// Insert new value to current possition
+        /// </summary>
+        /// <param name="row">number of row</param>
+        /// <param name="col">number of column</param>
+        /// <param name="val">value to insert</param>
+        public void At(int row, int col, float val)
+        {
+            storage.At(row, col, val);
+        }
+
+        /// <summary>
+        /// Gets a value at current possition
+        /// </summary>
+        /// <param name="row">number of row</param>
+        /// <param name="col">number of column</param>
+        public float At(int row, int col)
+        {
+            return storage.At(row, col);
+        }
+
+        /// <summary>
+        /// Check if elements are equals
+        /// </summary>
+        /// <param name="obj">opject to compare</param>
+        /// <returns>if they are equals</returns>
+        public override bool Equals(object obj)
+        {
+            return base.Equals(obj);
+        }
+
+        /// <summary>
+        /// Gets hash code of element
+        /// </summary>
+        /// <returns>hash code</returns>
+        public override int GetHashCode()
+        {
+            return base.GetHashCode();
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="_storage">
+        /// Saved at KnuthSparseMatrixStorage matrix.
+        /// </param>
+        internal KnuthSparseMatrix(KnuthSparseMatrixStorage<float> _storage)
+            : base(_storage)
+        {
+            storage = _storage;
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="_storage">
+        /// Saved at KnuthSparseMatrixStorage matrix.
+        /// </param>
+        internal KnuthSparseMatrix(KnuthSparseMatrix mtx)
+            : base(mtx.storage)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class.
+        /// </summary>
+        /// <param name="rows">
+        /// The number of rows.
+        /// </param>
+        /// <param name="columns">
+        /// The number of columns.
+        /// </param>
+        public KnuthSparseMatrix(int rows, int columns)
+            : this(new KnuthSparseMatrixStorage<float>(rows, columns, 0))
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="KnuthSparseMatrix"/> class. This matrix is square with a given size.
+        /// </summary>
+        /// <param name="order">the size of the square matrix.</param>
+        /// <exception cref="ArgumentException">
+        /// If <paramref name="order"/> is less than one.
+        /// </exception>
+        public KnuthSparseMatrix(int order)
+            : this(order, order)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="SparseMatrix"/> class from a 2D array. 
+        /// </summary>
+        /// <param name="array">The 2D array to create this matrix from.</param>
+        public KnuthSparseMatrix(float[,] array)
+            : this(array.GetLength(0), array.GetLength(1))
+        {
+            for (var i = 0; i < storage.RowCount; i++)
+            {
+                for (var j = 0; j < storage.ColumnCount; j++)
+                {
+                    storage.At(i, j, array[i, j]);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Create new instanse of sparse matrix 
+        /// </summary>
+        /// <param name="mtx">matrix in list<list> presentation to set as sparse</param>
+        public KnuthSparseMatrix(List<List<float>> array)
+            : this(array.Count, array[0].Count)
+        {
+            for (var i = 0; i < storage.RowCount; i++)
+            {
+                for (var j = 0; j < storage.ColumnCount; j++)
+                {
+                    storage.At(i, j, array[i][j]);
+                }
+            }
+        }
+        /// <summary>
+        /// Change signs of all values to opposite
+        /// </summary>
+        public void Negate()
+        {
+            KnuthNode<float> rw;
+            for (int i = 0; i < storage.Rows.Count; ++i)
+            {
+                rw = storage.Rows[i];
+                while (rw.Left != storage.Rows[i])
+                {
+                    rw.Val *= -1;
+                    rw = rw.Left;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Multiply two sparse matrix
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>result of multiplying this matrix</returns>
+        public static KnuthSparseMatrix operator *(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            if (mtx1.storage.Cols.Count != mtx2.storage.Rows.Count)
+                throw new Exception("Wrong size of matrix to multiple");
+            KnuthSparseMatrix Res = new KnuthSparseMatrix(mtx1.storage.Rows.Count, mtx2.storage.Cols.Count);
+            float tmp;
+            KnuthNode<float> Iter1, Iter2;
+            for (int i = 0; i < mtx1.storage.Rows.Count; ++i)
+            {
+                for (int j = 0; j < mtx2.storage.Cols.Count; ++j)
+                {
+                    tmp = 0; //initalize mult of i-th row and j-th column
+                    Iter1 = mtx1.storage.Rows[i].Left;
+                    Iter2 = mtx2.storage.Cols[j].Up;
+                    while (Iter1 != mtx1.storage.Rows[i] && Iter2 != mtx2.storage.Cols[j])
+                    {
+                        if (Iter1.Col == Iter2.Row)
+                        {
+                            tmp += Iter1.Val * Iter2.Val;
+                            Iter1 = Iter1.Left;
+                            Iter2 = Iter2.Up;
+                        }
+                        else if (Iter1.Col < Iter2.Row)
+                        {
+                            Iter2 = Iter2.Up;
+                        }
+                        else //(Iter2.Col < Iter1.Col)
+                        {
+                            Iter1 = Iter1.Left;
+                        }
+                    }
+                    if (tmp != 0) // paste not 0 element
+                    {
+                        Res.At(i, j, tmp);
+                    }
+
+                }
+            }
+            return Res;
+        }
+
+        /// <summary>
+        /// Adds two matrixes
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>result of summing this matrixes</returns>
+        public static KnuthSparseMatrix operator +(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            if (mtx1.storage.Rows.Count != mtx2.storage.Rows.Count || mtx1.storage.Cols.Count != mtx2.storage.Cols.Count)
+                throw new Exception("Wrong size of matrix to summ");
+            KnuthSparseMatrix Res = new KnuthSparseMatrix(mtx1.storage.Rows.Count, mtx1.storage.Cols.Count);
+            KnuthNode<float> Iter1, Iter2;
+            for (int i = 0; i < mtx1.storage.Rows.Count; ++i)
+            {
+                Iter1 = mtx1.storage.Rows[i].Left;
+                Iter2 = mtx2.storage.Rows[i].Left;
+                while (Iter1 != mtx1.storage.Rows[i] || Iter2 != mtx2.storage.Rows[i])
+                {
+                    if (Iter1 == mtx1.storage.Rows[i])
+                    {
+                        Res.At(Iter2.Row, Iter2.Col, Iter2.Val);
+                        Iter2 = Iter2.Left;
+                    }
+                    else if (Iter2 == mtx2.storage.Rows[i])
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val);
+                        Iter1 = Iter1.Left;
+                    }
+                    else if (Iter1.Col == Iter2.Col)
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val + Iter2.Val);
+                        Iter1 = Iter1.Left;
+                        Iter2 = Iter2.Left;
+                    }
+                    else if (Iter1.Col < Iter2.Col)
+                    {
+                        Res.At(Iter2.Row, Iter2.Col, Iter2.Val);
+                        Iter2 = Iter2.Left;
+                    }
+                    else //(Iter2.Col < Iter1.Col)
+                    {
+                        Res.At(Iter1.Row, Iter1.Col, Iter1.Val);
+                        Iter1 = Iter1.Left;
+                    }
+                }
+            }
+            return Res;
+        }
+
+        /// <summary>
+        /// Change signs of all values to opposite
+        /// </summary>
+        /// <param name="rightSide"></param>
+        /// <returns></returns>
+        public static KnuthSparseMatrix operator -(KnuthSparseMatrix rightSide)
+        {
+            if (rightSide == null)
+            {
+                throw new ArgumentNullException("rightSide");
+            }
+            KnuthSparseMatrix res = new KnuthSparseMatrix(rightSide);
+            res.Negate();
+            return res;
+        }
+
+        /// <summary>
+        /// Check if two matrixes is equals
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>true if matrixes are equal, else false</returns>
+        public static bool operator ==(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            return mtx1.Equals(mtx2);
+        }
+        /// <summary>
+        /// Check if two matrixes aren`t equals
+        /// </summary>
+        /// <param name="mtx1">first matrix</param>
+        /// <param name="mtx2">second matrix</param>
+        /// <returns>true if matrixes are not equals, else - false</returns>
+        public static bool operator !=(KnuthSparseMatrix mtx1, KnuthSparseMatrix mtx2)
+        {
+            return !mtx1.Equals(mtx2);
+        }
+
+        /// <summary>
+        /// Get current matrix in full standart form
+        /// </summary>
+        /// <returns>matrix in full form</returns>
+        public float[][] getFullMatrix()
+        {
+            return storage.getFullMatrix();
+        }
+    }
+}

src/Numerics/LinearAlgebra/Storage/KnuthNode.cs

@@ -0,0 +1,176 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace MathNet.Numerics.LinearAlgebra.Storage
+{
+    /// <summary>
+    /// Node to store elements in Knuth sparse matrix
+    /// </summary>
+     public partial class KnuthNode<T>
+         where T : struct, IEquatable<T>, IFormattable
+    {
+        /// <summary>
+        /// Define zero element for current type
+        /// </summary>
+        readonly T zero;
+
+        /// <summary>
+        /// Row position of element
+        /// </summary>
+        public int Row
+        {
+            set;
+            get;
+        }
+
+        /// <summary>
+        /// Column possition of element
+        /// </summary>
+        public int Col
+        {
+            set;
+            get;
+        }
+
+        /// <summary>
+        /// Stored value
+        /// </summary>
+        public T Val
+        {
+            set;
+            get;
+        }
+
+        /// <summary>
+        /// Pointer to lefter element in row 
+        /// </summary>
+        public KnuthNode<T> Left
+        {
+            set;
+            get;
+        }
+
+        /// <summary>
+        /// Pointer to upper element in column
+        /// </summary>
+        public KnuthNode<T> Up
+        {
+            set;
+            get;
+        }
+
+        /// <summary>
+        /// Initialize new default instanse of KnuthNode class
+        /// </summary>
+        public KnuthNode(T _zero)
+        {
+            zero = _zero;
+            Row = -1;
+            Col = -1;
+            Val = _zero;
+            Left = this;
+            Up = this;
+        }
+        /// <summary>
+        /// Initialize new instanse of KnuthNode class
+        /// </summary>
+        /// <param name="x">row number </param>
+        /// <param name="y">column number</param>
+        /// <param name="val">value</param>
+        public KnuthNode(int x, int y, T val, T _zero)
+        {
+            zero = _zero;
+            Row = x;
+            Col = y;
+            Val = val;
+            Left = this;
+            Up = this;
+        }
+
+        /// <summary>
+        /// Initialize new instanse of KnuthNode class
+        /// </summary>
+        /// <param name="x">row number </param>
+        /// <param name="y">column number</param>
+        /// <param name="val">value</param>
+        /// <param name="left">pointer to lefter element in row</param>
+        /// <param name="up">pointer to upper element in column</param>
+        public KnuthNode(int x, int y, T val, KnuthNode<T> left, KnuthNode<T> up, T _zero)
+        {
+            zero = _zero;
+            Row = x;
+            Col = y;
+            Val = val;
+            Left = left;
+            Up = up;
+            Left = this;
+            Up = this;
+        }
+
+        /// <summary>
+        /// Check if elements are equal
+        /// </summary>
+        /// <param name="nd2">second element</param>
+        /// <returns>if elements are equal</returns>
+        public bool Equals(KnuthNode<T> nd2)
+        {
+            if (nd2 == null)
+                return false;
+            return Val.Equals(nd2.Val) && Row == nd2.Row && Col == nd2.Col;
+        }
+
+        /// <summary>
+        /// Oparetor of equality
+        /// </summary>
+        /// <param name="nd1">first element</param>
+        /// <param name="nd2">second element</param>
+        /// <returns>if elements are equals</returns>
+        public static bool operator ==(KnuthNode<T> nd1, KnuthNode<T> nd2)
+        {
+            return nd1.Equals(nd2);
+        }
+
+        /// <summary>
+        /// Oparetor of equality
+        /// </summary>
+        /// <param name="nd1">first element</param>
+        /// <param name="nd2">second element</param>
+        /// <returns>if elements aren`t equals</returns>
+        public static bool operator !=(KnuthNode<T> nd1, KnuthNode<T> nd2)
+        {
+            return !nd1.Equals(nd2);
+        }
+
+        /// <summary>
+        /// Check if elements are equal
+        /// </summary>
+        /// <param name="nd2">second element</param>
+        /// <returns>if elements are equal</returns>
+        public override bool Equals(object obj)
+        {
+            return base.Equals(obj);
+        }
+
+         /// <summary>
+         /// Gets hash code of node
+         /// </summary>
+         /// <returns></returns>
+        public override int GetHashCode()
+        {
+            var hashNum = Math.Min(Row * Col, 25);
+            int hash = 17;
+            unchecked
+            {
+                for (var i = 0; i < hashNum; i++)
+                {
+                    var col = i % Row;
+                    var row = (i - col) / Col;
+                    hash = hash * 31 + Val.GetHashCode();
+                }
+            }
+            return hash;
+        }
+    }
+}

src/Numerics/LinearAlgebra/Storage/KnuthSparseMatrixStorage.cs

@@ -0,0 +1,453 @@
+
+namespace MathNet.Numerics.LinearAlgebra.Storage
+{
+    using System;
+    using System.Collections.Generic;
+
+    /// <summary>
+    /// class to storage matrix using algorithm of Knuth
+    /// </summary>
+    public class KnuthSparseMatrixStorage<T> : MatrixStorage<T>
+        where T : struct, IEquatable<T>, IFormattable
+    {
+        /// <summary>
+        /// Define zero element for cuurnt type
+        /// </summary>
+        readonly T zero;
+
+        /// <summary>
+        /// list to begin rows elements
+        /// </summary>
+        public List<KnuthNode<T>> Rows;
+
+        /// <summary>
+        /// list to begin cols elements
+        /// </summary>
+        public List<KnuthNode<T>> Cols;
+
+        /// <summary>
+        /// Compare to elements if they equals
+        /// </summary>
+        /// <param name="mtx2">element to compare</param>
+        /// <returns>if there quals</returns>
+        public bool Equals(KnuthSparseMatrixStorage<T> mtx2)
+        {
+            if (this.Rows.Count != mtx2.Rows.Count || this.Cols.Count != mtx2.Cols.Count)
+                return false;
+            KnuthNode<T> Iter1, Iter2;
+            //check rows 
+            for (int i = 0; i < this.Rows.Count; ++i)
+            {
+                Iter1 = this.Rows[i];
+                Iter2 = mtx2.Rows[i];
+                //empty and not empty row
+                if (Iter1.Left == Iter1 && Iter2.Left != Iter2)
+                    return false;
+                if (Iter2.Left == Iter2 && Iter1.Left != Iter1)
+                    return false;
+                while (Iter1.Left != this.Rows[i])
+                {
+                    //row of second mtx allready ends
+                    if (Iter2.Left == mtx2.Rows[i])
+                        return false;
+                    //not equal elements
+                    if (Iter1 != Iter2)
+                        return false;
+                    Iter1 = Iter1.Left;
+                    Iter2 = Iter2.Left;
+                }
+                //row of second mtx not allready ends
+                if (Iter2.Left != mtx2.Rows[i])
+                    return false;
+            }
+            //check cols 
+            for (int i = 0; i < this.Cols.Count; ++i)
+            {
+                Iter1 = this.Cols[i];
+                Iter2 = mtx2.Cols[i];
+                //empty and not empty coll
+                if (Iter1.Up == Iter1 && Iter2.Up != Iter2)
+                    return false;
+                if (Iter2.Up == Iter2 && Iter1.Up != Iter1)
+                    return false;
+                while (Iter1.Up != this.Cols[i])
+                {
+                    //coll of second mtx allready ends
+                    if (Iter2.Up == mtx2.Cols[i])
+                        return false;
+                    //not equal elements
+                    if (Iter1 != Iter2)
+                        return false;
+                    Iter1 = Iter1.Up;
+                    Iter2 = Iter2.Up;
+                }
+                //row of second mtx not allready ends
+                if (Iter2.Up != mtx2.Cols[i])
+                    return false;
+            }
+            //all elements are equals
+            return true;
+        }
+
+        public override bool Equals(MatrixStorage<T> other)
+        {
+            return base.Equals(other);
+        }
+
+        /// <summary>
+        /// Gets element at current position
+        /// </summary>
+        /// <param name="row">number of row</param>
+        /// <param name="column">number of column</param>
+        /// <returns>value on this position</returns>
+        public override T At(int row, int column)
+        {
+            return Get(row, column);
+        }
+
+        /// <summary>
+        /// Sets element to current position
+        /// </summary>
+        /// <param name="row"></param>
+        /// <param name="row">number of row</param>
+        /// <param name="column">number of column</param>
+        public override void At(int row, int column, T val)
+        {
+            Insert(row, column, val);
+        }
+
+        /// <summary>
+        /// Insert new value to current possition
+        /// </summary>
+        /// <param name="x">number of row</param>
+        /// <param name="y">number of column</param>
+        /// <param name="val">value to insert</param>
+        private void Insert(int x, int y, T val)
+        {
+            if (val.Equals(zero)) //var == 0
+            {
+                Delete(x, y);
+                return;
+            }
+            KnuthNode<T> rw = Rows[x];
+            KnuthNode<T> cl = Cols[y];
+            KnuthNode<T> toInsert = new KnuthNode<T>(x, y, val, zero);
+            //search needed row
+            while (rw.Left != Rows[x] && rw.Left.Col >= y)
+            {
+                rw = rw.Left;
+            }
+            //insert to row
+            if (rw.Col == y) // element allready exist
+            {
+                rw.Val = val;
+            }
+            else
+            {
+                toInsert.Left = rw.Left;
+                rw.Left = toInsert;
+            }
+            //search needed column
+            while (cl.Up != Cols[y] && cl.Up.Row >= x)
+            {
+                cl = cl.Up;
+            }
+            //insert to column
+            if (cl.Row == x) // element allready exist
+            {
+                cl.Val = val;
+            }
+            else
+            {
+                toInsert.Up = cl.Up;
+                cl.Up = toInsert;
+            }
+        }
+
+        /// <summary>
+        /// Delete current element
+        /// </summary>
+        /// <param name="x">number of row</param>
+        /// <param name="y">number of column</param>
+        private void Delete(int x, int y)
+        {
+            KnuthNode<T> rw = Rows[x].Left, prewRw = Rows[x];
+            KnuthNode<T> cl = Cols[y].Up, prewCl = Cols[y];
+            bool haver = false, havec = false;
+            for (int i = 0; i < Rows.Count; ++i)
+            {
+                if (rw.Col == y)
+                {
+                    haver = true;
+                    break;
+                }
+                prewRw = rw;
+                rw = rw.Left;
+            }
+            for (int i = 0; i < Cols.Count; ++i)
+            {
+                if (cl.Row == x)
+                {
+                    havec = true;
+                    break;
+                }
+                prewCl = cl;
+                cl = cl.Up;
+            }
+            if (havec && haver)
+            {
+                prewCl.Up = cl.Up;
+                prewRw.Left = rw.Left;
+            }
+        }
+
+        /// <summary>
+        /// Set bounds of new sparse matrix
+        /// </summary>
+        /// <param name="row">count of rows</param>
+        /// <param name="col">count of cols</param>
+        private void SetBounds(int row, int col)
+        {
+            Rows = new List<KnuthNode<T>>(row);
+            for (int i = 0; i < row; ++i)
+            {
+                Rows.Add(new KnuthNode<T>(zero));
+            }
+            Cols = new List<KnuthNode<T>>(col);
+            for (int i = 0; i < col; ++i)
+            {
+                Cols.Add(new KnuthNode<T>(zero));
+            }
+        }
+
+        /// <summary>
+        /// Create new instanse of sparse matrix 
+        /// </summary>
+        /// <param name="row">number of rows</param>
+        /// <param name="col">number of columns</param>
+        internal KnuthSparseMatrixStorage(int rows, int columns, T _zero)
+            : base(rows, columns)
+        {
+            zero = _zero;
+            SetBounds(rows, columns);
+        }
+
+        /// <summary>
+        /// Create new instanse of sparse matrix 
+        /// </summary>
+        /// <param name="mtx">matrix in standart presentation to set as sparse</param>
+        internal KnuthSparseMatrixStorage(int rows, int columns, T[][] mtx, T _zero)
+            : base(rows, columns)
+        {
+            zero = _zero;
+            SetBounds(mtx.Length, mtx[0].Length);
+            for (int i = 0; i < mtx.Length; ++i)
+            {
+                for (int j = 0; j < mtx[i].Length; ++j)
+                {
+                    if (!mtx[i][j].Equals(zero))
+                    {
+                        Insert(i, j, mtx[i][j]);
+                    }
+                }
+            }
+        }
+
+        /// <summary>
+        /// Create new instanse of sparse matrix 
+        /// </summary>
+        /// <param name="mtx">matrix in list<list> presentation to set as sparse</param>
+         internal KnuthSparseMatrixStorage(int rows, int columns, List<List<T>> mtx, T _zero)
+            : base(rows, columns)
+        {
+             zero = _zero;
+            SetBounds(mtx.Count, mtx[0].Count);
+            for (int i = 0; i < mtx.Count; ++i)
+            {
+                for (int j = 0; j < mtx[i].Count; ++j)
+                {
+                    if (!mtx[i][j].Equals(zero))
+                    {
+                        Insert(i, j, mtx[i][j]);
+                    }
+                }
+            }
+        }
+
+        /// <summary>
+        /// Copy existing instanse of matrix to new matrix
+        /// </summary>
+        /// <param name="mtx">matrix to copy</param>
+        public KnuthSparseMatrixStorage<T> Copy()
+        {
+            KnuthSparseMatrixStorage<T> res = new KnuthSparseMatrixStorage<T>(this.Rows.Count, this.Cols.Count, zero);
+            KnuthNode<T> rw;
+            for (int i = 0; i < Rows.Count; ++i)
+            {
+                rw = Rows[i].Left;
+                while (rw != Rows[i])
+                {
+                    res.Insert(rw.Row, rw.Col, rw.Val);
+                    rw = rw.Left;
+                }
+            }
+            return res;
+        }
+
+        /// <summary>
+        /// Gets element at curent position
+        /// </summary>
+        /// <param name="row"></param>
+        /// <param name="column"></param>
+        /// <returns></returns>
+        public T Get(int row, int column)
+        {
+            KnuthNode<T> rw = Rows[row].Left, prewRw = Rows[row];
+            KnuthNode<T> cl = Cols[column].Up, prewCl = Cols[column];
+            bool haver = false, havec = false;
+            for (int i = 0; i < Rows.Count; ++i)
+            {
+                if (rw.Col == column)
+                {
+                    haver = true;
+                    break;
+                }
+                prewRw = rw;
+                rw = rw.Left;
+            }
+            for (int i = 0; i < Cols.Count; ++i)
+            {
+                if (cl.Row == row)
+                {
+                    havec = true;
+                    break;
+                }
+                prewCl = cl;
+                cl = cl.Up;
+            }
+            if (havec && haver)
+            {
+                return cl.Val;
+            }
+            else
+            {
+                return zero;
+            }
+        }
+
+        /// <summary>
+        /// Gets or sets the value at the given row and column.
+        /// </summary>
+        /// <param name="row">
+        /// The row of the element.
+        /// </param>
+        /// <param name="column">
+        /// The column of the element.
+        /// </param>
+        /// <value>The value to get or set.</value>
+        //public T this[int row, int column]
+        //{
+        //    get
+        //    {
+        //        if(row< 0 || row >= Rows.Count)
+        //            throw new ArgumentOutOfRangeException("row");
+        //        if (column < 0 || column >= Cols.Count)
+        //            throw new ArgumentOutOfRangeException("column");
+        //        return Get(row, column);
+        //    }
+
+        //    set
+        //    {
+        //        if (row < 0 || row >= Rows.Count)
+        //            throw new ArgumentOutOfRangeException("row");
+        //        if (column < 0 || column >= Cols.Count)
+        //            throw new ArgumentOutOfRangeException("column");
+        //        Insert(row, column, value);
+        //    }
+        //}
+
+        public override void CopyTo(MatrixStorage<T> target, bool skipClearing = false)
+        {
+            base.CopyTo(target, skipClearing);
+        }
+        ///// <summary>
+        ///// Check if two matrixes is equals
+        ///// </summary>
+        ///// <param name="mtx1">first matrix</param>
+        ///// <param name="mtx2">second matrix</param>
+        ///// <returns>true if matrixes are equal, else false</returns>
+        //public static bool operator ==(KnuthSparseMatrixStorage<T> mtx1, KnuthSparseMatrixStorage<T> mtx2)
+        //{
+        //    return mtx1.Equals(mtx2);
+        //}
+        ///// <summary>
+        ///// Check if two matrixes aren`t equals
+        ///// </summary>
+        ///// <param name="mtx1">first matrix</param>
+        ///// <param name="mtx2">second matrix</param>
+        ///// <returns>true if matrixes are not equals, else - false</returns>
+        //public static bool operator !=(KnuthSparseMatrixStorage<T> mtx1, KnuthSparseMatrixStorage<T> mtx2)
+        //{
+        //    return !mtx1.Equals(mtx2);
+        //}
+
+        /// <summary>
+        /// Get current matrix in full standart form
+        /// </summary>
+        /// <returns>matrix in full form</returns>
+        public T[][] getFullMatrix()
+        {
+            T[][] res = new T[Rows.Count][];
+            for (int i = 0; i < Rows.Count; ++i)
+            {
+                res[i] = new T[Cols.Count];
+            }
+            //initialize
+            KnuthNode<T> iterator;
+            for (int i = 0; i < Rows.Count; ++i)
+            {
+                iterator = Rows[i].Left;
+                while (iterator != Rows[i])
+                {
+                    res[iterator.Row][iterator.Col] = iterator.Val;
+                    iterator = iterator.Left;
+                }
+            }
+            return res;
+        }
+
+        ///// <summary>
+        ///// Compare to elements if they equals
+        ///// </summary>
+        ///// <param name="mtx2">element to compare</param>
+        ///// <returns>if there quals</returns>
+        //public override bool Equals(object obj)
+        //{
+        //    return base.Equals(obj);
+        //}
+
+        /// <summary>
+        /// Serves as a hash function for a particular type.
+        /// </summary>
+        /// <returns>
+        /// A hash code for the current <see cref="T:System.Object"/>.
+        /// </returns>
+        public override int GetHashCode()
+        {
+            int RowCount = Rows.Count;
+            int ColumnCount = Cols.Count;
+            var hashNum = Math.Min(RowCount * ColumnCount, 25);
+            int hash = 17;
+            unchecked
+            {
+                for (var i = 0; i < hashNum; i++)
+                {
+                    var col = i % ColumnCount;
+                    var row = (i - col) / RowCount;
+                    hash = hash * 31 + Get(row, col).GetHashCode();
+                }
+            }
+            return hash;
+        }
+    }
+}

src/Numerics/Numerics.csproj

@@ -125,7 +125,13 @@
     <Compile Include="Distributions\Multivariate\MatrixNormal.cs" />
     <Compile Include="Distributions\Multivariate\Wishart.cs" />
     <Compile Include="LinearAlgebra\Storage\SparseVectorStorage.cs" />
+    <Compile Include="LinearAlgebra\Complex32\SparseMatrix.cs" />
+    <Compile Include="LinearAlgebra\Complex\KnuthSparseMatrix.cs" />
+    <Compile Include="LinearAlgebra\Double\KnuthSparseMatrix.cs" />
+    <Compile Include="LinearAlgebra\Single\KnuthSparseMatrix.cs" />
     <Compile Include="LinearAlgebra\Storage\DenseVectorStorage.cs" />
+    <Compile Include="LinearAlgebra\Storage\KnuthNode.cs" />
+    <Compile Include="LinearAlgebra\Storage\KnuthSparseMatrixStorage.cs" />
     <Compile Include="LinearAlgebra\Storage\MatrixStorage.Validation.cs" />
     <Compile Include="LinearAlgebra\Complex32\ExtensionMethods.cs" />
     <Compile Include="LinearAlgebra\Complex32\Factorization\Cholesky.cs" />
@@ -142,7 +148,6 @@
     <Compile Include="LinearAlgebra\Complex32\Solvers\IIterator.cs" />
     <Compile Include="LinearAlgebra\Complex32\Solvers\Preconditioners\IPreConditioner.cs" />
     <Compile Include="LinearAlgebra\Complex32\Solvers\StopCriterium\IIterationStopCriterium.cs" />
-    <Compile Include="LinearAlgebra\Complex32\SparseMatrix.cs" />
     <Compile Include="LinearAlgebra\Complex32\Vector.cs" />
     <Compile Include="LinearAlgebra\Complex\DenseMatrix.cs" />
     <Compile Include="LinearAlgebra\Complex\DenseVector.cs">