LA: more sensible choice of return storage type, with tests #146

src/Numerics/LinearAlgebra/Builder.cs

@@ -395,7 +395,7 @@ namespace MathNet.Numerics.LinearAlgebra
         {
             var storage = example.Storage;
             if (storage is DenseColumnMajorMatrixStorage<T>) return Dense(rows, columns);
-            if (!fullyMutable && storage is DiagonalMatrixStorage<T>) return Diagonal(rows, columns);
+            if (storage is DiagonalMatrixStorage<T>) return fullyMutable ? Sparse(rows, columns) : Diagonal(rows, columns);
             if (storage is SparseCompressedRowMatrixStorage<T>) return Sparse(rows, columns);
             return Dense(rows, columns);
         }
@@ -424,7 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra
             var storage1 = example.Storage;
             var storage2 = otherExample.Storage;
             if (storage1 is DenseColumnMajorMatrixStorage<T> || storage2 is DenseColumnMajorMatrixStorage<T>) return Dense(rows, columns);
-            if (!fullyMutable && storage1 is DiagonalMatrixStorage<T> && storage2 is DiagonalMatrixStorage<T>) return Diagonal(rows, columns);
+            if (storage1 is DiagonalMatrixStorage<T> && storage2 is DiagonalMatrixStorage<T>) return fullyMutable ? Sparse(rows, columns) : Diagonal(rows, columns);
             if (storage1 is SparseCompressedRowMatrixStorage<T> || storage2 is SparseCompressedRowMatrixStorage<T>) return Sparse(rows, columns);
             return Dense(rows, columns);
         }
@@ -1331,6 +1331,14 @@ namespace MathNet.Numerics.LinearAlgebra
             return example.Storage.IsDense || otherExample.Storage.IsDense ? Dense(example.Count) : Sparse(example.Count);
         }
 
+        /// <summary>
+        /// Create a new vector with a type that can represent and is closest to both provided samples.
+        /// </summary>
+        public Vector<T> SameAs(Matrix<T> matrix, Vector<T> vector, int length)
+        {
+            return matrix.Storage.IsDense || vector.Storage.IsDense ? Dense(length) : Sparse(length);
+        }
+
         /// <summary>
         /// Create a new dense vector with values sampled from the provided random distribution.
         /// </summary>

src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs

@@ -531,7 +531,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw DimensionsDontMatch<ArgumentException>(this, rightSide, "rightSide");
             }
 
-            var ret = Vector<T>.Build.SameAs(this, RowCount);
+            var ret = Vector<T>.Build.SameAs(this, rightSide, RowCount);
             DoMultiply(rightSide, ret);
             return ret;
         }
@@ -580,7 +580,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw DimensionsDontMatch<ArgumentException>(this, leftSide, "leftSide");
             }
 
-            var ret = Vector<T>.Build.SameAs(this, ColumnCount);
+            var ret = Vector<T>.Build.SameAs(this, leftSide, ColumnCount);
             DoLeftMultiply(leftSide, ret);
             return ret;
         }
@@ -727,7 +727,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw DimensionsDontMatch<ArgumentException>(this, rightSide, "rightSide");
             }
 
-            var result = Vector<T>.Build.SameAs(this, ColumnCount);
+            var result = Vector<T>.Build.SameAs(this, rightSide, ColumnCount);
             DoTransposeThisAndMultiply(rightSide, result);
             return result;
         }
@@ -926,7 +926,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw DimensionsDontMatch<ArgumentException>(this, other, "other");
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, other);
             DoPointwiseMultiply(other, result);
             return result;
         }
@@ -961,7 +961,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw DimensionsDontMatch<ArgumentException>(this, divisor);
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, divisor);
             DoPointwiseDivide(divisor, result);
             return result;
         }
@@ -996,7 +996,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw DimensionsDontMatch<ArgumentException>(this, divisor);
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, divisor);
             DoPointwiseModulus(divisor, result);
             return result;
         }
@@ -1074,7 +1074,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <returns>The kronecker product of the two matrices.</returns>
         public Matrix<T> KroneckerProduct(Matrix<T> other)
         {
-            var result = Build.SameAs(this, RowCount*other.RowCount, ColumnCount*other.ColumnCount);
+            var result = Build.SameAs(this, other, RowCount*other.RowCount, ColumnCount*other.ColumnCount);
             KroneckerProduct(other, result);
             return result;
         }

src/Numerics/LinearAlgebra/Matrix.cs

@@ -1176,7 +1176,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentNullException("lower");
             }
 
-            var result = Build.SameAs(this, lower, RowCount + lower.RowCount, ColumnCount + lower.ColumnCount);
+            var result = Build.SameAs(this, lower, RowCount + lower.RowCount, ColumnCount + lower.ColumnCount, RowCount != 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/Vector.Arithmetic.cs

@@ -218,7 +218,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, other);
             DoAdd(other, result);
             return result;
         }
@@ -346,7 +346,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, other);
             DoSubtract(other, result);
             return result;
         }
@@ -606,7 +606,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, other);
             DoPointwiseMultiply(other, result);
             return result;
         }
@@ -646,7 +646,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "divisor");
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, divisor);
             DoPointwiseDivide(divisor, result);
             return result;
         }
@@ -686,7 +686,7 @@ namespace MathNet.Numerics.LinearAlgebra
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "divisor");
             }
 
-            var result = Build.SameAs(this);
+            var result = Build.SameAs(this, divisor);
             DoPointwiseModulus(divisor, result);
             return result;
         }

src/UnitTests/LinearAlgebraTests/Complex/ReturnTypeTests.cs

@@ -0,0 +1,226 @@
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Complex;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
+{
+#if NOSYSNUMERICS
+    using Complex = Numerics.Complex;
+#else
+    using Complex = System.Numerics.Complex;
+#endif
+
+    [TestFixture]
+    public class ReturnTypeTests
+    {
+        readonly Vector<Complex> _vectorDense = Vector<Complex>.Build.Dense(3);
+        readonly Vector<Complex> _vectorSparse = Vector<Complex>.Build.Sparse(3);
+        readonly Matrix<Complex> _matrixDense = Matrix<Complex>.Build.Dense(3, 3);
+        readonly Matrix<Complex> _matrixSparse = Matrix<Complex>.Build.Sparse(3, 3);
+        readonly Matrix<Complex> _matrixDiagonal = Matrix<Complex>.Build.Diagonal(3, 3);
+
+        [Test]
+        public void VerifyExamples()
+        {
+            Assert.That(_vectorDense, Is.TypeOf<DenseVector>());
+            Assert.That(_vectorSparse, Is.TypeOf<SparseVector>());
+            Assert.That(_matrixDense, Is.TypeOf<DenseMatrix>());
+            Assert.That(_matrixSparse, Is.TypeOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal, Is.TypeOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Negate()
+        {
+            Assert.That(_vectorDense.Negate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Negate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Negate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Negate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Negate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Conjugate()
+        {
+            Assert.That(_vectorDense.Conjugate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Conjugate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Conjugate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Conjugate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Conjugate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Transpose()
+        {
+            Assert.That(_matrixDense.Transpose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Transpose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Transpose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void ConjugateTranspose()
+        {
+            Assert.That(_matrixDense.ConjugateTranspose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.ConjugateTranspose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.ConjugateTranspose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Add()
+        {
+            Assert.That(_vectorDense + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense + _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse + _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Subtract()
+        {
+            Assert.That(_vectorDense - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense - _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse - _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void PointwiseMultiply()
+        {
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixMultiply()
+        {
+            Assert.That(_matrixDense*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse*_matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixVectorMultiply()
+        {
+            Assert.That(_matrixDense*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense*_vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal*_vectorSparse, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_vectorDense*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense*_matrixSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_vectorDense*_matrixDiagonal, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDiagonal, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+        }
+
+        [Test]
+        public void Append()
+        {
+            Assert.That(_matrixDense.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void Stack()
+        {
+            Assert.That(_matrixDense.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void DiagonalStack()
+        {
+            Assert.That(_matrixDense.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            // Special Case
+            Assert.That(Matrix<Complex>.Build.DiagonalIdentity(2, 4).DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+    }
+}

src/UnitTests/LinearAlgebraTests/Complex32/ReturnTypeTests.cs

@@ -0,0 +1,222 @@
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Complex32;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
+{
+    using Numerics;
+
+    [TestFixture]
+    public class ReturnTypeTests
+    {
+        readonly Vector<Complex32> _vectorDense = Vector<Complex32>.Build.Dense(3);
+        readonly Vector<Complex32> _vectorSparse = Vector<Complex32>.Build.Sparse(3);
+        readonly Matrix<Complex32> _matrixDense = Matrix<Complex32>.Build.Dense(3, 3);
+        readonly Matrix<Complex32> _matrixSparse = Matrix<Complex32>.Build.Sparse(3, 3);
+        readonly Matrix<Complex32> _matrixDiagonal = Matrix<Complex32>.Build.Diagonal(3, 3);
+
+        [Test]
+        public void VerifyExamples()
+        {
+            Assert.That(_vectorDense, Is.TypeOf<DenseVector>());
+            Assert.That(_vectorSparse, Is.TypeOf<SparseVector>());
+            Assert.That(_matrixDense, Is.TypeOf<DenseMatrix>());
+            Assert.That(_matrixSparse, Is.TypeOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal, Is.TypeOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Negate()
+        {
+            Assert.That(_vectorDense.Negate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Negate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Negate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Negate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Negate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Conjugate()
+        {
+            Assert.That(_vectorDense.Conjugate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Conjugate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Conjugate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Conjugate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Conjugate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Transpose()
+        {
+            Assert.That(_matrixDense.Transpose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Transpose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Transpose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void ConjugateTranspose()
+        {
+            Assert.That(_matrixDense.ConjugateTranspose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.ConjugateTranspose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.ConjugateTranspose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Add()
+        {
+            Assert.That(_vectorDense + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense + _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse + _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Subtract()
+        {
+            Assert.That(_vectorDense - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense - _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse - _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void PointwiseMultiply()
+        {
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixMultiply()
+        {
+            Assert.That(_matrixDense*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse*_matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixVectorMultiply()
+        {
+            Assert.That(_matrixDense*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense*_vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal*_vectorSparse, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_vectorDense*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense*_matrixSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_vectorDense*_matrixDiagonal, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDiagonal, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+        }
+
+        [Test]
+        public void Append()
+        {
+            Assert.That(_matrixDense.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void Stack()
+        {
+            Assert.That(_matrixDense.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void DiagonalStack()
+        {
+            Assert.That(_matrixDense.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            // Special Case
+            Assert.That(Matrix<Complex32>.Build.DiagonalIdentity(2, 4).DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+    }
+}

src/UnitTests/LinearAlgebraTests/Double/ReturnTypeTests.cs

@@ -0,0 +1,220 @@
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
+{
+    [TestFixture]
+    public class ReturnTypeTests
+    {
+        readonly Vector<double> _vectorDense = Vector<double>.Build.Dense(3);
+        readonly Vector<double> _vectorSparse = Vector<double>.Build.Sparse(3);
+        readonly Matrix<double> _matrixDense = Matrix<double>.Build.Dense(3, 3);
+        readonly Matrix<double> _matrixSparse = Matrix<double>.Build.Sparse(3, 3);
+        readonly Matrix<double> _matrixDiagonal = Matrix<double>.Build.Diagonal(3, 3);
+
+        [Test]
+        public void VerifyExamples()
+        {
+            Assert.That(_vectorDense, Is.TypeOf<DenseVector>());
+            Assert.That(_vectorSparse, Is.TypeOf<SparseVector>());
+            Assert.That(_matrixDense, Is.TypeOf<DenseMatrix>());
+            Assert.That(_matrixSparse, Is.TypeOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal, Is.TypeOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Negate()
+        {
+            Assert.That(_vectorDense.Negate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Negate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Negate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Negate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Negate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Conjugate()
+        {
+            Assert.That(_vectorDense.Conjugate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Conjugate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Conjugate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Conjugate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Conjugate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Transpose()
+        {
+            Assert.That(_matrixDense.Transpose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Transpose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Transpose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void ConjugateTranspose()
+        {
+            Assert.That(_matrixDense.ConjugateTranspose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.ConjugateTranspose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.ConjugateTranspose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Add()
+        {
+            Assert.That(_vectorDense + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense + _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse + _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Subtract()
+        {
+            Assert.That(_vectorDense - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense - _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse - _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void PointwiseMultiply()
+        {
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixMultiply()
+        {
+            Assert.That(_matrixDense*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse*_matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixVectorMultiply()
+        {
+            Assert.That(_matrixDense*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense*_vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal*_vectorSparse, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_vectorDense*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense*_matrixSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_vectorDense*_matrixDiagonal, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDiagonal, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+        }
+
+        [Test]
+        public void Append()
+        {
+            Assert.That(_matrixDense.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void Stack()
+        {
+            Assert.That(_matrixDense.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void DiagonalStack()
+        {
+            Assert.That(_matrixDense.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            // Special Case
+            Assert.That(Matrix<double>.Build.DiagonalIdentity(2, 4).DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+    }
+}

src/UnitTests/LinearAlgebraTests/Single/ReturnTypeTests.cs

@@ -0,0 +1,220 @@
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Single;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
+{
+    [TestFixture]
+    public class ReturnTypeTests
+    {
+        readonly Vector<float> _vectorDense = Vector<float>.Build.Dense(3);
+        readonly Vector<float> _vectorSparse = Vector<float>.Build.Sparse(3);
+        readonly Matrix<float> _matrixDense = Matrix<float>.Build.Dense(3, 3);
+        readonly Matrix<float> _matrixSparse = Matrix<float>.Build.Sparse(3, 3);
+        readonly Matrix<float> _matrixDiagonal = Matrix<float>.Build.Diagonal(3, 3);
+
+        [Test]
+        public void VerifyExamples()
+        {
+            Assert.That(_vectorDense, Is.TypeOf<DenseVector>());
+            Assert.That(_vectorSparse, Is.TypeOf<SparseVector>());
+            Assert.That(_matrixDense, Is.TypeOf<DenseMatrix>());
+            Assert.That(_matrixSparse, Is.TypeOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal, Is.TypeOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Negate()
+        {
+            Assert.That(_vectorDense.Negate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Negate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Negate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Negate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Negate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Conjugate()
+        {
+            Assert.That(_vectorDense.Conjugate(), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.Conjugate(), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.Conjugate(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Conjugate(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Conjugate(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Transpose()
+        {
+            Assert.That(_matrixDense.Transpose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Transpose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Transpose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void ConjugateTranspose()
+        {
+            Assert.That(_matrixDense.ConjugateTranspose(), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.ConjugateTranspose(), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.ConjugateTranspose(), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Add()
+        {
+            Assert.That(_vectorDense + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense + _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse + _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense + _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse + _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal + _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void Subtract()
+        {
+            Assert.That(_vectorDense - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense - _vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse - _vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense - _matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse - _matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal - _matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void PointwiseMultiply()
+        {
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense.PointwiseMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse.PointwiseMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.PointwiseMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixMultiply()
+        {
+            Assert.That(_matrixDense*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixSparse, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense*_matrixDiagonal, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse*_matrixDiagonal, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDense, Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixSparse, Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal*_matrixDiagonal, Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.TransposeAndMultiply(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+        }
+
+        [Test]
+        public void MatrixVectorMultiply()
+        {
+            Assert.That(_matrixDense*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense*_vectorSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse*_vectorSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal*_vectorDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal*_vectorSparse, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_vectorDense*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDense, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorDense*_matrixSparse, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixSparse, Is.InstanceOf<SparseVector>());
+            Assert.That(_vectorDense*_matrixDiagonal, Is.InstanceOf<DenseVector>());
+            Assert.That(_vectorSparse*_matrixDiagonal, Is.InstanceOf<SparseVector>());
+
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDense.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixSparse.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorDense), Is.InstanceOf<DenseVector>());
+            Assert.That(_matrixDiagonal.TransposeThisAndMultiply(_vectorSparse), Is.InstanceOf<SparseVector>());
+        }
+
+        [Test]
+        public void Append()
+        {
+            Assert.That(_matrixDense.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Append(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Append(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void Stack()
+        {
+            Assert.That(_matrixDense.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.Stack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.Stack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+
+        [Test]
+        public void DiagonalStack()
+        {
+            Assert.That(_matrixDense.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixSparse), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDense.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixSparse.DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDense), Is.InstanceOf<DenseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixSparse), Is.InstanceOf<SparseMatrix>());
+            Assert.That(_matrixDiagonal.DiagonalStack(_matrixDiagonal), Is.InstanceOf<DiagonalMatrix>());
+
+            // Special Case
+            Assert.That(Matrix<float>.Build.DiagonalIdentity(2, 4).DiagonalStack(_matrixDiagonal), Is.InstanceOf<SparseMatrix>());
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -174,6 +174,7 @@
     <Compile Include="LinearAlgebraTests\Complex32\MatrixStructureTheory.cs" />
     <Compile Include="LinearAlgebraTests\Complex32\MatrixTests.Arithmetic.cs" />
     <Compile Include="LinearAlgebraTests\Complex32\MatrixTests.cs" />
+    <Compile Include="LinearAlgebraTests\Complex32\ReturnTypeTests.cs" />
     <Compile Include="LinearAlgebraTests\Complex32\Solvers\Iterative\BiCgStabTest.cs" />
     <Compile Include="LinearAlgebraTests\Complex32\Solvers\Iterative\GpBiCgTest.cs" />
     <Compile Include="LinearAlgebraTests\Complex32\Solvers\Iterative\MlkBiCgStabTest.cs" />
@@ -222,6 +223,7 @@
     <Compile Include="LinearAlgebraTests\Complex\MatrixStructureTheory.cs" />
     <Compile Include="LinearAlgebraTests\Complex\MatrixTests.Arithmetic.cs" />
     <Compile Include="LinearAlgebraTests\Complex\MatrixTests.cs" />
+    <Compile Include="LinearAlgebraTests\Complex\ReturnTypeTests.cs" />
     <Compile Include="LinearAlgebraTests\Complex\Solvers\Iterative\BiCgStabTest.cs" />
     <Compile Include="LinearAlgebraTests\Complex\Solvers\Iterative\GpBiCgTest.cs" />
     <Compile Include="LinearAlgebraTests\Complex\Solvers\Iterative\MlkBiCgStabTest.cs" />
@@ -250,6 +252,7 @@
     <Compile Include="LinearAlgebraTests\Complex\VectorTests.cs" />
     <Compile Include="LinearAlgebraTests\Complex\VectorTests.Norm.cs" />
     <Compile Include="LinearAlgebraTests\Double\MatrixStructureTheory.cs" />
+    <Compile Include="LinearAlgebraTests\Double\ReturnTypeTests.cs" />
     <Compile Include="LinearAlgebraTests\Double\SparseVectorArithmeticTheory.cs" />
     <Compile Include="LinearAlgebraTests\Double\DenseMatrixTests.cs" />
     <Compile Include="LinearAlgebraTests\Double\DenseVectorTest.TextHandling.cs" />
@@ -320,6 +323,7 @@
     <Compile Include="LinearAlgebraTests\Single\MatrixStructureTheory.cs" />
     <Compile Include="LinearAlgebraTests\Single\MatrixTests.Arithmetic.cs" />
     <Compile Include="LinearAlgebraTests\Single\MatrixTests.cs" />
+    <Compile Include="LinearAlgebraTests\Single\ReturnTypeTests.cs" />
     <Compile Include="LinearAlgebraTests\Single\Solvers\Iterative\BiCgStabTest.cs" />
     <Compile Include="LinearAlgebraTests\Single\Solvers\Iterative\GpBiCgTest.cs" />
     <Compile Include="LinearAlgebraTests\Single\Solvers\Iterative\MlkBiCgStabTest.cs" />