LA: Refactoring towards the generic matrix/vector types, test simplification

src/Numerics/LinearAlgebra/Builder.cs

@@ -0,0 +1,268 @@
+// <copyright file="Builder.cs" 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-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+#if NOSYSNUMERICS
+    using Complex64 = Numerics.Complex;
+#else
+    using Complex64 = System.Numerics.Complex;
+#endif
+
+namespace MathNet.Numerics.LinearAlgebra.Double
+{
+    internal class GenericBuilder : IGenericBuilder<double>
+    {
+        public double Zero
+        {
+            get { return 0d; }
+        }
+
+        public double One
+        {
+            get { return 1d; }
+        }
+
+        public Matrix<double> DenseMatrix(int rows, int columns)
+        {
+            return new DenseMatrix(rows, columns);
+        }
+
+        public Matrix<double> SparseMatrix(int rows, int columns)
+        {
+            return new SparseMatrix(rows, columns);
+        }
+
+        public Vector<double> DenseVector(int size)
+        {
+            return new DenseVector(size);
+        }
+
+        public Vector<double> SparseVector(int size)
+        {
+            return new SparseVector(size);
+        }
+    }
+}
+
+namespace MathNet.Numerics.LinearAlgebra.Single
+{
+    internal class GenericBuilder : IGenericBuilder<float>
+    {
+        public float Zero
+        {
+            get { return 0f; }
+        }
+
+        public float One
+        {
+            get { return 1f; }
+        }
+
+        public Matrix<float> DenseMatrix(int rows, int columns)
+        {
+            return new DenseMatrix(rows, columns);
+        }
+
+        public Matrix<float> SparseMatrix(int rows, int columns)
+        {
+            return new SparseMatrix(rows, columns);
+        }
+
+        public Vector<float> DenseVector(int size)
+        {
+            return new DenseVector(size);
+        }
+
+        public Vector<float> SparseVector(int size)
+        {
+            return new SparseVector(size);
+        }
+    }
+}
+
+namespace MathNet.Numerics.LinearAlgebra.Complex
+{
+    internal class GenericBuilder : IGenericBuilder<Complex64>
+    {
+        public Complex64 Zero
+        {
+            get { return Complex64.Zero; }
+        }
+
+        public Complex64 One
+        {
+            get { return Complex64.One; }
+        }
+
+        public Matrix<Complex64> DenseMatrix(int rows, int columns)
+        {
+            return new DenseMatrix(rows, columns);
+        }
+
+        public Matrix<Complex64> SparseMatrix(int rows, int columns)
+        {
+            return new SparseMatrix(rows, columns);
+        }
+
+        public Vector<Complex64> DenseVector(int size)
+        {
+            return new DenseVector(size);
+        }
+
+        public Vector<Complex64> SparseVector(int size)
+        {
+            return new SparseVector(size);
+        }
+    }
+}
+
+namespace MathNet.Numerics.LinearAlgebra.Complex32
+{
+    internal class GenericBuilder : IGenericBuilder<Numerics.Complex32>
+    {
+        public Numerics.Complex32 Zero
+        {
+            get { return Numerics.Complex32.Zero; }
+        }
+
+        public Numerics.Complex32 One
+        {
+            get { return Numerics.Complex32.One; }
+        }
+
+        public Matrix<Numerics.Complex32> DenseMatrix(int rows, int columns)
+        {
+            return new DenseMatrix(rows, columns);
+        }
+
+        public Matrix<Numerics.Complex32> SparseMatrix(int rows, int columns)
+        {
+            return new SparseMatrix(rows, columns);
+        }
+
+        public Vector<Numerics.Complex32> DenseVector(int size)
+        {
+            return new DenseVector(size);
+        }
+
+        public Vector<Numerics.Complex32> SparseVector(int size)
+        {
+            return new SparseVector(size);
+        }
+    }
+}
+
+namespace MathNet.Numerics.LinearAlgebra
+{
+    /// <summary>
+    /// Generic linear algebra type builder, for situations where a matrix or vector
+    /// must be created in a generic way. Usage of generic builders should not be
+    /// required in normal user code.
+    /// </summary>
+    public interface IGenericBuilder<T> where T : struct, IEquatable<T>, IFormattable
+    {
+        /// <summary>
+        /// Gets the value of <c>0.0</c> for type T.
+        /// </summary>
+        T Zero { get; }
+
+        /// <summary>
+        /// Gets the value of <c>1.0</c> for type T.
+        /// </summary>
+        T One { get; }
+
+        /// <summary>
+        /// Create a dense matrix of T with the given number of rows and columns.
+        /// </summary>
+        /// <param name="rows">The number of rows.</param>
+        /// <param name="columns">The number of columns.</param>
+        Matrix<T> DenseMatrix(int rows, int columns);
+
+        /// <summary>
+        /// Create a sparse matrix of T with the given number of rows and columns.
+        /// </summary>
+        /// <param name="rows">The number of rows.</param>
+        /// <param name="columns">The number of columns.</param>
+        Matrix<T> SparseMatrix(int rows, int columns);
+
+        /// <summary>
+        /// Create a dense vector of T with the given size.
+        /// </summary>
+        /// <param name="size">The size of the vector.</param>
+        Vector<T> DenseVector(int size);
+
+        /// <summary>
+        /// Create a sparse vector of T with the given size.
+        /// </summary>
+        /// <param name="size">The size of the vector.</param>
+        Vector<T> SparseVector(int size);
+    }
+
+    public static class Builder<T> where T : struct, IEquatable<T>, IFormattable
+    {
+        static Lazy<IGenericBuilder<T>> _singleton = new Lazy<IGenericBuilder<T>>(Create);
+
+        static IGenericBuilder<T> Create()
+        {
+            if (typeof(T) == typeof(Complex64))
+            {
+                return (IGenericBuilder<T>)(object)new Complex.GenericBuilder();
+            }
+
+            if (typeof(T) == typeof(Numerics.Complex32))
+            {
+                return (IGenericBuilder<T>)(object)new Complex32.GenericBuilder();
+            }
+
+            if (typeof(T) == typeof(double))
+            {
+                return (IGenericBuilder<T>)(object)new Double.GenericBuilder();
+            }
+
+            if (typeof(T) == typeof(float))
+            {
+                return (IGenericBuilder<T>)(object)new Single.GenericBuilder();
+            }
+
+            throw new NotSupportedException();
+        }
+
+        public static void Register(IGenericBuilder<T> builder)
+        {
+            _singleton = new Lazy<IGenericBuilder<T>>(() => builder);
+        }
+
+        public static IGenericBuilder<T> Instance
+        {
+            get { return _singleton.Value; }
+        }
+    }
+}

src/Numerics/LinearAlgebra/Common.cs

@@ -36,7 +36,6 @@ namespace MathNet.Numerics.LinearAlgebra
     using Complex64 = Numerics.Complex;
 #else
     using Complex64 = System.Numerics.Complex;
-
 #endif
 
     /// <summary>

src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@@ -357,33 +357,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             get { return _values; }
         }
 
-        /// <summary>
-        /// Creates a <c>DenseMatrix</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>DenseMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return new DenseMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<Complex> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Complex/DenseVector.cs

@@ -197,39 +197,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new DenseVector(array);
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex> CreateMatrix(int rows, int columns)
-        {
-            return new DenseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<Complex> CreateVector(int size)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Adds a scalar to each element of the vector and stores the result in the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@@ -186,35 +186,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 i => new Complex(distribution.Sample(), distribution.Sample())));
         }
 
-        /// <summary>
-        /// Creates a <c>DiagonalMatrix</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>DiagonalMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return fullyMutable
-                ? (Matrix<Complex>) new SparseMatrix(numberOfRows, numberOfColumns)
-                : new DiagonalMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<Complex> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new SparseVector(size);
-        }
-
         #region Elementary operations
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex/Solvers/IIterativeSolver.cs

@@ -32,6 +32,13 @@ using MathNet.Numerics.LinearAlgebra.Solvers.Status;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// Defines the interface for <see cref="IIterativeSolver"/> classes that solve the matrix equation Ax = b in
     /// an iterative manner.
@@ -64,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        Vector Solve(Matrix matrix, Vector vector);
+        Vector<Complex> Solve(Matrix<Complex> matrix, Vector<Complex> vector);
 
         /// <summary>
         /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
@@ -73,7 +80,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        void Solve(Matrix matrix, Vector input, Vector result);
+        void Solve(Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -82,7 +89,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        Matrix Solve(Matrix matrix, Matrix input);
+        Matrix<Complex> Solve(Matrix<Complex> matrix, Matrix<Complex> input);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -91,6 +98,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        void Solve(Matrix matrix, Matrix input, Matrix result);
+        void Solve(Matrix<Complex> matrix, Matrix<Complex> input, Matrix<Complex> result);
     }
 }

src/Numerics/LinearAlgebra/Complex/Solvers/IIterator.cs

@@ -34,6 +34,13 @@ using System;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// Defines the base interface for iterators that help control an iterative calculation.
     /// </summary>
@@ -80,7 +87,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<Complex> solutionVector, Vector<Complex> sourceVector, Vector<Complex> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs

@@ -212,14 +212,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="vector">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <returns>The result <see cref="Vector"/>, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex> Solve(Matrix<Complex> matrix, Vector<Complex> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -231,7 +231,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <param name="result">The result <see cref="Vector"/>, <c>x</c>.</param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -439,7 +439,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex> matrix, Vector<Complex> residual, Vector<Complex> x, Vector<Complex> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -459,7 +459,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex> result, Vector<Complex> source, Vector<Complex> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -483,7 +483,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <returns>The result <see cref="Matrix"/>, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex> Solve(Matrix<Complex> matrix, Matrix<Complex> input)
         {
             if (matrix == null)
             {
@@ -495,7 +495,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -507,7 +507,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <param name="result">The result <see cref="Matrix"/>, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex> matrix, Matrix<Complex> input, Matrix<Complex> result)
         {
             if (matrix == null)
             {
@@ -531,7 +531,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/CompositeSolver.cs

@@ -38,6 +38,13 @@ using System.Reflection;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// A composite matrix solver. The actual solver is made by a sequence of
     /// matrix solvers. 
@@ -413,14 +420,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex> Solve(Matrix<Complex> matrix, Vector<Complex> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -432,7 +439,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -482,8 +489,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
 
             // Create a copy of the solution and result vectors so we can use them
             // later on
-            var internalInput = (Vector)input.Clone();
-            var internalResult = (Vector)result.Clone();
+            var internalInput = input.Clone();
+            var internalResult = result.Clone();
 
             foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped))
             {
@@ -574,7 +581,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex> Solve(Matrix<Complex> matrix, Matrix<Complex> input)
         {
             if (matrix == null)
             {
@@ -586,7 +593,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -598,7 +605,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex> matrix, Matrix<Complex> input, Matrix<Complex> result)
         {
             if (matrix == null)
             {
@@ -622,7 +629,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs

@@ -261,14 +261,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex> Solve(Matrix<Complex> matrix, Vector<Complex> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -280,7 +280,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -524,7 +524,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex> matrix, Vector<Complex> residual, Vector<Complex> x, Vector<Complex> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -542,7 +542,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex> result, Vector<Complex> source, Vector<Complex> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -583,7 +583,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex> Solve(Matrix<Complex> matrix, Matrix<Complex> input)
         {
             if (matrix == null)
             {
@@ -595,7 +595,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -607,7 +607,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex> matrix, Matrix<Complex> input, Matrix<Complex> result)
         {
             if (matrix == null)
             {
@@ -631,7 +631,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs

@@ -97,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <summary>
         /// The collection of starting vectors which are used as the basis for the Krylov sub-space.
         /// </summary>
-        IList<Vector> _startingVectors;
+        IList<Vector<Complex>> _startingVectors;
 
         /// <summary>
         /// The number of starting vectors used by the algorithm
@@ -235,7 +235,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// Gets or sets a series of orthonormal vectors which will be used as basis for the 
         /// Krylov sub-space.
         /// </summary>
-        public IList<Vector> StartingVectors
+        public IList<Vector<Complex>> StartingVectors
         {
             [DebuggerStepThrough]
             get { return _startingVectors; }
@@ -281,14 +281,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex> Solve(Matrix<Complex> matrix, Vector<Complex> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -300,7 +300,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -639,7 +639,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         ///  the <paramref name="numberOfVariables"/> is smaller than 
         ///  the <paramref name="maximumNumberOfStartingVectors"/>.
         /// </returns>
-        static IList<Vector> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
+        static IList<Vector<Complex>> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
         {
             // Create no more starting vectors than the size of the problem - 1
             // Get random values and then orthogonalize them with
@@ -650,7 +650,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
             // mean = 0 and sd = 1
             var distribution = new Normal();
 
-            Matrix matrix = new DenseMatrix(numberOfVariables, count);
+            var matrix = new DenseMatrix(numberOfVariables, count);
             for (var i = 0; i < matrix.ColumnCount; i++)
             {
                 var samples = new Complex[matrix.RowCount];
@@ -670,10 +670,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
             var orthogonalMatrix = gs.Q;
 
             // Now transfer this to vectors
-            var result = new List<Vector>();
+            var result = new List<Vector<Complex>>();
             for (var i = 0; i < orthogonalMatrix.ColumnCount; i++)
             {
-                result.Add((Vector) orthogonalMatrix.Column(i));
+                result.Add(orthogonalMatrix.Column(i));
 
                 // Normalize the result vector
                 result[i].Multiply(1/result[i].L2Norm(), result[i]);
@@ -688,9 +688,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="arraySize">Number of vectors</param>
         /// <param name="vectorSize">Size of each vector</param>
         /// <returns>Array of random vectors</returns>
-        static Vector[] CreateVectorArray(int arraySize, int vectorSize)
+        static Vector<Complex>[] CreateVectorArray(int arraySize, int vectorSize)
         {
-            var result = new Vector[arraySize];
+            var result = new Vector<Complex>[arraySize];
             for (var i = 0; i < result.Length; i++)
             {
                 result[i] = new DenseVector(vectorSize);
@@ -706,7 +706,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="residual">Residual <see cref="Vector"/> data.</param>
         /// <param name="x">x <see cref="Vector"/> data.</param>
         /// <param name="b">b <see cref="Vector"/> data.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex> matrix, Vector<Complex> residual, Vector<Complex> x, Vector<Complex> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -724,7 +724,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex> result, Vector<Complex> source, Vector<Complex> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -748,7 +748,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex> Solve(Matrix<Complex> matrix, Matrix<Complex> input)
         {
             if (matrix == null)
             {
@@ -760,7 +760,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -772,7 +772,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex> matrix, Matrix<Complex> input, Matrix<Complex> result)
         {
             if (matrix == null)
             {
@@ -796,7 +796,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs

@@ -200,14 +200,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex> Solve(Matrix<Complex> matrix, Vector<Complex> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -219,7 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex> matrix, Vector<Complex> input, Vector<Complex> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -430,7 +430,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex> matrix, Vector<Complex> residual, Vector<Complex> x, Vector<Complex> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -448,7 +448,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex> result, Vector<Complex> source, Vector<Complex> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -482,7 +482,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex> Solve(Matrix<Complex> matrix, Matrix<Complex> input)
         {
             if (matrix == null)
             {
@@ -494,7 +494,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -506,7 +506,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex> matrix, Matrix<Complex> input, Matrix<Complex> result)
         {
             if (matrix == null)
             {
@@ -530,7 +530,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex/Solvers/Iterator.cs

@@ -37,6 +37,13 @@ using System.Linq;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// An iterator that is used to check if an iterative calculation should continue or stop.
     /// </summary>
@@ -212,7 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex> solutionVector, Vector<Complex> sourceVector, Vector<Complex> residualVector)
         {
             if (_stopCriterias.Count == 0)
             {

src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/Diagonal.cs

@@ -72,7 +72,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// The <see cref="Matrix"/> upon which this preconditioner is based.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex> matrix)
         {
             if (matrix == null)
             {
@@ -96,7 +96,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex> Approximate(Vector<Complex> rhs)
         {
             if (rhs == null)
             {
@@ -113,7 +113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -123,7 +123,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex> rhs, Vector<Complex> lhs)
         {
             if (rhs == null)
             {

src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/IPreConditioner.cs

@@ -30,6 +30,13 @@
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// The base interface for preconditioner classes.
     /// </summary>
@@ -54,20 +61,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// Initializes the preconditioner and loads the internal data structures.
         /// </summary>
         /// <param name="matrix">The matrix on which the preconditioner is based.</param>
-        void Initialize(Matrix matrix);
+        void Initialize(Matrix<Complex> matrix);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        Vector Approximate(Vector rhs);
+        Vector<Complex> Approximate(Vector<Complex> rhs);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        void Approximate(Vector rhs, Vector lhs);
+        void Approximate(Vector<Complex> rhs, Vector<Complex> lhs);
     }
 }

src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/Ilutp.cs

@@ -261,9 +261,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<Complex> UpperTriangle()
         {
-            return (Matrix)_upper.Clone();
+            return _upper.Clone();
         }
 
         /// <summary>
@@ -273,9 +273,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<Complex> LowerTriangle()
         {
-            return (Matrix)_lower.Clone();
+            return _lower.Clone();
         }
 
         /// <summary>
@@ -307,7 +307,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException"> If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex> matrix)
         {
             if (matrix == null)
             {
@@ -379,8 +379,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
                 _pivots[i] = i;
             }
 
-            Vector workVector = new DenseVector(sparseMatrix.RowCount);
-            Vector rowVector = new DenseVector(sparseMatrix.ColumnCount);
+            var workVector = new DenseVector(sparseMatrix.RowCount);
+            var rowVector = new DenseVector(sparseMatrix.ColumnCount);
             var indexSorting = new int[sparseMatrix.RowCount];
             
             // spaceLeft = lfilNnz * nnz(A)
@@ -536,7 +536,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// Pivot elements in the <paramref name="row"/> according to internal pivot array
         /// </summary>
         /// <param name="row">Row <see cref="Vector"/> to pivot in</param>
-        private void PivotRow(Vector row)
+        private void PivotRow(Vector<Complex> row)
         {
             var knownPivots = new Dictionary<int, int>();
             
@@ -588,7 +588,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// <param name="matrix">Source <see cref="Matrix"/>.</param>
         /// <param name="firstColumn">First column index to swap</param>
         /// <param name="secondColumn">Second column index to swap</param>
-        private static void SwapColumns(Matrix matrix, int firstColumn, int secondColumn)
+        private static void SwapColumns(Matrix<Complex> matrix, int firstColumn, int secondColumn)
         {
             for (var i = 0; i < matrix.RowCount; i++)
             {
@@ -605,7 +605,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// <param name="upperBound">Sort till upper bound</param>
         /// <param name="sortedIndices">Array with sorted vector indicies</param>
         /// <param name="values">Source <see cref="Vector"/></param>
-        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector<Complex> values)
         {
             // Copy the indices for the values into the array
             for (var i = 0; i < upperBound + 1 - lowerBound; i++)
@@ -630,7 +630,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex> Approximate(Vector<Complex> rhs)
         {
             if (rhs == null)
             {
@@ -647,7 +647,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -657,7 +657,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex> rhs, Vector<Complex> lhs)
         {
             if (rhs == null)
             {
@@ -682,7 +682,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
             // Solve equation here
             // Pivot(vector, result);
             // Solve L*Y = B(piv,:)
-            Vector rowValues = new DenseVector(_lower.RowCount);
+            var rowValues = new DenseVector(_lower.RowCount);
             for (var i = 0; i < _lower.RowCount; i++)
             {
                 _lower.Row(i, rowValues);
@@ -712,7 +712,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
 
             // We have a column pivot so we only need to pivot the
             // end result not the incoming right hand side vector
-            var temp = (Vector)lhs.Clone();
+            var temp = lhs.Clone();
             
             Pivot(temp, lhs);
         }
@@ -722,7 +722,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </summary>
         /// <param name="vector">Source <see cref="Vector"/>.</param>
         /// <param name="result">Result <see cref="Vector"/> after pivoting.</param>
-        private void Pivot(Vector vector, Vector result)
+        private void Pivot(Vector<Complex> vector, Vector<Complex> result)
         {
             for (var i = 0; i < _pivots.Length; i++)
             {

src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/IncompleteLU.cs

@@ -60,7 +60,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// Returns the upper triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<Complex> UpperTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -78,7 +78,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// Returns the lower triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<Complex> LowerTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -105,7 +105,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// <param name="matrix">The matrix upon which the preconditioner is based. </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex> matrix)
         {
             if (matrix == null)
             {
@@ -167,7 +167,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex> Approximate(Vector<Complex> rhs)
         {
             if (rhs == null)
             {
@@ -184,7 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -194,7 +194,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex> rhs, Vector<Complex> lhs)
         {
             if (rhs == null)
             {
@@ -224,7 +224,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
             //     z_i = l_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
             // }
             // NOTE: l_ii should be 1 because u_ii has to be the value
-            Vector rowValues = new DenseVector(_decompositionLU.RowCount);
+            var rowValues = new DenseVector(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
             {
                 // Clear the rowValues 

src/Numerics/LinearAlgebra/Complex/Solvers/Preconditioners/UnitPreconditioner.cs

@@ -33,6 +33,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
     using System;
     using Properties;
 
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// A unit preconditioner. This preconditioner does not actually do anything
     /// it is only used when running an <see cref="IIterativeSolver"/> without
@@ -54,7 +60,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null"/>. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex> matrix)
         {
             if (matrix == null)
             {
@@ -87,7 +93,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         ///     If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         ///   </para>
         /// </exception>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex> rhs, Vector<Complex> lhs)
         {
             if (rhs == null)
             {
@@ -116,7 +122,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
         /// <exception cref="ArgumentException">
         /// If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         /// </exception>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex> Approximate(Vector<Complex> rhs)
         {
             if (rhs == null)
             {
@@ -128,7 +134,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }

src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/DivergenceStopCriterium.cs

@@ -35,6 +35,13 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// Monitors an iterative calculation for signs of divergence.
     /// </summary>
@@ -217,7 +224,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex> solutionVector, Vector<Complex> sourceVector, Vector<Complex> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/FailureStopCriterium.cs

@@ -36,6 +36,13 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// Defines an <see cref="IIterationStopCriterium"/> that monitors residuals for NaN's.
     /// </summary>
@@ -75,7 +82,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex> solutionVector, Vector<Complex> sourceVector, Vector<Complex> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/IIterationStopCriterium.cs

@@ -29,6 +29,13 @@ using MathNet.Numerics.LinearAlgebra.Solvers.StopCriterium;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// The base interface for classes that provide stop criteria for iterative calculations. 
     /// </summary>
@@ -47,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<Complex> solutionVector, Vector<Complex> sourceVector, Vector<Complex> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/IterationCountStopCriterium.cs

@@ -35,6 +35,13 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// Defines an <see cref="IIterationStopCriterium"/> that monitors the numbers of iteration 
     /// steps as stop criterium.
@@ -130,7 +137,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex> solutionVector, Vector<Complex> sourceVector, Vector<Complex> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex/Solvers/StopCriterium/ResidualStopCriterium.cs

@@ -36,6 +36,13 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
 {
+
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
     /// <summary>
     /// Defines an <see cref="IIterationStopCriterium"/> that monitors residuals as stop criterium.
     /// </summary>
@@ -222,7 +229,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex> solutionVector, Vector<Complex> sourceVector, Vector<Complex> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -328,33 +328,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfDiagonalInit(rows, columns, init));
         }
 
-        /// <summary>
-        /// Creates a <c>SparseMatrix</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 SparseMatrix(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>
         /// Returns a new matrix containing the lower triangle of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Complex/SparseVector.cs

@@ -133,39 +133,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new SparseVector(SparseVectorStorage<Complex>.OfInit(length, init));
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex> CreateMatrix(int rows, int columns)
-        {
-            return new SparseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<Complex> CreateVector(int size)
-        {
-            return new SparseVector(size);
-        }
-
         /// <summary>
         /// Conjugates vector and save result to <paramref name="result"/>
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@@ -352,33 +352,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             get { return _values; }
         }
 
-        /// <summary>
-        /// Creates a <c>DenseMatrix</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>DenseMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return new DenseMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<Complex32> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DenseVector.cs

@@ -192,39 +192,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new DenseVector(array);
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex32> CreateMatrix(int rows, int columns)
-        {
-            return new DenseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<Complex32> CreateVector(int size)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Adds a scalar to each element of the vector and stores the result in the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@@ -181,35 +181,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 i => new Complex32((float) distribution.Sample(), (float) distribution.Sample())));
         }
 
-        /// <summary>
-        /// Creates a <c>DiagonalMatrix</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>DiagonalMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return fullyMutable
-                ? (Matrix<Complex32>) new SparseMatrix(numberOfRows, numberOfColumns)
-                : new DiagonalMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<Complex32> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new SparseVector(size);
-        }
-
         #region Elementary operations
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex32/Solvers/IIterativeSolver.cs

@@ -32,6 +32,8 @@ using MathNet.Numerics.LinearAlgebra.Solvers.Status;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
 {
+    using Numerics;
+
     /// <summary>
     /// Defines the interface for <see cref="IIterativeSolver"/> classes that solve the matrix equation Ax = b in
     /// an iterative manner.
@@ -64,7 +66,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        Vector Solve(Matrix matrix, Vector vector);
+        Vector<Complex32> Solve(Matrix<Complex32> matrix, Vector<Complex32> vector);
 
         /// <summary>
         /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
@@ -73,7 +75,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        void Solve(Matrix matrix, Vector input, Vector result);
+        void Solve(Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -82,7 +84,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        Matrix Solve(Matrix matrix, Matrix input);
+        Matrix<Complex32> Solve(Matrix<Complex32> matrix, Matrix<Complex32> input);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -91,6 +93,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        void Solve(Matrix matrix, Matrix input, Matrix result);
+        void Solve(Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result);
     }
 }

src/Numerics/LinearAlgebra/Complex32/Solvers/IIterator.cs

@@ -34,6 +34,8 @@ using System;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
 {
+    using Numerics;
+
     /// <summary>
     /// Defines the base interface for iterators that help control an iterative calculation.
     /// </summary>
@@ -80,7 +82,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<Complex32> solutionVector, Vector<Complex32> sourceVector, Vector<Complex32> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs

@@ -206,14 +206,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="vector">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <returns>The result <see cref="Vector"/>, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex32> Solve(Matrix<Complex32> matrix, Vector<Complex32> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -225,7 +225,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <param name="result">The result <see cref="Vector"/>, <c>x</c>.</param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -433,7 +433,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex32> matrix, Vector<Complex32> residual, Vector<Complex32> x, Vector<Complex32> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -453,7 +453,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex32> result, Vector<Complex32> source, Vector<Complex32> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -477,7 +477,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <returns>The result <see cref="Matrix"/>, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex32> Solve(Matrix<Complex32> matrix, Matrix<Complex32> input)
         {
             if (matrix == null)
             {
@@ -489,7 +489,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -501,7 +501,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <param name="result">The result <see cref="Matrix"/>, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result)
         {
             if (matrix == null)
             {
@@ -525,7 +525,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/CompositeSolver.cs

@@ -38,6 +38,8 @@ using System.Reflection;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
 {
+    using Numerics;
+
     /// <summary>
     /// A composite matrix solver. The actual solver is made by a sequence of
     /// matrix solvers. 
@@ -413,14 +415,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex32> Solve(Matrix<Complex32> matrix, Vector<Complex32> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -432,7 +434,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -482,8 +484,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
 
             // Create a copy of the solution and result vectors so we can use them
             // later on
-            var internalInput = (Vector)input.Clone();
-            var internalResult = (Vector)result.Clone();
+            var internalInput = input.Clone();
+            var internalResult = result.Clone();
 
             foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped))
             {
@@ -574,7 +576,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex32> Solve(Matrix<Complex32> matrix, Matrix<Complex32> input)
         {
             if (matrix == null)
             {
@@ -586,7 +588,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -598,7 +600,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result)
         {
             if (matrix == null)
             {
@@ -622,7 +624,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs

@@ -255,14 +255,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex32> Solve(Matrix<Complex32> matrix, Vector<Complex32> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -274,7 +274,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -523,7 +523,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex32> matrix, Vector<Complex32> residual, Vector<Complex32> x, Vector<Complex32> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -541,7 +541,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex32> result, Vector<Complex32> source, Vector<Complex32> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -582,7 +582,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex32> Solve(Matrix<Complex32> matrix, Matrix<Complex32> input)
         {
             if (matrix == null)
             {
@@ -594,7 +594,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -606,7 +606,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result)
         {
             if (matrix == null)
             {
@@ -630,7 +630,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs

@@ -91,7 +91,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <summary>
         /// The collection of starting vectors which are used as the basis for the Krylov sub-space.
         /// </summary>
-        IList<Vector> _startingVectors;
+        IList<Vector<Complex32>> _startingVectors;
 
         /// <summary>
         /// The number of starting vectors used by the algorithm
@@ -229,7 +229,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// Gets or sets a series of orthonormal vectors which will be used as basis for the 
         /// Krylov sub-space.
         /// </summary>
-        public IList<Vector> StartingVectors
+        public IList<Vector<Complex32>> StartingVectors
         {
             [DebuggerStepThrough]
             get { return _startingVectors; }
@@ -275,14 +275,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex32> Solve(Matrix<Complex32> matrix, Vector<Complex32> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -294,7 +294,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -638,7 +638,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         ///  the <paramref name="numberOfVariables"/> is smaller than 
         ///  the <paramref name="maximumNumberOfStartingVectors"/>.
         /// </returns>
-        static IList<Vector> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
+        static IList<Vector<Complex32>> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
         {
             // Create no more starting vectors than the size of the problem - 1
             // Get random values and then orthogonalize them with
@@ -649,7 +649,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
             // mean = 0 and sd = 1
             var distribution = new Normal();
 
-            Matrix matrix = new DenseMatrix(numberOfVariables, count);
+            var matrix = new DenseMatrix(numberOfVariables, count);
             for (var i = 0; i < matrix.ColumnCount; i++)
             {
                 var samples = new Complex32[matrix.RowCount];
@@ -669,10 +669,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
             var orthogonalMatrix = gs.Q;
 
             // Now transfer this to vectors
-            var result = new List<Vector>();
+            var result = new List<Vector<Complex32>>();
             for (var i = 0; i < orthogonalMatrix.ColumnCount; i++)
             {
-                result.Add((Vector) orthogonalMatrix.Column(i));
+                result.Add(orthogonalMatrix.Column(i));
 
                 // Normalize the result vector
                 result[i].Multiply(1/result[i].L2Norm().Real, result[i]);
@@ -687,9 +687,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="arraySize">Number of vectors</param>
         /// <param name="vectorSize">Size of each vector</param>
         /// <returns>Array of random vectors</returns>
-        static Vector[] CreateVectorArray(int arraySize, int vectorSize)
+        static Vector<Complex32>[] CreateVectorArray(int arraySize, int vectorSize)
         {
-            var result = new Vector[arraySize];
+            var result = new Vector<Complex32>[arraySize];
             for (var i = 0; i < result.Length; i++)
             {
                 result[i] = new DenseVector(vectorSize);
@@ -705,7 +705,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="residual">Residual <see cref="Vector"/> data.</param>
         /// <param name="x">x <see cref="Vector"/> data.</param>
         /// <param name="b">b <see cref="Vector"/> data.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex32> matrix, Vector<Complex32> residual, Vector<Complex32> x, Vector<Complex32> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -723,7 +723,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex32> result, Vector<Complex32> source, Vector<Complex32> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -747,7 +747,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex32> Solve(Matrix<Complex32> matrix, Matrix<Complex32> input)
         {
             if (matrix == null)
             {
@@ -759,7 +759,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -771,7 +771,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result)
         {
             if (matrix == null)
             {
@@ -795,7 +795,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs

@@ -194,14 +194,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<Complex32> Solve(Matrix<Complex32> matrix, Vector<Complex32> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -213,7 +213,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<Complex32> matrix, Vector<Complex32> input, Vector<Complex32> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -429,7 +429,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<Complex32> matrix, Vector<Complex32> residual, Vector<Complex32> x, Vector<Complex32> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -447,7 +447,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<Complex32> result, Vector<Complex32> source, Vector<Complex32> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -481,7 +481,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<Complex32> Solve(Matrix<Complex32> matrix, Matrix<Complex32> input)
         {
             if (matrix == null)
             {
@@ -493,7 +493,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -505,7 +505,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<Complex32> matrix, Matrix<Complex32> input, Matrix<Complex32> result)
         {
             if (matrix == null)
             {
@@ -529,7 +529,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Complex32/Solvers/Iterator.cs

@@ -37,6 +37,8 @@ using System.Linq;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
 {
+    using Numerics;
+
     /// <summary>
     /// An iterator that is used to check if an iterative calculation should continue or stop.
     /// </summary>
@@ -212,7 +214,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex32> solutionVector, Vector<Complex32> sourceVector, Vector<Complex32> residualVector)
         {
             if (_stopCriterias.Count == 0)
             {

src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/Diagonal.cs

@@ -67,7 +67,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// The <see cref="Matrix"/> upon which this preconditioner is based.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex32> matrix)
         {
             if (matrix == null)
             {
@@ -91,7 +91,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex32> Approximate(Vector<Complex32> rhs)
         {
             if (rhs == null)
             {
@@ -108,7 +108,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -118,7 +118,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex32> rhs, Vector<Complex32> lhs)
         {
             if (rhs == null)
             {

src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/IPreConditioner.cs

@@ -30,6 +30,8 @@
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
 {
+    using Numerics;
+
     /// <summary>
     /// The base interface for preconditioner classes.
     /// </summary>
@@ -54,20 +56,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// Initializes the preconditioner and loads the internal data structures.
         /// </summary>
         /// <param name="matrix">The matrix on which the preconditioner is based.</param>
-        void Initialize(Matrix matrix);
+        void Initialize(Matrix<Complex32> matrix);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        Vector Approximate(Vector rhs);
+        Vector<Complex32> Approximate(Vector<Complex32> rhs);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        void Approximate(Vector rhs, Vector lhs);
+        void Approximate(Vector<Complex32> rhs, Vector<Complex32> lhs);
     }
 }

src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/Ilutp.cs

@@ -256,9 +256,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<Complex32> UpperTriangle()
         {
-            return (Matrix)_upper.Clone();
+            return _upper.Clone();
         }
 
         /// <summary>
@@ -268,9 +268,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<Complex32> LowerTriangle()
         {
-            return (Matrix)_lower.Clone();
+            return _lower.Clone();
         }
 
         /// <summary>
@@ -302,7 +302,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException"> If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex32> matrix)
         {
             if (matrix == null)
             {
@@ -374,8 +374,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
                 _pivots[i] = i;
             }
 
-            Vector workVector = new DenseVector(sparseMatrix.RowCount);
-            Vector rowVector = new DenseVector(sparseMatrix.ColumnCount);
+            var workVector = new DenseVector(sparseMatrix.RowCount);
+            var rowVector = new DenseVector(sparseMatrix.ColumnCount);
             var indexSorting = new int[sparseMatrix.RowCount];
             
             // spaceLeft = lfilNnz * nnz(A)
@@ -531,7 +531,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// Pivot elements in the <paramref name="row"/> according to internal pivot array
         /// </summary>
         /// <param name="row">Row <see cref="Vector"/> to pivot in</param>
-        private void PivotRow(Vector row)
+        private void PivotRow(Vector<Complex32> row)
         {
             var knownPivots = new Dictionary<int, int>();
             
@@ -583,7 +583,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <param name="matrix">Source <see cref="Matrix"/>.</param>
         /// <param name="firstColumn">First column index to swap</param>
         /// <param name="secondColumn">Second column index to swap</param>
-        private static void SwapColumns(Matrix matrix, int firstColumn, int secondColumn)
+        private static void SwapColumns(Matrix<Complex32> matrix, int firstColumn, int secondColumn)
         {
             for (var i = 0; i < matrix.RowCount; i++)
             {
@@ -600,7 +600,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <param name="upperBound">Sort till upper bound</param>
         /// <param name="sortedIndices">Array with sorted vector indicies</param>
         /// <param name="values">Source <see cref="Vector"/></param>
-        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector<Complex32> values)
         {
             // Copy the indices for the values into the array
             for (var i = 0; i < upperBound + 1 - lowerBound; i++)
@@ -625,7 +625,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex32> Approximate(Vector<Complex32> rhs)
         {
             if (rhs == null)
             {
@@ -642,7 +642,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -652,7 +652,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex32> rhs, Vector<Complex32> lhs)
         {
             if (rhs == null)
             {
@@ -677,7 +677,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
             // Solve equation here
             // Pivot(vector, result);
             // Solve L*Y = B(piv,:)
-            Vector rowValues = new DenseVector(_lower.RowCount);
+            var rowValues = new DenseVector(_lower.RowCount);
             for (var i = 0; i < _lower.RowCount; i++)
             {
                 _lower.Row(i, rowValues);
@@ -707,7 +707,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
 
             // We have a column pivot so we only need to pivot the
             // end result not the incoming right hand side vector
-            var temp = (Vector)lhs.Clone();
+            var temp = lhs.Clone();
             
             Pivot(temp, lhs);
         }
@@ -717,7 +717,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </summary>
         /// <param name="vector">Source <see cref="Vector"/>.</param>
         /// <param name="result">Result <see cref="Vector"/> after pivoting.</param>
-        private void Pivot(Vector vector, Vector result)
+        private void Pivot(Vector<Complex32> vector, Vector<Complex32> result)
         {
             for (var i = 0; i < _pivots.Length; i++)
             {

src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/IlutpElementSorter.cs

@@ -30,6 +30,8 @@
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
 {
+    using Numerics;
+
     /// <summary>
     /// An element sort algorithm for the <see cref="Ilutp"/> class.
     /// </summary>
@@ -47,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <param name="upperBound">The stopping index.</param>
         /// <param name="sortedIndices">An array that will contain the sorted indices once the algorithm finishes.</param>
         /// <param name="values">The <see cref="Vector"/> that contains the values that need to be sorted.</param>
-        public static void SortDoubleIndicesDecreasing(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        public static void SortDoubleIndicesDecreasing(int lowerBound, int upperBound, int[] sortedIndices, Vector<Complex32> values)
         {
             // Move all the indices that we're interested in to the beginning of the
             // array. Ignore the rest of the indices.
@@ -73,7 +75,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <param name="upperBound">The stopping index.</param>
         /// <param name="sortedIndices">An array that will contain the sorted indices once the algorithm finishes.</param>
         /// <param name="values">The <see cref="Vector"/> that contains the values that need to be sorted.</param>
-        private static void HeapSortDoublesIndices(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        private static void HeapSortDoublesIndices(int lowerBound, int upperBound, int[] sortedIndices, Vector<Complex32> values)
         {
             var start = ((upperBound - lowerBound + 1) / 2) - 1 + lowerBound;
             var end = (upperBound - lowerBound + 1) - 1 + lowerBound;
@@ -95,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <param name="count">Length of <paramref name="values"/></param>
         /// <param name="sortedIndices">Indicies of <paramref name="values"/></param>
         /// <param name="values">Target <see cref="Vector"/></param>
-        private static void BuildDoubleIndexHeap(int start, int count, int[] sortedIndices, Vector values)
+        private static void BuildDoubleIndexHeap(int start, int count, int[] sortedIndices, Vector<Complex32> values)
         {
             while (start >= 0)
             {
@@ -111,7 +113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <param name="values">Target <see cref="Vector"/></param>
         /// <param name="begin">Root position</param>
         /// <param name="count">Length of <paramref name="values"/></param>
-        private static void SiftDoubleIndices(int[] sortedIndices, Vector values, int begin, int count)
+        private static void SiftDoubleIndices(int[] sortedIndices, Vector<Complex32> values, int begin, int count)
         {
             var root = begin;
 

src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/IncompleteLU.cs

@@ -55,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// Returns the upper triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<Complex32> UpperTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -73,7 +73,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// Returns the lower triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<Complex32> LowerTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -100,7 +100,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <param name="matrix">The matrix upon which the preconditioner is based. </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex32> matrix)
         {
             if (matrix == null)
             {
@@ -162,7 +162,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex32> Approximate(Vector<Complex32> rhs)
         {
             if (rhs == null)
             {
@@ -179,7 +179,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -189,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex32> rhs, Vector<Complex32> lhs)
         {
             if (rhs == null)
             {
@@ -219,7 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
             //     z_i = l_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
             // }
             // NOTE: l_ii should be 1 because u_ii has to be the value
-            Vector rowValues = new DenseVector(_decompositionLU.RowCount);
+            var rowValues = new DenseVector(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
             {
                 // Clear the rowValues 

src/Numerics/LinearAlgebra/Complex32/Solvers/Preconditioners/UnitPreconditioner.cs

@@ -32,6 +32,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
 {
     using System;
     using Properties;
+    using Numerics;
 
     /// <summary>
     /// A unit preconditioner. This preconditioner does not actually do anything
@@ -54,7 +55,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null"/>. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<Complex32> matrix)
         {
             if (matrix == null)
             {
@@ -87,7 +88,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         ///     If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         ///   </para>
         /// </exception>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<Complex32> rhs, Vector<Complex32> lhs)
         {
             if (rhs == null)
             {
@@ -116,7 +117,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
         /// <exception cref="ArgumentException">
         /// If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         /// </exception>
-        public Vector Approximate(Vector rhs)
+        public Vector<Complex32> Approximate(Vector<Complex32> rhs)
         {
             if (rhs == null)
             {
@@ -128,7 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }

src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/DivergenceStopCriterium.cs

@@ -35,6 +35,8 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
 {
+    using Numerics;
+
     /// <summary>
     /// Monitors an iterative calculation for signs of divergence.
     /// </summary>
@@ -217,7 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex32> solutionVector, Vector<Complex32> sourceVector, Vector<Complex32> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/FailureStopCriterium.cs

@@ -36,6 +36,8 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
 {
+    using Numerics;
+
     /// <summary>
     /// Defines an <see cref="IIterationStopCriterium"/> that monitors residuals for NaN's.
     /// </summary>
@@ -75,7 +77,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex32> solutionVector, Vector<Complex32> sourceVector, Vector<Complex32> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/IIterationStopCriterium.cs

@@ -29,6 +29,8 @@ using MathNet.Numerics.LinearAlgebra.Solvers.StopCriterium;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
 {
+    using Numerics;
+
     /// <summary>
     /// The base interface for classes that provide stop criteria for iterative calculations. 
     /// </summary>
@@ -47,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<Complex32> solutionVector, Vector<Complex32> sourceVector, Vector<Complex32> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/IterationCountStopCriterium.cs

@@ -35,6 +35,8 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
 {
+    using Numerics;
+
     /// <summary>
     /// Defines an <see cref="IIterationStopCriterium"/> that monitors the numbers of iteration 
     /// steps as stop criterium.
@@ -130,7 +132,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex32> solutionVector, Vector<Complex32> sourceVector, Vector<Complex32> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex32/Solvers/StopCriterium/ResidualStopCriterium.cs

@@ -36,6 +36,8 @@ using System.Diagnostics;
 
 namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
 {
+    using Numerics;
+
     /// <summary>
     /// Defines an <see cref="IIterationStopCriterium"/> that monitors residuals as stop criterium.
     /// </summary>
@@ -222,7 +224,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<Complex32> solutionVector, Vector<Complex32> sourceVector, Vector<Complex32> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -323,33 +323,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfDiagonalInit(rows, columns, init));
         }
 
-        /// <summary>
-        /// Creates a <c>SparseMatrix</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<Complex32> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return new SparseMatrix(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<Complex32> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new SparseVector(size);
-        }
-
         /// <summary>
         /// Returns a new matrix containing the lower triangle of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/SparseVector.cs

@@ -128,39 +128,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new SparseVector(SparseVectorStorage<Complex32>.OfInit(length, init));
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<Complex32> CreateMatrix(int rows, int columns)
-        {
-            return new SparseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<Complex32> CreateVector(int size)
-        {
-            return new SparseVector(size);
-        }
-
         /// <summary>
         /// Conjugates vector and save result to <paramref name="result"/>
         /// </summary>

src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@@ -349,33 +349,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             get { return _values; }
         }
 
-        /// <summary>
-        /// Creates a <c>DenseMatrix</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>DenseMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return new DenseMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<double> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Double/DenseVector.cs

@@ -192,39 +192,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new DenseVector(array);
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<double> CreateMatrix(int rows, int columns)
-        {
-            return new DenseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<double> CreateVector(int size)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Adds a scalar to each element of the vector and stores the result in the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@@ -180,35 +180,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 i => distribution.Sample()));
         }
 
-        /// <summary>
-        /// Creates a <c>DiagonalMatrix</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>DiagonalMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<double> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return fullyMutable
-                ? (Matrix<double>) new SparseMatrix(numberOfRows, numberOfColumns)
-                : new DiagonalMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<double> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new SparseVector(size);
-        }
-
         #region Elementary operations
 
         /// <summary>

src/Numerics/LinearAlgebra/Double/Solvers/IIterativeSolver.cs

@@ -64,7 +64,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        Vector Solve(Matrix matrix, Vector vector);
+        Vector<double> Solve(Matrix<double> matrix, Vector<double> vector);
 
         /// <summary>
         /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
@@ -73,7 +73,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        void Solve(Matrix matrix, Vector input, Vector result);
+        void Solve(Matrix<double> matrix, Vector<double> input, Vector<double> result);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -82,7 +82,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        Matrix Solve(Matrix matrix, Matrix input);
+        Matrix<double> Solve(Matrix<double> matrix, Matrix<double> input);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -91,6 +91,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        void Solve(Matrix matrix, Matrix input, Matrix result);
+        void Solve(Matrix<double> matrix, Matrix<double> input, Matrix<double> result);
     }
 }

src/Numerics/LinearAlgebra/Double/Solvers/IIterator.cs

@@ -80,7 +80,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<double> solutionVector, Vector<double> sourceVector, Vector<double> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs

@@ -204,7 +204,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="vector">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <returns>The result <see cref="Vector"/>, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<double> Solve(Matrix<double> matrix, Vector<double> vector)
         {
             if (vector == null)
             {
@@ -223,7 +223,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <param name="result">The result <see cref="Vector"/>, <c>x</c>.</param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<double> matrix, Vector<double> input, Vector<double> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -278,7 +278,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
             // Compute r_0 = b - Ax_0 for some initial guess x_0
             // In this case we take x_0 = vector
             // This is basically a SAXPY so it could be made a lot faster
-            Vector residuals = new DenseVector(matrix.RowCount);
+            var residuals = new DenseVector(matrix.RowCount);
             CalculateTrueResidual(matrix, residuals, result, input);
 
             // Choose r~ (for example, r~ = r_0)
@@ -287,13 +287,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
             // create seven temporary vectors needed to hold temporary
             // coefficients. All vectors are mangled in each iteration.
             // These are defined here to prevent stressing the garbage collector
-            Vector vecP = new DenseVector(residuals.Count);
-            Vector vecPdash = new DenseVector(residuals.Count);
-            Vector nu = new DenseVector(residuals.Count);
-            Vector vecS = new DenseVector(residuals.Count);
-            Vector vecSdash = new DenseVector(residuals.Count);
-            Vector temp = new DenseVector(residuals.Count);
-            Vector temp2 = new DenseVector(residuals.Count);
+            var vecP = new DenseVector(residuals.Count);
+            var vecPdash = new DenseVector(residuals.Count);
+            var nu = new DenseVector(residuals.Count);
+            var vecS = new DenseVector(residuals.Count);
+            var vecSdash = new DenseVector(residuals.Count);
+            var temp = new DenseVector(residuals.Count);
+            var temp2 = new DenseVector(residuals.Count);
 
             // create some temporary double variables that are needed
             // to hold values in between iterations
@@ -431,7 +431,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        private static void CalculateTrueResidual(Matrix<double> matrix, Vector<double> residual, Vector<double> x, Vector<double> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -451,7 +451,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        private bool ShouldContinue(int iterationNumber, Vector<double> result, Vector<double> source, Vector<double> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -475,7 +475,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <returns>The result <see cref="Matrix"/>, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<double> Solve(Matrix<double> matrix, Matrix<double> input)
         {
             if (matrix == null)
             {
@@ -487,7 +487,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -499,7 +499,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <param name="result">The result <see cref="Matrix"/>, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<double> matrix, Matrix<double> input, Matrix<double> result)
         {
             if (matrix == null)
             {
@@ -523,7 +523,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Double/Solvers/Iterative/CompositeSolver.cs

@@ -410,14 +410,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<double> Solve(Matrix<double> matrix, Vector<double> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -429,7 +429,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<double> matrix, Vector<double> input, Vector<double> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -479,8 +479,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
 
             // Create a copy of the solution and result vectors so we can use them
             // later on
-            var internalInput = (Vector)input.Clone();
-            var internalResult = (Vector)result.Clone();
+            var internalInput = input.Clone();
+            var internalResult = result.Clone();
 
             foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped))
             {
@@ -571,7 +571,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<double> Solve(Matrix<double> matrix, Matrix<double> input)
         {
             if (matrix == null)
             {
@@ -583,7 +583,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -595,7 +595,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<double> matrix, Matrix<double> input, Matrix<double> result)
         {
             if (matrix == null)
             {
@@ -619,7 +619,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs

@@ -259,14 +259,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<double> Solve(Matrix<double> matrix, Vector<double> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -278,7 +278,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<double> matrix, Vector<double> input, Vector<double> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -332,11 +332,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
 
             // x_0 is initial guess
             // Take x_0 = 0
-            Vector xtemp = new DenseVector(input.Count);
+            var xtemp = new DenseVector(input.Count);
 
             // r_0 = b - Ax_0
             // This is basically a SAXPY so it could be made a lot faster
-            Vector residuals = new DenseVector(matrix.RowCount);
+            var residuals = new DenseVector(matrix.RowCount);
             CalculateTrueResidual(matrix, residuals, xtemp, input);
 
             // Define the temporary scalars
@@ -344,26 +344,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
 
             // Define the temporary vectors
             // rDash_0 = r_0
-            Vector rdash = DenseVector.OfVector(residuals);
+            var rdash = DenseVector.OfVector(residuals);
 
             // t_-1 = 0
-            Vector t = new DenseVector(residuals.Count);
-            Vector t0 = new DenseVector(residuals.Count);
+            var t = new DenseVector(residuals.Count);
+            var t0 = new DenseVector(residuals.Count);
 
             // w_-1 = 0
-            Vector w = new DenseVector(residuals.Count);
+            var w = new DenseVector(residuals.Count);
 
             // Define the remaining temporary vectors
-            Vector c = new DenseVector(residuals.Count);
-            Vector p = new DenseVector(residuals.Count);
-            Vector s = new DenseVector(residuals.Count);
-            Vector u = new DenseVector(residuals.Count);
-            Vector y = new DenseVector(residuals.Count);
-            Vector z = new DenseVector(residuals.Count);
+            var c = new DenseVector(residuals.Count);
+            var p = new DenseVector(residuals.Count);
+            var s = new DenseVector(residuals.Count);
+            var u = new DenseVector(residuals.Count);
+            var y = new DenseVector(residuals.Count);
+            var z = new DenseVector(residuals.Count);
 
-            Vector temp = new DenseVector(residuals.Count);
-            Vector temp2 = new DenseVector(residuals.Count);
-            Vector temp3 = new DenseVector(residuals.Count);
+            var temp = new DenseVector(residuals.Count);
+            var temp2 = new DenseVector(residuals.Count);
+            var temp3 = new DenseVector(residuals.Count);
 
             // for (k = 0, 1, .... )
             var iterationNumber = 0;
@@ -527,7 +527,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        private static void CalculateTrueResidual(Matrix<double> matrix, Vector<double> residual, Vector<double> x, Vector<double> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -545,7 +545,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        private bool ShouldContinue(int iterationNumber, Vector<double> result, Vector<double> source, Vector<double> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -586,7 +586,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<double> Solve(Matrix<double> matrix, Matrix<double> input)
         {
             if (matrix == null)
             {
@@ -598,7 +598,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -610,7 +610,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<double> matrix, Matrix<double> input, Matrix<double> result)
         {
             if (matrix == null)
             {
@@ -634,7 +634,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs

@@ -89,7 +89,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <summary>
         /// The collection of starting vectors which are used as the basis for the Krylov sub-space.
         /// </summary>
-        private IList<Vector> _startingVectors;
+        private IList<Vector<double>> _startingVectors;
 
         /// <summary>
         /// The number of starting vectors used by the algorithm
@@ -227,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// Gets or sets a series of orthonormal vectors which will be used as basis for the 
         /// Krylov sub-space.
         /// </summary>
-        public IList<Vector> StartingVectors
+        public IList<Vector<double>> StartingVectors
         {
             [DebuggerStepThrough]
             get
@@ -279,14 +279,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<double> Solve(Matrix<double> matrix, Vector<double> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -298,7 +298,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<double> matrix, Vector<double> input, Vector<double> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -642,7 +642,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         ///  the <paramref name="numberOfVariables"/> is smaller than 
         ///  the <paramref name="maximumNumberOfStartingVectors"/>.
         /// </returns>
-        private static IList<Vector> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
+        private static IList<Vector<double>> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
         {
             // Create no more starting vectors than the size of the problem - 1
             // Get random values and then orthogonalize them with
@@ -653,7 +653,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
             // mean = 0 and sd = 1
             var distribution = new Normal();
 
-            Matrix matrix = new DenseMatrix(numberOfVariables, count);
+            var matrix = new DenseMatrix(numberOfVariables, count);
             for (var i = 0; i < matrix.ColumnCount; i++)
             {
                 var samples = distribution.Samples().Take(matrix.RowCount).ToArray();
@@ -667,10 +667,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
             var orthogonalMatrix = gs.Q;
 
             // Now transfer this to vectors
-            var result = new List<Vector>();
+            var result = new List<Vector<double>>();
             for (var i = 0; i < orthogonalMatrix.ColumnCount; i++)
             {
-                result.Add((Vector)orthogonalMatrix.Column(i));
+                result.Add(orthogonalMatrix.Column(i));
                 
                 // Normalize the result vector
                 result[i].Multiply(1 / result[i].L2Norm(), result[i]);
@@ -685,9 +685,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="arraySize">Number of vectors</param>
         /// <param name="vectorSize">Size of each vector</param>
         /// <returns>Array of random vectors</returns>
-        private static Vector[] CreateVectorArray(int arraySize, int vectorSize)
+        private static Vector<double>[] CreateVectorArray(int arraySize, int vectorSize)
         {
-            var result = new Vector[arraySize];
+            var result = new Vector<double>[arraySize];
             for (var i = 0; i < result.Length; i++)
             {
                 result[i] = new DenseVector(vectorSize);
@@ -703,7 +703,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="residual">Residual <see cref="Vector"/> data.</param>
         /// <param name="x">x <see cref="Vector"/> data.</param>
         /// <param name="b">b <see cref="Vector"/> data.</param>
-        private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        private static void CalculateTrueResidual(Matrix<double> matrix, Vector<double> residual, Vector<double> x, Vector<double> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -721,7 +721,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        private bool ShouldContinue(int iterationNumber, Vector<double> result, Vector<double> source, Vector<double> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -745,7 +745,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<double> Solve(Matrix<double> matrix, Matrix<double> input)
         {
             if (matrix == null)
             {
@@ -757,7 +757,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -769,7 +769,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<double> matrix, Matrix<double> input, Matrix<double> result)
         {
             if (matrix == null)
             {

src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs

@@ -192,14 +192,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<double> Solve(Matrix<double> matrix, Vector<double> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -211,7 +211,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<double> matrix, Vector<double> input, Vector<double> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -427,7 +427,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        private static void CalculateTrueResidual(Matrix<double> matrix, Vector<double> residual, Vector<double> x, Vector<double> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -445,7 +445,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        private bool ShouldContinue(int iterationNumber, Vector<double> result, Vector<double> source, Vector<double> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -479,7 +479,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<double> Solve(Matrix<double> matrix, Matrix<double> input)
         {
             if (matrix == null)
             {
@@ -491,7 +491,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -503,7 +503,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<double> matrix, Matrix<double> input, Matrix<double> result)
         {
             if (matrix == null)
             {
@@ -527,7 +527,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Double/Solvers/Iterator.cs

@@ -212,7 +212,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<double> solutionVector, Vector<double> sourceVector, Vector<double> residualVector)
         {
             if (_stopCriterias.Count == 0)
             {

src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Diagonal.cs

@@ -66,7 +66,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// The <see cref="Matrix"/> upon which this preconditioner is based.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<double> matrix)
         {
             if (matrix == null)
             {
@@ -90,7 +90,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<double> Approximate(Vector<double> rhs)
         {
             if (rhs == null)
             {
@@ -107,7 +107,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -117,7 +117,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<double> rhs, Vector<double> lhs)
         {
             if (rhs == null)
             {

src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/IPreConditioner.cs

@@ -54,20 +54,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// Initializes the preconditioner and loads the internal data structures.
         /// </summary>
         /// <param name="matrix">The matrix on which the preconditioner is based.</param>
-        void Initialize(Matrix matrix);
+        void Initialize(Matrix<double> matrix);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        Vector Approximate(Vector rhs);
+        Vector<double> Approximate(Vector<double> rhs);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        void Approximate(Vector rhs, Vector lhs);
+        void Approximate(Vector<double> rhs, Vector<double> lhs);
     }
 }

src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Ilutp.cs

@@ -255,9 +255,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<double> UpperTriangle()
         {
-            return (Matrix)_upper.Clone();
+            return _upper.Clone();
         }
 
         /// <summary>
@@ -267,9 +267,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<double> LowerTriangle()
         {
-            return (Matrix)_lower.Clone();
+            return _lower.Clone();
         }
 
         /// <summary>
@@ -301,7 +301,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException"> If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<double> matrix)
         {
             if (matrix == null)
             {
@@ -373,8 +373,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
                 _pivots[i] = i;
             }
 
-            Vector workVector = new DenseVector(sparseMatrix.RowCount);
-            Vector rowVector = new DenseVector(sparseMatrix.ColumnCount);
+            var workVector = new DenseVector(sparseMatrix.RowCount);
+            var rowVector = new DenseVector(sparseMatrix.ColumnCount);
             var indexSorting = new int[sparseMatrix.RowCount];
             
             // spaceLeft = lfilNnz * nnz(A)
@@ -530,7 +530,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// Pivot elements in the <paramref name="row"/> according to internal pivot array
         /// </summary>
         /// <param name="row">Row <see cref="Vector"/> to pivot in</param>
-        private void PivotRow(Vector row)
+        private void PivotRow(Vector<double> row)
         {
             var knownPivots = new Dictionary<int, int>();
             
@@ -582,7 +582,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <param name="matrix">Source <see cref="Matrix"/>.</param>
         /// <param name="firstColumn">First column index to swap</param>
         /// <param name="secondColumn">Second column index to swap</param>
-        private static void SwapColumns(Matrix matrix, int firstColumn, int secondColumn)
+        private static void SwapColumns(Matrix<double> matrix, int firstColumn, int secondColumn)
         {
             for (var i = 0; i < matrix.RowCount; i++)
             {
@@ -599,7 +599,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <param name="upperBound">Sort till upper bound</param>
         /// <param name="sortedIndices">Array with sorted vector indicies</param>
         /// <param name="values">Source <see cref="Vector"/></param>
-        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector<double> values)
         {
             // Copy the indices for the values into the array
             for (var i = 0; i < upperBound + 1 - lowerBound; i++)
@@ -624,7 +624,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<double> Approximate(Vector<double> rhs)
         {
             if (rhs == null)
             {
@@ -641,7 +641,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -651,7 +651,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<double> rhs, Vector<double> lhs)
         {
             if (rhs == null)
             {
@@ -676,7 +676,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
             // Solve equation here
             // Pivot(vector, result);
             // Solve L*Y = B(piv,:)
-            Vector rowValues = new DenseVector(_lower.RowCount);
+            var rowValues = new DenseVector(_lower.RowCount);
             for (var i = 0; i < _lower.RowCount; i++)
             {
                 _lower.Row(i, rowValues);
@@ -706,7 +706,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
 
             // We have a column pivot so we only need to pivot the
             // end result not the incoming right hand side vector
-            var temp = (Vector)lhs.Clone();
+            var temp = lhs.Clone();
             
             Pivot(temp, lhs);
         }
@@ -716,7 +716,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </summary>
         /// <param name="vector">Source <see cref="Vector"/>.</param>
         /// <param name="result">Result <see cref="Vector"/> after pivoting.</param>
-        private void Pivot(Vector vector, Vector result)
+        private void Pivot(Vector<double> vector, Vector<double> result)
         {
             for (var i = 0; i < _pivots.Length; i++)
             {

src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/IlutpElementSorter.cs

@@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <param name="upperBound">The stopping index.</param>
         /// <param name="sortedIndices">An array that will contain the sorted indices once the algorithm finishes.</param>
         /// <param name="values">The <see cref="Vector"/> that contains the values that need to be sorted.</param>
-        public static void SortDoubleIndicesDecreasing(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        public static void SortDoubleIndicesDecreasing(int lowerBound, int upperBound, int[] sortedIndices, Vector<double> values)
         {
             // Move all the indices that we're interested in to the beginning of the
             // array. Ignore the rest of the indices.
@@ -73,7 +73,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <param name="upperBound">The stopping index.</param>
         /// <param name="sortedIndices">An array that will contain the sorted indices once the algorithm finishes.</param>
         /// <param name="values">The <see cref="Vector"/> that contains the values that need to be sorted.</param>
-        private static void HeapSortDoublesIndices(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        private static void HeapSortDoublesIndices(int lowerBound, int upperBound, int[] sortedIndices, Vector<double> values)
         {
             var start = ((upperBound - lowerBound + 1) / 2) - 1 + lowerBound;
             var end = (upperBound - lowerBound + 1) - 1 + lowerBound;
@@ -95,7 +95,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <param name="count">Length of <paramref name="values"/></param>
         /// <param name="sortedIndices">Indicies of <paramref name="values"/></param>
         /// <param name="values">Target <see cref="Vector"/></param>
-        private static void BuildDoubleIndexHeap(int start, int count, int[] sortedIndices, Vector values)
+        private static void BuildDoubleIndexHeap(int start, int count, int[] sortedIndices, Vector<double> values)
         {
             while (start >= 0)
             {
@@ -111,7 +111,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <param name="values">Target <see cref="Vector"/></param>
         /// <param name="begin">Root position</param>
         /// <param name="count">Length of <paramref name="values"/></param>
-        private static void SiftDoubleIndices(int[] sortedIndices, Vector values, int begin, int count)
+        private static void SiftDoubleIndices(int[] sortedIndices, Vector<double> values, int begin, int count)
         {
             var root = begin;
 

src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/IncompleteLU.cs

@@ -54,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// Returns the upper triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<double> UpperTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -72,7 +72,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// Returns the lower triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<double> LowerTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -99,7 +99,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <param name="matrix">The matrix upon which the preconditioner is based. </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<double> matrix)
         {
             if (matrix == null)
             {
@@ -161,7 +161,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<double> Approximate(Vector<double> rhs)
         {
             if (rhs == null)
             {
@@ -178,7 +178,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -188,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<double> rhs, Vector<double> lhs)
         {
             if (rhs == null)
             {
@@ -218,7 +218,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
             //     z_i = l_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
             // }
             // NOTE: l_ii should be 1 because u_ii has to be the value
-            Vector rowValues = new DenseVector(_decompositionLU.RowCount);
+            var rowValues = new DenseVector(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
             {
                 // Clear the rowValues 

src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/UnitPreconditioner.cs

@@ -54,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null"/>. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<double> matrix)
         {
             if (matrix == null)
             {
@@ -87,7 +87,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         ///     If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         ///   </para>
         /// </exception>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<double> rhs, Vector<double> lhs)
         {
             if (rhs == null)
             {
@@ -116,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
         /// <exception cref="ArgumentException">
         /// If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         /// </exception>
-        public Vector Approximate(Vector rhs)
+        public Vector<double> Approximate(Vector<double> rhs)
         {
             if (rhs == null)
             {
@@ -128,7 +128,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }

src/Numerics/LinearAlgebra/Double/Solvers/StopCriterium/DivergenceStopCriterium.cs

@@ -217,7 +217,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<double> solutionVector, Vector<double> sourceVector, Vector<double> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Double/Solvers/StopCriterium/FailureStopCriterium.cs

@@ -75,7 +75,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<double> solutionVector, Vector<double> sourceVector, Vector<double> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Double/Solvers/StopCriterium/IIterationStopCriterium.cs

@@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<double> solutionVector, Vector<double> sourceVector, Vector<double> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Double/Solvers/StopCriterium/IterationCountStopCriterium.cs

@@ -130,7 +130,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<double> solutionVector, Vector<double> sourceVector, Vector<double> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Double/Solvers/StopCriterium/ResidualStopCriterium.cs

@@ -222,7 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<double> solutionVector, Vector<double> sourceVector, Vector<double> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -321,33 +321,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfDiagonalInit(rows, columns, init));
         }
 
-        /// <summary>
-        /// Creates a <c>SparseMatrix</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 SparseMatrix(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>
         /// Returns a new matrix containing the lower triangle of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Double/SparseVector.cs

@@ -128,39 +128,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new SparseVector(SparseVectorStorage<double>.OfInit(length, init));
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<double> CreateMatrix(int rows, int columns)
-        {
-            return new SparseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<double> CreateVector(int size)
-        {
-            return new SparseVector(size);
-        }
-
         /// <summary>
         /// Adds a scalar to each element of the vector and stores the result in the result vector.
         /// Warning, the new 'sparse vector' with a non-zero scalar added to it will be a 100% filled

src/Numerics/LinearAlgebra/IBuilder.cs

@@ -0,0 +1,64 @@
+// <copyright file="IBuilder.cs" 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-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+namespace MathNet.Numerics.LinearAlgebra
+{
+    public interface IBuilder<T> where T : struct, IEquatable<T>, IFormattable
+    {
+        /// <summary>
+        /// Creates a <strong>Matrix</strong> 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 <strong>Matrix</strong> with the given dimensions.
+        /// </returns>
+        /// <remarks>
+        /// Creates a matrix of the same matrix type as the current matrix.
+        /// </remarks>
+        Matrix<T> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false);
+
+        /// <summary>
+        /// Creates a Vector 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 Vector with the given dimension.
+        /// </returns>
+        /// <remarks>
+        /// Creates a vector of the same type as the current matrix.
+        /// </remarks>
+        Vector<T> CreateVector(int size, bool fullyMutable = false);
+    }
+}

src/Numerics/LinearAlgebra/Matrix.cs

@@ -59,6 +59,8 @@ namespace MathNet.Numerics.LinearAlgebra
             ColumnCount = storage.ColumnCount;
         }
 
+        static readonly IGenericBuilder<T> Builder = Builder<T>.Instance;
+
         /// <summary>
         /// Gets the raw matrix data storage.
         /// </summary>
@@ -235,31 +237,29 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Creates a <strong>Matrix</strong> for the given number of rows and columns.
+        /// Create a matrix of the same kind 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 <strong>Matrix</strong> with the given dimensions.
-        /// </returns>
-        /// <remarks>
-        /// Creates a matrix of the same matrix type as the current matrix.
-        /// </remarks>
-        public abstract Matrix<T> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false);
+        /// <param name="rows">The number of rows.</param>
+        /// <param name="columns">The number of columns.</param>
+        /// <remarks>Creates a matrix of the same matrix type as the current matrix.</remarks>
+        public Matrix<T> CreateMatrix(int rows, int columns)
+        {
+            return Storage.IsDense
+                ? Builder.DenseMatrix(rows, columns)
+                : Builder.SparseMatrix(rows, columns);
+        }
 
         /// <summary>
-        /// Creates a Vector with a the given dimension.
+        /// Create a vector of the same kind with the given size.
         /// </summary>
         /// <param name="size">The size of the vector.</param>
-        /// <param name="fullyMutable">True if all fields must be mutable.</param>
-        /// <returns>
-        /// A Vector with the given dimension.
-        /// </returns>
-        /// <remarks>
-        /// Creates a vector of the same type as the current matrix.
-        /// </remarks>
-        public abstract Vector<T> CreateVector(int size, bool fullyMutable = false);
+        /// <remarks>Creates a vector of the same type as the current matrix.</remarks>
+        public Vector<T> CreateVector(int size)
+        {
+            return Storage.IsDense
+                ? Builder.DenseVector(size)
+                : Builder.SparseVector(size);
+        }
 
         /// <summary>
         /// Copies a row into an Vector.
@@ -1060,7 +1060,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension);
             }
 
-            var result = CreateMatrix(RowCount, ColumnCount + right.ColumnCount, fullyMutable: true);
+            var result = CreateMatrix(RowCount, ColumnCount + right.ColumnCount);
             Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, skipClearing: true);
             right.Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, right.RowCount, 0, ColumnCount, right.ColumnCount, skipClearing: true);
             return result;
@@ -1116,7 +1116,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension, "lower");
             }
 
-            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount, fullyMutable: true);
+            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount);
             Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, skipClearing: true);
             lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, 0, lower.ColumnCount, skipClearing: true);
             return result;
@@ -1170,7 +1170,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentNullException("lower");
             }
 
-            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount, fullyMutable: true);
+            var result = CreateMatrix(RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
             Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount);
             lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, ColumnCount, lower.ColumnCount);
             return result;

src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@@ -349,33 +349,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             get { return _values; }
         }
 
-        /// <summary>
-        /// Creates a <c>DenseMatrix</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>DenseMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return new DenseMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<float> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Returns the transpose of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Single/DenseVector.cs

@@ -191,39 +191,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new DenseVector(array);
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<float> CreateMatrix(int rows, int columns)
-        {
-            return new DenseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<float> CreateVector(int size)
-        {
-            return new DenseVector(size);
-        }
-
         /// <summary>
         /// Adds a scalar to each element of the vector and stores the result in the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@@ -180,35 +180,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 i => (float) distribution.Sample()));
         }
 
-        /// <summary>
-        /// Creates a <c>DiagonalMatrix</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>DiagonalMatrix</c> with the given dimensions.
-        /// </returns>
-        public override Matrix<float> CreateMatrix(int numberOfRows, int numberOfColumns, bool fullyMutable = false)
-        {
-            return fullyMutable
-                ? (Matrix<float>) new SparseMatrix(numberOfRows, numberOfColumns)
-                : new DiagonalMatrix(numberOfRows, numberOfColumns);
-        }
-
-        /// <summary>
-        /// Creates a <see cref="Vector{T}"/> 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="Vector{T}"/> with the given dimension.
-        /// </returns>
-        public override Vector<float> CreateVector(int size, bool fullyMutable = false)
-        {
-            return new SparseVector(size);
-        }
-
         #region Elementary operations
 
         /// <summary>

src/Numerics/LinearAlgebra/Single/Solvers/IIterativeSolver.cs

@@ -64,7 +64,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        Vector Solve(Matrix matrix, Vector vector);
+        Vector<float> Solve(Matrix<float> matrix, Vector<float> vector);
 
         /// <summary>
         /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the
@@ -73,7 +73,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        void Solve(Matrix matrix, Vector input, Vector result);
+        void Solve(Matrix<float> matrix, Vector<float> input, Vector<float> result);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -82,7 +82,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        Matrix Solve(Matrix matrix, Matrix input);
+        Matrix<float> Solve(Matrix<float> matrix, Matrix<float> input);
 
         /// <summary>
         /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
@@ -91,6 +91,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        void Solve(Matrix matrix, Matrix input, Matrix result);
+        void Solve(Matrix<float> matrix, Matrix<float> input, Matrix<float> result);
     }
 }

src/Numerics/LinearAlgebra/Single/Solvers/IIterator.cs

@@ -80,7 +80,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<float> solutionVector, Vector<float> sourceVector, Vector<float> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Single/Solvers/Iterative/BiCgStab.cs

@@ -204,14 +204,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="vector">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <returns>The result <see cref="Vector"/>, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<float> Solve(Matrix<float> matrix, Vector<float> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -223,7 +223,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Vector"/>, <c>b</c>.</param>
         /// <param name="result">The result <see cref="Vector"/>, <c>x</c>.</param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<float> matrix, Vector<float> input, Vector<float> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -278,7 +278,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
             // Compute r_0 = b - Ax_0 for some initial guess x_0
             // In this case we take x_0 = vector
             // This is basically a SAXPY so it could be made a lot faster
-            Vector residuals = new DenseVector(matrix.RowCount);
+            var residuals = new DenseVector(matrix.RowCount);
             CalculateTrueResidual(matrix, residuals, result, input);
 
             // Choose r~ (for example, r~ = r_0)
@@ -287,13 +287,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
             // create seven temporary vectors needed to hold temporary
             // coefficients. All vectors are mangled in each iteration.
             // These are defined here to prevent stressing the garbage collector
-            Vector vecP = new DenseVector(residuals.Count);
-            Vector vecPdash = new DenseVector(residuals.Count);
-            Vector nu = new DenseVector(residuals.Count);
-            Vector vecS = new DenseVector(residuals.Count);
-            Vector vecSdash = new DenseVector(residuals.Count);
-            Vector temp = new DenseVector(residuals.Count);
-            Vector temp2 = new DenseVector(residuals.Count);
+            var vecP = new DenseVector(residuals.Count);
+            var vecPdash = new DenseVector(residuals.Count);
+            var nu = new DenseVector(residuals.Count);
+            var vecS = new DenseVector(residuals.Count);
+            var vecSdash = new DenseVector(residuals.Count);
+            var temp = new DenseVector(residuals.Count);
+            var temp2 = new DenseVector(residuals.Count);
 
             // create some temporary float variables that are needed
             // to hold values in between iterations
@@ -431,7 +431,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        private static void CalculateTrueResidual(Matrix<float> matrix, Vector<float> residual, Vector<float> x, Vector<float> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -451,7 +451,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        private bool ShouldContinue(int iterationNumber, Vector<float> result, Vector<float> source, Vector<float> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -475,7 +475,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <returns>The result <see cref="Matrix"/>, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<float> Solve(Matrix<float> matrix, Matrix<float> input)
         {
             if (matrix == null)
             {
@@ -487,7 +487,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -499,7 +499,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
         /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
         /// <param name="result">The result <see cref="Matrix"/>, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<float> matrix, Matrix<float> input, Matrix<float> result)
         {
             if (matrix == null)
             {
@@ -523,7 +523,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Single/Solvers/Iterative/CompositeSolver.cs

@@ -413,14 +413,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<float> Solve(Matrix<float> matrix, Vector<float> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -432,7 +432,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<float> matrix, Vector<float> input, Vector<float> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -482,8 +482,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
 
             // Create a copy of the solution and result vectors so we can use them
             // later on
-            var internalInput = (Vector)input.Clone();
-            var internalResult = (Vector)result.Clone();
+            var internalInput = input.Clone();
+            var internalResult = result.Clone();
 
             foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped))
             {
@@ -574,7 +574,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<float> Solve(Matrix<float> matrix, Matrix<float> input)
         {
             if (matrix == null)
             {
@@ -586,7 +586,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -598,7 +598,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<float> matrix, Matrix<float> input, Matrix<float> result)
         {
             if (matrix == null)
             {
@@ -622,7 +622,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Single/Solvers/Iterative/GpBiCg.cs

@@ -253,14 +253,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<float> Solve(Matrix<float> matrix, Vector<float> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -272,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<float> matrix, Vector<float> input, Vector<float> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -521,7 +521,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<float> matrix, Vector<float> residual, Vector<float> x, Vector<float> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -539,7 +539,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<float> result, Vector<float> source, Vector<float> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -580,7 +580,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<float> Solve(Matrix<float> matrix, Matrix<float> input)
         {
             if (matrix == null)
             {
@@ -592,7 +592,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -604,7 +604,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<float> matrix, Matrix<float> input, Matrix<float> result)
         {
             if (matrix == null)
             {

src/Numerics/LinearAlgebra/Single/Solvers/Iterative/MlkBiCgStab.cs

@@ -88,7 +88,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <summary>
         /// The collection of starting vectors which are used as the basis for the Krylov sub-space.
         /// </summary>
-        private IList<Vector> _startingVectors;
+        private IList<Vector<float>> _startingVectors;
 
         /// <summary>
         /// The number of starting vectors used by the algorithm
@@ -226,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// Gets or sets a series of orthonormal vectors which will be used as basis for the 
         /// Krylov sub-space.
         /// </summary>
-        public IList<Vector> StartingVectors
+        public IList<Vector<float>> StartingVectors
         {
             [DebuggerStepThrough]
             get
@@ -278,14 +278,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<float> Solve(Matrix<float> matrix, Vector<float> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -297,7 +297,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<float> matrix, Vector<float> input, Vector<float> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -641,7 +641,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         ///  the <paramref name="numberOfVariables"/> is smaller than 
         ///  the <paramref name="maximumNumberOfStartingVectors"/>.
         /// </returns>
-        private static IList<Vector> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
+        private static IList<Vector<float>> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables)
         {
             // Create no more starting vectors than the size of the problem - 1
             // Get random values and then orthogonalize them with
@@ -652,7 +652,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
             // mean = 0 and sd = 1
             var distribution = new Normal();
 
-            Matrix matrix = new DenseMatrix(numberOfVariables, count);
+            var matrix = new DenseMatrix(numberOfVariables, count);
             for (var i = 0; i < matrix.ColumnCount; i++)
             {
                 var samples = new float[matrix.RowCount];
@@ -670,10 +670,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
             var orthogonalMatrix = gs.Q;
 
             // Now transfer this to vectors
-            var result = new List<Vector>();
+            var result = new List<Vector<float>>();
             for (var i = 0; i < orthogonalMatrix.ColumnCount; i++)
             {
-                result.Add((Vector)orthogonalMatrix.Column(i));
+                result.Add(orthogonalMatrix.Column(i));
                 
                 // Normalize the result vector
                 result[i].Multiply(1 / result[i].L2Norm(), result[i]);
@@ -688,9 +688,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="arraySize">Number of vectors</param>
         /// <param name="vectorSize">Size of each vector</param>
         /// <returns>Array of random vectors</returns>
-        private static Vector[] CreateVectorArray(int arraySize, int vectorSize)
+        private static Vector<float>[] CreateVectorArray(int arraySize, int vectorSize)
         {
-            var result = new Vector[arraySize];
+            var result = new Vector<float>[arraySize];
             for (var i = 0; i < result.Length; i++)
             {
                 result[i] = new DenseVector(vectorSize);
@@ -706,7 +706,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="residual">Residual <see cref="Vector"/> data.</param>
         /// <param name="x">x <see cref="Vector"/> data.</param>
         /// <param name="b">b <see cref="Vector"/> data.</param>
-        private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        private static void CalculateTrueResidual(Matrix<float> matrix, Vector<float> residual, Vector<float> x, Vector<float> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -724,7 +724,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        private bool ShouldContinue(int iterationNumber, Vector<float> result, Vector<float> source, Vector<float> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -748,7 +748,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<float> Solve(Matrix<float> matrix, Matrix<float> input)
         {
             if (matrix == null)
             {
@@ -760,7 +760,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -772,7 +772,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<float> matrix, Matrix<float> input, Matrix<float> result)
         {
             if (matrix == null)
             {
@@ -796,7 +796,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector)input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs

@@ -192,14 +192,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="vector">The solution vector, <c>b</c>.</param>
         /// <returns>The result vector, <c>x</c>.</returns>
-        public Vector Solve(Matrix matrix, Vector vector)
+        public Vector<float> Solve(Matrix<float> matrix, Vector<float> vector)
         {
             if (vector == null)
             {
                 throw new ArgumentNullException();
             }
 
-            Vector result = new DenseVector(matrix.RowCount);
+            var result = new DenseVector(matrix.RowCount);
             Solve(matrix, vector, result);
             return result;
         }
@@ -211,7 +211,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution vector, <c>b</c></param>
         /// <param name="result">The result vector, <c>x</c></param>
-        public void Solve(Matrix matrix, Vector input, Vector result)
+        public void Solve(Matrix<float> matrix, Vector<float> input, Vector<float> result)
         {
             // If we were stopped before, we are no longer
             // We're doing this at the start of the method to ensure
@@ -427,7 +427,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="residual">Residual values in <see cref="Vector"/>.</param>
         /// <param name="x">Instance of the <see cref="Vector"/> x.</param>
         /// <param name="b">Instance of the <see cref="Vector"/> b.</param>
-        static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b)
+        static void CalculateTrueResidual(Matrix<float> matrix, Vector<float> residual, Vector<float> x, Vector<float> b)
         {
             // -Ax = residual
             matrix.Multiply(x, residual);
@@ -445,7 +445,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="source">Source <see cref="Vector"/>.</param>
         /// <param name="residuals">Residual <see cref="Vector"/>.</param>
         /// <returns><c>true</c> if continue, otherwise <c>false</c></returns>
-        bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals)
+        bool ShouldContinue(int iterationNumber, Vector<float> result, Vector<float> source, Vector<float> residuals)
         {
             if (_hasBeenStopped)
             {
@@ -479,7 +479,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <returns>The result matrix, <c>X</c>.</returns>
-        public Matrix Solve(Matrix matrix, Matrix input)
+        public Matrix<float> Solve(Matrix<float> matrix, Matrix<float> input)
         {
             if (matrix == null)
             {
@@ -491,7 +491,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
                 throw new ArgumentNullException("input");
             }
 
-            var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount);
+            var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount);
             Solve(matrix, input, result);
             return result;
         }
@@ -503,7 +503,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
         /// <param name="matrix">The coefficient matrix, <c>A</c>.</param>
         /// <param name="input">The solution matrix, <c>B</c>.</param>
         /// <param name="result">The result matrix, <c>X</c></param>
-        public void Solve(Matrix matrix, Matrix input, Matrix result)
+        public void Solve(Matrix<float> matrix, Matrix<float> input, Matrix<float> result)
         {
             if (matrix == null)
             {
@@ -527,7 +527,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative
 
             for (var column = 0; column < input.ColumnCount; column++)
             {
-                var solution = Solve(matrix, (Vector) input.Column(column));
+                var solution = Solve(matrix, input.Column(column));
                 foreach (var element in solution.EnumerateNonZeroIndexed())
                 {
                     result.At(element.Item1, column, element.Item2);

src/Numerics/LinearAlgebra/Single/Solvers/Iterator.cs

@@ -212,7 +212,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<float> solutionVector, Vector<float> sourceVector, Vector<float> residualVector)
         {
             if (_stopCriterias.Count == 0)
             {

src/Numerics/LinearAlgebra/Single/Solvers/Preconditioners/Diagonal.cs

@@ -66,7 +66,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// The <see cref="Matrix"/> upon which this preconditioner is based.</param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<float> matrix)
         {
             if (matrix == null)
             {
@@ -90,7 +90,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<float> Approximate(Vector<float> rhs)
         {
             if (rhs == null)
             {
@@ -107,7 +107,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -117,7 +117,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<float> rhs, Vector<float> lhs)
         {
             if (rhs == null)
             {

src/Numerics/LinearAlgebra/Single/Solvers/Preconditioners/IPreConditioner.cs

@@ -54,20 +54,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// Initializes the preconditioner and loads the internal data structures.
         /// </summary>
         /// <param name="matrix">The matrix on which the preconditioner is based.</param>
-        void Initialize(Matrix matrix);
+        void Initialize(Matrix<float> matrix);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        Vector Approximate(Vector rhs);
+        Vector<float> Approximate(Vector<float> rhs);
 
         /// <summary>
         /// Approximates the solution to the matrix equation <b>Mx = b</b>.
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        void Approximate(Vector rhs, Vector lhs);
+        void Approximate(Vector<float> rhs, Vector<float> lhs);
     }
 }

src/Numerics/LinearAlgebra/Single/Solvers/Preconditioners/Ilutp.cs

@@ -255,9 +255,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<float> UpperTriangle()
         {
-            return (Matrix)_upper.Clone();
+            return _upper.Clone();
         }
 
         /// <summary>
@@ -267,9 +267,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// This method is used for debugging purposes only and should normally not be used.
         /// </remarks>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<float> LowerTriangle()
         {
-            return (Matrix)_lower.Clone();
+            return _lower.Clone();
         }
 
         /// <summary>
@@ -301,7 +301,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException"> If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<float> matrix)
         {
             if (matrix == null)
             {
@@ -373,8 +373,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
                 _pivots[i] = i;
             }
 
-            Vector workVector = new DenseVector(sparseMatrix.RowCount);
-            Vector rowVector = new DenseVector(sparseMatrix.ColumnCount);
+            var workVector = new DenseVector(sparseMatrix.RowCount);
+            var rowVector = new DenseVector(sparseMatrix.ColumnCount);
             var indexSorting = new int[sparseMatrix.RowCount];
             
             // spaceLeft = lfilNnz * nnz(A)
@@ -530,7 +530,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// Pivot elements in the <paramref name="row"/> according to internal pivot array
         /// </summary>
         /// <param name="row">Row <see cref="Vector"/> to pivot in</param>
-        private void PivotRow(Vector row)
+        private void PivotRow(Vector<float> row)
         {
             var knownPivots = new Dictionary<int, int>();
             
@@ -582,7 +582,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <param name="matrix">Source <see cref="Matrix"/>.</param>
         /// <param name="firstColumn">First column index to swap</param>
         /// <param name="secondColumn">Second column index to swap</param>
-        private static void SwapColumns(Matrix matrix, int firstColumn, int secondColumn)
+        private static void SwapColumns(Matrix<float> matrix, int firstColumn, int secondColumn)
         {
             for (var i = 0; i < matrix.RowCount; i++)
             {
@@ -599,7 +599,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <param name="upperBound">Sort till upper bound</param>
         /// <param name="sortedIndices">Array with sorted vector indicies</param>
         /// <param name="values">Source <see cref="Vector"/></param>
-        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        private static void FindLargestItems(int lowerBound, int upperBound, int[] sortedIndices, Vector<float> values)
         {
             // Copy the indices for the values into the array
             for (var i = 0; i < upperBound + 1 - lowerBound; i++)
@@ -624,7 +624,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<float> Approximate(Vector<float> rhs)
         {
             if (rhs == null)
             {
@@ -641,7 +641,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -651,7 +651,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<float> rhs, Vector<float> lhs)
         {
             if (rhs == null)
             {
@@ -676,7 +676,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
             // Solve equation here
             // Pivot(vector, result);
             // Solve L*Y = B(piv,:)
-            Vector rowValues = new DenseVector(_lower.RowCount);
+            var rowValues = new DenseVector(_lower.RowCount);
             for (var i = 0; i < _lower.RowCount; i++)
             {
                 _lower.Row(i, rowValues);
@@ -706,7 +706,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
 
             // We have a column pivot so we only need to pivot the
             // end result not the incoming right hand side vector
-            var temp = (Vector)lhs.Clone();
+            var temp = lhs.Clone();
             
             Pivot(temp, lhs);
         }
@@ -716,7 +716,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </summary>
         /// <param name="vector">Source <see cref="Vector"/>.</param>
         /// <param name="result">Result <see cref="Vector"/> after pivoting.</param>
-        private void Pivot(Vector vector, Vector result)
+        private void Pivot(Vector<float> vector, Vector<float> result)
         {
             for (var i = 0; i < _pivots.Length; i++)
             {

src/Numerics/LinearAlgebra/Single/Solvers/Preconditioners/IlutpElementSorter.cs

@@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <param name="upperBound">The stopping index.</param>
         /// <param name="sortedIndices">An array that will contain the sorted indices once the algorithm finishes.</param>
         /// <param name="values">The <see cref="Vector"/> that contains the values that need to be sorted.</param>
-        public static void SortDoubleIndicesDecreasing(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        public static void SortDoubleIndicesDecreasing(int lowerBound, int upperBound, int[] sortedIndices, Vector<float> values)
         {
             // Move all the indices that we're interested in to the beginning of the
             // array. Ignore the rest of the indices.
@@ -73,7 +73,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <param name="upperBound">The stopping index.</param>
         /// <param name="sortedIndices">An array that will contain the sorted indices once the algorithm finishes.</param>
         /// <param name="values">The <see cref="Vector"/> that contains the values that need to be sorted.</param>
-        private static void HeapSortDoublesIndices(int lowerBound, int upperBound, int[] sortedIndices, Vector values)
+        private static void HeapSortDoublesIndices(int lowerBound, int upperBound, int[] sortedIndices, Vector<float> values)
         {
             var start = ((upperBound - lowerBound + 1) / 2) - 1 + lowerBound;
             var end = (upperBound - lowerBound + 1) - 1 + lowerBound;
@@ -95,7 +95,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <param name="count">Length of <paramref name="values"/></param>
         /// <param name="sortedIndices">Indicies of <paramref name="values"/></param>
         /// <param name="values">Target <see cref="Vector"/></param>
-        private static void BuildDoubleIndexHeap(int start, int count, int[] sortedIndices, Vector values)
+        private static void BuildDoubleIndexHeap(int start, int count, int[] sortedIndices, Vector<float> values)
         {
             while (start >= 0)
             {
@@ -111,7 +111,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <param name="values">Target <see cref="Vector"/></param>
         /// <param name="begin">Root position</param>
         /// <param name="count">Length of <paramref name="values"/></param>
-        private static void SiftDoubleIndices(int[] sortedIndices, Vector values, int begin, int count)
+        private static void SiftDoubleIndices(int[] sortedIndices, Vector<float> values, int begin, int count)
         {
             var root = begin;
 

src/Numerics/LinearAlgebra/Single/Solvers/Preconditioners/IncompleteLU.cs

@@ -54,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// Returns the upper triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the upper triagonal elements.</returns>
-        internal Matrix UpperTriangle()
+        internal Matrix<float> UpperTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -72,7 +72,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// Returns the lower triagonal matrix that was created during the LU decomposition.
         /// </summary>
         /// <returns>A new matrix containing the lower triagonal elements.</returns>
-        internal Matrix LowerTriangle()
+        internal Matrix<float> LowerTriangle()
         {
             var result = new SparseMatrix(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
@@ -99,7 +99,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <param name="matrix">The matrix upon which the preconditioner is based. </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null" />.</exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<float> matrix)
         {
             if (matrix == null)
             {
@@ -161,7 +161,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <returns>The left hand side vector.</returns>
-        public Vector Approximate(Vector rhs)
+        public Vector<float> Approximate(Vector<float> rhs)
         {
             if (rhs == null)
             {
@@ -178,7 +178,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "rhs");
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }
@@ -188,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </summary>
         /// <param name="rhs">The right hand side vector.</param>
         /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<float> rhs, Vector<float> lhs)
         {
             if (rhs == null)
             {
@@ -218,7 +218,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
             //     z_i = l_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
             // }
             // NOTE: l_ii should be 1 because u_ii has to be the value
-            Vector rowValues = new DenseVector(_decompositionLU.RowCount);
+            var rowValues = new DenseVector(_decompositionLU.RowCount);
             for (var i = 0; i < _decompositionLU.RowCount; i++)
             {
                 // Clear the rowValues 

src/Numerics/LinearAlgebra/Single/Solvers/Preconditioners/UnitPreconditioner.cs

@@ -54,7 +54,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// </param>
         /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <see langword="null"/>. </exception>
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
-        public void Initialize(Matrix matrix)
+        public void Initialize(Matrix<float> matrix)
         {
             if (matrix == null)
             {
@@ -87,7 +87,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         ///     If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         ///   </para>
         /// </exception>
-        public void Approximate(Vector rhs, Vector lhs)
+        public void Approximate(Vector<float> rhs, Vector<float> lhs)
         {
             if (rhs == null)
             {
@@ -116,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
         /// <exception cref="ArgumentException">
         /// If the size of <paramref name="rhs"/> is different the number of rows of the coefficient matrix.
         /// </exception>
-        public Vector Approximate(Vector rhs)
+        public Vector<float> Approximate(Vector<float> rhs)
         {
             if (rhs == null)
             {
@@ -128,7 +128,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Preconditioners
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
             }
 
-            Vector result = new DenseVector(rhs.Count);
+            var result = new DenseVector(rhs.Count);
             Approximate(rhs, result);
             return result;
         }

src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/DivergenceStopCriterium.cs

@@ -217,7 +217,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<float> solutionVector, Vector<float> sourceVector, Vector<float> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/FailureStopCriterium.cs

@@ -75,7 +75,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<float> solutionVector, Vector<float> sourceVector, Vector<float> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/IIterationStopCriterium.cs

@@ -47,7 +47,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector);
+        void DetermineStatus(int iterationNumber, Vector<float> solutionVector, Vector<float> sourceVector, Vector<float> residualVector);
 
         /// <summary>
         /// Gets the current calculation status.

src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/IterationCountStopCriterium.cs

@@ -130,7 +130,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<float> solutionVector, Vector<float> sourceVector, Vector<float> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Single/Solvers/StopCriterium/ResidualStopCriterium.cs

@@ -222,7 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.StopCriterium
         /// on the invocation of this method. Therefore this method should only be called if the 
         /// calculation has moved forwards at least one step.
         /// </remarks>
-        public void DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
+        public void DetermineStatus(int iterationNumber, Vector<float> solutionVector, Vector<float> sourceVector, Vector<float> residualVector)
         {
             if (iterationNumber < 0)
             {

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -321,33 +321,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfDiagonalInit(rows, columns, init));
         }
 
-        /// <summary>
-        /// Creates a <c>SparseMatrix</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 SparseMatrix(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>
         /// Returns a new matrix containing the lower triangle of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Single/SparseVector.cs

@@ -128,39 +128,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new SparseVector(SparseVectorStorage<float>.OfInit(length, init));
         }
 
-        /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="rows">
-        /// The number of rows.
-        /// </param>
-        /// <param name="columns">
-        /// The number of columns.
-        /// </param>
-        /// <returns>
-        /// A matrix with the given dimensions.
-        /// </returns>
-        public override Matrix<float> CreateMatrix(int rows, int columns)
-        {
-            return new SparseMatrix(rows, columns);
-        }
-
-        /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
-        /// </summary>
-        /// <param name="size">
-        /// The size of the <strong>Vector</strong> to create.
-        /// </param>
-        /// <returns>
-        /// The new <c>Vector</c>.
-        /// </returns>
-        public override Vector<float> CreateVector(int size)
-        {
-            return new SparseVector(size);
-        }
-
         /// <summary>
         /// Adds a scalar to each element of the vector and stores the result in the result vector.
         /// Warning, the new 'sparse vector' with a non-zero scalar added to it will be a 100% filled

src/Numerics/LinearAlgebra/Vector.cs

@@ -58,6 +58,8 @@ namespace MathNet.Numerics.LinearAlgebra
             Count = storage.Length;
         }
 
+        static readonly IGenericBuilder<T> Builder = Builder<T>.Instance;
+
         /// <summary>
         /// Gets the raw vector data storage.
         /// </summary>
@@ -131,21 +133,29 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Creates a matrix with the given dimensions using the same storage type
-        /// as this vector.
+        /// Create a matrix of the same kind with the provided number of rows and columns.
         /// </summary>
         /// <param name="rows">The number of rows.</param>
         /// <param name="columns">The number of columns.</param>
-        /// <returns>A matrix with the given dimensions.</returns>
-        public abstract Matrix<T> CreateMatrix(int rows, int columns);
+        /// <remarks>Creates a matrix of the same matrix type as the current matrix.</remarks>
+        public Matrix<T> CreateMatrix(int rows, int columns)
+        {
+            return Storage.IsDense
+                ? Builder.DenseMatrix(rows, columns)
+                : Builder.SparseMatrix(rows, columns);
+        }
 
         /// <summary>
-        /// Creates a <strong>Vector</strong> of the given size using the same storage type
-        /// as this vector.
+        /// Create a vector of the same kind with the provided dimension.
         /// </summary>
-        /// <param name="size">The size of the <strong>Vector</strong> to create.</param>
-        /// <returns>The new <c>Vector</c>.</returns>
-        public abstract Vector<T> CreateVector(int size);
+        /// <param name="size">The size of the vector.</param>
+        /// <remarks>Creates a vector of the same type as the current matrix.</remarks>
+        public Vector<T> CreateVector(int size)
+        {
+            return Storage.IsDense
+                ? Builder.DenseVector(size)
+                : Builder.SparseVector(size);
+        }
 
         /// <summary>
         /// Returns a deep-copy clone of the vector.