diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index e4aa4434..681787a3 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -341,7 +341,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </returns>
         public override Complex At(int row, int column)
         {
-            return _data[(column * _rowCount) + row];
+            return _storage[row, column];
         }
 
         /// <summary>
@@ -358,7 +358,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </param>
         public override void At(int row, int column, Complex value)
         {
-            _data[(column * _rowCount) + row] = value;
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index 98c6040b..9d4c16c6 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -63,7 +63,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </exception>
         public DiagonalMatrix(int order) : base(order)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex>(order, order);
+            _storage = new SparseDiagonalMatrixStorage<Complex>(order, order, Complex.Zero);
             _data = _storage.Data;
         }
 
@@ -78,7 +78,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </param>
         public DiagonalMatrix(int rows, int columns) : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns, Complex.Zero);
             _data = _storage.Data;
         }
 
@@ -94,7 +94,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <param name="value">The value which we assign to each diagonal element of the matrix.</param>
         public DiagonalMatrix(int rows, int columns, Complex value) : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns, Complex.Zero);
             _data = _storage.Data;
 
             for (var i = 0; i < _data.Length; i++)
@@ -112,7 +112,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <param name="diagonalArray">The one dimensional array which contain diagonal elements.</param>
         public DiagonalMatrix(int rows, int columns, Complex[] diagonalArray) : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns, diagonalArray);
+            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns, Complex.Zero, diagonalArray);
             _data = _storage.Data;
         }
 
@@ -125,7 +125,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public DiagonalMatrix(Complex[,] array) : this(array.GetLength(0), array.GetLength(1))
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex>(array.GetLength(0), array.GetLength(1));
+            _storage = new SparseDiagonalMatrixStorage<Complex>(array.GetLength(0), array.GetLength(1), Complex.Zero);
             _data = _storage.Data;
 
             for (var i = 0; i < RowCount; i++)
@@ -165,7 +165,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override Complex At(int row, int column)
         {
-            return row == column ? _data[row] : 0.0;
+            return _storage[row, column];
         }
 
         /// <summary>
@@ -185,14 +185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override void At(int row, int column, Complex value)
         {
-            if (row == column)
-            {
-                _data[row] = value;
-            }
-            else if (((value.Real != 0.0) && !double.IsNaN(value.Real)) || ((value.Imaginary != 0.0) && !double.IsNaN(value.Imaginary)))
-            {
-                throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
-            }
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 4a6a5ebc..686fad85 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -341,7 +341,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </returns>
         public override Complex32 At(int row, int column)
         {
-            return _data[(column * _rowCount) + row];
+            return _storage[row, column];
         }
 
         /// <summary>
@@ -358,7 +358,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </param>
         public override void At(int row, int column, Complex32 value)
         {
-            _data[(column * _rowCount) + row] = value;
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index 48909ea2..a9800b3f 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
@@ -64,7 +64,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(int order)
             : base(order)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex32>(order, order);
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(order, order, Complex32.Zero);
             _data = _storage.Data;
         }
 
@@ -80,7 +80,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(int rows, int columns)
             : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns, Complex32.Zero);
             _data = _storage.Data;
         }
 
@@ -97,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(int rows, int columns, Complex32 value)
             : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns, Complex32.Zero);
             _data = _storage.Data;
 
             for (var i = 0; i < _data.Length; i++)
@@ -116,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(int rows, int columns, Complex32[] diagonalArray)
             : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns, diagonalArray);
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns, Complex32.Zero, diagonalArray);
             _data = _storage.Data;
         }
 
@@ -130,7 +130,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(Complex32[,] array)
             : this(array.GetLength(0), array.GetLength(1))
         {
-            _storage = new SparseDiagonalMatrixStorage<Complex32>(array.GetLength(0), array.GetLength(1));
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(array.GetLength(0), array.GetLength(1), Complex32.Zero);
             _data = _storage.Data;
 
             for (var i = 0; i < RowCount; i++)
@@ -170,7 +170,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override Complex32 At(int row, int column)
         {
-            return row == column ? _data[row] : 0.0f;
+            return _storage[row, column];
         }
 
         /// <summary>
@@ -190,14 +190,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override void At(int row, int column, Complex32 value)
         {
-            if (row == column)
-            {
-                _data[row] = value;
-            }
-            else if (((value.Real != 0.0) && !double.IsNaN(value.Real)) || ((value.Imaginary != 0.0) && !double.IsNaN(value.Imaginary)))
-            {
-                throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
-            }
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index 47ad0963..a5262fbc 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -341,7 +341,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public override double At(int row, int column)
         {
-            return _data[(column * _rowCount) + row];
+            return _storage[row, column];
         }
         
         /// <summary>
@@ -358,7 +358,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </param>
         public override void At(int row, int column, double value)
         {
-            _data[(column * _rowCount) + row] = value;
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 7768bec2..2ab9e782 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -62,7 +62,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </exception>
         public DiagonalMatrix(int order) : base(order)
         {
-            _storage = new SparseDiagonalMatrixStorage<double>(order, order);
+            _storage = new SparseDiagonalMatrixStorage<double>(order, order, 0d);
             _data = _storage.Data;
         }
 
@@ -77,7 +77,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </param>
         public DiagonalMatrix(int rows, int columns) : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns, 0d);
             _data = _storage.Data;
         }
 
@@ -93,7 +93,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="value">The value which we assign to each diagonal element of the matrix.</param>
         public DiagonalMatrix(int rows, int columns, double value) : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns, 0d);
             _data = _storage.Data;
 
             for (var i = 0; i < _data.Length; i++)
@@ -111,7 +111,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="diagonalArray">The one dimensional array which contain diagonal elements.</param>
         public DiagonalMatrix(int rows, int columns, double[] diagonalArray) : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns, diagonalArray);
+            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns, 0d, diagonalArray);
             _data = _storage.Data;
         }
 
@@ -124,7 +124,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public DiagonalMatrix(double[,] array) : this(array.GetLength(0), array.GetLength(1))
         {
-            _storage = new SparseDiagonalMatrixStorage<double>(array.GetLength(0), array.GetLength(1));
+            _storage = new SparseDiagonalMatrixStorage<double>(array.GetLength(0), array.GetLength(1), 0d);
             _data = _storage.Data;
 
             for (var i = 0; i < RowCount; i++)
@@ -164,7 +164,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override double At(int row, int column)
         {
-            return row == column ? _data[row] : 0.0;
+            return _storage[row, column];
         }
 
         /// <summary>
@@ -184,14 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override void At(int row, int column, double value)
         {
-            if (row == column)
-            {
-                _data[row] = value;
-            }
-            else if (value != 0.0 && !Double.IsNaN(value))
-            {
-                throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
-            }
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 80739f0c..fd0bf68f 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -341,7 +341,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </returns>
         public override float At(int row, int column)
         {
-            return _data[(column * _rowCount) + row];
+            return _storage[row, column];
         }
 
         /// <summary>
@@ -358,7 +358,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </param>
         public override void At(int row, int column, float value)
         {
-            _data[(column * _rowCount) + row] = value;
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index b64c64b4..1c6b099c 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -63,7 +63,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(int order)
             : base(order)
         {
-            _storage = new SparseDiagonalMatrixStorage<float>(order, order);
+            _storage = new SparseDiagonalMatrixStorage<float>(order, order, 0f);
             _data = _storage.Data;
         }
 
@@ -79,7 +79,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(int rows, int columns)
             : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns, 0f);
             _data = _storage.Data;
         }
 
@@ -96,7 +96,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(int rows, int columns, float value)
             : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns);
+            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns, 0f);
             _data = _storage.Data;
 
             for (var i = 0; i < _data.Length; i++)
@@ -115,7 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(int rows, int columns, float[] diagonalArray)
             : base(rows, columns)
         {
-            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns, diagonalArray);
+            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns, 0f, diagonalArray);
             _data = _storage.Data;
         }
 
@@ -129,7 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(float[,] array)
             : this(array.GetLength(0), array.GetLength(1))
         {
-            _storage = new SparseDiagonalMatrixStorage<float>(array.GetLength(0), array.GetLength(1));
+            _storage = new SparseDiagonalMatrixStorage<float>(array.GetLength(0), array.GetLength(1), 0f);
             _data = _storage.Data;
 
             for (var i = 0; i < RowCount; i++)
@@ -169,7 +169,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override float At(int row, int column)
         {
-            return row == column ? _data[row] : 0.0f;
+            return _storage[row, column];
         }
 
         /// <summary>
@@ -189,14 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
         public override void At(int row, int column, float value)
         {
-            if (row == column)
-            {
-                _data[row] = value;
-            }
-            else if (value != 0.0 && !float.IsNaN(value))
-            {
-                throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
-            }
+            _storage[row, column] = value;
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
index 5149d628..70be6cf5 100644
--- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
@@ -24,6 +24,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Data = data;
         }
 
+        public T this[int row, int column]
+        {
+            get { return Data[(column*RowCount) + row]; }
+            set { Data[(column*RowCount) + row] = value; }
+        }
+
         public void Clear()
         {
             Array.Clear(Data, 0, Data.Length);
@@ -125,6 +131,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 throw new ArgumentOutOfRangeException("columnCount");
             }
 
+            // Copy
+
             for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnMax; j++, jj++)
             {
                 //Buffer.BlockCopy(Data, j*RowCount + sourceRowIndex, target.Data, jj*target.RowCount + targetRowIndex, rowCount * System.Runtime.InteropServices.Marshal.SizeOf(typeof(T)));
diff --git a/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs
index 14b0e836..29a4f727 100644
--- a/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs
@@ -9,19 +9,41 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
         public int RowCount { get; private set; }
         public int ColumnCount { get; private set; }
         public T[] Data { get; private set; }
+        readonly T _zero;
 
-        internal SparseDiagonalMatrixStorage(int rows, int columns)
+        internal SparseDiagonalMatrixStorage(int rows, int columns, T zero)
         {
             RowCount = rows;
             ColumnCount = columns;
             Data = new T[Math.Min(rows, columns)];
+            _zero = zero;
         }
 
-        internal SparseDiagonalMatrixStorage(int rows, int columns, T[] data)
+        internal SparseDiagonalMatrixStorage(int rows, int columns, T zero, T[] data)
         {
             RowCount = rows;
             ColumnCount = columns;
             Data = data;
+            _zero = zero;
+        }
+
+        public T this[int row, int column]
+        {
+            get
+            {
+                return row == column ? Data[row] : _zero;
+            }
+            set
+            {
+                if (row == column)
+                {
+                    Data[row] = value;
+                }
+                else if (!_zero.Equals(value))
+                {
+                    throw new IndexOutOfRangeException("Cannot set an off-diagonal element in a diagonal matrix.");
+                }
+            }
         }
 
         public void Clear()
@@ -51,7 +73,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Array.Copy(Data, 0, target.Data, 0, Data.Length);
         }
 
-        public void CopyTo(DenseColumnMajorMatrixStorage<T> target)
+        public void CopyTo(DenseColumnMajorMatrixStorage<T> target, bool targetKnownClear = false)
         {
             if (target == null)
             {
@@ -64,7 +86,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 throw new ArgumentException(message, "target");
             }
 
-            target.Clear();
+            if (!targetKnownClear)
+            {
+                target.Clear();
+            }
 
             for (int i = 0; i < Data.Length; i++)
             {